Physical mechanism of initial breakdown pulses and narrow bipolar events in lightning discharges

1.1 Phenomenology of IBPs

The initial breakdown stage of a lightning flash encompasses its first several to tens of milliseconds, and it is characterized by a sequence of pulses typically detected with electric field change sensors on the ground [e.g., Clarence and Malan, 1957; Kitagawa and Brook, 1960; Weidman and Krider, 1979; Bils et al., 1988; Villanueva et al., 1994; Shao and Krehbiel, 1996]. Recent results by Marshall et al. [2014a] suggest that these initial breakdown pulses (IBPs) should be observable in all lightning discharges. A “classical” IBP has a duration of tens of microseconds, and it is one of the largest pulses at the beginning of both cloud-to-ground (CG) and intracloud (IC) flashes [e.g., Nag et al., 2009]. Figure 1 provides a model-based depiction of a typical IBP waveform observed experimentally (further details are presented in section 10). Marshall et al. [2014a] performed a statistical analysis of IBPs in CG flashes in three different locations. They reported average amplitudes (range normalized to 100 km) of approximately ∼1–3 V/m, and these amplitudes are a fraction of ∼0.2–0.5 of the return stroke amplitude. In IC lightning, classic IBPs are usually the highest amplitude pulses in an entire flash [Bils et al., 1988; Villanueva et al., 1994]. Nonetheless, Nag et al. [2009] demonstrated that a significantly larger number of low-amplitude short-duration IBPs occur during the initial states of lightning discharges and that classical IBPs are only a small fraction of the total number of IBPs occurring in a given lightning flash. They call these lower-amplitude shorter-duration pulses narrow IBPs. They point out that ∼22–26% of the IBPs have duration shorter than 1 μs [Nag et al., 2009].

Comparing IBP trains of CG lighting that were and were not followed by a return stroke, Nag and Rakov, [2008, 2009] concluded that IBPs (in −CG lightning) are formed when the downward negative leader extends from its initiation region, just below the main negative charge region, and encounters an appreciable lower positive charge region in the thundercloud. These authors speculate that when the lower positive charge region is small IBPs are unlikely to be produced [Nag and Rakov, 2009, Figure 3]. The simultaneous observation, using both high-speed cameras and fast electric field antennas, of the initial lightning development just below the cloud base [Stolzenburg et al., 2013; Campos and Saba, 2013] has indicated that IBP trains are related to bursts of light that appear to retrace the existing lightning channel. Multistation records of IBPs [Karunarathne et al., 2013] have determined that as a flash evolves the location of IBP sources inside the cloud coincide with the position of negative leaders as determined by a VHF lightning mapping system.

1.2 Phenomenology of NBEs

Another class of electric field signatures that are not well understood at present is referred to as narrow bipolar events (NBEs). Legacy studies showed that NBEs are a rather unusual kind of intracloud discharge, typically occurring in complete isolation of any other kind of lightning process inside the thundercloud [Le Vine, 1980; Willett et al., 1989]. However, more recent investigations show that a considerable fraction occur as the initiation process in otherwise normal lightning discharges [Rison et al., 1999; Thomas et al., 2001; Wu et al., 2014]. For example, Wu et al. [2014] found that 103 out of 638 (∼16%) recorded positive NBEs were followed by lightning. NBEs are known to be the “strongest sources of radio frequency radiation from lightning” [Le Vine, 1980]. At VLF/LF frequencies (3–300 kHz) their amplitudes are comparable to typical (range normalized to 100 km) return stroke amplitudes, of ∼6 V/m [Rakov and Uman, 2003, page 154]. For instance, using an array of field change sensors [e.g., Kitagawa and Brook, 1960], Smith et al. [1999] recorded amplitudes of ∼10 V/m. However, it is at HF/VHF (3–300 MHz) that these discharges are particularly intense. Thomas et al. [2001], for example, reported a lightning mapping array observation of a NBE with source peak power of 300 kW in the 60–66 MHz passband of the receivers. This value contrasts with typical <10 kW source amplitudes associated to the stepping of negative lightning leaders inside the thundercloud [Rison et al., 1999; Thomas et al., 2001].

These electric field signatures have received their name because they are characterized by a short-duration bipolar waveform. A typical (far-field) NBE waveform “resembles a single full cycle of a distorted sine-wave, with a risetime of 1–2 μs, then falling and crossing zero in a few microseconds, followed by the second half cycle, which is lower amplitude but longer-duration compared to the first half cycle; the typical NBE full-width is 10–20 μs” [Jacobson and Light, 2012]. Figure 2a provides a model-based depiction of a typical (far-field) NBE waveform observed experimentally (further details are presented in section 11). Average values reported in the literature for the range-normalized peak amplitude, for the ratio of the initial electric field peak to the opposite polarity overshoot, for the total pulse duration, and for the duration of the first half cycle, are ∼8–21 V/m, ∼2–10, ∼24–26 μs, and ∼2–5 μs, respectively [Willett et al., 1989; Smith et al., 1999; Nag et al., 2010]. NBEs occur with both positive and negative polarity, with typically negative ones occurring at higher altitudes [Smith et al., 2004]. The most common ones are of positive polarity and occur between ∼7 and 16 km altitude, while negative between ∼15 and 20 km [Smith et al., 2004; Wu et al., 2012].

The intense VHF content of NBEs has made possible their detection by satellites, first by the Blackbeard VHF receiver onboard the ALEXIS satellite [Holden et al., 1995; Massey and Holden, 1995] and later by the FORTE satellite [Jacobson et al., 1999]. VHF NBE emissions are observed by satellites in the form of “transionospheric pulse pairs” or TIPPs [Holden et al., 1995]. TIPPs consist of two pulses, each a few microseconds long, separated by tens of microseconds [Jacobson et al., 1999]. The second pulse is caused by the reflection from the conducting Earth of the electromagnetic radiation emitted by the NBE source [Smith et al., 1999]. The FORTE satellite carried both a VHF and an optical detector. A puzzling finding provided by these instruments is that the strong NBE VHF radiation is accompanied by very little light emission, in comparison to other lightning discharge processes [Light and Jacobson, 2002; Jacobson et al., 2013].

Some publications refer to the source of a NBE as a “compact intracloud lightning discharge” [e.g., Nag et al., 2010; Rakov, 2013]. In the present work we use the terminology “narrow bipolar event” [e.g., Eack, 2004; Jacobson et al., 2013] (equivalent to narrow bipolar pulse [e.g., Smith et al., 1999; Thomas et al., 2001]), because it describes the electric field change waveform, as shown in Figure 2a. This choice is motivated by the analogy with the name “initial breakdown pulse,” which also refers to the electric field change phenomenology (Figure 1). Later in this paper we relate both IBPs and NBEs to the same phenomenon inside the thundercloud: the development of the initial lightning leader (more details are in section 7).

