A model for electric field enhancement in lightning leader tips to levels allowing X-ray and γ ray emissions

1 Introduction

The hypothesis by Wilson [1924] that postulated acceleration of electrons to high energies in thundercloud electric fields presently is worldwide recognized owing to the well-documented detections of terrestrial gamma ray flashes (TGFs) in the near space [Fishman et al., 1994; Smith et al., 2005; Briggs et al., 2010, 2011; Tavani et al., 2011] and X-ray and γ ray enhancements detected in situ experiments with balloons, aboard airplane, and at the Earth's surface at the sea level and at mountain altitudes [Shaw, 1967; McCarthy and Parks, 1985; Suszcynsky et al., 1996; Eack et al., 1996a, 1996b, 2000; Brunetti et al., 2000; Khaerdinov et al., 2005; Chilingarian et al., 2009, 2011; Chilingaryan et al., 2010; Torii et al., 2009, 2011; Tsuchiya et al., 2007, 2009, 2011, 2012]. These events, with photon spectra extending up to tens MeV, last from hundreds microsecond [Briggs et al., 2010] up to minutes [Tsuchiya et al., 2007]. They mainly precede the lightning activity, although there are evidences that some TGFs occur after the lightning electromagnetic pulses (EMPs) [Cummer et al., 2005; Cohen et al., 2006]. The origin of these high-energy phenomena is understood in the framework of experimentally proven numerical models [e.g., Lehtinen et al., 1999; Dwyer and Smith, 2005; Gurevich et al., 2007; Babich et al., 2008, 2010; Dwyer, 2012, and citations therein], assuming a development of relativistic runaway electron avalanches (RREAs) [Gurevich et al., 1992] seeded by cosmic rays in vast and rather long-living thundercloud electric fields whose strength is many times lower not only than the runaway threshold but even than the conventional self-breakdown strength at sea level.

Our paper is devoted to less documented high-energy processes in thunderstorm atmosphere, namely, flashes of penetrating photon emissions observed in direct correlation with EMPs of triggered lightning and in natural negative lightning [Moore et al., 2001; Dwyer, 2003a; Dwyer et al., 2004, 2005; Howard et al., 2008]. These flashes are many times narrower than cited above events and last for about ~1 µs. The photon energies in them are usually around 100 keV, only sometimes up to 1 MeV. The photons, obviously, are bremsstrahlung from microsecond pulses of energetic electrons produced in rather short-living and localized electric field. The observations are puzzling because energetic electrons are generally thought not to be present in conventional gas discharges. Electric discharges in the atmosphere are powered by electric fields that accelerate free electrons to energies around 15 eV, where secondary electrons are created through impact ionization of N₂ and O₂ or indirectly by photoionization of O₂ by photons emitted by excited N₂. The production and loss of free electrons depend on a magnitude of the electric field strength. The threshold electric field, E_k, for a discharge is where sources (ionization) equal losses (attachment of free electrons to O₂ molecules, recombination, etc.). At sea level pressures, the E_k is typically of 32 kV/cm. The bulk of free electrons in a discharge are usually at energies below a few eV, with the tail of the energy distribution function reaching above the ionization threshold and sustaining the discharge. Clearly, these energies are not sufficient to generate X-rays.

The electron energy losses generally increase with electron energy reaching in the atmosphere, a maximum at the electron ε~150 eV [Peterson and Green, 1968]. At higher energies the losses decrease until about 1 MeV where they start to increase energy again. It means that the frictional force, F(ε), exerted by the atmosphere through the electron collisions with air molecules, varies in the same manner. To accelerate the low-energy electrons to high energies, the electric field strength must be of 270 kV/cm [Babich and Stankevich, 1972], or ~8 E_k, to overcome the maximum of the frictional force. If such fields exist in discharges, electrons may achieve very high energies because they are accelerated past the maximum and into the so-called runaway regime, where the frictional force decreases with increasing electron energy. Such high fields are difficult to achieve, however, because the ionization creates an electrically good conductor that tends to limit the field. It has been interesting therefore to understand how lightning appears able to generate X-rays and γ rays because it opens to a new aspect of lightning discharges or more generally to the gas discharge.

