plasma-expert/reference/sources/plasma-nature-lightning-cha...


								=== PAGE 2 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								collision rates between charged and neutral particles, and low degree of ionization (see, e.g., Becker et al.,

								2004). Streamers, for instance, only partially meet the three standard conditions that traditionally define a

								plasma (Bittencourt, 2004; pp. 6–11). These criteria define the plasma's ability to shield short-range electro-

								static interactions between individual particles, remain quasi-neutral, and respond collectively to long-range

								electromagnetic forces. The three conditions can be estimated for typical streamer properties at atmospheric

								pressure, that is, electron temperature of 23,000 K, or ∼2 eV, and electron density of 1018–1020 m−3 (Raizer,

								1991; section 12.3). First, the Debye length is ∼1–10 μm, which is relatively smaller than the streamer radius,

								∼0.1–1 mm (Naidis, 2009). Second, there are many electrons in a Debye sphere, ∼500–5,000. Third, the

								electron-neutral collision frequency is ∼1012 s−1, which is higher than the frequency of relevant processes,

								including the plasma frequency. Therefore, it can be said that the first two conditions are approximately

								met, but not the third one. On the other hand, it is easy to show that all three conditions are met in the

								return stroke channel. Therefore, even though the formal definition of a plasma is not always met within the

								many elements of a lightning flash, we refer to its constituting ionized gas as a “plasma,” because it remains

								quasi-neutral and responds collectively to applied electric fields.

								The aforementioned collective behavior in lightning is evidenced in the many types of ionization waves

								(e.g., streamer front, leader front, dart leaders, and return strokes), its ability to shield itself from exter-

								nally applied electric fields, and its negative differential resistance, which in its turn map into several

								phenomenological features, including its fractal structure, the contrasting behavior of positively and neg-

								atively charged extremities, and the fact that leader channels are enveloped by streamer zones and corona

								sheaths. This manuscript focus on perhaps the most important feature attributed to the plasma nature of

								lightning—its nonlinear resistance. A correct description of the channel resistance is required to better

								characterize lightning electromagnetic emissions, to correctly predict its deleterious effects in man-made

								structures, to quantify the impacts of lightning in atmospheric chemistry, and to address fundamental open

								questions regarding lightning initiation, propagation, and polarity asymmetries. The nonlinear plasma resis-

								tance is in its turn dependent on the history of energy deposition and losses in the channel and cannot be

								accurately determined without properly tracking the evolution of all other channel properties, including

								electric field, electron density, temperature, and radius.

								Efforts to characterize the nonlinear resistance and overall plasma properties of the lightning channel can

								be classified into three categories: (1) LTE gas-dynamic models (Aleksandrov et al., 2000; Chemartin et al.,

								2009; Hill, 1971; Paxton et al., 1986; Plooster, 1971; Ripoll et al., 2014a), (2) streamer-to-leader transition

								models (Aleksandrov et al., 2001; Bazelyan et al., 2007; da Silva & Pasko, 2013; da Silva, 2015; Gallimberti,

								1979; Gallimberti et al., 2002; Popov, 2003; 2009), and (3) semiempirical resistance models (Baker, 1990;

								De Conti et al., 2008; Koshak et al., 2015; Mattos & Christopoulos, 1990; Theethayi & Cooray, 2005). The

								three categories are described in the upcoming paragraphs. Instructive discussions and additional references

								regarding each of the three categories can also be found in sections 2.5, 2.3, and 4.4, respectively, of Bazelyan

								and Raizer's (2000) textbook. On a separate note, the literature concerning the resistance of short spark

								discharges in the laboratory is very rich and has provided many insights into building the models cited above

								(see, e.g., Engel et al., 1989; Kushner et al., 1985; Marode et al., 1979; Naidis, 1999; Riousset et al., 2010;

								Takaki & Akiyama, 2001). It is outside of our scope to provide a detailed review of these investigations, but

								it can easily be found elsewhere (da Silva & Pasko, 2013; Engel et al., 1989; Montano et al., 2006).

								The first group of investigations evaluates the resistance of a lightning channel under the assumption that

								the plasma is in LTE. In this framework, the electrical conductivity is only a function of temperature, that

								is, 𝜎= 𝜎(T), which is valid for atmospheric-pressure arcs at temperatures higher than ∼10,000 K, where T

								or simply the word “temperature” here and in the remainder of this manuscript corresponds to the tem-

								perature of the neutral gas. (The 10,000-K threshold is a rough estimate; see section 2.2 for justifications.)

								Following the return stroke simulations performed by Plooster (1971), these models describe how Joule

								heating deposition in the channel core heats the air and causes rapid hydrodynamic expansion. They solve

								a system of three equations accounting for conservation of mass, momentum, and energy (or enthalpy) of

								the neutral gas (air). They are often solved in a 1-D radial domain, with the exception being the work of

								Chemartin et al. (2009) where efforts are made to capture the 3-D tortuosity of a plasma arc. A few of these

								models also present a detailed description of the plasma radiative transfer (see, e.g., Paxton et al., 1986;

								Ripoll et al., 2014a).

								DA SILVA ET AL.

								9443


								=== PAGE 3 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								The second class is dedicated to a detailed description of the streamer-to-leader transition process, which

								takes place during the discharge onset or at the tip of a growing channel. Streamer-to-leader transition is the

								name given to the sequence of processes converting cold and low-conductivity plasma channels (streamers)

								into hot and highly conducting ones (leaders), a condition required to allow lightning channels to propagate

								for several kilometers in the atmosphere before decaying (Bazelyan & Raizer, 2000, p. 59). These models

								account for the hydrodynamic expansion of the neutral gas, such as the ones described in the first category.

								However, following in the footsteps of the seminal monograph by Gallimberti (1979), they also account for

								a non-LTE plasma conductivity arising from the detailed kinetic balance of an air plasma. The more recent

								models describe in detail the energy exchange between charged and neutral particles accounting for the

								partitioning of electronic power between elastic collisions, and excitation of vibrational and electronic states,

								and also delayed vibrational energy relaxation of nitrogen molecules (see, e.g., da Silva & Pasko, 2013). The

								non-LTE conductivity regime encompasses temperatures lower than ∼10,000 K. The models in this category

								(cited in this paragraph) do not account for photoionization, which is important at the high temperatures

								present in the return stroke channel.

								The third category groups investigations where a semiempirical expression for the channel resistance (per

								unit length) as a function of time, R(t), has been employed in return stroke simulations. The reasoning

								behind such approach is that it is impractical to use the self-consistent gas-dynamic simulations to calculate

								the resistance of a channel that is 5 (or more) orders of magnitude longer than wider. Therefore, a paramet-

								ric dependence for R(t) facilitates the implementation of a height-dependent, transmission-line-like return

								stroke model. These investigations use expressions for R(t) derived by Barannik et al. (1975), Kushner et al.

								(1985), and others, as reviewed by De Conti et al. (2008). To the best of our knowledge, only Liang et al.

								(2014) present an effort to couple a self-consistent resistance calculation with a transmission-line-like return

								stroke model. These authors use a two-temperature plasma model to infer the electronic conductivity. The

								model does not account for channel expansion or plasma chemistry, and it is unclear how well it compares

								to the conventional gas-dynamic return stroke simulations. Nonetheless, investigations such as done by De

								Conti et al. (2008) and Liang et al. (2014) raise the need for accurate and computationally efficient models

								for R(t).

								The objective of this work is to fill a gap in the peer-reviewed literature by introducing a comprehensive—yet

								simple—model that can exemplify the plasma nature of lightning channels (section 2.1). We describe

								a series of parameterizations that allow us to capture both the low-temperature/non-LTE and the

								high-temperature/LTE regimes, account for radial expansion, and include negative-ion chemistry, at little

								computational cost (section 2.2). The model is first tested by calculating the time scale for streamer-to-leader

								transition (section 3.1), it is then validated against experimental data on the steady-state negative differential

								resistance of plasma arcs (section 3.2), and finally, compared to well-established gas-dynamic return stroke

								simulations (section 3.3). As an application of the model, we simulate optical emissions of rocket-triggered

								lightning and compare to the experimental findings of Quick and Krider (2017) (section 3.4).

								2. Model Formulation

								2.1. Basic Equations

								In this work we describe the minimal model to qualitatively capture the consequences of the plasma nature

								of lightning channels. The key simplification here is to solve a set of zero-dimensional equations (i.e., with

								zero spatial dimensions) that describe the temporal dynamics of the plasma in a given cross section of the

								channel. Starting from a general 3-D problem, we can progressively reduce the dimensionality of the system.

								A schematical representation of the model is given in Figure 1a. It can be assumed that the lightning channel

								is a long cylinder. The axial symmetry indicates that the plasma conditions do not depend on the polar

								coordinate. Furthermore, the 2-D long cylinder geometry can be reduced to a 1-D radial one, by noting that

								variations along the channel have significantly larger length scales than along the radial direction. Thus,

								the change in plasma properties are driven by the conduction current created by the overall lightning tree

								dynamics and merely imposed in that channel section. Finally, the 1-D radial dynamics can be averaged

								over to produce self-similar solutions of average channel properties. The minimal set of equations can be

								written as follows:

								E = RI =

								I

								𝜎𝜋r2

								c

								=

								I

								e𝜇ene𝜋r2

								c

								(1)

								DA SILVA ET AL.

								9444


								=== PAGE 4 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 1. (a) Schematical representation of how the model simulates a cross sectional area of the lightning channel,

								provided only the current passing through that region I(t) and the channel initial conditions. (b) Current waveforms

								adopted in this study: constant current versus four-parameter pulsed profile. (c) Radial temperature profile and

								corresponding channel expansion. Lightning leader channels are surrounded by streamer zones and corona sheaths,

								which are not depicted in panel (a).

								𝜌mcp

								dT

								dt = 𝜂T𝜎E2 −4𝜅T

								r2

								g

								(T −Tamb

								) −4𝜋𝜖

								(2)

								dne

								dt = (𝜈i −𝜈a2 −𝜈a3

								) ne + 𝜈dnn + kepn2

								LTE −kepne

								(ne + nn

								)

								(3)

								dnn

								dt =

								(

								𝜈a2 + 𝜈a3

								)

								ne −𝜈dnn −knpnn

								(

								ne + nn

								)

								(4)

								dr2

								c

								dt = 4Da

								(5)

								dr2

								g

								dt = 4𝜅T

								𝜌mcp

								(6)

								Equation (1) is the Ohm's law applied to the channel's cross section, which relates the axial electric field E

								to the electrical current I, via the resistance per unit channel length R = 1∕𝜎𝜋r2

								c, where 𝜎is the electrical

								conductivity and rc is the plasma channel or current-carrying radius. (For the remainder of this manuscript,

								we refer to the resistance per unit channel length R as simply the resistance.) The electrical conductivity

								is given by 𝜎= e𝜇ene under the assumption that only the electron contribution is important, where e is the

								electronic charge, 𝜇e is the electron mobility, and ne is the electron density. This is a reasonable approxima-

								tion because the ion mobility is of the order of 10−4 m2·V−1·s−1 (at 1 atm), while the electron mobility is 2–4

								orders of magnitude larger in the range of typical electric fields present in electrical discharges (see, e.g.,

								Figure 3a).

								Equation (2) describes the rate of change of air temperature T, where 𝜌m is the air mass density and cp is

								the specific heat at constant pressure. The first term on the right-hand side is the rate of Joule heating of air,

								where 𝜂T ≃10% is the fraction of electron Joule heating power contributing to air heating. The second term

								represents cooling due to heat conduction, where rg is the thermal radius (delimiting the hot air region), 𝜅T

								is the thermal conductivity, and Tamb = 300 K is the ambient air temperature. The third term corresponds to

								DA SILVA ET AL.

								9445


								=== PAGE 5 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 2. The solid lines show the local-thermodynamic equilibrium (LTE) properties of air as a function of temperature used in the present paper. (a) Mass

								density 𝜌m, (b) specific heat at constant pressure cp, (c) the product 𝜌mcp, (d) thermal conductivity 𝜅T, (e) electrical conductivity 𝜎LTE, and (f) net emission

								coefficient 𝜖for an optically thin plasma. The red dashed line in panels (a) and (c) show the ideal gas law trend 𝜌m ∝1∕T. In the original references, data are

								only available for temperatures to the left of the vertical dash-dotted line. For higher temperatures, we perform an analytical extrapolation using the data in the

								range highlighted in green. The air-plasma properties shown in the figure are taken from Boulos et al. (1994, pp. 413–417), unless otherwise noted.

								energy loss due to radiative emission, where 𝜖is the net radiation emission coefficient. Equation (2) assumes

								isobaric air heating and neglects cooling by convection.

								Equation (3) describes the change in electron density ne. The first term on the right-hand side describes

								the rate of change due to field-induced, electron-impact processes, where 𝜈i, 𝜈a2, and 𝜈a3 are the ionization,

								two-, and three-body attachment frequencies, respectively. The second term describes electron detachment

								from negative ions, where 𝜈d is the detachment frequency and nn is the negative-ion density. The third term

								describes the effective rate of thermal ionization, where kep is the rate coefficient for electron-positive ion

								recombination, and nLTE is the electron density in local thermodynamical equilibrium (LTE), defined as

								nLTE = 𝜎LTE∕e𝜇e. The LTE conductivity 𝜎LTE is only a function of temperature (see, e.g., Figure 2e). The fourth

								DA SILVA ET AL.

								9446


								=== PAGE 6 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								term represents plasma decay due to electron-positive ion recombination. Charge neutrality is assumed;

								thus, the positive-ion density is equal to ne + nn.

								Equation (4) describes the evolution of an effective or generic negative-ion density nn. This quantity rep-

								resents O−and O−

								2 , the dominant negative ions in atmospheric discharges. These species are created by

								two- (𝜈a2) three-body attachment (𝜈a3), respectively. The last term in equation (4) represents a sink of

								negative ions due to negative-positive ion recombination, where knp is the corresponding rate coefficient.

								Equations (5) and (6) describe the rate of expansion of the current-carrying radius rc and of the thermal

								radius rg, respectively, where Da is the ambipolar diffusion coefficient. For all purposes, rc represents the dis-

								charge channel radius, because it enters in the calculation of Joule heating power deposited in the channel

								via equation (1). The parameter rg is best interpreted as a measure of the curvature of the radial temper-

								ature profile, and its only contribution in the system of equations is in the thermal conduction cooling in

								equation (2).

								The set of six equations (1)–(6) is solved to obtain the temporal dynamics of six unknowns E, T, ne, nn, rg,

								and rc, respectively. The input parameters are the source current dynamics I(t) and the initial conditions for

								the five state variables (T, ne, nn, rg, and rc), as shown in Figure 1a. The initial value of the electric field is

								given directly from equation (1).

								In order to solve equations (1)–(6), several coefficients are required. These coefficients are a function of E∕𝛿,

								T, or both. The quantity E∕𝛿is the so-called reduced electric field, where 𝛿is the reduction of air density

								in comparison to the sea level, room temperature value, defined precisely as 𝛿= 𝜌m(h, T)∕𝜌m(h = 0km, T =

								300 K); h here corresponds to the altitude above mean sea level. Figure 2 shows all LTE plasma coefficients

								used: (a) 𝜌m, (b) cp, (c) 𝜌mcp, (d) 𝜅T, (e) 𝜎LTE, and (f) 𝜖. The LTE parameters are, by definition, only function

								of temperature. Note that the assumption of isobaric heating combined with the ideal gas law would lead to

								a dependence 𝜌m ∝1∕T between mass density and temperature. This trend is shown in Figures 2a and 2c

								as a red dashed line. However, in the present work, we use the full equilibrium calculations given by Boulos

								et al. (1994), shown as blue solid lines in the figure.

								Figure 3 shows the field-dependent coefficients: (a) 𝜇e, (b) effective frequencies of electron production and

								loss processes, (c, d) recombination coefficients, and (e, f) Da. The conventional breakdown threshold is

								defined by the equality between electron-impact ionization (𝜈i) and two-body attachment (𝜈a2) in Figure 3b.

								For the coefficients used here its numerical value is Ek∕𝛿= 28.4 kV/cm. Figures 3c and 3d show both

								the electron-positive ion (kep) and negative-positive ion (knp) recombination coefficients, as a function of

								the reduced electric field and temperature, respectively. Similarly, Figures 3e and 3f show the ambipolar

								diffusion as a function of electric field and temperature, respectively.

								The coefficients have been obtained from the following references: 𝜌m, cp, and 𝜅T (Boulos et al., 1994); 𝜎LTE

								(Boulos et al., 1994; Yos, 1963); 𝜖(Naghizadeh-Kashani et al., 2002); 𝜇e (Cho & Rycroft, 1998); 𝜈i and 𝜈a2

								(Benilov & Naidis, 2003); 𝜈a3 (Morrow & Lowke, 1997); 𝜈d (Luque & Gordillo-Vázquez, 2012); kep and knp

								(Kossyi et al., 1992); and Da is defined by the Einstein relation (Raizer, 1991, p. 20). Both kep and Da effectively

								depend on the electron temperature Te. The expression for Te(E∕𝛿, T) is taken from Vidal et al. (2002). The

								rate coefficients are given for an air composition of 80% N2 and 20% O2. All rate coefficients used in this

								manuscript have been summarized in the form of two Matlab functions and made publicly available online

								(da Silva, 2019a).

								2.2. Key Assumptions

								1. Externally driven electrical current. A key assumption of the model is that the electrical current is gener-

								ated by the overall lightning discharge electrodynamics and merely imposed to the channel cross section

								of interest. This allows one to calculate the channel properties for a given constant or pulsed current

								waveform. Here we use two types of waveforms: a constant current (in sections 3.1 and 3.2) and a

								four-parameter pulsed current waveform (in sections 3.3 and 3.4). The pulsed current waveform quali-

								tatively captures most impulsive processes taking place in the lightning channel, and it is given by the

								following mathematical expression:

								I(t) =

								{ Ip t∕𝜏r

								if

								t ≤𝜏r

								(Ip −Icc) exp(−t∕𝜏f) + Icc if

								t > 𝜏r

								(7)

								DA SILVA ET AL.

								9447


								=== PAGE 7 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 3. Electric-field-dependent coefficients used in this investigation. (a) Electron mobility 𝜇e. (b) Effective frequencies of electron production and loss

								processes 𝜈i, 𝜈a2, 𝜈a3, and 𝜈d, from equation (3). (c, d) Recombination coefficients kep and knp. (e, f) Ambipolar diffusion coefficient Da. Panel (c) shows the

								recombination coefficients as a function of E∕𝛿for two different temperature values. Contrastingly, panel (d) shows the same coefficients as a function of T for

								two values of E∕𝛿. The same strategy is used to display Da in panels (e) and (f). Panel (d) also shows the rate coefficient for three-body electron-positive-ion

								recombination (electrons are the third body), or more precisely kep3ne, with ne = 1020 m−3. This process is not included in the model, and the coefficient is just

								shown for comparison with the two-body rate. Expressions for the rate coefficients shown in this figure are given by da Silva and Pasko (2013); see text for

								references.

								The four parameters in the waveform are peak current Ip, rise time 𝜏r, fall time 𝜏f, and continuing current

								Icc. These four parameters can be adjusted to represent a first or subsequent return stroke with or without

								continuing current. They can also be adjusted to allow the model to simulate the surge current injected

								in the leader channel following the stepping process (see, e.g., Winn et al., 2011), a dart leader reioniza-

								tion wave, or ICC pulses happening during the initial continuous current (ICC) stage of a rocket-triggered

								lightning flash. A schematical representation of this waveform is given in Figure 1b. It should be noted

								that several different analytical functions have been used to simulate the current waveform propagating

								through the lightning channel, such as the Heidler function (Heidler, 1985; Rakov & Uman, 1998), the dou-

								ble exponential (Bruce & Golde, 1941), or the asymmetric Gaussian (e.g., da Silva et al., 2016). The model

								DA SILVA ET AL.

								9448


								=== PAGE 8 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								can handle any of them as input; equation (7) is chosen for its simplicity and to facilitate the comparison

								with the work of Plooster (1971) and Paxton et al. (1986) in section 3.3 below.

								The overall strategy of prescribing I(t) and calculating the channel properties has been success-

								fully employed by a number of researchers to investigate the dynamics of streamer-to-leader and

								streamer-to-spark transition (Aleksandrov et al., 2001; da Silva & Pasko, 2012; Gallimberti et al., 2002;

								Popov, 2003) and to simulate the channel decay following a return stroke (Aleksandrov et al., 2000; Hill,

								1971; Paxton et al., 1986; Plooster, 1971). Although insightful, this strategy does not reveal the full lightning

								electrodynamics, because changes in the plasma conductivity should feedback into how much current is

								flowing in the channel. However, the approach used here allows us to provide a detailed characterization

								of the plasma-channel nonlinear resistance R(t) for a given current I(t). This manuscript should be seen

								as an initial effort toward quantifying the effects of the nonlinear plasma resistance into the overall elec-

								trodynamics of lightning leaders. Future investigations can leverage this model by replacing equation (1)

								with lumped or distributed circuit equations that describe the lightning discharge tree.

								2. Averaged radial dynamics. The radial profile of temperature is assumed to follow a step function so that

								T(r) = T for r ≤rg and T(r) = Tamb for r > rg. The radial expansion is given by an increase of rg at a rate

								given by equation (6). It is assumed here that the expansion rate is determined by thermal conduction

								or, in other words, the radial temperature profile follows the equation 𝜕T∕𝜕t = k∇2T, where k = 𝜅T∕𝜌mcp.

								The solution for this equation under a delta function initial condition is T(r, t) = exp(−r2∕4kt)∕

								√

								4𝜋kt. The

								solution is a Gaussian function with half-width rg =

								√

								4kt. Taking the time derivative of this expression,

								one obtains the expansion rate of the thermal radius in equation (6).

								The second term in the right-hand side (rhs) of equation (2) is the spatially averaged Laplacian of temper-

								ature, that is, the rhs of the heat conduction equation. The method for evaluating that term is illustrated

								in Figure 1c. It is assumed that the thermal conduction-driven expansion conserves the area under the

								curve in Figure 1c, or the quantity A = (T −Tamb)𝜋r2

								g. Therefore, 𝜕T∕𝜕t|thermal

								conduction is determined from set-

								ting 𝜕A∕𝜕t = 0. This is a rather robust assumption since it is virtually equivalent to enforcing energy

								conservation. However, in reality, the shape of the profile is not preserved as assumed here.

								Similar results are obtained by assuming that the plasma distribution expands with ambipolar diffusion,

								leading to the expansion rate given in equation (5). In this case, the conserved quantity is A = ne𝜋r2

								c,

								or simply the number of electrons per unit channel length. Conservation of A in this case is equivalent

								to conservation of mass. This analysis also yields a radially averaged ambipolar diffusion sink term in

								equation (3). However, this loss process is negligible in comparison to chemically driven losses and, there-

								fore, it is not included in equation (3). Our considerations here are similar to Braginskii's (1958), where

								the plasma channel boundary is assumed to behave as a moving piston that “snowplows” the ambient gas.

								Both models yield a channel radius expansion as rc ∝

								√

								t, but Braginskii's expansion rate is not deter-

								mined by ambipolar diffusion. In a comparison between several semiempirical models of the lightning

								return stroke resistance, De Conti et al. (2008) concluded that the model accounting for channel expan-

								sion rc ∝

								√

								t effects in the resistance yielded the most robust return stroke radiated electromagnetic field

								signatures.

								3. Thermal ionization rate. At temperatures of several thousand Kelvin, the plasma-channel composition is

								roughly made of equal parts electrons and NO+ ions (Aleksandrov et al., 1997; da Silva & Pasko, 2013;

								Popov, 2003). The NO+ ions are formed by associative ionization of N and O atoms at a rate F = kassocnOnN.

								The plasma density is dictated by a balance between associative ionization and electron-positive ion recom-

								bination, that is, by F = kepnenNO+ ≈kepn2

								e. Without knowing the precise rate F, we know that at high

								temperatures this equation should yield the LTE conductivity given in Figure 2e, or the corresponding elec-

								tron density nLTE = 𝜎LTE∕e𝜇e. This can be achieved by setting the rate of thermal (associative) ionization

								to be equal to F = kepn2

								LTE, as done in equation (3).

								Therefore, equation (3) is designed to essentially have two different modes of operation. At low

								(near-ambient) temperatures, the plasma population balance is driven by electron-impact ionization,

								attachment, and detachment, that is, the typical chemistry considered in the streamer breakdown of short

								air gaps (da Silva & Pasko, 2013; Flitti & Pancheshnyi, 2009; Liu & Pasko, 2004; Naidis, 2005; Pancheshnyi

								et al., 2005). However, at high temperatures (≳10,000 K) the equation yields the LTE conductivity 𝜎LTE(T),

								in alignment with the typical approach used for the simulation of free-burning arcs (Chemartin et al.,

								2009; Lowke et al., 1992) or used in gas-dynamic return stroke simulations (Aleksandrov et al., 2000;

								Paxton et al., 1986; Plooster, 1971). It is not possible to state exactly what is the minimum temperature at

								which the assumption of LTE regime yields accurate calculations. Both T and Te depend on the history of

								DA SILVA ET AL.

								9449


								=== PAGE 9 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								energy deposition and losses in the channel, which in its turn depend on the electric field and the elec-

								tron density. In this manuscript, we loosely give the value of 10,000 K as an estimate. This is the value at

								which the electron temperature is only 5% and 50% larger than the neutral gas one for electric fields of 10

								and 1000 V/m, respectively. In the present work, the electron temperature is obtained under the assump-

								tion that the electron energy balance equation is in steady state. Therefore, yielding the simple relation

								Te = T + f(E∕𝛿), where the function f(E∕𝛿) ∝(E∕𝛿)0.46 is taken from Vidal et al. (2002). Essentially, this

								equation asserts that the non-equilibrium results from the presence of an electric field in the discharge

								plasma and that equilibrium is only achieved when E = 0.

								In some types of plasmas the high-temperature density is given by a balance between electron impact ion-

								ization (driven by high T and not high E∕𝛿) and three-body electron-positive ion recombination (electrons

								are the third body). One such example are Argon arc discharges at atmospheric pressure (see, e.g., Sanson-

								nens et al., 2000; Tanaka et al., 2003). In this case, plasma losses would happen at a rate ≈kep3n3

								e, and using

								the assumptions discussed in the last two paragraphs, the plasma production rate would be ≈kep3n3

								LTE,

								where kep3 is the three-body electron-positive ion recombination rate coefficient given in units of m6/s.

								Owing to the cubic power law dependence, three-body electron-positive ion recombination is important

								when the plasma density is high. In this work, we assume that the high-temperature balance is given by the

								two-body processes, because they are the dominant ones in the temperature range between 2,000–9,000 K

								(i.e., in the transition to LTE regime), as discussed by Bazelyan and Raizer (2000, pp. 75–80) and Aleksan-

								drov et al. (2001). To verify that this assumption is true, we first plot the rate coefficients kep ∝T−1.5

								e

								and

								kep3 ∝T−4.5

								e

								in Figure 3d with rate coefficients taken from Kossyi et al. (1992) for an air plasma. Figure 3d

								actually shows kep3ne so that the units match, with ne = 1020 m−3, a typically large value in our simula-

								tions. It can be seen that due to the weaker fall off with temperature, two-body recombination increasingly

								dominates over three-body in the temperature range of interest. Second, we show later in section 3.3 quan-

								titative comparisons between the two rates for specific simulation results obtained with our model, further

								justifying our use of two-body process rates.

								4. Negative ion chemistry. Equation (4) describes the evolution of an effective or generic negative-ion density

								nn, representing O−(created by two-body attachment) and O−

								2 (created by three-body attachment), the

								dominant negative ions in ambient-temperature discharges. In the hot lightning channel, negative ions

								disappear, and the plasma composition is given by a balance of positive ions and electrons. By comparing

								equations (3) and (4), we can see that attachment works as a sink in the former, but as a source in the lat-

								ter. Detachment plays the opposite role. Therefore, the attachment-detachment cycle does not represent

								a true plasma loss. Effectively, electrons can be thought to be temporarily stored in negative ions to be

								released at a later time, after substantial accumulation. It is assumed here that O−

								2 created by three-body

								attachment quickly converts into O−in collisions with atomic oxygen favored by elevated temperatures

								in the lightning channel (da Silva & Pasko, 2013, Figure 11a). Therefore, detachment is dominantly

								driven by collisions between O−and N2 (Luque & Gordillo-Vázquez, 2012; Rayment & Moruzzi, 1978).

								These assumptions allow us to account for effects of negative-ion chemistry in a simple yet reasonably

								accurate manner.

								5. Fast air heating. The coefficient 𝜂T in the first term on the rhs of equation (2) is the fraction of elec-

								tronic power (or Joule heating rate 𝜎E2) that is directly transferred into random translational kinetic

								energy of neutrals and, thus, contributes to air heating. This quantity has been calculated to be 𝜂T ≃0.1

								at near-ambient temperatures (da Silva & Pasko, 2013; da Silva, 2015), largely arising from surplus energy

								from the quenching of excited electronic states and molecular (electron-impact) dissociation, which

								consist the so-called fast air heating mechanism (Popov, 2001, 2011; da Silva & Pasko, 2014).

								Most of the remainder electronic power is spent into the excitation of vibrational energy levels of nitrogen

								molecules. However, as temperature increases, rates of vibrational-translational energy relaxation quickly

								accelerate, effectively making 𝜂T

								≈1 for temperatures of 2,000 K and above (provided that radiative

								losses are treated in a separate sink term in the rhs of the energy balance equation). This delayed vibra-

								tional energy relaxation is typically described with an extra equation for the total vibrational energy of N2

								molecules. In the present work, we capture this phenomenology, without the need for an extra equation,

								by adopting a parametric dependence of 𝜂T on temperature, given by 𝜂T = 0.1+0.9[tanh(T∕Tamb−4)+1]∕2.

								The added second term in this expression simulates the acceleration of vibrational energy relaxation,

								yielding 𝜂T = 1 for T >2,000 K with a smooth ramp transition between 1,000–1,500 K.

								DA SILVA ET AL.

								9450


								=== PAGE 10 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								3. Results and Discussion

								3.1. Streamer-to-Leader Transition

								The most fundamental step in the formation of a lightning channel is the streamer-to-leader transition.

								Streamers are the precursor stage. They are thin filamentary discharge channels that propagate as a non-

								linear electron-impact ionization wave, self-enhancing the electric field at its tips. Their conductivity is of

								the order of 0.1–1 S/m. They require electric fields higher than 17% of the conventional breakdown thresh-

								old for stable propagation. Streamer lifetimes are rather short, approximately tens of microseconds, limited

								by attachment to oxygen molecules. Leaders are a necessity for the breakdown of air gaps longer than one

								meter (Bazelyan & Raizer, 2000, p. 59). It takes several milliseconds for a leader to come from the cloud to the

								ground. The only way to keep the leader channel conductive for so long is by substantially heating the air. In

								the hot air plasma, attachment loses its importance; instead, the electron density decays via electron-positive

								ion recombination, which is substantially slower. The transition between streamer and leader happens in

								a region in space called stem, a converging point where several streamers in a streamer corona are rooted.

								In this region the small current carried by individual streamers can add up to values ≳1 A to produce air

								heating and create a leader channel.

								da Silva and Pasko (2013) developed a first-principles model to investigate the dynamics of streamer-to-

								leader transition. It consists of four main blocks: (1) a set of fully nonlinear gas-dynamic equations that

								described the heating and radial expansion of the neutral gas; (2) a detailed kinetic scheme accounting for

								the most important processes in an air discharge plasma; (3) energy exchange between charged and neu-

								tral particles accounting for the partitioning of electronic power between elastic collisions, and excitation

								of vibrational and electronic states; and (4) delayed vibrational energy relaxation of nitrogen molecules. da

								Silva and Pasko's (2013) model was validated against streamer-to-spark transition time scales measured in

								centimeter-long laboratory discharges ( ˇCernák et al., 1995; Larsson, 1998). That model was also applied to

								simulation of leader speeds at reduced air densities and for interpretation of the phenomenology of gigan-

								tic jets (da Silva & Pasko, 2012), as well as to study the mechanism of infrasound emissions in sprites

								(da Silva & Pasko, 2014). Figure 4a shows, as discontinuous traces, the air heating rate calculated with da

								Silva and Pasko's (2013) model with an assumed Gaussian initial distribution of electron density in the

								streamer channel. The peak ne value is 2 ×1020 m−3 and the e-folding spatial scale is rc = 0.3, 0.5, and 1

								mm, respectively. The streamer-to-leader transition time scale 𝜏h is defined as the time required to heat the

								channel up to 2000 K; the heating rate shown in the figure is simply 1∕𝜏h. The 2000-K threshold is chosen

								because when the channel reaches this temperature level a thermal-ionizational plasma instability is trig-

								gered: vibrational relaxation is accelerated, temperature raises very sharply, N + O associative ionization

								starts to take place, and transition to leader mode is unavoidable.

								The present work's goal is to propose the minimal physical model to describe the dynamics of the leader

								plasma. As discussed in section 2.2, the model uses a simplified plasma chemistry and parameterized radial

								dynamics. As a means of validation, in Figure 4 we compare the present model with the simulations of

								da Silva and Pasko (2013). Figure 4a uses the same initial conditions as the previous work and an ini-

								tial current-carrying radius rc = 0.5 mm. The figure shows order-of-magnitude agreement between the two

								models. However, there is an inherently different slope between the two curves, attributed to the multiple

								parameterizations and simplifications introduced in this paper. The other three panels in the figure show

								the effects of the initial conditions in the air heating rate: ne (b), rc (c), and rg (d). The current-carrying

								radius is the parameter that has the largest influence on the heating rate (Figure 4c). The thermal radius rg

								has no effect on the heating rate at all (Figure 4d), because this quantity is exclusively related to the cooling

								rate of the channel (see equation (2)), which is negligible in submicrosecond time scales. The dependence

								on initial electron density is slightly more complicated. The heating rate is ∝∫

								𝜏h

								0

								𝜎E2dt which, according

								to equation (1) is also ∝∫

								𝜏h

								0

								I2∕nedt. The inverse 1∕ne dependence can be qualitatively seen when com-

								paring the 1020- and 1022-m−3 cases. But reducing the initial electron density tends to increase the electric

								field according to Ohm's law. If the electric field goes beyond Ek, ionization increases ne until the field drops

								down to the Ek level. This self-regulatory mechanism imposes a maximum heating rate given by the 1018-

								to 1020-m−3 curves in Figure 4b.

								For the sake of comparison, we have repeated the calculations shown here with a full LTE version of the

								model. This is done by replacing equations (3) an (4) with 𝜎= 𝜎LTE and by setting 𝜂T = 1. The calculated

								air heating rate is in the range of 1012–1015 s−1 for currents between 1 and 100 A. They are not shown in

								Figure 4 because they lie completely outside of the vertical-axis limits. This result indicates that a full-LTE

								DA SILVA ET AL.

								9451


								=== PAGE 11 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 4. Calculated heating rate (1∕𝜏h) leading to the conversion of a streamer into a leader channel. The title in the

								four panels list the initial conditions for electron density ne (ne in the figure), current-carrying radius rc (rc), and

								thermal radius rg (rg) used in the simulations. The ambient neutral temperature is 300 K, and there are no negative

								ions initially. Panel (a) shows as discontinuous traces the calculation of da Silva and Pasko (2013) for the same initial

								conditions, but three different values of rc. The gray shaded area delimiting the calculations of da Silva and Pasko

								(2013) is repeated in all four panels for comparison. Panels (b)–(d) emphasize the effect of changing the initial

								conditions for ne (b), rc (c), and rg (d).

								model completely overestimates the air heating rate, and cannot capture the finite streamer-to-leader (or

								to-spark) transition time scale, well known from laboratory studies to be a fraction of 1 μs ( ˇCernák et al.,

								1995; Larsson, 1998). The reason for the unreasonably high air heating rate of a full-LTE model lies in the

								fact that the LTE conductivity at 300 K is substantially lower than the typical conductivity in a streamer

								channel (see Figure 2e). Since conductivity is lower, the resistance per unit length R is larger, and so is the

								Joule heating rate RI2, which is the same argument presented when discussing Figure 4b.

								In summary, the present model compares very well to a first-principles theoretical simulation that has been

								validated with spark data from laboratory discharges. The proposed computer-simulation tool is able to

								account for the finite time scale of streamer-to-leader transition, something that a full-LTE model cannot.

								The following input parameters are used as initial conditions in all simulations below, unless otherwise

								noted: ne = 1020 m−3, rc = 0.5 mm, rg = 5 mm, nn = 0, T = 300 K.

								DA SILVA ET AL.

								9452


								=== PAGE 12 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 5. (a) Temporal dynamics of resistance in a discharge channel for several current values. Solid and dashed lines

								show the contrast between full model versus suppressed channel expansion, respectively. The figure also shows the

								data by Tanaka et al. (2000) as a solid black line, with the gray shaded area marking ±50% variability. (b–d) Resistance

								value at 10 ms as a function of current. Panel (b) also shows the data from Tanaka et al. (2000) at 10 ms (square with

								±50% error bar), as well as, the steady-state arc resistance measured by King (1961) (black solid line with ±50% gray

								shaded band). Panels (b)–(d) emphasize the effect of changing the initial conditions for ne (b), rc (c), and rg (d), with

								the initial conditions being listed in the panel titles and legends. The gray shaded band marking the results from King

								(1961) are repeated in panels (b-d) for comparison with our simulations.

								3.2. Steady-State Negative Differential Resistance

								The behavior of the steady-state resistance of arc channels has been used to discuss the phenomenology

								of lightning channels (Hare et al., 2019; Heckman, 1992; Krehbiel et al., 1979; Mazur & Ruhnke, 2014;

								Williams, 2006; Williams & Heckman, 2012; Williams & Montanyà, 2019). Steady-state plasma arcs exhibit

								the so-called negative differential resistance, that is, the resistance decreases with increasing electrical cur-

								rent. Such behavior is reproduced in our simulations and shown in Figure 5. Figure 5a shows the temporal

								evolution of resistance in the discharge channel for several values of electrical current between 1 A and

								10 kA. It is easy to see that, owing to channel expansion, there is no true steady-state resistance. A con-

								stant value for the steady-state resistance can only be obtained if channel expansion is suppressed (compare

								the solid and dashed lines). At low currents (see the 1-A curve), one can start to see the channel recovery

								DA SILVA ET AL.

								9453


								=== PAGE 13 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Table 1

								Fit Parameters for the Resistance per Unit Channel Length Formula R = A∕Ib

								Reference

								Current range (A)

								Time scale (s)

								A (Ω Ab/m)

								b

								Mean fit error (%)

								This Work

								100–104

								10−2

								4.27×103

								1.18

								35

								This Work

								100–104

								1

								4.81×103

								1.37

								74

								This Work: Region I

								100–101

								10−2

								1.24×104

								1.84

								9

								This Work: Region II

								101–103

								10−2

								2.82×103

								1.16

								4

								This Work: Region III

								103–104

								10−2

								0.18×103

								0.75

								1

								King (1961)

								100–104

								—

								2.87×103

								1.16

								25

								Bazelyan and Raizer (1998)

								—

								—

								3×104

								2

								—

								starting as early as 0.1 ms. The recovery in this case is due to the fact that the channel cools down to a suffi-

								cient level that three-body attachment becomes important, accelerating the rate of plasma density depletion.

								For currents higher than 10 A, the resistance is still decreasing at the 0.1 s mark; in some cases after a partial

								recovery. In Figure 5a we also show data from Tanaka et al. (2000) used by Chemartin et al. (2009) to validate

								their 3-D free burning arc simulations. Tanaka et al. (2000) report on 1.6-m-long arcs with 100-A current.

								Their measurements are shown in Figure 17 of Chemartin et al. (2009). We obtain a good agreement between

								our simulations and the measurements despite the fact that the 3-D tortuous nature of the arc channel is

								neglected in the present work.

								For the purpose of evaluating the negative differential resistance behavior predicted by our simulations, we

								evaluate the resistance (per unit length) at 10 ms for several different values of electrical current. The results

								are shown in Figure 5b alongside measurements from King (1961). We chose to compare our simulations to

								King's measurements because this work has been featured in a number of manuscripts in lightning-research

								literature (e.g., Heckman, 1992; Mazur & Ruhnke, 2014; Williams, 2006; Williams & Heckman, 2012). The

								data from King (1961) is shown as a black solid line with a ±50% variability gray shaded band. The gray

								band is repeated in panels (b)–(d) for comparison with our simulations. The time instant of 10 ms is chosen

								because it is when the time-dependent data from Tanaka et al. (2000) (shown as a square with ±50% error

								bar) best aligns with King's curve. Our calculations in Figure 5a show good agreement with King's curve;

								the average difference between the two is 40%. Figures 5b–5d show the effects of the initial conditions in the

								steady-state resistance: ne (b), rc (c), and rg (d). It can be seen that changes in the initial conditions have very

								little impact on the resistance in the 10-ms time scale. It is as if the channel “forgets” the initial conditions

								(Aleksandrov et al., 2001). Given the uncertainty in determining the initial conditions of the channel, this

								result lends robustness to the resistance calculations shown hereafter. However, in shorter time scales the

								resistance R does depend on the initial conditions. Similarly to the discussion in section 3.1, the dependence

								on ne and rg is weak, but the dependence on rc can be more noticeable. The dependence on the initial channel

								radius becomes weaker and weaker at higher currents. As an example, at the 10-μs mark, we find that the

								ratio R(rc=2 mm)/R(rc=0.5 mm) is of the order of 700 for a constant current of 10 A. The same ratio is only

								0.63 for a current of 1,000 A.

								The dependence of the resistance on electrical current can be approximated by the analytical formula

								R = A∕Ib, where A and b are positive constants. It is easy to see that with this dependence dR∕dI < 0 always,

								in accordance with the terminology “negative differential resistance.” The limiting case b = 1 corresponds

								to a constant steady-state electric field inside the channel (with numerical value equal to A). We have eval-

								uated the fit parameters that best match our model for the standard set of initial conditions (shown in the

								title of Figure 5a). The results are shown in Table 1 alongside the fit parameters for the King's curve and also

								values given by Bazelyan and Raizer (1998). It can be seen that the exponent b that best fits both the present

								work and King (1961) are very close to each other (b = 1.16–1.18). The empirical trend given by Bazelyan

								and Raizer (1998) has a substantially steeper slope (b = 2). If we run the simulation for a longer time, up to

								1 s, the power law index increases from 1.18 to 1.37 (see second row in Table 1). However, the mean fit error

								doubles indicating that the curve deviates further from the power law approximation.

								It can be seen from Table 1 that fitting the power law dependence to a four-decade current range produces

								errors of 35–74%. A better fit can be produced by braking down the current range in three regions: (I) 100–101

								A, (II) 101–103 A, and (III) 103–104 A. The three regions are marked in Figure 5c. It can be seen in Table 1

								DA SILVA ET AL.

								9454


								=== PAGE 14 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								that the three regions have different power law indexes, progressively lower as current increases. Detailed

								analysis of the temporal evolution of energy deposition in the channel reveals that the steady state is given by

								different mechanisms in the three regions. In Region I the steady state is given by a balance of Joule heating

								and heat conduction, that is, between the first and second terms in the right-hand side of equation (2).

								Meanwhile, In Region III the steady state is given by a balance with radiative emission, that is, between

								the first and third terms in the right-hand side of equation (2). Region II is marked by a comparable role

								between the two loss processes; radiative emission is important in the submillisecond time scale, while heat

								conduction is significant at later stages.

								3.3. Energy Deposition in Return Strokes

								The return stroke follows the attachment of lightning leader channels to ground structures. In the case of a

								negative cloud-to-ground discharge, the return stroke effectively lowers several coulombs of negative charge

								originally deposited along the downward propagating stepped leader. The high-current return stroke wave

								(with typically tens of kiloamperes) rapidly heats the channel to peak temperatures of the order of 30,000

								K, emitting intense optical radiation, and creating a channel expansion shock wave (that produces audible

								thunder). According to Rakov and Uman (1998), models that describe the lightning return stroke can be

								divided into four categories: gas-dynamic or physical, electromagnetic, distributed-circuit, and engineering

								models. The basic set of equations described in this manuscript fits into the first category, where the cur-

								rent flowing through the channel is an input parameter and all other channel properties can be calculated

								from first principles. Some of the most well-accepted investigations within this framework are the papers

								by Plooster (1971) and Paxton et al. (1986). These authors solve the hydrodynamic equations of motion for

								atmospheric-pressure air in a Lagrangian frame of reference. A description of this simulation approach,

								which shows contemporary versions of the pertinent equations, is given by Aleksandrov et al. (2000). The

								model resolves the 1-D radial profiles of all state variables and captures the shock wave expansion as driven

								by ohmic heating. The plasma is assumed to be in LTE and the conductivity is simply 𝜎= 𝜎LTE(T). These

								models also describe the radial transport of radiation, and primarily differ by its implementation and com-

								prehensiveness. Plooster (1971) used a single temperature-independent opacity to obtain radiation loss and

								absorption in each radial grid point, while Paxton et al. (1986) used a detailed multigroup radiative transport

								algorithm using a diffusion approximation. A detailed discussion on plasma radiative transport is given by

								Ripoll et al. (2014a).

								In Figure 6a we present a comparison between our model's results and the seminal works of Plooster (1971)

								and Paxton et al. (1986). The current waveform has the qualitative shape depicted in Figure 1b, with a rise

								time of 5 μs and a fall time of 50 μs (or simply written as 5/50 μs). The peak current is 20 kA, a typical value

								for first return strokes, and no continuing current is incorporated. The current waveform is the same one

								used in the two papers for the simulation case shown in Figure 8 of Paxton et al. (1986). We generate initial

								conditions by starting the simulation with the standard streamer-like channel parameters used in section 3.1

								and running a constant 10-A current through the channel during 4 μs. This strategy ensures that the channel

								has the properties of a leader discharge prior to the return stroke. These initial conditions are rc = 1 mm, rg

								= 1 cm, ne = 9×1017 m−3, and T = 5000 K. Additionally, instead of using the value of 𝜌m(T = 5000K) for

								the air mass density, the ambient value 𝜌m(T = 300K) = 0.7 kg/m3 is used. These initial conditions are very

								similar to the ones used in the aforementioned references. Note that even a steady current as low as 10 A can

								produce a leader with temperature of ∼5,000 K. This value is within the estimate for the predart and postdart

								leader channel temperatures provided by Rakov (1998), which are 3,000 K and 20,000 K, respectively. It

								can be seen from Figure 6a that our model compares very well with simulation results of Plooster (1971),

								predicting a peak temperature of 36,000 K. The mean difference between the two curves is 3%.

								Both curves (Plooster's and ours) deviate from the results of Paxton et al. (1986). It can be seen from Figure 6a

								that a better agreement with Paxton et al. (1986) can be found by simply multiplying the radiative emis-

								sion coefficient (last term in equation (2)) by a factor of 10. This fact can be better understood by looking

								at the energy deposition in the return stroke channel, depicted in Figure 6b. The figure shows (in order)

								the four terms in the energy equation (2): the internal energy is given by 𝜌mcpT𝜋r2

								c, the Joule heating by

								∫𝜂T𝜎E2𝜋r2

								cdt, the thermal conduction by ∫

								4𝜅T

								r2

								g

								(

								T −Tamb

								)

								𝜋r2

								cdt, and the radiative emission by ∫4𝜋𝜖𝜋r2

								cdt.

								It can be seen that the channel's temperature is dictated by a balance between Joule heating and cooling by

								radiative emission. Therefore, simply increasing the rate of channel cooling by radiation can lower the peak

								temperature and provide a better agreement with Paxton's results. As mentioned above, the models pre-

								DA SILVA ET AL.

								9455


								=== PAGE 15 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 6. (a) Evolution of temperature in a 20-kA return stroke channel: comparison between the present investigation and established results (Paxton et al.,

								1986; Plooster, 1971). (b) Energy deposition in the return stroke channel. The four lines, in the order listed in the figure legend, correspond to the four terms in

								the energy balance equation (2). Panel (c) is a zoom-in into the gray shaded rectangle in panel (b). Panels (d)–(f) show the radius, resistance per unit length,

								and rates of electron-positive ion recombination, respectively. Panel (f) justifies a posteriori neglecting the three-body process in equation (3).

								sented by Plooster (1971) and Paxton et al. (1986) are essentially the same and only differ by the treatment

								of radiative emission, lending further credence to the idea that peak temperatures are dictated by radiative

								emission.

								An important conclusion to be drawn here is that the effective representation of the radiative emission

								through a net emission coefficient (𝜖in Figure 2f) produces a proper description of the channel temperature

								dynamics, especially because all four curves in Figure 6a have similar qualitative shape and rate of cooling

								after the peak. Moreover, at 35 μs the total deposited energy in our simulations of 5.6 kJ/m compares well

								to the estimates of 2 and 3.8 kJ/m by Plooster (1971) and Paxton et al. (1986), respectively (see also Rakov &

								Uman, 1998, Table I). The state of the art in lightning spectroscopy is the recent investigations by Walker and

								Christian (2017, 2019). From the ratio of several atomic spectral lines recorded with 1-μs temporal resolution,

								these authors report peak temperatures ranging between 32 and 42 kK for five rocket-triggered lightning

								strikes with peak currents varying between 8.1 and 17.3 kA (Walker & Christian, 2019, Figure 4). There is

								not a clear linear correlation between peak current and peak temperature in their dataset and the average

								peak temperature between the five strikes is ≈36±4 kK. Remarkably, our work and Plooster's do a better job

								reproducing the measured peak temperatures than Paxton's. Further work is required to explain the highest

								value registered by Walker and Christian (2019), in excess of 42,000 K.

								DA SILVA ET AL.

								9456


								=== PAGE 16 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Figure 6d shows the channel radius as a function of time. We have verified that the proposed averaged

								radial dynamics qualitatively captures the radial expansion and also provides order-of-magnitude quantita-

								tive agreement with previous investigations alike (Braginskii, 1958; Koshak et al., 2015; Plooster, 1971). All

								of these models (including ours) predict an initial rapid channel expansion rate, leveling off when the chan-

								nel is cooling down. During the initial return stroke stage (0.5–5 μs), our calculated radius is 8–42% smaller

								than the results obtained by Braginskii (1958) and Plooster (1971), shown in Table II of Plooster (1971).

								Koshak et al. (2015) improved on the channel radial expansion rate derived by Braginskii (1958) and found a

								good agreement with Plooster (1971) at the 35-μs mark. Both investigations yielded a 1.5-cm radius at 35 μs,

								while our simulations yielded a value 57% lower. Generally, the results are in good agreement with previous

								investigations. However, it should be noted that our peak channel expansion rate is ∼500 m/s, which is a

								factor of 4 lower than in Koshak et al. (2015).

								Figure 6e presents the resistance (per unit length) as a function of time. It can be seen that the resistance

								drops by more than two orders of magnitude while the current is rising, illustrating how negative differential

								resistance works for a current changing over time. After that, while the current is decreasing exponentially

								in time, the resistance achieves a stable value between 0.6–1 Ω/m. This leveling off is in agreement with

								the trend seen in measurements (Jayakumar et al., 2006, Figure 4). Jayakumar et al. (2006) measured the

								electrical current to ground and the vertical electric field in close vicinity to a series of rocket-triggered

								lightning strikes in Florida. At the instant of peak power, these authors found resistance values between

								0.67 and 31 Ω/m in eight different strikes. In our calculations, we obtain R = 0.6 Ω/m, which is close to the

								lowest resistance value in their dataset. This value is closer to the measurements than the early estimate

								of 0.035 Ω/m by Rakov (1998). Additionally, Jayakumar et al. (2006) registered input electrical energies

								between 0.9–6.4 kJ/m, also in range with our calculations.

								Figure 6f shows the rates of electron-positive ion recombination. The figure shows a comparison between

								the rate of two- and three-body recombination with coefficients taken from Kossyi et al. (1992). The

								figure is included here to justify the model design assumptions discussed in section 2.2 (item #3). In

								the regime studied here and with the rate coefficients for an air plasma available in the literature, the

								three-body recombination rate is substantially slower than the two-body counterpart, justifying neglecting

								it in equation (3).

								3.4. Behavior of Light Emission in Return Strokes

								The net emission coefficient 𝜖describes the radiative emission in all bands of the optical spectrum, encom-

								passing the infrared, visible, and ultraviolet (Naghizadeh-Kashani et al., 2002). Most of the radiation

								escaping the plasma is in the vacuum ultraviolet range (wavelengths lower than 200 nm) and is caused by

								atomic emissions. However, this band is not easily detected because the radiation is readily absorbed by

								atmospheric-pressure air surrounding the plasma discharge (Cressault et al., 2015). Spectroscopic measure-

								ments of rocket-triggered lightning strikes show characteristic line emissions associated with neutral, singly,

								and doubly ionized nitrogen and oxygen, neutral argon, neutral hydrogen, and neutral copper (from the

								triggering wire) and present no detected molecular emissions (Walker & Christian, 2017).

								For the purposes of comparing our simulations with observations, we estimate the power (per unit channel

								length) emitted in the visible range as 𝜂vis4𝜋𝜖𝜋r2

								c, where 𝜂vis is the fraction of optical radiation emitted in

								the visible range (380–780 nm). We use a constant fraction 𝜂vis = 3% for the sake of simplicity. In reality 𝜂vis

								depends on the radial distribution of the plasma temperature and the cumulative balance of emission and

								absorption. Table 2 shows seven estimates of 𝜂vis based on different references and techniques. Perhaps the

								most pertinent is estimate #2, which is calculated by taking the ratio of 𝜖vis in the visible range calculated

								by Cressault et al. (2011, Figure 2) to 𝜖in the total optical range calculated by (Naghizadeh-Kashani et al.,

								2002, Figure 13) for an optically thin plasma. This strategy places 𝜂vis between 0.1% and 10% in the tempera-

								ture range between 3,000 and 30,000 K. Within this range, we adopt the value of 3% because it yields a good

								agreement with experimental data from Quick and Krider (2017) discussed below.

								Figure 7 shows properties of return stroke light emission and comparison to rocket-triggered lightning data

								collected by Quick & Krider (2017, Figures 15 and 16). From a 200-m distance to the lightning striking point,

								Quick and Krider (2017) recorded the luminosity of a 62-m-long channel segment near the ground. The

								radiometers used had an approximately flat spectral response in the 400- to 1,000-nm range. Figure 7a shows

								the simulated temporal dynamics of visible power and electrical current in the channel, for conditions that

								resemble the aforementioned observations. The waveform is 0.5/50 μs with a 12-kA peak current, similar to

								DA SILVA ET AL.

								9457


								=== PAGE 17 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Table 2

								Fraction of Optical Power Radiated in the Visible Range by an Air Plasma

								#

								Estimation method and reference

								𝜂vis (%)

								1

								Black-body spectral radiance (Siegel, 2001, p. 22) (3,000–30,000 K)

								5.3–49

								2

								Visible 𝜖vis calculated by Cressault et al. (2011) (3,000–30,000 K)

								0.1–10

								3

								Visible radiance calculated by Cressault et al. (2015) (8,000–30,000 K)

								0.2–0.6

								4

								20-kJ/m hot air shock (Ripoll et al., 2014a, Figure 9 and section 3.1 )

								14.3

								5

								Several simulations in Table 1 of Ripoll et al. (2014a)

								4–30

								6

								Section 4.2 of Ripoll et al. (2014a)

								5.3–21.7

								7

								A 12-kA discharge (Ripoll et al., 2014b, Figures 9b and 10b)

								30

								Empirical (this work)

								3

								the median case in the data set (Quick & Krider, 2017, Table 1). Figure 7b shows the first 3 μs of light emis-

								sion, evidencing a 0.1-μs delay between the rise of current and optical emissions in the channel. Figure 7c

								shows the effects of increasing peak current, which lead to higher emitted power and longer duration of the

								light emission.

								The delay shown in Figure 7b is evaluated at the 20% of peak level. The 0.1-μs value is in excellent agreement

								with experimental results by Carvalho et al. (2014, 2015) and Quick and Krider (2017) who found delays of

								Figure 7. (a) Temporal evolution of power per unit channel length emitted by a return stroke in the visible range (left-hand side axis) and electrical current

								(right-hand side). Panel (b) is a zoom-in into the gray shaded rectangle in panel (a). (c) Visible power emitted for several different peak current values.

								(d) Visible peak power versus peak current for four different current waveforms. (e) Energy emitted in the visible range versus charge transferred to the ground

								(the integration time is 2 ms). Panels (d) and (e) show a comparison with the experimental data from Quick and Krider (2017). The big crosses indicate the

								average ± standard deviation in the dataset. The data were collected during a study conducted by the University of Arizona at the International Center for

								Lightning Research and Testing, in Camp Blanding, FL, in 2012.

								DA SILVA ET AL.

								9458


								=== PAGE 18 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								0.09 ± 0.05 and 0.09 ± 0.06 μs, respectively. Differently than Quick and Krider (2017), Carvalho et al. (2014)

								recorded luminosity from a 3-m-long channel segment near the ground. From such a short segment, the

								luminosity rise time is not masked by the geometrical growth of the return stroke in the field of view. The

								fact that both experimental investigations observing different channel lengths (62 and 3 m) yielded similar

								results lends robustness to the ∼0.1 μs measured delay. Furthermore, analysis of different types of pulses

								occurring in the return stroke channel (Zhou et al., 2014) and of several channel segments at different heights

								(Carvalho et al., 2015) have led to the general conclusion that current and luminosity have similar rise times

								and the delay between the two has the same order of magnitude as such time scales. More precisely, Carvalho

								et al. (2015) found that the delay is approximately linearly dependent on the current rise time according to

								the following fit formula: delay = 0.35 𝜏1.03

								rise , where 𝜏rise is the 10–90% current rise time given in microseconds.

								The fit comprises rise times between ∼0.1 μs (for return strokes) and ∼100 μs (for M components). Using

								this formula, we obtain a delay of 0.14 μs for the simulation shown in Figure 7b, once more indicating good

								agreement between simulation and measurements.

								In our simulations the delay between the rise of current and optical emissions highlighted in Figure 7b has

								a clear interpretation. It is attributed to the finite time scale of channel heating and expansion. Since the

								initial channel temperature for the simulations shown in this section is 5000 K, non-LTE effects play a minor

								role here. From equations (1) and (2), the air heating rate thus is 𝜕T∕𝜕t ≃(I2∕𝜎LTE𝜋2r4

								c −4𝜋𝜖)∕𝜌mcp. What

								determines the finite 0.1-μs delay, in a return stroke with 0.5-μs rise time, are the coefficients 𝜌m, cp, 𝜎LTE,

								and 𝜖, as well as the channel expansion rc(t). A comparison with a full-LTE version of the simulation code

								yielded a similar time delay between current and optical emissions, but the full-LTE model overestimated

								the peak optical power by a factor of 3–4.

								Figures 7d and 7e show the peak visible power versus peak current and total energy versus charge, respec-

								tively. The integration time for the charge and energy is 2 ms. The figures show simulations for different

								current waveforms and comparison with light emitted by rocket-triggered lightning. The data correspond

								to optical irradiance from 55 rocket-triggered lightning strikes (with currents and charges ranging between

								3–20 kA and 0.3–3 C, respectively) observed in Florida by Quick and Krider (2017) in 2012. The irradi-

								ance data is converted to power per unit channel length according to equation (2) in the original reference.

								The simulations use the same initial conditions as in Figure 6, and the results indicate a direct relationship

								between current and power and between total energy and charge. Additionally, the calculations (under the

								𝜂vis = 3% assumption) present good agreement with the observational data, especially near the average val-

								ues (the big crosses in the figures). The peak visible power shows little dependence on the current waveform

								parameters in the range used (𝜏r = 0.5 and 5 μs, and 𝜏f = 50 and 150 μs). The rise time also does not affect the

								relationship between energy emitted and charge transferred to the ground, shown in Figure 7e. The same

								figure also shows that strokes with a narrower current pulse (i.e., with shorter fall time) are more efficient

								in converting electrical energy into optical.

								There are two important issues that must be noted about the comparison made in Figures 7d and 7e.

								First, the radiometers used by Quick and Krider (2017) have a flat spectral response in the 400-1,000

								nm range. According to Ripoll et al. (2014a, 2014b), about twice as much energy is emitted in this range

								than in the visible, because it includes part of the infrared spectrum. Second, Quick and Krider (2017)

								state that rocket-triggered lightning strikes radiate around half as much energy as first strikes in natural

								cloud-to-ground flashes. But the simulations use initial conditions that best resemble first strikes in natu-

								ral lightning, similarly to the works by Plooster (1971) and Paxton et al. (1986). Therefore, if we attempt

								to scale the numerical results to correspond to optical power emitted in the 400-1,000 nm range (×2) by

								rocket-triggered lightning (×1/2), the factors of two cancel and the curves would stay in the same place in

								Figures 7d and 7e, which lends further credence to the comparison. Nonetheless, it should be noted that

								our numerical investigations did not capture the approximate quadratic scaling between peak luminosity

								and peak current, that is, luminosity ∝I2

								p, seen in observations (Carvalho et al., 2015; Quick & Krider, 2017;

								Zhou et al., 2014). Further work is required to explain all experimentally inferred relationships between

								current and luminosity derived from close-by observations of rocket-triggered lightning.

								When analyzing the light emission of return strokes, two additional factors must be noted. First, in

								rocket-triggered lightning there is a nonnegligible amount of copper emission within the visible spectrum,

								DA SILVA ET AL.

								9459


								=== PAGE 19 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								arising from the vaporization of the copper wire that connects the rocket to the ground (Walker &

								Christian, 2017). Second, there is a geometric growth effect of the optical emission within the field of view

								of the detector. For the sake of simplicity, these two effects are neglected in the simulations by assuming that

								the fraction of total energy radiated by neutral copper is small in comparison to all other emissions from

								the air plasma, and by assuming that within the narrow field of view of the detector (only 62 m of chan-

								nel length) the return stroke current amplitude does not change considerably. All these uncertainties are

								encapsulated within the parameter 𝜂vis, adjusted within reason to fit the measurements.

								In all simulations shown in Figure 7, the total energy deposited in the channel by Joule heating ranges

								between 10 J/m and 18 kJ/m. At the instant of peak electrical power, the channel resistance varies between

								0.6–130 Ω/m within all simulation cases presented in this section. For peak currents larger than 5 kA, this

								quantity shows little dependence on the current rise time and fall time values used, and can be fitted by

								the following formula R = A∕Ip, where A = 13 kA Ω/m (the mean error between fit and simulation results

								is lower than 3%). From this formula it is easy to see that in the range of peak currents between 10 and 20

								kA, the channel resistance per unit length at the instant of peak electrical power reduces from 1.3 to 0.65

								Ω/m. Once again these values are in good agreement with the experimental findings of Jayakumar et al.

								(2006, Table 2).

								4. Summary and Conclusions

								In summary, in this manuscript we introduced, validated, and used a physics-based computational tool to

								calculate the lightning channel's nonlinear plasma resistance. A model that bridges an existing gap in the

								literature, by providing a self-consistent evaluation of the plasma properties at little computational cost

								(i.e., at the cost of solving five ordinary differential equations). In this paper, we showed how the proposed

								computer-simulation tool can perform well in a wide range of current values, from 1 to 104 A. It can capture

								well the non-LTE plasma regime, by reproducing the finite time scale for streamer-to-leader transition with

								reasonable accuracy. Furthermore, in the high-current/full-LTE regime, the model can capture well the

								temporal evolution of the neutral-gas temperature and the estimated energy deposition by a return stroke,

								in good agreement with the work of Plooster (1971) and Paxton et al. (1986).

								The model also describes well the negative differential resistance behavior of steady-state arc discharges, in

								good agreement with the experimental findings of King (1961) and Tanaka et al. (2000). The steady-state

								resistance in the millisecond time scale has an inverse power law dependence on the current, that is,

								R = A∕Ib, where A and b are fitting constants. We found that the power law index b decreases with increas-

								ing current, because at different current regimes the steady state is dictated by distinct physical processes.

								At low currents (I < 10 A) the steady state is given by a balance of Joule heating and heat conduction, while

								at high currents (I > 1 kA) the steady state is given by a balance with radiative losses. The intermediate cur-

								rent range is marked by a comparable role between the two loss processes, with radiative emission being

								important in the submillisecond time scale, while heat conduction being significant at later stages.

								We presented a detailed description of the light emission in a return stroke. We showed that the proposed

								model can reproduce the experimentally inferred direct relationship between peak current and peak radi-

								ated power and between charge transferred to ground and total energy radiated, as experimentally inferred

								by Quick and Krider (2017). The caveat is that the quadratic power law relationship between the two remains

								unexplained. The model also captures the 0.1-μs delay between the rise of current and optical emissions in

								rocket-triggered lightning return strokes, as measured with high precision by Carvalho et al. (2014, 2015).

								It has been suggested that the negative differential resistance behavior of lightning channels plays an impor-

								tant role in the mechanism of current cutoff, which in its turn makes some flashes transfer charge to ground

								by a series of (discrete) return strokes, while others by a single stroke followed by a long continuing current

								(Krehbiel et al., 1979; Hare et al., 2019; Heckman, 1992; Mazur et al., 1995). Recent review articles argue that

								the role of negative differential resistance in the channel cutoff remains to be quantified (Mazur & Ruhnke,

								2014; Williams, 2006; Williams & Heckman, 2012; Williams & Montanyà, 2019). The model described in

								this manuscript can be applied for simulating multiple return strokes in a flash and other types of processes

								taking place in the lightning channel, such as dart-leader ionization waves and M components, provided

								that the current waveform is given (see Figure 1b). Suggestions of future work include coupling this tool to

								DA SILVA ET AL.

								9460


								=== PAGE 20 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								distributed circuit models of the lightning return stroke, or to fractal models of the growing lightning-leader

								network. We speculate that this strategy will provide important insights into the physics of lightning channel

								cutoff.

								References

								Aleksandrov, N. L., Bazelyan, E. M., Kochetov, I. V., & Dyatko, N. A. (1997). The ionization kinetics and electric field in the leader channel

								in long air gaps. Journal of Physics D: Applied Physics, 30, 1616–1624. https://doi.org/10.1088/0022-3727/30/11/011

								Aleksandrov, N. L., Bazelyan, E. M., & Konchakov, A. M. (2001). Plasma parameters in the channel of a long leader in air. Plasma Physics

								Reports, 27, 875–885. https://doi.org/10.1134/1.1409721

								Aleksandrov, N. L., Bazelyan, E. M., & Shneider, M. N. (2000). Effect of continuous current during pauses between successive strokes on

								the decay of the lightning channel. Plasma Physics Reports, 26(10), 893–901. https://doi.org/10.1134/1.1316830

								Baker, L. (1990). Return-stroke transmission line model. In R. L. Gardner (Ed.), Lightning electromagnetics (pp. 63–74). Philadelphia, PA:

								Taylor and Francis.

								Barannik, S. I., Vasserman, S. B., & Lukin, A. N. (1975). Resistance and inductance of a gas arc. Soviet physics. Technical physics, 19(11),

								1499–1453.

								Bazelyan, E. M., & Raizer, Y. P. (1998). Spark discharge. Boca Raton/Florida: CRC Press.

								Bazelyan, E. M., & Raizer, Y. P. (2000). Lightning physics and lightning protection. Bristol/Philadelphia: Institute of Physics Publishing.

								Bazelyan, E. M., Raizer, Y. P., & Aleksandrov, N. L. (2007). The effect of reduced air density on streamer-to-leader transition and on

								properties of long positive leader. Journal of Physics D: Applied Physics, 40(14), 4133–4144. https://doi.org/10.1088/0022-3727/40/14/007

								Becker, K. H., Kogelschatz, U., Schoenbach, K. H., & Barker, R. J. (2004). Non-equilibrium air plasmas at atmospheric pressure (pp. 700).

								Boca Raton, FL: CRC Press. https://doi.org/10.1201/9781482269123

								Benilov, M. S., & Naidis, G. V. (2003). Modelling of low-current discharges in atmospheric-pressure air taking account of non-equilibrium

								effects. Journal of Physics D: Applied Physics, 36, 1834–1841. https://doi.org/10.1088/0022-3727/36/15/314

								Bittencourt, J. A. (2004). Fundamentals of plasma physics (3rd ed.). New York: Springer.

								Boulos, M. I., Fauchais, P., & Pfender, E. (1994). Thermal plasmas, fundamentals and applications (Vol. 1). New York: Plenum.

								Braginskii, S. I. (1958). Theory of the development of a spark channel. Soviet Physics-JETP, 34(6), 1068–1074.

								Bruce, C. E. R., & Golde, R. H. (1941). The lightning discharge. Journal of the Institution of Electrical Engineers—Part II: Power Engineering,

								88(6), 487–505. https://doi.org/10.1049/ji-2.1941.0065

								Carvalho, F. L., Jordan, D. M., Uman, M. A., Ngin, T., Gamerota, W. R., & Pilkey, J. T. (2014). Simultaneously measured lightning return

								stroke channel-base current and luminosity. Geophysical Research Letters, 41, 7799–7805. https://doi.org/10.1002/2014GL062190

								Carvalho, F. L., Uman, M. A., Jordan, D. M., & Ngin, T. (2015). Lightning current and luminosity at and above channel bottom for return

								strokes and M-components. Journal Geophysical Research: Atmospheres, 120, 10,645–10,663. https://doi.org/10.1002/2015JD023814

								ˇCernák, M., van Veldhuizen, E. M., Morva, I., & Rutgers, W. R. (1995). Effect of cathode surface properties on glow-to-arc transition in a

								short positive corona gap in ambient air. Journal of Physics D: Applied Physics, 28, 1126–1132. https://doi.org/10.1088/0022-3727/28/6/

								015

								Chemartin, L., Lalande, P., Montreuil, E., Delalondre, C., Cheron, B. G., & Lago, F. (2009). Three dimensional simulation of a DC free

								burning arc. Application to lightning physics. Atmospheric Research, 91(2-4), 371–380. https://doi.org/10.1016/j.atmosres.2008.07.009

								Cho, M., & Rycroft, M. J. (1998). Computer simulation of the electric field structure and optical emission from cloud-top to the ionosphere.

								Journal of Atmospheric and Solar-Terrestrial Physics, 60(7-9), 871–888. https://doi.org/10.1016/S1364-6826(98)00017-0

								Cressault, Y., Bauchire, J. M., Hong, D., Rabat, H., Riquel, G., Sanchez, F., & Gleizes, A. (2015). Radiation of long and high power arcs.

								Journal of Physics D: Applied Physics, 48(41), 415201. https://doi.org/10.1088/0022-3727/48/41/415201

								Cressault, Y., Gleizes, A., & Riquel, G. (2011). Calculation of the radiation emitted by isothermal arc plasmas in air and air-metal mixtures.

								In XX International Symposium on Plasma Chemistry (ISPC), Philadelphia, USA.

								da Silva, C. L. (2015). Numerical modeling of leader discharge mechanisms in lightning, blue jets, gigantic jets, and sprites (Ph.D. thesis),

								Penn State University University Park, PA.

								da Silva, C. L. (2019a). Rate coefficients used in the plasma nature of lightning channels and the resulting nonlinear resistance. Zenodo,

								Software https://doi.org/10.5281/zenodo.2597562

								da Silva, C. L. (2019b). Simulation data output used in the plasma nature of lightning channels and the resulting nonlinear resistance.

								Zenodo, Dataset https://doi.org/10.5281/zenodo.2597552

								da Silva, C. L., Merrill, R. A., & Pasko, V. P. (2016). Mathematical constraints on the use of transmission line models to investigate the

								preliminary breakdown stage of lightning flashes. Radio Science, 51, 367–380. https://doi.org/10.1002/2015RS005853

								da Silva, C. L., & Pasko, V. P. (2012). Simulation of leader speeds at gigantic jet altitudes. Geophysical Research Letters, 39, L13805. https://

								doi.org/10.1029/2012GL052251

								da Silva, C. L., & Pasko, V. P. (2013). Dynamics of streamer-to-leader transition at reduced air densities and its implications for propaga-

								tion of lightning leaders and gigantic jets. Journal of Geophysical Research: Atmospheres, 118, 13,561–13,590. https://doi.org/10.1002/

								2013JD020618

								da Silva, C. L., & Pasko, V. P. (2014). Infrasonic acoustic waves generated by fast air heating in sprite cores. Geophysical Research Letters,

								41, 1789–1795. https://doi.org/10.1002/2013GL059164

								De Conti, A., Visacro, S., Theethayi, N., & Cooray, V. (2008). A comparison of different approaches to simulate a nonlinear channel

								resistance in lightning return stroke models. Journal of Geophysical Research, 113, D14129. https://doi.org/10.1029/2007JD009395

								Engel, T. G., Donaldson, A. L., & Kristiansen, M. (1989). The pulsed discharge arc resistance and its functional behavior. IEEE Transactions

								on Plasma Science, 17(2), 323–329. https://doi.org/10.1109/27.24643

								Flitti, A., & Pancheshnyi, S. (2009). Gas heating in fast pulsed discharges in N2–O2 mixtures. The European Physical Journal Applied

								Physics, 45(2), 21001. https://doi.org/10.1051/epjap/2009011

								Gallimberti, I. (1979). The mechanism of the long spark formation. Journal de Physique Colloques, 40(C7), 193–250.

								Gallimberti, I., Bacchiega, G., Bondiou-Clergerie, A., & Lalande, P. (2002). Fundamental processes in long air gap discharges. Comptes

								Rendus Physique, 3, 1335–1359. https://doi.org/10.1016/S1631-0705(02)01414-7

								Hare, B. M., Scholten, O., Dwyer, J., Trinh, T. N. G., Buitink, S., ter Veen, S., et al. (2019). Needle-like structures discovered on positively

								charged lightning branches. Nature, 568(7752), 360–363. https://doi.org/10.1038/s41586-019-1086-6

								Heckman, S. (1992). Why does a lightning flash have multiple strokes? (Ph.D. thesis), Massachusetts Institute of Technology, Cambridge.

								Acknowledgments

								This research was supported by a

								Research Infrastructure Development

								award from New Mexico NASA

								EPSCoR. We have made the simulation

								data output (da Silva, 2019b) and rate

								coefficients (da Silva, 2019a) shown in

								this manuscript publicly available

								online.

								DA SILVA ET AL.

								9461


								=== PAGE 21 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Heidler, F. (1985). Traveling current source model for LEMP calculation, Proc. 6th Int. Zurich Symp. Electromagnetic Compatibility

								(pp. 157–162). Switzerland: Zurich.

								Hill, R. . D. (1971). Channel heating in return-stroke lightning. Journal of Geophysical Research, 76(3), 637–645. https://doi.org/10.1029/

								JC076i003p00637

								Jayakumar, V., Rakov, V. A., Miki, M., Uman, M. A., Schnetzer, G. H., & Rambo, K. J. (2006). Estimation of input energy in rocket-triggered

								lightning. Geophysical Research Letters, 33, L05702. https://doi.org/10.1029/2005GL025141

								King, L. A. (1961). The voltage gradient of the free-burning arc in air or nitrogen. British Electrical and Allied Industries Research

								Association, Reference G/XT172.

								Koshak, W. J., Solakiewicz, R. J., & Peterson, H. S. (2015). A return stroke NOx production model. Journal of the Atmospheric Sciences,

								72(2), 943–954. https://doi.org/10.1175/JAS-D-14-0121.1

								Kossyi, I. A., Kostinsky, A. Y., Matveyev, A. A., & Silakov, V. P. (1992). Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen

								mixtures. Plasma Sources Science and Technology, 1, 207–220. https://doi.org/10.1088/0963-0252/1/3/011

								Krehbiel, P. R., Brook, M., & McCrory, R. A. (1979). An analysis of the charge structure of lightning discharges to ground. Journal of

								Geophysical Research, 84(C5), 2432–2456. https://doi.org/10.1029/JC084iC05p02432

								Kushner, M. J., Kimura, W. D., & Byron, S. R. (1985). Arc resistance of laser-triggered spark gaps. Journal of Applied Physics, 58(5),

								1744–1751. https://doi.org/10.1063/1.336023

								Larsson, A. (1998). The effect of a large series resistance on the streamer-to-spark transition in dry air. Journal of Physics D: Applied Physics,

								31, 1100–1108. https://doi.org/10.1088/0022-3727/31/9/011

								Liang, C., Carlson, B., Lehtinen, N., Cohen, M., Marshall, R. A., & Inan, U. S. (2014). Differing current and optical return stroke speeds in

								lightning. Geophysical Research Letters, 41, 2561–2567. https://doi.org/10.1002/2014GL059703

								Liu, N. Y., & Pasko, V. P. (2004). Effects of photoionization on propagation and branching of positive and negative streamers in sprites.

								Journal of Geophysical Research, 109, A04301. https://doi.org/10.1029/2003JA010064

								Lowke, J. J., Kovitya, P., & Schmidt, H. P. (1992). Theory of free-burning arc columns including the influence of the cathode. Journal of

								Physics D: Applied Physics, 25(11), 1600–1606. https://doi.org/10.1088/0022-3727/25/11/006

								Luque, A., & Gordillo-Vázquez, F. J. (2012). Mesospheric electric breakdown and delayed sprite ignition caused by electron detachment.

								Nature Geoscience, 5(1), 22–25. https://doi.org/10.1038/NGEO1314

								Marode, E., Bastien, F., & Bakker, M. (1979). A model of the streamer-induced spark formation based on neutral dynamics. Journal of

								Applied Physics, 50, 140–146. https://doi.org/10.1063/1.325697

								Mattos, M. A. F., & Christopoulos, C. (1990). A model of the lightning channel, including corona, and prediction of the generated

								electromagnetic fields. Journal of Physics D: Applied Physics, 23(1), 40. https://doi.org/10.1088/0022-3727/23/1/007

								Mazur, V., Krehbiel, P. R., & Shao, X.-M. (1995). Correlated high-speed video and radio interferometric observations of a cloud-to-ground

								lightning flash. Journal Geophysical Research, 100(D12), 25,731–25,753. https://doi.org/10.1029/95JD02364

								Mazur, V., & Ruhnke, L. H. (2014). The physical processes of current cutoff in lightning leaders. Journal Geophysical Research: Atmospheres,

								119, 2796–2810. https://doi.org/10.1002/2013JD020494

								Montano, R., Becerra, M., Cooray, V., Rahman, M., & Liyanage, P. (2006). Resistance of spark channels. IEEE Transactions on Plasma

								Science, 34(5), 1610–1619. https://doi.org/10.1109/TPS.2006.883350

								Morrow, R., & Lowke, J. J. (1997). Streamer propagation in air. Journal of Physics D: Applied Physics, 30, 614–627. https://doi.org/10.1088/

								0022-3727/30/4/017

								Naghizadeh-Kashani, Y., Cressault, Y., & Gleizes, A. (2002). Net emission coefficient of air thermal plasmas. Journal of Physics D: Applied

								Physics, 35(22), 2925. https://doi.org/10.1088/0022-3727/35/22/306

								Naidis, G. (1999). Simulation of streamer-to-spark transition in short non-uniform air gaps. Journal of Physics D: Applied Physics, 32(20),

								2649–2654. https://doi.org/10.1088/0022-3727/32/20/311

								Naidis, G. V. (2005). Dynamics of streamer breakdown of short non-uniform air gaps. Journal of Physics D: Applied Physics, 38(21),

								3889–3893. https://doi.org/10.1088/0022-3727/38/21/009

								Naidis, G. V. (2009). Positive and negative streamers in air: Velocity-diameter relation. Physical Review E, 79(5), 57401. https://doi.org/10.

								1103/PhysRevE.79.057401

								Pancheshnyi, S., Nudnova, M., & Starikovskii, A. (2005). Development of a cathode-directed streamer discharge in air at different

								pressures: Experiment and comparison with direct numerical simulation. Physical Review E, 71, 16407. https://doi.org/10.1103/

								PhysRevE.71.016407

								Paxton, A. H., Gardner, R. L., & Baker, L. (1986). Lightning return stroke. A numerical calculation of the optical radiation. Physics Fluids,

								29(8), 2736–2741. https://doi.org/10.1063/1.865514

								Plooster, M. N. (1971). Numerical model of the return stroke stroke of the lightning discharge. Physics Fluids, 14(10), 2124–2133. https://

								doi.org/10.1063/1.1693303

								Popov, N. A. (2001). Investigation of the mechanism for rapid heating of nitrogen and air in gas discharges. Plasma Physics Reports, 27(10),

								886–896. https://doi.org/10.1134/1.1409722

								Popov, N. A. (2003). Formation and development of a leader channel in air. Plasma Physics Reports, 29(8), 695–708. https://doi.org/10.

								1134/1.1601648

								Popov, N. A. (2009). Study of the formation and propagation of a leader channel in air. Plasma Physics Reports, 35, 785–793. https://doi.

								org/10.1134/S1063780X09090074

								Popov, N. A. (2011). Fast gas heating in a nitrogen-oxygen discharge plasma: I. Kinetic mechanism. Journal of Physics D: Applied Physics,

								44, 285201. https://doi.org/10.1088/0022-3727/44/28/285201

								Quick, M. G., & Krider, P. E. (2017). Optical power and energy radiated by return strokes in rocket-triggered lightning. Journal Geophysics

								Research: Atmospheres, 122, 8816–8832. https://doi.org/10.1002/2017JD027363

								Raizer, Y. P. (1991). Gas discharge physics. New York: Springer-Verlag.

								Rakov, V. A. (1998). Some inferences on the propagation mechanisms of dart leaders and return strokes. Journal Geophysics Research,

								103(D2), 1879–1887. https://doi.org/10.1029/97JD03116

								Rakov, V. A., & Uman, M. A. (1998). Review and evaluation of lightning return stroke models including some aspects of their application.

								IEEE Transactions on Electromagnetic Compatibility, 40(4), 403–426. https://doi.org/10.1109/15.736202

								Rakov, V. A., & Uman, M. A. (2003). Lightning: Physics and effects (pp. 687). Cambridge: University Press.

								Rayment, S. W., & Moruzzi, J. L. (1978). Electron detachment studies between O−ions and nitrogen. International Journal of Mass

								Spectrometry and Ion Processes, 26(3), 321–326. https://doi.org/10.1016/0020-7381(78)80033-3

								Riousset, J. A., Pasko, V. P., & Bourdon, A. (2010). Air-density-dependent model for analysis of air heating associated with streamers,

								leaders, and transient luminous events. Journal Geophysical Research, 115, A12321. https://doi.org/10.1029/2010JA015918

								DA SILVA ET AL.

								9462


								=== PAGE 22 ===

								Journal of Geophysical Research: Atmospheres

								10.1029/2019JD030693

								Ripoll, J.-F., Zinn, J., Colestock, P. L., & Jeffery, C. A. (2014b). On the dynamics of hot air plasmas related to lightning discharges: 2.

								Electrodynamics. Journal Geophysical Research: Atmospheres, 119, 9218–9235. https://doi.org/10.1002/2013JD020068

								Ripoll, J.-F., Zinn, J., Jeffery, C. A., & Colestock, P. L. (2014a). On the dynamics of hot air plasmas related to lightning discharges: 1. Gas

								dynamics. Journal Geophysical Research: Atmospheres, 119, 9196–9217. https://doi.org/10.1002/2013JD020067

								Sansonnens, L., Haidar, J., & Lowke, J. J. (2000). Prediction of properties of free burning arcs including effects of ambipolar diffusion.

								Journal of Physics D: Applied Physics, 33(2), 148. https://doi.org/10.1088/0022-3727/33/2/309

								Siegel, R. (2001). Thermal radiation heat transfer (4th ed., Vol. 1, pp. 864). New York: CRC Press.

								Takaki, K., & Akiyama, H. (2001). The resistance of a high-current pulsed discharge in nitrogen. Japanese Journal of Applied Physics,

								40(Part 1, No. 2B), 979–983. https://doi.org/10.1143/jjap.40.979

								Tanaka, S., Sunabe, K., & Goda, Y. (2000). Three dimensional behaviour analysis of D.C. free arc column by image processing technique.

								In XIII Int'l. Conf. on Gas Discharges and Their Applications, Glasgow.

								Tanaka, M., Terasaki, H., Ushio, M., & Lowke, J. J. (2003). Numerical study of a free-burning argon arc with anode melting. Plasma

								Chemistry and Plasma Processing, 23(3), 585–606. https://doi.org/10.1023/A:1023272007864

								Theethayi, N., & Cooray, V. (2005). On the representation of the lightning return stroke process as a current pulse propagating along a

								transmission line. IEEE Transactions on Power Delivery, 20(2), 823–837. https://doi.org/10.1109/TPWRD.2004.839188

								Vidal, F., Gallimberti, I., Rizk, F. A. M., Johnston, T. W., Bondiou-Clergerie, A., Comtois, D., et al. (2002). Modeling of the air plasma near

								the tip of the positive leader. IEEE Transactions on Plasma Science, 30(3), 1339–1349. https://doi.org/10.1009/TPS.2002.801538

								Walker, T. D., & Christian, H. J. (2017). Triggered lightning spectroscopy: Part 1. A qualitative analysis. Journal of Geophysical Research:

								Atmospheres, 122, 8000–8011. https://doi.org/10.1002/2016JD026419

								Walker, T. D., & Christian, H. J. (2019). Triggered lightning spectroscopy: 2. A quantitative analysis. Journal of Geophysical Research:

								Atmospheres, 124, 3930–3942. https://doi.org/10.1029/2018JD029901

								Williams, E. R. (2006). Problems in lightning physics—The role of polarity asymmetry. Plasma Sources Science and Technology, 15,

								S91—S108. https://doi.org/10.1088/0963-0252/15/2/S12

								Williams, E. R., & Heckman, S. (2012). Polarity asymmetry in lightning leaders: The evolution of ideas on lightning behavior from strikes

								to aircraft. Journal Aerospace Lab, 5, AL05–04.

								Williams, E., & Montanyà, J. (2019). A closer look at lightning reveals needle-like structures. Nature, 568, 319–320. https://doi.org/10.1038/

								d41586-019-01178-7

								Winn, W. P., Aulich, G. D., Hunyady, S. J., Eack, K. B., Edens, H., Krehbiel, P. R., et al. (2011). Lightning leader stepping, K changes, and

								other observations near an intracloud flash. Journal Geophysical Research, 116, D23115. https://doi.org/10.1029/2011JD015998

								Yos, J. M. (1963). Transport properties of nitrogen, hydrogen, oxygen, and air to 30,000◦K (Tech. Rep.) Wilmington, MA: AVCO Corp.

								Tech. Memo. RAD-TM-63-7.

								Zhou, M., Wang, D., Wang, J., Takagi, N., Gamerota, W. R., Uman, M. A., et al. (2014). Correlation between the channel-bottom light

								intensity and channel-base current of a rocket-triggered lightning flash. Journal Geophysical Research: Atmospheres, 119, 13,457–13,473.

								https://doi.org/10.1002/2014JD022367

								DA SILVA ET AL.

								9463