1.3 Previous Modeling of IBP and NBE Sources

There are two classes of theoretical models for the source of electrical current, I(z,t), that generates IBP and NBE pulses. The first one is referred to as transmission line model and it was applied to make inferences about the sources of both IBP and NBE signatures [e.g., Shao and Heavner, 2006; Watson and Marshall, 2007; Nag and Rakov, 2010; Karunarathne et al., 2014]. The second one is based on the assumption that relativistic runaway electron avalanches seeded by cosmic ray air showers are the source of NBEs [e.g., Gurevich et al., 2002, 2004; Gurevich and Zybin, 2004; Tierney et al., 2005; Arabshahi et al., 2014]. Having obtained the current distribution I(z,t) using either modeling strategy, the electric field generated by these sources can be calculated and compared to measurements. The vertical electric field just above the surface of the perfectly conducting ground, at a distance D from the source, can be conveniently expressed as a sum of three components [Uman et al., 1975]:

(1)

where , , and c₀ is the speed of light in vacuum. In equation 1 it is assumed that the source current is distributed between heights h₁ (lower) and h₂ (upper). The three terms on the right-hand side of equation 1 are commonly referred to as electrostatic, induction, and radiation fields, respectively. For distances far away from the source, the amplitude of these components decrease with distance as 1/R³, 1/R², and 1/R, respectively.

Transmission line (TL) models for IC discharges are an extension of the TL models developed for simulation of lightning return strokes [Uman and McLain, 1969] (detailed reviews on TL return stroke models are available in the literature [e.g., Baba and Rakov, 2007]). In this framework the existence of a conducting channel inside the thunderstorm is assumed and a current pulse is injected at one of its extremities. For instance, consider the existence of a conducting channel between heights h₁ and h₂ [e.g., Nag and Rakov, 2010, Figure 5]. Consider also that a current pulse I(h₁,t) is injected at the channel's lower extremity at height h₁. Disregarding effects of reflection at the upper end, the solution for the current distribution on the channel is

(2)

where f(z) is a function that prescribes the attenuation (or amplification) of the current wave and v is its speed. There are several variations of this model depending on the functional dependence of f(z). Watson and Marshall [2007] and Karunarathne et al. [2014] employed these models to match multistation measurements of NBEs and IBPs, respectively. They have used several variations of the TL model choosing f(z) to be constant, linearly decreasing, exponentially decreasing/increasing, and following the Kumaraswamy's distribution (which is zero at the channel's extremities and peaks somewhere in the middle of the TL). These authors demonstrated that, for a specific pulse measured, different choices of f(z) lead to (order of magnitude) similar inferred parameters for the current source. For example, they inferred that the source of both NBEs and IBPs should be ∼100–1000 m long [Karunarathne et al., 2014, Table 4].

Equation 2 is a general solution for the wave equation for an infinite TL [e.g., Bazelyan and Raizer, 2000, pages 175–176]. The limitation of this approach is that it neglects reflections at the channel's termination. For instance, for electromagnetic waves propagating on the TL at approximately the speed of light, reflections will occur in a time scales of ∼0.3–3 μs, for channels ∼100–1000 m long, respectively. This time scale is comparable or shorter than the above mentioned duration of both IBPs (section 1) and NBEs (section 2). Nag and Rakov [2010] introduced the effects of reflections by using effective reflection coefficients at the TL extremities and writing the total current wave as a summation of an infinite number of waves that are injected in the TL every time the original wave reaches one of the TL terminations. A detailed description of this technique can be found in engineering electromagnetics textbooks [e.g., Inan and Inan, 1998, section 2.3]. A weakness of TL models is that they do not explain what is the physical source of the current pulse injected in the channel and what process creates the channel in the first place. It should be noted that both Watson and Marshall [2007] and Nag and Rakov [2010] suggest, but do not substantiate using physics-based modeling, that runaway breakdown might have created the initial channel in which the NBE current pulse was injected. This assumption is based on a different modeling strategy that is described next below.

The second class of NBE models are based on the idea that high-energy electrons, from cosmic ray extensive air showers (EASs), reaching thundercloud altitudes, may act as seed particles for relativistic runaway electron avalanches (RREAs) [Gurevich and Zybin, 2004]. A requirement for RREA multiplication is that the electric field within the thunderstorm is larger than the runaway breakdown threshold in kilometer-long scales. This threshold is ∼3 kV/cm at sea level and it decreases exponentially with altitude following the ambient air density [e.g., Gurevich et al., 2002; Dwyer, 2003]. Depending on the energy of the seed particles and the thunderstorm field strength the runaway electrons can rapidly multiply and produce a large number of low-energy electrons. These low-energy electrons would drift in the thunderstorm field and produce a current distribution that radiates in the LF/VLF frequency range. Its electromagnetic radiation would have a characteristic bipolar waveform, similar to NBEs. However, as suggested by Dwyer and Uman [2014, section 6.5] and further quantified by Arabshahi et al. [2014], the required thunderstorm and EAS properties to match observed NBE waveforms appear to be unrealistic. For instance, Arabshahi et al. [2014] demonstrated that the required minimum cosmic ray energy to reproduce measured NBEs through the EAS–RREA model is ∼10²¹ eV. These authors estimate that the frequency of occurrence of an EAS with such high energy, over a typical thunderstorm area of 100 km², is 1 every 67 years. These authors also concluded that, required for this mechanism, thunderstorm electric field magnitudes should be comparable to the conventional breakdown threshold [Arabshahi et al., 2014, Figure 2]. The conventional breakdown threshold is defined by the equality between electron-impact ionization and attachments rates, it has a value of ∼30 kV/cm at ground level, and it decreases exponentially with altitude, following air density. It should be noted that such high electric fields have never been measured inside thunderstorms [e.g., Dwyer and Uman, 2014, Table 3.1].

1.4 Purpose of This Work

In view of the above listed properties of IBPs and NBEs, and strengths and weaknesses of existing theoretical models, we propose in this article a new physics-based model approach. In section 6 we summarize principal kinds/polarities of pulses and put them in context of the large-scale thunderstorm spatial charge distribution. In section 7 we propose that both IBPs and NBEs may be viewed as the electromagnetic transient associated to the sudden (i.e., stepwise) elongation of the negative leader tip in the ambient thunderstorm field. In section 8 we introduce a new numerical model to simulate the initial stage of the lightning leader tree development. The model is built on a bidirectional (zero-net-charge) lightning leader concept, first envisioned by H. W. Kasemir in 1950 [Kasemir, 2013]. We simulate a finite-length finite-conductivity leader elongating in the thunderstorm electric field, and we solve a set of integrodifferential equations to retrieve the full dynamics of charges and currents induced in it. Hence, our proposed approach consists of a generalization of the transmission line [e.g., Nag and Rakov, 2010] and electrostatic [e.g., Pasko, 2014] approximations used for analysis of in-cloud discharge processes. In sections 10–11 we present the full dynamics of charges and currents in the elongating leader that emit electromagnetic radiation observed on the ground in the form of IBPs and NBEs, and we show how IBP and NBE waveforms are related to the leader properties. Finally, in sections 13–16 we discuss the implications of our results to a number of different phenomena including: conventional lightning, transient luminous events, and terrestrial gamma ray flashes.