The phenomenon of high-field runaway in dense gases was first observed in the pioneering series of laboratory experiments about 50 years ago in atmospheric air [Stankevich and Kalinin, 1967; Tarasova and Khudyakova, 1969; Tarasova et al., 1975] and helium at 1 atmosphere pressure [Noggle et al., 1968] using high-voltage pulses with subnanosecond risetimes, which allow overcoming the field limitation imposed by the increasing conductivity. Since then, numerous laboratory experiments have been conducted [e.g., Babich et al., 1990; Babich, 2003, 2005; Babich and Loĭko, 2009; Babich et al., 2014a, 2014b].

The first analytical estimate of the fluency of this process in weakly ionized gases was published more than half a century ago by Gurevich [1961] and was followed by more accurate studies, for instance by Babich and Stankevich [1972] which for the first time included space charge effects, Bakhov et al. [2000] and Chanrion and Neubert [2010] adopting a full Monte Carlo and particle-in-cell approach, respectively. However, the mechanism capable of accounting for the impulsive electric field enhancement at the tip of a lightning leader to a level sufficient for generation of intensive fluxes of runaway electrons (REs) is not fully understood.

One process that circumvents the need for high electric fields was proposed by Gurevich et al. [1992]. They considered seed electrons, created with energies in the runaway range by cosmic rays and further accelerated and multiplied in lightning—and thunderstorm electric fields, thereby forming a RREA. While this process may possibly power avalanches and associated emissions in clouds, more recent observations suggest that the mechanism primarily is “thermal runaway” [Gurevich et al., 2007] or, more precisely, “high-field runaway” [Dwyer, 2008], with the source region close to the lightning leader tip [Dwyer et al., 2011]. As noted above, this process is not straightforward because of the high electric fields needed.

A model of electron acceleration at lightning leader tips has been proposed by Moss et al. [2006] and further elaborated by Celestin and Pasko [2011]. It involves the two basic modes of gas discharges: the streamer and the leader. A streamer is an ionization wave that in the atmosphere propagates at velocities approaching magnitudes on the order of 10⁹ cm/s (cf., simulations by Celestin and Pasko [2011] with E = 50 kV/cm and 1 atm) and where the space charge in the wavefront can lead to magnitudes of electric field strength of several E_k. In the wake, it creates a small-scale filament (of 1 mm or less in radius in the atmosphere at sea level pressure) of enhanced ionization and electric conductivity. Many streamers are thought to be injected from a lightning leader tip and to assist its propagation. A lightning leader propagates slower, has a high degree of ionization, a high electric conductivity, and can carry hundreds of amperes of current. It is created by thermal ionization of the neutral gas in the leader channel. When negative lightning propagates (in a direction opposite to the vector of the electric field strength), it often occurs in steps. In between the steps, a negative leader segment, a so-called space leader, may be formed meters ahead [e.g., Gorin et al., 1976; Reess et al., 1995]. The stepped propagation occurs when the main leader channel connects to the space leader [e.g., Biagi et al., 2010].

In Moss et al. [2006] and Celestin and Pasko [2011], a two-step process is proposed. The first step occurs in the streamer fields. The authors believed that “while during the most of expansion dynamics of the streamer the maximum field in the streamer head stays close to ~5 E_k, just prior to the branching of the streamer, this peak field in the streamer head can reach very high values…” sufficient to accelerate thermal electrons to the runaway energy range. In the second step, the electrons are further accelerated to 10–100 keV and higher when the lightning channel connects to the space stem and transfers its potential to the stem. It then creates an impulsive high electric field at the tip of the stem, which is now the new tip of the lightning channel. The number of REs generated in a streamer corona flash was estimated to ~10¹²–10¹⁵ [Moss et al., 2006] and 10¹⁷ [Celestin and Pasko, 2011], the latter number being in the range thought to be required to explain the γ ray bursts observed from satellites, the terrestrial gamma ray flashes (TGFs) [Dwyer, 2008]. For comparison, Gurevich et al. [2007] estimated a production of ~3 × 10¹² REs in one leader step.

Because the leader-streamer interaction is a very complex process that cannot yet be simulated realistically, estimates are made that are sensitive to assumptions on parameters. Celestin and Pasko [2011], for instance, assumed that the RREA length l_av is ~7 cm and the avalanche region is L~1 m giving L/l_av~15 avalanche lengths, whereas Gurevich et al. [2007] assume l_av~1 cm and L~10 cm (for the head length), which give a total of L/l_av~10 avalanche lengths with a the field intensity >200 kV/cm (~6 E_k). However, according to Monte Carlo simulations [Babich and Bochkov, 2011; Kutsyk et al., 2012a, 2012b], the avalanche length is ~16 cm for 218 kV/cm, which decreases the number of avalanche lengths of the assumed region and thereby the number of REs.