2.1 Thunderstorm Charge Distribution and Source Heights of IBPs and NBEs

Figure 3 shows a typical representation of a normal-polarity thunderstorm charge structure. The simplified picture is built on the general view inferred by Stolzenburg et al. [1998, Figure 3] from nearly 50 electric field soundings in different kinds of thunderstorms. The charge configuration consists of the well-known tripolar charge distribution [Williams, 1989; Rakov and Uman, 2003, page 69] supplemented with a screening charge layer on the top [MacGorman and Rust, 1998, pages 51 and 72–74]. According to Stolzenburg et al. [1998], this thunderstorm charge structure is representative of the updraft region of either an isolated thunderstorm or of a convective element of a multicell thunderstorm (e.g., a mesoscale convective system).

Each charge region in Figure 3 is assumed to be a disk of uniform charge density and its size is represented by the colored rectangles in the figure. The total charge in each disk is also displayed. The electrostatic potential U_amb (the subscript “amb” stands for ambient), and its associated electrostatic field is also shown in Figure 3. The size, height, and magnitude of each charge region is consistent with work by Krehbiel et al. [2008]. The chosen charge magnitudes are such that produce maximum large-scale in-cloud electric fields of ∼50–100 kV/m, consistent with observations [e.g., Marshall et al., 1995; Stolzenburg et al., 2007]. Similar to the vertical electric field profiles obtained from balloon soundings [Marshall et al., 1995, Figures 2–6], we also show in Figure 3 the breakeven threshold for reference. The breakeven threshold can be linked to acceleration of relativistic (∼1 MeV, i.e., cosmic ray secondaries) electrons and is “the electric field strength at which the rate of energy gain from the field is equal to the rate of energy loss (mostly due to ionization of the air) for an electron moving directly along the field line” [e.g., Dwyer and Uman, 2014, section 3.7]. The breakeven threshold field is commonly associated with the maximum large-scale electric field observed in thunderstorms [e.g., Marshall et al., 1995]. The depicted charge structure forms three electric field regions, A/B/C, with upward/downward/upward directed ambient electric field, respectively. The direction of the ambient electric field defines the polarity of the electric field pulses observed on the ground, as discussed in section 7.

To the right of Figure 3 we present a summary of the different types/polarities of electric field pulses and their region of occurrence inside the thundercloud. Three-dimensional VHF mapping of lightning sources [e.g., Proctor, 1991; Rison et al., 1999; Thomas et al., 2001] combined with balloon electric field soundings [e.g., Coleman et al., 2003] have shown clearly that −CG lightning initiates just below the main negative charge regions, while +IC lightning starts in the gap between the main dipole charge centers in the storm. Recent studies have also shown that the source heights for IBPs in −CG and +IC lightning (approximately) coincide with the location of initial VHF sources as measured by a lightning mapping array [Karunarathne et al., 2013; Marshall et al., 2013]. Hence, for the thunderstorm charge configuration depicted in Figure 3, IBPs in −CGs and in +IC lightning occur in Regions A and B, respectively.

The NBE source heights have been estimated by measuring the direct NBE electromagnetic wave and its reflection in the ionosphere [Smith et al., 2004; Wu et al., 2012]. It has been found that −NBE nearly always occur above +NBEs in the same storm. Statistical analysis of a large pool of data lead previous investigators to draw the qualitative picture that +NBEs are generated between the main negative and positive charge regions, while −NBEs occur between the main positive charge region and the screening layer on top of the storm [Wu et al., 2012, Figure 7]. Following the nomenclature used in Figure 3 of this work, positive/negative NBEs occur in Regions B/C of the thunderstorm, respectively. Note that the actual source heights of charge layers vary from storm to storm depending on storm type, geographic location, etc. They may even vary within the same storm over time. Hence, the definition of Regions A, B, and C, in Figure 3, is only relative to the charge configuration adopted here. For instance, Region B is located between 7 and 10 km altitude. In contrast, Smith et al. [2004] found that most +NBE sources were located between 7 and 15 km altitude, thus suggesting that “Region B in their thunderstorms” extended to higher altitudes. Also note that in our schematics the truncation of Region C at 13 km altitude does not imply that −NBEs cannot occur at higher altitudes, since −NBEs have been observed to occur up to ∼20 km [Wu et al., 2012, Figure 6] and ∼30 km altitude [Smith et al., 2004, Figure 5].

In our model, the upper negative screening charge layer plays an important role in enhancing the electric field at the thunderstorm top, creating a favorable scenario for the generation of high-altitude −NBEs. There is a long history of studies on thunderstorm screening charge layers, as reviewed in detail by MacGorman and Rust [1998, pages 72–74]. Over 60 years ago, Grenet [1947] and Vonnegut [1953] independently suggested that screening layers would form at cloud tops and were important in thunderstorm electricity. Brown et al. [1971] were the first to model screening layer formation. Winn et al. [1978] provided the first direct measurements of screening layers. Many more measurements of screening charge layers were later performed by Stolzenburg et al. [1998], which contributed to the conceptual model of thunderstorm charge structure shown in Figure 3. For a detailed description of self-consistent modeling of thundercloud screening charge layers we guide the readers to the work of Riousset et al. [2010].

2.2 Conceptual Model for IBP and NBE Generation

In contrast to existing models of IBP and NBE sources (section 3), we propose that they are generated by a common physical mechanism and that they are linked to the initial development of the bidirectional lightning leader tree. We propose that IBPs and NBEs are the electromagnetic response to the stepping process of the initial lightning leader. In this paper we use the word “stepping” to describe the intermittent propagation of negative polarity tips of IC lightning leaders. Recent studies have advanced our understanding on how negative stepped leaders in CG lightning propagate [e.g., Rakov, 2013, Section 13]. For instance, high-speed camera recordings by Hill et al. [2011] (at 300,000 frames per second) showed that negative lightning leaders below the cloud base have average interstep intervals of ∼16 μs and step length of ∼5 m, as they approach the ground. Additionally, high-speed camera observations of stepped leaders in natural [Hill et al., 2011] and triggered [Biagi et al., 2010] lightning have showed that the stepping process in negative lightning leaders proceeds via the formation of a space leader ahead of the main leader channel. The space leader is analogous to the pilot systems observed in meter-scale negative polarity discharges in laboratory [e.g., Gorin et al., 1976; Ortega et al., 1994; Reess et al., 1995; Kochkin et al., 2015].