On the other hand, the RREA can continue developing in the weaker electric fields at larger distances away from a leader tip. Therefore, it is more reasonable to avoid addressing the number of RREA lengths and directly estimate the RREA enhancement as ~ exp(U_L/< ε_re >) [Dwyer, 2004]. Here U_L is the available potential drop in the cloud, and < ε_re > is the characteristic energy of electrons in the RREA. With U_L ≈ 20 MV, which Gurevich et al. [2007] used in the region where RREA is being developed, and < ε_re > ≈ 7 MeV (valid in the range of small overvoltages δ = eE/(F_minP) below the self-breakdown limit in dry air, inherent for the weak thundercloud fields [Babich et al., 2004; Dwyer, 2003b, 2004, 2007; Kutsyk et al., 2012a, 2012b]), the above formula gives the RREA enhancement of 17. Note, here we do not use the conventional definition of the overvoltage relative to the self-breakdown magnitude E_k/P ≈ 32 kV/(cm · atm) valid for the homogeneous field in dry air; the overvoltage is defined relative to the minimum F_min ≈ 2.18 keV/(cm · atm) of the electron drag force F(ε) achieved in air at electron energy ε of 1 MeV. Such definition is convenient for describing the runaway process [Gurevich et al., 1992; Roussel-Dupre et al., 1994; Dwyer, 2003b].

Any proposed model of the increase of the electric field at the leader tip should be consistent with observations, such as those by Howard et al. [2008] taken in the nature and Kochkin et al. [2012] in the laboratory. Most likely, the EMP registered by Howard et al. [2008] is generated as a result of the current jump during the new step formation. The pulse is bipolar, i.e., consists of two peaks. The first peak is positive, lasting ~100 ns, and the second peak is negative and ~1 µs. In Figure 4 of Howard et al. [2008], the X-rays appear 0.5 µs after the positive peak maximum. In our opinion, the positive peak is connected with the growth of the leader current during the new step formation when the potential is being carried along the space leader channel, and the longer negative peak is connected with the fading current in the new leader segment. Maybe it is possible to reconcile the data with the model by Celestin and Pasko [2011] if one assumes that the time lag between the EMP and the X-ray onsets is exactly the formation time of the streamer corona. Note, however, that Howard et al. [2008] observed that the X-ray source was located below that of the EMP. According to the model by Celestin and Pasko, the positions of the sources should coincide. On the other hand, Howard et al. [2008] mention some uncertainty in the X-ray source position; therefore, one cannot be confident while interpreting their observations in the framework of any mechanism. Kochkin et al. [2012], who investigated the generation of X-rays in long spark laboratory discharges, observed that the radiation appears not at the stage of the streamer corona development but during a contact of positive and negative streamers.

A consistent model, capable of describing the high-field runaway in front of the lightning leader, should include a development of powerful RREA producing the high-energy bremsstrahlung detected in correlation with lightning leaders. It is well known that “tip of long leader possesses a potential of ~1 MV in laboratory sparks and of ~10 MV or even higher in lightning” [Bazelyan and Raizer, 2000, 2001, page 75]. Hence, the voltage in front of the leaders is sufficient for the electron energizing up to the MeV energies; however, to start the RREA development, sufficiently strong electric field is required to produce seed REs. Our goal is more limited than complete describing the high-field runaway process including the RREA development: we only plan to consider the initial stage of the process by demonstrating that electric field at the leader front locally can achieve magnitudes allowing low-energy electrons to overcome the runaway energy threshold.

2 A New Mechanism of Runaway Electron Generation

According to the deterministic definition, the electric field in the atmosphere at sea level pressures must be E > E_max ≈ 260–270 kV/cm ≈ (8.1–8.4) E_k for the low-energy electrons be capable of overcoming the maximum frictional force at ~150 eV and efficiently populate the runaway energy range [e.g., Babich and Stankevich, 1972; Babich et al., 1990; Babich, 2003, 2005]. At lower E fields, a deterministic runaway threshold can instead be introduced ε_th,run(E) as the minimum energy required for electrons to be in the runaway regime at given E. It is usually defined simply from the equation F(ε) = eE [e.g., Gurevich, 1961; Babich and Stankevich, 1972; Babich et al., 1990; Babich, 1995]. More precise definitions of ε_th,run are available in Bakhov et al. [2000], Babich and Bochkov [2011], and Kutsyk et al. [2012a, 2012b].

In the stochastic approach, a probability can be calculated for low-energy electrons to be accelerated into the runaway regime [Bakhov et al., 2000; Chanrion and Neubert, 2010]. As E decreases from E_max (where most electrons will be—or are—in the runaway regime), the combined effect of the increase in ε_th,run (further out in the tail of the thermal electron energy distribution) and increasing difficulty of gaining sufficient energy, the probability for low-energy electrons to be accelerated past ε_th,run dramatically decreases. For instance, in air at 1 atm, the probability is 10⁻⁷ and 4 · 10⁻¹², respectively, for E = 6 E_k (~190 kV/cm) and E = 5 E_k (~160 kV/cm), in spite of a rather low corresponding energy threshold of ε_th,run ~277 eV and 470 eV [Chanrion and Neubert, 2010]. Similarly, studying the runaway process in pure nitrogen at 1 atm [Bakhov et al., 2000], a rate was computed of thermal electrons transition to the runaway mode, i.e., a rate at which thermal electrons are accelerated past so called “stochastic runaway threshold ε_th,run,stoch” defined as an energy, above which the probability of an electron to return to the low-energy range is extremely low. Specifically, with ε_th,run,stoch = 8 keV, the computed rate decreases from 8 · 10¹⁰ s⁻¹ at 400 kV/cm to 10⁴ s⁻¹ at 240 kV/cm [Bakhov et al., 2000]. At weaker fields, the dependence on E is even stronger. Therefore, for the thermal electron acceleration into the runaway mode to be significant in the framework of the stochastic description of the process, strong fields with intensities of 200 kV/cm (6.3 E_k) or higher are required [Bakhov et al., 2000; Chanrion and Neubert, 2010].

The mechanism for a generation of such field magnitudes in the vicinity of the lightning channel tip is at present speculative. According to Bazelyan and Raizer [2000, 2001, page 75], the field at the head of the leader channel does not significantly exceed the self-breakdown limit of ~30 kV/cm and is even less inside the channel. It is well known that strong fields are generated in heads of the streamers constituting the streamer corona of the leader; however, it is not clear what the lifetime of such fields is and therefore their efficiency in acceleration, although estimates have been made [Celestin and Pasko, 2011].

Furthermore, assuming a hemispherical approximation of the ionization wavefront D'yakonov and Kocharovsky [1988, 1989] showed that the maximum value of the normal component of the field strength at the streamer front is of the order of a magnitude E_inf, at which the ionization rate by electron impacts ν_ion(E), as a function of E, experiences an inflection (E_inf~5 E_k in air), where sharp ν_ion(E) increase with E changes to extremely weak growth with E. The maximum is first met in the region on the hemispherical space charge sheath that is normal to the background field (the direction of propagation). A further increase in the streamer sheath field leads to an increase of the ionization in the transverse direction relative to the direction of the streamer propagation such that E_m,tip~E_inf is met over a larger region of the sheath including the lateral sides. The result is that the velocities of the ionization wave in the longitudinal and transversal directions become almost equal and that the streamer experiences strong transversal expansion, eventually limiting the field intensity at the leader front [Raizer and Simakov, 1998]. This effect has also been demonstrated in numerical simulations [Vitello et al., 1994; Babaeva and Naidis, 1996; Kulikovsky, 1997, 1998; Celestin and Pasko, 2011]. The overall conclusion from these works is that the steady state magnitude of the field intensity during streamer development in a homogeneous medium does not exceed 150 − 170 kV/cm ≈ (4.7–5.3)E_k at sea level pressures, for which the runaway probability is as low as 4 ⋅ 10^− 15 − 10^− 10 [Chanrion and Neubert, 2010].

In the model that apparently allows the required fields to develop was proposed by Celestin and Pasko [2011]. They assume that REs are generated during a new leader step formation, where propagation of the electrical potential along the leader channel leads to a formation of streamer corona at a lateral channel surface. During branching of the streamers constituting the corona, the field is being increased up to ~240 kV/cm at the time of branching.