Similar investigations of IC lightning propagation mechanisms with high-speed cameras may not be feasible, since the cloud is opaque to optical radiation. For this reason, we have to rely on phenomenology inferred from studies at radio frequencies. Combining ground-based and balloon-borne electric field change measurements with VHF mapping of IC lightning leader locations, Winn et al. [2011] have provided a comprehensive description on the stepping process of IC lightning leaders. These authors have shown that the time interval between steps of IC lightning is ∼0.5–7 ms, while the estimated step length is ∼50–600 m. Note that these numbers are significantly larger than in negative leaders in CG lightning. According to Winn et al. [2011], to date there is no clear experimental evidence of whether space leaders are involved in the stepping process of IC lightning leaders. Therefore, in our model, we introduce a simplified representation of the complex sequence of processes that might be involved in the stepping of an IC negative leader. We simulate the stepping process by allowing the leader to “lurch forward” during a step [Winn et al., 2011, section 7].

Our proposed conceptual picture of how IBPs and NBEs are generated is schematically shown in Figure 4. The figure shows specifically how IBPs in +IC lightning and +NBEs are created in Region B of the thunderstorm, but a completely analogous process can be assumed to occur in Regions A and C to generate the other kinds/polarities of pulses. Consider the existence of an initial lightning leader of length ℓ₀ between heights h₁ and h₂, as shown in Figure 4c. Consider now the lengthening of its negative extremity through one step of length ℓ_step. The upper (negatively charged) tip moves from height h₂ to . Figures 4a and 4b show the electrostatic potential and linear charge density distributions, respectively, before and after the stepping process, as inferred from electrostatic models of the bidirectional lightning leader development [e.g., Mazur and Ruhnke, 1998; Pasko, 2014]. In Figure 4a the ambient electric field is assumed constant and directed downward; hence, the potential linearly increases with altitude, resembling the potential distribution in Region B of the storm (Figure 3).

The stepping process creates a new conducting segment ahead of the main leader on a time scale τ_step. Currents are induced in order to redistribute the charge on the leader toward the new steady state. The existing leader has a finite line conductivity G₀, expressed in units of S m. The line conductivity G₀ is simply the reciprocal of the channel's linear resistance, i.e., 1/G₀ has units of Ω/m. The line conductivity can be evaluated from detailed plasma physics models [e.g., Chemartin et al., 2009; da Silva and Pasko, 2013a] by integrating the leader conductivity, σ(r,t), over the channel's cross section, i.e., , where a is the channel radius.

Between the two electrostatic steady states shown in Figures 4a and 4b there is a transient current surge that radiates in the LF/VLF frequency range. The induced current will have the same direction of the ambient electric field, which is downward in Figure 3a, i.e., E_amb<0. The radiated electric field will have a positive polarity first peak, since the radiation component of the total electric field is (see equation 1). Therefore, the location h₀ of the existing leader inside the thundercloud determines the polarity of the observed electric field waveform. In the proposed framework, the characteristics of the emitted waveform depend on a minimal number of four parameters: ℓ₀ and G₀ (that describe the existing leader) and ℓ_step and τ_step (that describe the stepping process). In this paper we demonstrate that both IBP and NBE waveforms can be reproduced with realistic values for these four parameters (see likely range of values below). This represents an advance in comparison to existing TL and EAS–RREA models (section 3) that either have no physics-based explanation for the input parameters or require unrealistic ones.

We simulate the initial stages of the lightning tree development, where it can be assumed that the discharge is predominantly vertical and has a length scale, ℓ₀, between tens of meters up to ∼2 km. Beyond this length the leader network starts to branch horizontally inside the thunderstorm charge centers [e.g., Rison et al., 1999; Thomas et al., 2001; Behnke et al., 2005]. Winn et al. [2011, Figure 9] used electric field change data recorded at distances of ∼200 m from the lightning channels combined with VHF mapping to estimate leader step lengths of ∼50–600 m. Marshall et al. [2013, Table 2] used an array of fast antennas to estimate vertical displacements in IBP bursts of ∼200–1800 m. For this reason, we assume that ℓ_step can vary between tens of meters up to ∼1.5 km. Nonetheless, it should be kept in mind that IC lightning leaders steps longer than 600 m have not being unambiguously measured yet.

The stepping time scale, τ_step, is bounded by (1) the streamer-to-leader transition time scale and (2) the observed interstep interval during the in-cloud development of lightning flashes. The former is the physical time to convert a streamer corona to a hot and highly conducting leader channel and it is estimated, from a first-principles heating model [da Silva and Pasko, 2013a], to be between one and several microseconds. The latter is measured to be between ∼0.01 and 7 ms [e.g., Winn et al., 2011; Marshall et al., 2013]. It is likely that only a fraction of the interstep time interval is used to form the new leader segment. A considerable fraction of this interval may be used for building up of streamer zones and space leaders (if present). For this reason we estimate that τ_step should be between several to tens of microseconds. The leader's line conductivity can be estimated from experiments in long laboratory sparks [Frind, 1960; King, 1961]. The experimental results of King [1961] have been widely used in lightning literature to discuss the negative differential resistance of lightning leader channels [e.g., Heckman, 1992; Williams and Heckman, 2012; Mazur and Ruhnke, 2014]. From King's Volt-Ampere curve [e.g., Mazur and Ruhnke, 2014, Figure 4], one can estimate that between currents of ∼10–1000 A the leader (steady state) line conductivity lies between 0.01 and 1 S m. This range of possible values for G₀ is adopted in the present work. This range is also consistent with the values used in return stroke simulations. For instance, Rakov [1998] estimates that ahead of the return stroke wave the lightning channel line conductivity is 0.29 S m. However, in some works, values of up to 14 S m are adopted [Baba and Rakov, 2007, Table 1]. The likely range of values for the four input parameters of the model are listed in the bottom of Table 1.

Simulating the VHF radiation from the source of IBPs and NBEs is out of the scope of the present work. Nonetheless, in the proposed framework it can be clearly seen that charge redistribution leads to the induction of electric current/potential waves in the conducting leader. In regions where the potential difference ΔU = U − U_amb exceeds the threshold for streamer corona formation, streamer corona flashes will occur around the conducting leader, especially at its tips. This process will copiously emit VHF radiation. The strength and location of streamer flashes around the conducting leader are modulated by the electric current/potential waves propagating on it. Hence, the envelope of HF/VHF radiation should be directly related to the LF/VLF radiation waveform.

2.3 Numerical Model of Charges and Currents in a Developing Lightning Leader