However, we consider another approach by noting that if the medium is nonuniform from preionized channels formed by previous streamers, the limitation of the electric field magnitude is removed, because a preionized channel tends to focus the ionization wave. A reason of the focusing due to the ionization nonuniformity is as follows. It is known that the lower is the electron density, in which the ionization wave propagates, the lower is the ionization wave velocity. Obviously, the electron density decreases in the direction from the inner areas of the preionized channel to its periphery. Hence, the further is the locus of considered point from the channel axis, the slower develops the ionization process in it. Thereby, the transverse expansion of the ionization wavefront occurs much slower than the longitudinal propagation; therefore, the radius of the front remains to be limited almost to that of the already existing preionized channel, and as a result, the field intensity at the wavefront can reach magnitudes of 240 kV/cm or higher required for the efficient RE generation [Bakhov et al., 2000].

We recall that negative lightning usually propagates in steps with a streamer corona at the tip of the leader as observed in laboratory gas discharges [Bazelyan and Raizer, 2000, 2001, page 77]. At the external boundary of the streamer corona, a plasma is formed, extending in the direction along the electric field vector, from both ends of which streamers start in the directions to and from the leader head (bi-leader). Further, from this plasma, the so-called space leader is formed as described earlier. The space leader with streamer coronas is illustrated in Figure 1a. The positive head of the space leader moves toward the negative head of the main leader channel and at contact, a united conducting channel is being formed where charges are redistributed and the potential of the old leader U_L is transferred to the space leader. The same conceptual framework was described earlier and adopted by Moss et al. [2006] and Celestin and Pasko [2011]. It is shown in Figure 1b. The potential jump in the body of the space leader will launch new ionization waves propagating along old streamer channels in the negative streamer corona region of the new leader tip, shown in Figure 1c.

We propose therefore that the electric fields can reach amplitudes required to generate bursts of energetic electrons and associated bremsstrahlung in the streamer channels of the space leader. We note that this may happen by the same physical mechanism in the streamers at both ends of the space leader, that is, (1) when opposite polarity streamers, growing toward each other from the main and space leaders, connect and (2) when the potential of the main leader reaches negative pole of the space leader. In the following, we demonstrate the latter case by modeling an ionization wave in a preformed conducting streamer channel.

3 Mathematical Formulation of the Problem

Given the present lack of knowledge of the physical conditions in the leader channels of lightning and of leader-streamer interactions as mentioned earlier, we limited ourselves to numerical simulations of a generalized channel model relevant for the stage illustrated in Figure 1c, namely, the propagation of a secondary ionization wave along old streamer channels at the time of the potential jump in the body of the space leader. To avoid misunderstanding, it is noted that the earlier stage of collision of the leader and space leader, as in Mallios et al. [2013], is not considered. The simulations start at the termination of the motion of the leader step which generates a sharp increase of the external field, Е₀, and a new ionization wave along old streamer channels is initiated.

As according to estimations by Bazelyan and Raizer [2000, 2001, page 76], a radius of the streamer zone (corona) of laboratory long leaders is limited to the magnitudes of ~1.5 m at the leader tip potential 1.5 MV, the radius of the streamer zone of the lightning leader, likely, can reach tens of meters. Unfortunately, computers still only allow simulations of streamers a few centimeters in length (at sea level pressures); hence, only the initial stage of the ionization wave propagation can be simulated. According to Bazelyan and Raizer [2000, 2001, page 75], the upper limit of the field at the head of the leader channel does not exceed significantly the self-breakdown limit ~30 kV/cm; we assume that voltage jump up to U_L generates at the front of new leader step (Figure 1c) an electric field of 50 kV/cm. This magnitude is taken to be the initial condition on the external field strength throughout the simulation domain.

The problem is formulated in the cylindrical coordinates (z, r) as shown in Figure 2. The simulation domain has the longitudinal and transversal dimensions L₀ × R₀ relative to the vector of static and homogeneous external electric field given by , where is a unit vector along the symmetry axis z of the leader. At time t = 0, this domain is bridged by an ionized channel, where most of the electrons are attached to oxygen molecules. In view of the adopted cylindrical symmetry, the problem is two dimensional.

The kinetics of low-energy electrons (e), positive (p), and negative (n) ions in the discharge is described by the following set of the equations:

(1)

where n_e, n_p, and n_n are the densities; , , and are the drift velocities, μ_e and μ_p,n are the electron and ion mobilities; D_e is the electron diffusion coefficient; ν_ion is the ionization rate of the air molecules by impacts of low-energy electrons; β_ep and β_pn are, respectively, coefficients of recombination of electrons with positive ions and positive and negative ions; is the rate of the electron attachment to the oxygen molecules; k_diss and k_thr are the rates of dissociative and three-body attachment; N and are the local densities of air and oxygen molecules at point z; and S_ph is the photoionization source.

The temperatures of the negative ions and the gas as a hole are assumed not to increase during the development time of the simulated process. Therefore, the rate of collisions capable of destroying negative ions is too low for the electron detachment to have significant effect on the electron density (see Neubert and Chanrion [2013] for further discussion on electron detachment).

The initial densities of the charged particles in the old streamer channel are set be as follows:

(2)

Here and are the initial densities of electrons and positive ions at the streamer channel symmetry axis, and l_r is the characteristic width (“radius”) of the radial distribution.

To start the ionization wave at the initial moment of time, a cloud of neutral plasma (disturbance) was set in the vicinity of the coordinate onset (Figure 2), consisting of electrons and positive ions with Gaussian spatial distribution of the densities:

(3)

where the amplitude is set equal to 10²⁰ m⁻³ throughout the simulations. Merging equations 2 and 3 gives the initial conditions for the set of equation 1:

(4)

The set in equation 1 was closed by equations for the self-consistent electric field:

(5)

where ϕ_int is the potential of the space charge field produced in the process of the plasma neutrality violation, ρ_q = e ⋅ (n_p − n_n − n_e) is the charge density, ε₀ is the dielectric permittivity of free space, and e is the elementary charge.

For solving the first equation in 5, the magnitudes of the potential at the border of the computational area are required. The potential in a boundary point was computed using the general solution of Poisson equation:

(6)

The integration is conducted over a region V_ρ where an equality |ρ_q| ≥ 0.001 ⋅ |ρ_q|^max is satisfied. Here |ρ_q|^max is the maximum value of the charge density within the computational domain [Bourdon et al., 2007].

To calculate the photoionization source S_ph, the model by Luque et al. [2007], developed further by Bourdon et al. [2007], was used. In this model, the computation of the integral in the classical model by Zheleznyak et al. [1982] is replaced by solving three Helmholtz equations:

(7)

where , p = 1 atm is the air pressure, is the oxygen partial pressure, and p_q is the quenching pressure [Zheleznyak et al., 1982]. The other symbols are defined in Bourdon et al. [2007].

At the boundary of the computational area, the photoionization source was calculated as follows:

(8)

where

The expression 8 was derived using the general solution of Helmholtz equation:

(9)

When deriving equation 8, we decomposed the function in equation 9 over the Taylor series in the vicinity of point and retained only the zero term. This procedure is the same used in electrostatics while decomposing the electric potential over the multipoles. The integration is conducted over the entire volume of the simulation domain V_sim.

The number of electrons in the runaway mode N_re(t) at the moment of time t was computed as follows:

(10)

where ν_run(E) is the rate of transition of low-energy electrons to the runaway mode. We use a dependence of the runaway rate on the field intensity calculated by Bakhov et al. [2000] in nitrogen of atmospheric density in the range of the field from 240 to 400 kV/cm for various magnitudes of the stochastic runaway threshold ε_th,run,stoch.

Computations were executed assuming STP conditions with the oxygen partial pressure  = 0.2 atm, the quenching pressure p_q = 0.0395 atm [Bourdon et al., 2007], L₀ = 4 cm, R₀ = 6 mm, E₀ = 50 kV/cm, and l_z = l_r = 500 µm. The magnitude of the initial electron density at the channel symmetry axis was varied between simulations.

Initially, we need to calculate the magnitude. With this goal, a preliminary problem is solved using the equation 1 without the flux terms; we let the field inside the old streamer channel be E = 10 kV/cm and densities at the channel symmetry axis  = 10²⁰ m⁻³. The simulations are conducted until the instant of time, , at which the electron density at the old streamer channel symmetry axis achieves the given magnitude. Typical values of were found to be ~1 µs with . Further, while solving the full set of equation 1, we consider  = 1 µs as the initial moment of time ( ).