Consider a lightning leader as a long (with length ℓ changing with time) and thin (with radius a ≪ ℓ) conductor formed inside the thunderstorm and, initially, situated at an altitude h₀, as schematically shown in Figure 4. The leader is a good (but not perfect) conductor with an effective line conductivity G. As the leader develops in the thunderstorm ambient electric field , charge and current are induced on the leader's surface (linear charge density q and current I are expressed in units of C/m and A, respectively). Charge redistribution on the leader takes place through surface electromagnetic waves, and the full dynamics of these waves can be obtained by the solution of the complete set of Maxwell's equations for the electric and magnetic fields. The properties of electromagnetic waves guided by a single-conductor transmission line were first investigated by Sommerfeld [1899] over a century ago. These “Sommerfeld” waves were investigated with purpose of engineering applications [Goubau, 1950] and were also used to simulate lightning dart leaders and return strokes [Borovsky, 1995]. Their dominant mode, as described in the aforementioned works, is a transverse magnetic one, with the only nonzero components being E_z, E_r, and B_φ, in cylindrical coordinates. The number of unknown variables can be further reduced by introducing the electric U and magnetic potentials. Consequently, we can describe the leader electrodynamics by a set of four equations and four unknowns:

(3)

(4)

(5)

(6)

Equations 3 and 4 define the retarded electric and magnetic potentials, respectively. The integration limits are from the lower h₁ to the upper h₂ leader tips (that may vary over time). The potentials at a location z and time t are generated by charges and currents at a retarded time , where is the distance between source and observation points and c is the speed of light in the medium. Equation 5 is obtained from the relationship between electric field and potentials (E_z=−∂U/∂z − ∂A/∂t) and by ensuring Ohm's law along the leader (I = GE_z). Equation 6 is the continuity equation. The system 3–6, although already complete, can be supplemented with the Lorenz gauge condition ∂U/∂t + c²∂A/∂z = 0.

In order to solve the system 3–6, we put forward a numerical method based on the time domain electric field integral equation [e.g., Miller et al., 1973; Davies and Duncan, 1994]. A similar approach has been used to simulate lightning return strokes [Moini et al., 2000; Liang et al., 2014] and to estimate lightning electromagnetic fields capable of generating terrestrial gamma ray flashes [Carlson et al., 2010]. We divide the leader in N equal segments of length Δz, as schematically shown in Figure 4d. Charges and currents are defined by piecewise constant functions on a staggered grid, where the position of charges z_q,i=(i − 1/2)Δz and currents z_I,i=iΔz are offset by Δz/2. The time step is defined as Δt = Δz/c, in such a way that the electric potential, at position z_q,i and time t^k, can be conveniently discretized as , where [K_i,j]⁻¹ is the capacitance (per unit length) matrix. Expressions for K_i,j are provided, for example, by Balanis [1989, page 674]. A similar expression can be defined for the magnetic potential: , where [L_i,j] is the inductance (per unit length) matrix. Equations 5–6 are discretized with an implicit first-order finite-difference scheme. Equation 5, for example, becomes . Finally, introducing the above-defined expressions for and , into equations 5–6, one can obtain a tridiagonal system of equations for currents that can be easily solved with Thomas algorithm [e.g., Hockney and Eastwood, 1988, page 185]. The proposed time-marching algorithm is more efficient than the original method proposed by Miller et al. [1973] in the sense that it requires the inversion of a tridiagonal (sparse) matrix at every time step instead of a fully populated (dense) matrix. The choice of Δt made here also makes the numerical expressions for the potentials more intuitive, since cumbersome temporal interpolations are not necessary. The implicit time-marching scheme is stable as long as Δz > 4a (see detailed numerical analysis by Davies and Duncan [1994]). In all simulated cases in this paper we adopt Δz = 5 m and a = 1 cm.

In equations 3–4 ε = 5.3ε₀ and μ = μ₀ [Moini et al., 2000]. The increase of ε with respect to its vacuum value is performed to account for the fact that the leader has a larger capacitance than a thin conductor of 1 cm radius. The linear charge density q(z) distributed along the leader channel generates an electric field component perpendicular to it that creates an streamer corona envelope. The corona sheath redistributes the leader charge around it, increasing its capacitance [Baum and Baker, 1990; Maslowski and Rakov, 2006]. Instead of including a conductor with larger radius in our simulations, we increase the leader's capacitance by increasing ε. It introduces a completely analogous effect (see explanation below) without compromising the stability of the numerical scheme. The total capacitance of our simulated leader is [Jackson, 2000]. Now, let us hypothetically replace the conductor of radius a embedded in the medium with permittivity ε, by another one with larger radius R_c embedded in the atmosphere, with permittivity ε₀. The capacitance of the new object is . Equalizing these two expressions for capacitance, and assuming ℓ≈1 km, one can find that the increase in permittivity (ε = 5.3ε₀), proposed by Moini et al. [2000] and adopted in the present work, produces the same capacitance as a leader with an effective corona sheath radius R_c≃90 m. The increase in the leader's capacitance also reduces the group velocity of waves propagating on it. The speed of light in the medium becomes  m/s, which better reflects observations [Moini et al., 2000].

In the present paper we specifically simulate one of the initial steps of the leader's negative extremity. On these short time scales, we can assume that the leader properties such as its capacitance matrix and line conductivity are constant over time. Moreover, we assume that the leader line conductivity is uniform and equal to a constant value G₀. The stepping process is simulated by raising the conductivity of the new leader segment from a very small number, 10⁻⁹G₀, to the leader steady state value, G₀, on a time scale τ_step. This process is parameterized with a Gaussian function for 0≤t≤τ_step, where .

The numerical model described in this section can be used to simulate the complete dynamics of the asymmetric bidirectional lightning leader tree. The model can be extended to account for three-dimensional features of the leader network, such as branching, for example [e.g., Lalande and Mazur, 2012]. It also allows for implementation of nonlinear effects in the resistance (the so-called negative differential resistance [e.g., Mazur and Ruhnke, 2014]) and capacitance (through modification of the corona sheath region [e.g., Maslowski and Rakov, 2006]). In the present work, we study the first elongation of the leader's negative tip as an isolated process and show how its dynamics produce the observed IBP and NBE electric field waveforms.

3.1 Dynamics of Charges and Currents Generating IBPs in CG Lightning

In this section we present the full dynamics of charges and currents during the first stepwise elongation of the bidirectional leader negative extremity during the initial stages of a CG lightning discharge. The radiated electric field has the waveform of a classic IBP [Nag et al., 2009]. We empirically choose the four input parameters to match IBP-4 recorded by Karunarathne et al. [2014] in 2011 near the NASA Kennedy Space Flight Center. The input parameters used are summarized in the first line of Table 1. The simulated multistation electric field waveforms are shown in Figure 1 in the same format as presented by Karunarathne et al. [2014, Figure 4]. The full dynamics of charges (a), potential (b), current (c), and line conductivity (d) are shown in Figure 5. At t = 0, one can see in the figure the existence of a short (ℓ₀=700 m) leader. The leader is evident by the equipotential region in Figure 5a, centered at h₀=6 km height, in Region A of the thunderstorm. The stepping process occurs in an analogous form to the schematics shown in Figure 4, but in the downward direction (the step length is ℓ_step=480 m). This is evident from Figure 5d, where the line conductivity in the newly formed leader segment, below the main leader, is raised to the leader's value, G₀=0.85 S m, on a time scale τ_step=28 μs. As conductivity is raised in the new leader portion, charge flows into it (Figure 5b). The current associated with this charge redistribution is shown in Figure 5c, where it can be seen to peak just below the connection point between the old leader and its new segment. Note that the only information prescribed in the model is the leader elongation through the raise in conductivity in the new section. The current pulse is not injected in the channel, as in TL models. It is actually a self-consistent response of the stepping process, as calculated from equations 3–6. Finally, we use the obtained current distribution, I(z,t), in equation 1 to obtain the multistation IBP waveforms shown in Figure 1.

The “bidirectional” aspect of the simulated lightning leader is an intrinsic property of the model as was first envisioned by H. W. Kasemir in 1950 [Kasemir, 2013]. At all times the leader is bidirectional, meaning it has simultaneously positive and negative polarity extremities going in opposite directions. In the results shown in this paper we do not follow the long leader evolution during later stages, as done, for example, by Mazur and Ruhnke [1998] and Pasko [2014] with purely electrostatic models. We only simulate the first ∼40 μs (see Figure 5) around the first negative leader step. On this short time scale the positive leader propagating (continuously) with a uniform speed of 10⁵ m/s only extends 4 m. This distance is negligible in comparison to other length scales in the model and, therefore, can be neglected.

Figure 6 shows a more detailed view on the leader's current distribution. Figures 6a–6h show snapshots of the simulated I(z,t) (solid line) and comparison to TL model results by Karunarathne et al. [2014] (dashed line). Figures 6i–6j present 2-D distributions of I(z,t) comparing our model results (i) to the previous work (j). The triangles in Figures 6i–6j mark the time instants of the snapshots in Figures 6a–6h. As discussed in section 3, Karunarathne et al. [2014] model IBP sources by injecting a current pulse I(h₁,t) at the bottom of a TL and obtain I(z,t) from equation 2. These authors have tried several dependences for f(z) (see equation 2) and have found that the best match for IBP-4 follows the Kumaraswamy distribution, f(x) = abx^{a − 1}(1 − x^a)^{b − 1}, where x = (z − h₁)/(h₂−h₁), and a and b are fitting parameters (the authors refer to this model as MTLK). The fitting parameters that Karunarathne et al. [2014] found to best mimic IBP-4 waveform are listed in their Table 2. Note that the TL fitting approach requires knowledge of around twice as many parameters. Despite the different nature of the two modeling approaches, the agreement between the two models, shown in Figure 6, is excellent; see, for instance, the vertical profiles near the time instant of peak current in Figure 6d. The dashed lines in Figures 6i and 6j illustrate the boundaries of the conductor in both models. Note the leader stepping in our model represented by a conductor elongation at t = 0 (in Figure 6i). The total length of the leader in present modeling, ℓ = ℓ₀+ℓ_step=1180 m, is comparable to fitting value 1016 m found by Karunarathne et al. [2014].

We found that the source peak current is 39.4 kA, while Karunarathne et al. [2014] estimate it to be 46.4 kA. Despite the difference, the calculated multistation E field amplitudes have the same value (compare our Figure 1 to their Figure 4). This happens because the electric field far away from the source depends not simply on the current but on the current moment. In the limit that , equation 1 can be simplified and becomes

(7)

where M_Q is the charge moment, M_I is the current moment, and ∂M_I/∂t is its derivative. In equation 7 these physical quantities are all evaluated at retarded time . The current moment is the parameter that can be unambiguously inferred from the remote sensing of the low-frequency electromagnetic fields, but not the current itself [e.g., Cummer, 2003]. From the model-generated I(z,t) we can calculate the charge moment, current moment, and its derivative as

(8)

(9)

(10)

Figure 7 shows the time dynamics of charge moment change, current moment, and its derivative. The charge moment change is defined as ΔM_Q(t) = M_Q(t) − M_Q(t = 0). The figure shows an excellent agreement between our results (solid line) and best fit TL model by Karunarathne et al. [2014] (dashed line). Equation 7 is useful because it clearly illustrates the 1/R dependence of the different electric field components. Far away from the source the third (radiation) term is dominant. Hence, . The current moment has a positive sign, i.e., upward. The first peak of the current moment derivative is positive (Figure 7c) because of the current raise from zero to its peak. Therefore, the polarity of the IBP pulse is negative. A summary of the different pulse polarities and its region of occurrence in the thundercloud is shown in Figure 3. It can be seen from equation 7 that for the radiation component has opposite sign with respect to that of electrostatic and induction ones, while for all three components have the same sign [e.g., Nag and Rakov, 2010, Figure A3]. The horizontal distance is equivalent to the well-known reversal distance for the electrostatic field due to an elevated finite-length vertical dipole [Ogawa and Brook, 1964; Rakov and Uman, 2003, page 71].

3.2 Effects of the Leader Step Characteristics on NBE Waveforms

In this section, we demonstrate that the proposed model can reproduce both IBP and NBE observations. Figure 2 shows simulated two-station E field waveform as observed by Eack [2004]. This (referred hereafter and in Table 1 as NBE-1) observation was used in a number of papers in existing literature for model validation purposes [Watson and Marshall, 2007; Nag and Rakov, 2010; Arabshahi et al., 2014]. A +NBE in our model is generated in Region B of the thunderstorm (see Figure 3) by the first upward step of the negative leader extremity, as schematically shown in Figure 4. The four leader characterizing parameters used in this simulation are listed in the second line of Table 1. For the sake of brevity, we do not show full dynamics of currents such as done in section 10. Nonetheless, a completely analogous process to the one illustrated in Figure 5 takes place in this case as well, with the only difference that the stepping is upward. Figure 8 shows the calculated source charge moment change, current moment, and its derivative. The figure shows good agreement between our simulations results and best transmission line model fit by Watson and Marshall [2007]. The input parameters used by Watson and Marshall [2007] are listed in their Table 2 and the resulting E field waveform is shown in their Figure 3b. We note that the total leader length in our simulation is 1630 m, which is 2–3 times larger than their estimate (of 630 m). Nonetheless, the similarity between ΔM_Q, M_I, and ∂M_I/∂t, calculated with the different models (Figure 8), leads to agreement in the simulated E field waveform. Nag and Rakov [2010, Table C1] also simulated NBE-1. They estimated that the TL model-based channel length is 650 m.