For the dependencies of μ_e, μ_p,n, D_e, ν_ion, β_ep, β_pn, k_diss, and k_thr on the field intensity and pressure, we follow Bochkov et al. [2013]. To calculate the photoionization rate, we find the parameters λ_j and A_j (j = 1–3) in Bourdon et al. [2007] and the dependence of the quantity ξν_u/ν_ion on the field intensity in Liu and Pasko [2004].

The RE number N_re(t) was computed with the stochastic runaway threshold ε_th,run,stoch = 8 keV using the Monte Carlo data on the runaway rate ν_run(E) available in Bakhov et al. [2000], which with an error, not exceeding 13%, is approximated as follows:

(11)

4 Results and Discussion

The set in equation 1 was solved using the finite differences technique. To get rid of numerical diffusion while calculating the flux terms, the technique of flux-corrected transport [Zalesak, 1979] was used. The QUICKEST 3 [Leonard, 1991] and “donor cell” schemes were used, respectively, as the high-order and low-order schemes.

The time step chosen was Δt = 5 · 10⁻¹³ s. Adaptive mesh was accepted for the axis z: in the domain of the ionization wavefront with size of 1 cm, the spatial step was Δz = 10⁻⁵ m, and in the remaining area away from the front, the step was Δz = 5 · 10⁻⁵ m. The domain with step Δz = 10⁻⁵ m moved consistently with the ionization wavefront propagation. Along the r axis, the mesh was uniform with constant step Δr = 10⁻⁵ m. The simulations were terminated at the moment of time , such that a position of the maximum of the field intensity appeared at the point with coordinate z_f = 3 cm.

The results for  = 10¹⁰–10¹⁵ m⁻³ are presented in Table 1, where E_f is the maximum value of the field intensity at the ionization wavefront, υ_f is the propagation velocity of the front at , and k_RE is the number of REs produced by the ionization wave per unit length at the final section of its path from 2.9 to 3.0 cm:

(12)

Table 1. Results of Computations

, m−3	, ns	Ef, kV/cm	υf, m/s	kRE, 1/m
1010	3.86	205	1.3 · 107	2 · 104
1011	3.80	267	2 · 107	5 · 107
1012	3.66	296	2 · 107	2 · 109
1013	3.50	286	3 · 107	9 · 108
1014	3.31	253	3 · 107	2 · 107
1015	3.14	204	3 · 107	3 · 104

The field intensity and electron density on the symmetry axis at various moments of time are shown in Figure 3 for  = 10¹² m⁻³. Two-dimensional distributions at  = 3.6 ns are illustrated in Figure 4. It is seen that, although the primary channel is significantly expended as a result of the ionization wave development, the transverse size of the front domain with the field intensity close to the runaway threshold ~240 kV/cm remains almost the same as the characteristic size of the primary channel l_r = 500 µm (cf., Figure 2).

It is seen from the data in Table 1 that with increasing preionization level, given by , the maximum value of the field intensity E_f at the ionization wavefront increases as we expected. However, when gets too high, the initial electron conductivity tends to reduce the accumulated space charge density at the wavefront; as a result, E_f decreases with increase as seen for  > 10¹³ m⁻³. This decrease has been observed in other works that look at the effects of the increase of the photoionization level with altitude [Liu and Pasko, 2004; Luque et al., 2007; Chanrion and Neubert, 2008] or at the effect of the ionization from high-energy electrons ahead of the streamer [Chanrion et al., 2014]. We also point out that the photoionization in air at ground pressure produces electron density of 10¹³ m⁻³ at the vicinity of the wavefront [Bazelyan and Raizer, 1998, 1997, page 92] which is within the range of the used magnitudes of the preionization level, . However, the ionization by electron impacts in the front dominates because we analyze the proforce ionization wave, i.e., propagating along .

It is seen from Table 1 that with  = 10¹¹–10¹⁴ m⁻³, the field intensity is ~250–300 kV/cm and within the range required for efficient RE generation; on the order of ~10⁷–10⁹ electrons/m are estimated with energies above the stochastic runaway threshold ε_th,run,stoch = 8 keV [Bakhov et al., 2000]. From these values, it is possible to estimate a number of REs generated by one leader step. Note, however, that two uncertainties are inevitable, for such estimation. First, the number of streamers with front field sufficiently strong for the RE generation a priori is not known—only rather uncertain estimates are possible. Second, there is an uncertainty connected with the extremely strong dependence of the number of REs on the field intensity, which, obviously, at each instant of time, strongly varies in different streamers constituting the corona. In fact, the runaway probability varies 7 orders of magnitude for E/E_k increasing from 6 to 9 [Chanrion and Neubert, 2010].