The four input parameters used in our model affect the characteristics of the simulated LF/VLF E field waveform, as shown in Figures 9 and 10. To fully demonstrate this effect we run a series of three cases, each one changing ℓ_step in the range between 50 and 1500 m. The remaining input parameters for each case are listed in Table 1. In all calculations below we assume that the observer is at a distance D = 100 km from the source and that the source is at h₀=8 km from the ground. The obtained waveforms for all simulations within Case 1 are shown in Figure 9. It is evident that larger steps produce higher peaks and longer total durations. To quantify these effects, we define in Figure 9b four waveform parameters: E₁ is the amplitude of first peak, E₂ is the amplitude of the opposite polarity overshoot, Δt₁ is the duration of the initial half cycle, and Δt₂ is the total pulse duration. The total duration is defined by the interval where E field modulus is greater than 10% of the peak-to-peak amplitude. These quantities are completely analogous to the data-derived ones shown in Figures 9 and 10 of Nag et al. [2010]. Figure 10 shows the dependence of E₁ (a), |E₁/E₂| (b), Δt₁ (c), and Δt₂ (d) on leader step length for all three simulated cases. Case 2 has τ_step twice as long as Case 1, while Case 3 has both τ_step and ℓ₀ twice as long as Case 1 (see Table 1). For comparison purposes, the shaded bands in Figure 10 show average ± one standard deviation values for 48 NBEs observed by Nag et al. [2010]. The best match within all simulated cases is marked with a circle, which is the waveform shown in Figure 9b.

In Figure 10a the increase of E₁ with ℓ_step can be clearly seen. This trend does not continue for all values because for larger steps the leader tends to escape from the high field region in the thunderstorm (see Figure 3). This concept can be understood by noticing that [Pasko, 2014]. The characteristic time scale for charge moment change is ∼τ_step. Hence, M_I∼ΔM_Q/τ_step and . According to equation 7 , and therefore, . This last expression qualitatively explain the numerically derived results in Figure 10a. It predicts the increase of E₁ with: ℓ_step (shown in all three cases), ℓ₀ (compare Case 3 with Case 2), and with 1/τ_step (compare Case 2 with Case 1). It also predicts that E₁ increases with the average ambient electric field in the region where the leader is developing, E_amb. This is evident from the numerical results by comparing long steps of Case 3 with Case 2. In the latter the upper tip (at height h₂) is located at a region of higher electric field than the former (see Figure 3). The longer model leader in Case 3 is compensated by the weaker model ambient electric field and the resulting E₁ amplitudes of both cases (for ℓ_step>1000 m) are comparable.

Figure 10 shows that amplitude E₁ (a) and duration Δt₂ (d) increase with ℓ_step. Waveforms recorded by electric field change sensors show that the amplitude and duration of IBPs tend to increase concurrently when comparing narrow IBPs with classic IBPs [Nag et al., 2009]. Our modeling results indicate that this trend may be also extended to NBE pulses as well. Therefore, Figure 10 illustrates a key conclusion from our work: that the difference between narrow IBPs, classic IBPs, and NBEs might be explained by the increase in a single parameter: the length of the initial negative leader step. Figure 10d shows that the total duration Δt₂ has similar increase with both τ_step (compare Case 2 with Case 1) and with ℓ₀ (compare Case 3 with Case 2). On the other hand, Figure 10c shows that the duration of first half cycle Δt₁ has stronger dependence on τ_step rather than ℓ₀. One of the criteria used by Wu et al. [2014] to automatically identify NBEs in their data is the abrupt rise in the waveforms' first peak (i.e., very short risetime, as shown in their Figure 2). From the four input parameters of our model, τ_step is the one that best reproduces this waveform property. A faster stepping rate (i.e., large 1/τ_step) leads to shorter risetimes.

4.1 On the Origin of the Initial Leader Segment

In the proposed model we assume the existence of an initial conducting segment of length ℓ₀, without explicitly defining time evolution or physical scenario of how it was created. There are currently two major hypothesis for the lightning initiation mechanism [e.g., Solomon et al., 2001; Rakov, 2013, section 2; Dwyer and Uman, 2014, section 3]. The first one, sometimes referred to as conventional breakdown mechanism, asserts that lightning initiates with streamer corona emissions from thundercloud hydrometeors [Griffiths and Phelps, 1976a, 1976b]. The second one, sometimes referred to as runaway breakdown, states that lightning initiation is facilitated by runaway electron avalanches seeded by cosmic ray secondary electrons [Gurevich et al., 1992, 1999]. It should be emphasized that the model proposed here does not favor or require any of these two mechanisms to explain lightning initiation.

Although we simulate only the electric field pulse associated to the first stepwise elongation of the negative leader extremity, the formation of the first leader segment should produce electromagnetic radiation itself. Recently, Marshall et al. [2014b] have identified slow electric field changes preceding the first IBP in both IC and CG lightning. These initial electric field changes (IECs) have average durations of 0.18 and 1.53 ms for CG and IC flashes, respectively. Hence, IECs are significantly slower than IBPs and NBEs, which have durations of ∼20 μs. Marshall et al. [2014b] also estimate that the peak current moment for IEC is ∼0.2 kA km, which is significantly smaller than the current moment of IBPs and NBEs (compare to Figures 7b and 8b). As discussed by Marshall et al. [2014b], the IEC's small amplitude and slow rate of change might have made these signatures overlooked in previous measurements.

In the framework of our model, IECs can be interpreted as the electromagnetic signature associated with the formation of the initial leader segment of length ℓ₀. Therefore, the measurements made by Marshall et al. [2014b] represent an experimental validation of our proposed mechanism for IBPs. We note that our unified perspective of how IBPs and NBEs are generated inside the thundercloud postulates that an E field change analogous to IECs in CG and IC lightning should precede NBEs. The existence of such E field change prior to the occurrence of an NBE awaits experimental verification. This change may be hindered depending on specifics of dynamics and duration of formation of the initial conducting segment.

Figure 2b shows a model-generated close range E field change signature associated with an NBE. The waveform looks substantially different from its far-field counterpart, shown in Figure 2a. The reason for this difference is that at close range the electrostatic and induction components of the total electric field are dominant over the radiation one. Another feature that can be seen in the close range waveform is the electrostatic field offset, of about −4.8 V/m, prior to the electric field pulse. This electrostatic offset is produced by the existing initial leader segment. Using only the first term on the right-hand side of equation 7, one can estimate that the creation of the initial leader segment corresponds to a charge moment change of about −76 C m. This amount is well within the range of IECs in IC flashes reported by Marshall et al. [2014b, Table 1].

4.2 Bouncing-Wave Mechanism of IBP/NBE Source Currents

The time scale for round trip of electromagnetic waves on the leader channel, simulated in section 10, is  μs, while peak current is established at ∼24 μs. Hence, wave reflections in the remote leader end (at height h₂=6.175 km) play a key role in establishing a standing wave on the channel. This phenomenon was referred to as bouncing-wave mechanism in the context of TL models [Nag and Rakov, 2010]. In the framework of present modeling, reflections are a natural response of the system given by the fact that the leader's conductivity terminates at its tips. Reflections are fully accounted for without the need to introduce artificial reflection coefficients. Nag and Rakov [2010] discuss that secondary peaks in NBE waveforms are evidence of current reflections at the source. From the perspective of our modeling framework the relationship between the source dynamics and secondary peaks in the radiated electric field is quite intuitive. Strong radiation from reflections will only occur if the peak current of incident wave (at the remote leader end) is comparable to the overall peak current during the whole stepping process. This is not the case of IBP-4 in Figure 1. One can see from snapshots in Figures 6a–6h that the current disturbance is generated at height ∼h₁=5.825 km and propagates both upward and downward. Due to the finite conductivity of the leader channel, the waveform arriving at height h₂ has significantly lower amplitude than at the source. Its reflection does not produce a peak in the radiated electric field (see Figure 1).

In the limit of an infinitely long TL, the electric and magnetic potentials can be simplified and rewritten as ΔU = q/C₀ and A = L₀I, respectively. The quantities and are the average linear capacitance and inductance, respectively [e.g., Bazelyan and Raizer, 2000, pages 39 and 176]. Using the above assumption, equations 5 and 6 can be simplified and become L₀∂I/∂t + ∂U/∂z + I/G₀=0 and C₀∂U/∂t + ∂I/∂z = 0, respectively. These expressions are the widely known (resistive) TL equations contained in engineering electromagnetics textbooks [e.g., Inan and Inan, 1998, page 25]. From a simple inspection of the above equations one can see that the solution has the dependence [e.g., Bazelyan and Raizer, 2000, page 176], meaning that waves propagating on a resistive transmission line attenuate on a time scale 2L₀G₀. For the parameters used in this work L₀≈2.3 μH/m. Comparing the two time scales defined in this section, one can see that a secondary peak, associated to the current wave reflection on the remote leader end, will occur if or  Ω. In other words, a secondary peak will occur if the leader resistance ℓ₀/G₀ is lower than the wave impedance . Assuming that the maximum value of G₀ is 1 S m, as discussed in section 7, one can find that only very short leader channels  m will generate IBPs and NBEs with secondary peaks, in agreement with findings by Nag and Rakov [2010].

4.3 NBEs, Conventional Lightning, Blue Jets, and Gigantic Jets

As reviewed in section 2, one of the most peculiar features of NBEs is that they are quite infrequently followed by conventional lightning [e.g., Smith et al., 1999]. However, some observations report that +NBEs appear to be the first event in a +IC lightning flash [e.g., Rison et al., 1999; Thomas et al., 2001]. The physical mechanism proposed in this paper states that NBEs, similarly to IBPs, are generated by (one of) the first elongation (or stepping) of the (upper) negative leader extremity in a +IC lightning, thus supporting the latter experimental results. At the current stage of theoretical development we can only speculate why most +NBEs reported in the literature are not followed by conventional lightning signals. A plausible scenario is that, due to a specific thundercloud charge configuration, after a large step (that generates an NBE pulse), the negative leader tip does not have significant potential difference with respect to the ambient potential and cannot ionize the air ahead of it to further propagate. This speculation is supported by recent experimental results that show that +NBEs that occurred as the initial event in +IC lightning occur mostly below 10 km altitude [Wu et al., 2014]. These results can be interpreted as indicative that +NBEs occurring at higher altitudes, i.e., further away from the high electric field region (Region B in Figure 3), cannot support further leader growth to become a +IC lightning.

Krehbiel et al. [2008] proposed an unified perspective on how the different kinds of lightning-related discharges, including blue jets [Wescott et al., 1995] and gigantic jets [Pasko et al., 2002; Su et al., 2003], are formed inside the thundercloud. They showed that gigantic jets are initiated in the same way as +IC lighting and owing to a charge imbalance in the thundercloud the upward leaders can penetrate through the upper positive charge region and escape the thundercloud to be observed in the stratosphere and mesosphere [Krehbiel et al., 2008; Riousset et al., 2010; da Silva and Pasko, 2012, 2013b]. Similarly, blue jets are discharges initiated between the upper positive charge region and screening charge above it (Region C in Figure 3) [Krehbiel et al., 2008; Edens, 2011]. Therefore, from the perspective of the mechanism proposed here, it is quite natural that +NBEs might also happen as the first event during the formation of a gigantic jet. Similarly, −NBEs might occur as the first discharge process in a blue jet. In context of gigantic jets, this idea has been experimentally proved to be correct: Lu et al. [2011a] have observed two gigantic jets that initiated with a +NBE signal. In context of blue jets, this theoretical prediction awaits experimental verification.

4.4 IBPs and Terrestrial Gamma Ray Flashes

Terrestrial gamma ray flashes (TGFs) are gamma ray bursts produced by thunderstorms and observed by low-Earth orbit satellites [Fishman et al., 1994]. TGFs were originally thought to be produced in association with sprites in the mesosphere [e.g., Taranenko and Roussel-Dupré, 1996; Inan et al., 1996], but later studies have shown that TGFs are generated within the first few milliseconds of a +IC lightning discharge development and have constrained their source altitudes to be inside thunderstorms [e.g., Stanley et al., 2006; Shao et al., 2010; Lu et al., 2010]. Additional studies point out that TGFs are coincident with a ULF (30–300 Hz) pulse with duration of 2–6 ms, which corresponds to charge moment changes of tens of C km and current moment of tens of kA km [Cummer et al., 2005; Lu et al., 2010, 2011b]. The TGF-associated lightning radio signals contain single- or multiple-VLF impulses (with durations <100 μs) superimposed on the slow ULF pulse [Lu et al., 2011b]. Østgaard et al. [2013] presented for the first time simultaneous detection of a TGF and optical lightning signals. The conclusion to be drawn from the aforementioned investigations is that TGFs are generated during the initial breakdown stage of +IC lightning [e.g., Marshall et al., 2013].

Pasko [2014] has employed an electrostatic model to simulate the development of the bidirectional IC lightning discharge [e.g., Mazur and Ruhnke, 1998] to demonstrate that the slow ULF pulse measured in association with TGFs may be actually due to the growth of the IC lightning leader itself. Marshall et al. [2013] combined electric field change [e.g., Karunarathne et al., 2013] sensors with azimuthal magnetic field ones (same as that used by Lu et al. [2011b]) to demonstrate that the VLF pulses measured synchronously with TGFs are indeed IBPs (or IBP trains). A similar idea has also been discussed by Shao et al. [2010]. The modeling results presented in our paper show that (positive) IBPs are related to the individual steps of the negative leader tip of the +IC lightning across the gap between the main negative and positive charge regions in a thunderstorm (Region B in Figure 3). Hence, our results corroborate to the idea that both (slow) ULF and (fast) VLF radio signals measured in association with TGFs are byproducts of the IC lightning leader growth itself. We note that TGFs have also been found to occur closely in time with NBEs [Stanley et al., 2006; Shao et al., 2010; Lu et al., 2011b].