Proceeding, with these cautions in mind, to estimate the number of REs expected from a lighting leader step, we note that according to Bazelyan and Raizer [2000, 2001, page 77], about N_str ≈ 1.2 · 10⁵ streamers develop in a streamer zone within the radius of R_s ≈ 1.5 m of a leader head at a potential U_L = 1.5 MV. The lifetime t_life of one streamer propagating with a minimum velocity υ_s,min~10⁵ m/s is of t_life~R_s/υ_s,min~10 µs. Only in channels of relatively young streamers, developing during time in the range  1.4–2.3 µs (corresponding to  = 10¹⁴–10¹¹ m⁻³ at which E_f > 240 kV/cm) up to a length  ≈ 10–20 cm, the electron density is sufficient to start new ionization wave-generating REs. Note that we let streamer radius changes weakly relative to the initial radius l_r = 500 µm. If the streamer becomes much wider, we can assume that the streamer branches. The radii of the branches are less than that of the parent streamer; therefore, the processes of field enhancement and electron acceleration continue at the branch heads.

With the above considerations, we estimate the number of REs generated in a leader step N_RE ≈ l_s·k_RE ·  N_str~10¹², where we have assumed that the ionization waves propagate about 10% of the streamer length,  ≈ 0.1. The leader potential is probably considerably higher than the 1.5 MV used by Bazelyan and Raizer [2000, 2001, page 76]; therefore, the number of REs generated by lightning discharges is likely larger. Note that the above estimation, N_RE ≈ 10¹², is rather close to estimated by Babich et al. [2013] numbers 10¹⁰ − 4 ⋅ 10¹¹ of high-energy REs generated by one leader step, which are required for reproducing the X-ray energy of ~1–2 MeV measured by Dwyer et al. [2005].

We only simulated a single streamer, with ambient empty space around it. In reality, the corona consists of multiple streamers developing simultaneously. There is a chance that they affect each other and modify the electric field. Ignoring this fact is a shortcoming inherent for the available models [Moss et al., 2006; Celestin and Pasko, 2011; Bazelyan and Raizer, 2000, 2001, chap. 2], considering a development of a single streamer in the corona. The single-streamer approximation is substantiated by Bazelyan and Raizer [2000, 2001, page 78], who estimated that N_str ≈ 1.2 · 10⁵ streamers developing in a streamer zone with size R ≈ 1.5 m are sufficiently distant from each other and can be considered as lonely, developing in the self-consistent electric field combining the front leader field and the field of the considered streamer. Nevertheless, accurate analyses are required of multistreamer problem.

5 Conclusion

We have shown that the inhomogeneous environment of a negative lightning leader may lead to impulsive electric field enhancement to magnitudes allowing acceleration of low-energy electrons to runaway energies. The high fields are created because the ionization wave speed depends not only on the ionization rate ν_ion(E) but also on the ionization degree of the ambient plasma, being generally lower at lower plasma densities. If preionized streamer channels are present in front of the leader, this effect tends to focus the ionization wave to the channel, thereby enhancing the peak field in the ionization wave. The inhomogeneity is important for transverse scales of the order of the radius of a streamer head dimension.

Numerical simulations were executed in a self-consistent electric field allowing for the drift motion of electrons and positive and negative ions, and for electron diffusion, impact ionization of air molecules by low-energy electrons, recombination, attachment of low-energy electrons to the oxygen molecules, and photoionization of the air. We demonstrated that with characteristic size 500 µm of the initial spatial distribution of the electron density with maximum magnitudes in the range from 10¹¹ to 10¹⁴ m⁻³, an electric field with intensity 250–300 kV/cm, corresponding to E/E_k~8–10, is produced in an ionization wave launched by an impulsive increase in the electric potential. The enhancement is sufficient for electrons to reach the runaway mode and, consequently, for producing the X-ray and γ ray pulses observed in correlation with the lightning leader steps [Dwyer et al., 2005]. The estimated number of runaway electrons is consistent with the numbers required [Babich et al., 2013] for reproducing the X-ray energy generated by stepped leader [Dwyer et al., 2005]. We predict that the RE numbers generated by lightning exceed 10¹² by a large factor.