Remote Sensing of Source Currents of Narrow Bipolar Events Using Measured Electric Fields

In this work the electric field of narrow bipolar events (NBEs) measured at a remote location is used to extract the current waveform of the source discharge. All calculations correspond to a vertical linear current source above a perfectly conducting ground plane. The current study uses the well established formulation of electromagnetic fields in the frequency domain, and develops a deconvolution based technique to obtain exact reconstruction of the source current, improving upon previous modeling of NBEs, which often require tuning several inter‐dependent parameters to determine the current that best reproduces the observed electric field. Our proposed solution, although readily available in standard electromagnetic textbooks, has never been employed in the context of lightning related discharges, and offers a simple and efficient alternative to previous conventional time domain calculations.


Introduction
Various aspects of lightning physics are not well understood.Primary among them is the question of how lightning is initiated inside thunderclouds.The study of radio signals emitted by the underlying processes involved in the initial stages of lightning development provide a powerful means of gaining insight into said processes.There are various types of bipolar pulses radiated during the early stages of lightning development.These include initial breakdown pulses (IBPs), which are relatively slow-rising, wide bipolar waveforms, which often have several small pulses superimposed on the initial half cycle (Weidman & Krider, 1979;Karunarathne et al., 2014), narrow bipolar events (NBEs), also known as compact intracloud discharges, which are bipolar waveforms with a short risetime (1-2 μs), pulse width of 10-30 μs, and opposite polarity half-cycle of lower amplitude and longer duration (Le Vine, 1980;Willett et al., 1989;Smith et al., 1999), and energetic in-cloud pulses (EIPs), which are intracloud discharges with high peak currents (>200 kA) and significantly longer duration than NBEs (Lyu et al., 2015).
Most NBEs occur in isolation of any other discharge activity in thunderclouds (Smith et al., 1999), but some are the initiating event of lightning flashes (Rison et al., 1999;Nag et al., 2010;Wu et al., 2014;Lyu et al., 2015;Rison et al., 2016;Lyu et al., 2019).In a study of 157 NBEs, Nag et al. (2010) showed that only a fourth were associated with lightning flashes, and the rest occurred in isolation.However, not all lightning flashes are initiated by NBEs.The rest are initiated by an impulsive discharge with a VHF signature of duration ∼1 μs (Marshall et al., 2014(Marshall et al., , 2019;;Lyu et al., 2019).NBEs are accompanied by strong VHF radiation and are recognized as the strongest natural source of terrestrial VHF radiation (Le Vine, 1980;Willett et al., 1989).The majority of NBEs are of positive polarity (positive charge moving downwards) and occur between the main negative and upper positive charge regions of a thundercloud, while negative polarity NBEs occur between the upper positive and negative screening charge layers (Rison et al., 1999;da Silva & Pasko, 2015;Tilles et al., 2019).NBEs are believed to be caused by a type of discharge known as fast breakdown which is a system of propagating streamers or a discharge wave driven by streamers (Rison et al., 2016;Tilles et al., 2019;Huang et al., 2021).This view of a system of streamers is also supported in (Soler et al., 2020;Liu et al., 2019;Liu & Dwyer, 2020).Because of the short duration of NBEs, the assumed length of the discharge channel is several hundred meters or less (Watson & Marshall, 2007;Dwyer & Uman, 2014;Rison et al., 2016).It has been suggested recently that relativistic runaway discharges driven by photoelectric feedback can be sources of NBEs on these compact scales (Pasko et al., 2023).Significantly longer channel lengths were estimated in (Eack, 2004;da Silva & Pasko, 2015).In addition to ground based observations, NBEs have also been observed by satellites, and in this context are referred to as transionospheric pulse pairs (TIPPs) since they are observed as distinct double pulses with the second pulse being the reflection from the conducting ground (Massey & Holden, 1995;Light & Jacobson, 2002).It has been recently proposed that orbital observations of 10s of μs blue flashes originating from thunderclouds represent an optical equivalent of radio waves observations of NBEs (Neubert et al., 2021;Li et al., 2023, and references therein).
A class of models known as transmission line (TL) or modified transmission line (MTL) models are commonly used to model the source current for NBEs and IBPs (Shao & Heavner, 2006;Watson & Marshall, 2007;Nag et al., 2010;Karunarathne et al., 2014;Rison et al., 2016).In addition, da Silva and Pasko (2015) developed a numerical model to explain the initial stage of the lightning leader tree development, and demonstrated that these pulses can be explained by the stepping of a negative lightning leader.To the best of our knowledge, the most recent work on remote sensing of currents in compact intra-cloud discharges is due to Karunarathne et al. (2019) where the authors develop a matrix inversion method to extract source current of IBPs from far field electric field measurements.This method is illustrated in further studies of IBPs (Karunarathne et al., 2020(Karunarathne et al., , 2021)).Once the source current is obtained from any of these models, the radiated electric or magnetic field pulses can be readily calculated.We advocate here the use of a frequency domain solution of Maxwell's equations which can be found in standard electromagnetic textbooks (for example, Balanis, 2012) in the form of simple equations.This solution is computationally less expensive than the commonly used time domain solution.Using a simple deconvolution process one may extract the source current from electric or magnetic field data using the response of the sourcepropagation-observer system to an impulsive current source.Unlike TL or MTL models, which require the optimization of several model parameters, such as channel length, peak current, current pulse width, propagation speed of the current wave, and current attenuation factor, this method only requires the length of the discharge channel, which can often be directly inferred when VHF source data accompanies the measurements of radiated electric or magnetic fields (for example, Rison et al., 2016;Liu et al., 2022).

Model Formulation
We model the source of NBEs as a straight vertical channel of length l located at an altitude h, spanning h l/2 to h + l/2, above a perfectly conducting ground plane.The point of observation is placed on the ground at a horizontal distance d away from the source.The modeling problem is then reduced to determining the current I (z, t) flowing in the channel that can best describe the NBE under consideration.Once this current is known, the electric ⃗ E or magnetic ⃗ B fields are readily calculated.Note that for an observer on the perfectly conducting ground, the boundary conditions necessitate that only the normal component of ⃗ E and the tangential component of ⃗ B survive.In our case these are E z and B ϕ components, where z is the vertical coordinate, and ϕ is the azimuthal coordinate of the cylindrical coordinate system (r, ϕ, z).The conducting ground plane can be replaced by an image source beneath the plane, carrying current in the same direction as the original source (Balanis, 2012, p. 315).This simply doubles the magnitude of the fields on the ground.
The TL or MTL models were originally developed to model the lightning return stroke.In these models a current pulse is injected at one end of the channel, for example, the bottom end at h 1 = h l/2, which then travels upwards along the channel with a propagation velocity v.The current at any time t and at location z on the channel is given by where the function f(z) specifies the spatial distribution (or attenuation) of the current pulse along the channel.Different transmission line models are distinguished by different f(z).Watson and Marshall (2007) and Rison et al. (2016) used MTL models to best reproduce E-field change measurements for NBEs.Similarly, Karunarathne et al. ( 2014) used several types of MTL models to model E-field measurements of IBPs.A detailed review of these models can be found in (Rakov & Uman, 2003, p. 400).

Field Equations
The classic equations given by Uman et al. (1975) are the most commonly used method to calculate the electric and magnetic fields for a given source current I (z, t).The electric field, as a function of time t is given by where ϵ 0 is the permittivity of free space, c is the speed of light, R z′) is the distance of an element dz′ along the channel from the point of observation, and t r = t R/c (or τ r = τ R/c) is the retarded time.This equation is referred to as Uman's equation in this work.The corresponding equation for the magnetic field B ⟶ can also be found in (Uman et al., 1975).These equations are a special case of the more general Jefimenko's equations (see e.g., Griffiths, 2013, p. 449).The three terms on the right hand side of Equation 2 from left to right are known as the electrostatic, induction, and radiation terms, respectively.
For the purpose of this work, we use the frequency domain solution of the wave equation presented in (Balanis, 2012, p. 284).For our problem geometry, the electric field spectrum is given by where η = ̅̅̅̅̅̅̅̅̅̅ μ 0 /ϵ 0 √ is the characteristic impedance of free space, β = ω/c is the phase constant (also called wavenumber), I(z, ω) is the spectrum of I (z, t) in the frequency domain, that is, I(z,ω) = F[I(z, t)], where F denotes Fourier transform, and , and Equation 3 is referred to as Balanis' equation in this work.It can be shown that Equations 2 and 3 provide solutions for the field that are exactly equivalent (see Section 1 in Supporting Information S1).
More generally, for any arbitrary current density ⃗ J in free space where ⃗ R = ⃗r ⃗r′, and ⃗r and ⃗r′ are the locations of the observer and source, respectively.In our cylindrical (r, ϕ, z) coordinates ⃗r = r îr and ⃗r′ = z′ îz , where îr and îz are unit vectors along r and z coordinates, respectively.Once ⃗ E( ⃗r,ω) is obtained, its inverse Fourier transform gives the time domain electric field, that is, ⃗r,ω) ] .The corresponding equation for the magnetic field can be found in (Balanis, 2012, p. 285).
Equation 3 is applicable to an arbitrary current source I (z, t), which is not necessarily in the form of Equation 1.In a discrete domain, F and F 1 are easily implemented using standard functions provided for most programming languages, that is, fft() and ifft(), respectively, in MATLAB.Further, Equation 3 is computationally less expensive than Equation 2 since the electrostatic term in Equation 2 involves an additional integration in time, and the radiation term includes an additional differentiation, analytical expressions for which may either not be present or difficult to obtain.In a computational domain with a spatial grid with n points and n samples in time, Equation 2 requires the evaluation of n 2 integrations and n differentiations, while Equation 3 requires the evaluation of only n integrations.
A similar frequency domain treatment is presented in (LeVine & Meneghini, 1976;Le Vine & Meneghini, 1978) for an arbitrarily oriented linear current source driven by a traveling current wave and located above a perfectly conducting ground plane.It should be noted that LeVine and Meneghini (1976) simply used their frequency domain formulation to obtain and study field spectrums, without using these to calculate the time domain field.Where time domain fields needed to be calculated for a source current in the form of a traveling current wave I(t z/v), Le Vine and Meneghini (1978) used analytical time domain expressions.Here the time domain result is obtained by an inverse Fourier transform only after E z (ω) is calculated in frequency domain using Equation 3.
We note that due to the β in the denominator, Equation 3 results in a singularity for ω = 0. Two approaches to calculate E z (ω = 0) are given in Section 2 in Supporting Information S1 and are shown to be equivalent.For any current I (z, t), we have In practical implementations, to ensure that Equation 6does not diverge, we need to add the same current pulse as I(z, t), but with reversed polarity placed at some time after the original pulse, that is, I(z, t T ).Care needs to be taken such that there is enough separation T between the original pulse and the reverse polarity pulse so that there is no overlap between them.An example is schematically shown in the upper half panel of Figure 1a.Since the original pulse I(z, t) moves some charge Q along the channel, and I(z, t T ) moves the same amount of charge but with opposite polarity, that is, Q, the net charge transfer is zero at some late time 2T.The reason behind this step is detailed in Section 2 in Supporting Information S1.Lastly, since I(z, t) is real and causal, only one half of its spectrum is enough to specify the entire spectrum.Hence, we obtain an analytic signal for I(z, t) using the Hilbert transform, for which the negative half of its spectrum vanishes (Oppenheim & Schafer, 1999, p. 775).Following this we can compute the Fourier transform to obtain I(z, ω) and use it in Equation 3.

Extraction of Current From Electric Field Measurements
The first transmission line model with a current pulse propagating with speed v was given by Uman and McLain (1969), who adapted it from an earlier model for lightning return stroke current by Dennis and Pierce (1964).The propagation speed v was first introduced to these models by Dennis and Pierce (1964) who defined it as velocity of charge moving in the channel.We want to emphasize that v should only be interpreted as the velocity of propagating current wave, and not as the velocity of moving charge.Charges are created and neutralized along the channel, and the propagation is simply the propagation of phase.The front of the return stroke is itself considered to be propagating with speed v f .Depending on the model, v may or may not be equal to v f .The model by Dennis and Pierce (1964) was itself adapted from the return stroke model of Bruce and Golde (1941), which is the earliest of such models.This form of a propagating current pulse, with the additional parameter v, may be suited to model return strokes, but does not have any physical justification for situations of fixed length channels/regions pertinent to isolated NBE and IBP sources.In this work we develop a new model in which the current I (z, t) has the form where i(t) is the current waveform, and g(z) specifies its spatial distribution on the channel.This eliminates the need for a traveling current pulse with propagation speed v.The current at any moment in time and location on the channel is completely determined by Equation 7. da Silva and Pasko (2015, Figure 6) show the current distribution along the channel at several instances of time for two different models of an IBP: the transmission line model of Karunarathne et al. (2014), and the physics-based model of da Silva and Pasko (2015).For both models the figure shows that the current approaches zero at the channel ends, and the peak current occurs close to the middle of the channel.A zero current at all times at the channel extremities also ensures overall charge neutrality of the channel.Consequently, to model any bipolar lightning ⃗ E or ⃗ B pulse, a half-wave sinusoidal distribution of the current along the source channel can be used as a simple model that captures principal physics based distribution including zero current values at both ends of the channel and maximum of the current in the middle of the channel (da Silva & Pasko, 2015).The employed current model is expected to be valid in situations when isolated elongated conductor is polarized upon application of an external electric field (for example, Pasko et al., 1998;Pasko, 2014;da Silva & Pasko, 2015;Cummer, 2020).We emphasize that our approach does not specify the physical mechanism that leads to creation of the channel/region with enhanced conductivity.We define a fixed g (z) for modeling all pulses considered here as: It should be highlighted that a uniform current distribution along the channel, that is, g(z) = 1, would return essentially the same current waveform as given by the distribution in Equation 8, only scaled by a factor of 2/π.The half-wave cosine distribution of Equation 8 is chosen here for reasons mentioned above, and because it results in a physically realizable charge distribution along the channel.A simple application of the continuity equation to the current model given by Equations 7 and 8 will result in a continuous linear charge density ρ l (z) on the channel that is proportional to ∂g(z)/∂z, that is, ρ l (z) ∝ sin (π(z h)/l), for an arbitrary i(t).For the g(z) = 1 case, ρ l (z) is expressed using two delta functions of opposite signs at the channel ends, h l/2 and h + l/2.In both cases the conducting channel is polarized similar to a dipole, with one half of the channel being positively charged, and the other half negatively charged, and the E-field waveform is completely defined by the time dependent dipole and current moments of the polarized conducting channel.
Similar to the technique developed by Cummer and Inan (2000) to extract the current from the ELF sferics of lightning discharges using a deconvolution method, here we use E-field waveforms measured at an arbitrary but known distance from the source to determine the NBE current i(t).Cummer and Inan (2000) assume a current that only varies with time and is constant along the channel length at any instance of time, that is, g(z) = 1.Instead of i(t), they equivalently determine the current moment M I = i(t)l (in the present case, M I = (2/π)i(t)l).Our choice of g(z) given by Equation 8gives a more realistic estimation of the current based on the discussion in the preceding paragraph.
Taking the Fourier transform of Equation 7and substituting resultant current I (z, ω) = g(z)i(ω) in Equation 3, we arrive at where E IR z (ω) can be obtained from Equation 3 by setting i(ω) = 1 corresponding to i(t) = δ(t), where δ(t) is the Dirac delta function.Hence, the source-propagation-observer system is linear time-invariant with E IR z (ω) being the impulse frequency response.If E z (t) is known, then the current i(t) can be recovered using deconvolution, that is, ]. (10) Note that in the limit lim ω→0 E IR z (ω = 0) = ∞, and hence, the DC current component i (ω = 0) = 0, which is as expected for a physical i(t) representing a current pulse.
Similar to the modification to current described earlier and depicted in Figure 1a to calculate the electric field using Equation 3, in the inverse process of calculating the current using Equation 10, the E-field waveform needs to be modified in the manner depicted in the lower half panel of Figure 1a.To ensure that F[E z (t)] does not diverge, E z (t) needs to be appended with E z (t T).This ensures that the field value is zero at t = 2T.The example shown in Figure 1a forms an input-output pair for Equation 3 and vice-versa for Equation 10.

Results and Discussion
We use the transmission line model generated data to demonstrate application of our proposed method to NBEs reported in (Eack, 2004).We use actual data of their NBE1 and NBE3 reported in (Rison et al., 2016).For Eack's NBE, the modified transmission line exponential increasing (MTLEI) model was developed by Watson and Marshall (2007), where f(z) is an exponential factor as in MTLE, with exponential growth of the current pulse rather than exponential decay.As opposed to these MTL models, which require optimizing several model parameters, the method developed in this work requires no such process, and only needs the height of channel center h, length l, and horizontal distance to observer d as model inputs.In (Rison et al., 2016) the source location can be inferred from location of VHF sources.
Eack's NBE was measured both in the near field (d = 2.8 km from source) and far field (d = 200 km from source).Since it is the same NBE observed at different locations, as was expected, the current recovered from both the cases are nearly identical, as can be seen in Figure 1b.The E-field waveforms, for both near and far field, used to obtain the current are shown in panels (c) and (d) of Figure 1, respectively.The electric fields obtained using Balanis' equation are identical to Uman's equation for both the near and far field cases.Eack (2004) noted that, for this NBE, the initial discharge occurred at an altitude of h = 7 km.They calculated the vertical extent of the discharge by multiplying the NBE pulse width by the discharge propagation velocity and arrived at an estimate of 4 km.Recent interferometer observations of VHF sources estimate vertical extents of NBE discharges to be between 100 m and 1 km such that prior presence of 4 km long conducting channel might be unlikely.da Silva and Pasko (2015) assumed a pre-existing channel of 800 m, and a step length of another 800 m for the negative leader stepping.This 1.6 km length also seems larger than the possible range mentioned above.Hence, we use a channel length l = 600 m, that is the length estimated by the MTLEI model of Watson and Marshall (2007), and that is more inline with recent measurements of NBEs at close range (for example, Rison et al., 2016).Since Rison's NBEs were measured in the near field, where the electrostatic and induction terms had significant contributions, they form an excellent test scenario.In their modified transmission line with an exponential attenuation factor f(z) (MTLE) model, Rison et al. (2016) used the first two terms of Uman's equation for the electrostatic and induction terms, and Equation 11 from (Shao et al., 2005) for the radiation term.The (Shao et al., 2005) equation is for the case of source in free space.When a perfectly conducting ground is present, it needs to be modified to account for the ground.To do this we follow an analysis similar to (Krider, 1992;Shao et al., 2004).Here, we refer to the modified version as Shao's equation.This gives a result different than when Uman's equation is used for all the three E-field terms.
Figure 2 shows the results for the E-field waveform for both NBE1 and NBE3 using the MTLE model with Uman's equation, Balanis' equation, and Shao's equation.Additionally, the measured waveforms are also shown.The measured data are openly provided by Li et al. (2022a).All the same parameters were used for the MTLE model as used in (Rison et al., 2016).We note that for both NBEs Uman's and Balanis' equation give identical results as expected.Li et al. (2022a) obtained similar results using Uman's equation and the full wave Finite Difference Time Domain solution in (Li et al., 2020).However, these results do not agree with results from Shao's equation.This is on account of the opposite polarity second peak of the radiation term which is almost completely suppressed by Shao's equation (see Figure 3 of (Li et al., 2022a)).The secondary peak is caused by a discontinuity at the far end of the channel as the current pulse leaves the channel.This discontinuity is ignored in Shao's equation.The difference in results would not arise for a current that decays completely before reaching the far end, as there would be no current discontinuity here.In the derivation of this equation for a TL model with the far end of the channel extending at the same speed as the propagation speed of the current pulse v, Shao et al. (2004) argue that the current at this end will always be zero, since during the time the current pulse takes to travel from the channel base to the far end, the far end will extend further away at the same speed.This seems reasonable for a return stroke where the channel extends to lengths of several kilometres, and the current pulse decays to zero.However, for short channels of lengths a few hundred meters, similar to the present case of NBEs, it is unclear if this assumption can be safely made (Li et al., 2022a).
Additional explanation of Shao's derivation was given in (Shao et al., 2005), where the instantaneous extension of the channel with the propagating front of the current pulse was said to induce a current source at the top of the channel, which essentially cancels out the traveling current discontinuity.Giving this a physical explanation, Li et al. (2022a) developed a rebounding MTLE model, where they consider a secondary current pulse traveling in the direction opposite to the main current pulse.Essentially, once the current reaches the far end, instead of Following Rison et al. (2016) we assume d = 5.5 km and h = 6.0 km for their NBE1, and d = 3.3 km and h = 6.6 km for NBE3, and l = 500 m for both NBEs.Figures 3a and 3c show the current pulse i(t) extracted for NBE1 and NBE3, respectively.The currents were noisy, as expected from the data, and were smoothed using a lowpass filter with normalized passband frequency 0.01π rad/sample.The sampling frequency was ∼16 MHz.The smoothed currents are also shown in panels (a) and (c) for each NBE.These currents can then be used to calculate the E-field using either Uman's or Balanis' equations, or Equation 9.This is shown in panels (b) and (d) of Figure 3 for NBE1 and NBE3, respectively.We further show 10 point moving average of the original data in these panels.The fields calculated using the extracted currents show excellent agreement with observed data, and are also able to reproduce the finer details in the NBE waveforms, which is not possible with TL or MTL models.
A major concern in any deconvolution process is the lack of uniqueness of the solution.From Equation 9 it is clearly seen that zeros in the impulse response spectrum can be the cause of potential non-uniqueness.Cummer and Inan (2000) address this issue by approximating the time domain counterpart of the convolution in Equation 9 as a linear system, and adding a linear regularization term.This regularization term ensures uniqueness and the system can then be solved for the current using any well known linear solver such as LU decomposition or conjugate gradient.In our model we do not explicitly address the non-uniqueness of deconvolution.The only test for the uniqueness of the solution, is to check whether the obtained solution (the extracted current in our case) produces the observed/measured data when convolved with the impulse response (Wertheim, 1975).We performed this test for all NBEs considered here as well as for several other pulses used in our testing.In every single case, the extracted current successfully reproduced the field waveform which was used as input to the model.

Conclusions
We present here a method for the remote sensing of discharge currents producing NBEs, using measured E-field waveforms.The results of this technique were demonstrated for several NBEs.In addition, we introduce the use of a simple full wave solution of Maxwell's equations in the frequency domain to calculate electromagnetic fields from lightning discharges.This is a well established solution which can be found in textbooks, but has never been used to obtain time domain waveforms using frequency domain calculations.This solution further lends itself naturally to the deconvolution technique for current extraction.
applied for the first time to analysis of radiation from NBE source currents • The proposed new approach is validated by detailed comparisons with previously published NBE studies • The methodology is developed that allows exact reconstruction of source currents that radiate the observed wideband NBE waveforms Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.(a) A sample pair of i(t) and E z (t).This i(t) can be used as input to Equation 3 to obtain E z (t).Similarly, this E z (t) when used as input to Equation 10 will return i(t) as shown.(b) Source currents extracted for Eack's NBE using both near field (d = 2.8 km) and far field (d = 200 km) E-field waveforms.Also, h = 7 km, and l = 600 m.(c) Eack's NBE E-field waveform for the near field and (d) far field pulses calculated using Uman's (solid line) and Balanis' (open circles) equations.

Figure 2 .
Figure 2. NBE E-field waveforms calculated using Uman's, Shao's, and Balanis' equations for (Rison et al., 2016)'s (a) NBE1 and (b) NBE3.The current profile, as well as all other model parameters, obtained from their MTLE model were used.Rison's results, which used Shao's equation, are successfully reproduced.

Figure 3 .
Figure 3. (a) NBE source current extracted for Rison's NBE1 (h = 6.0 km, l = 500 m, and d = 5.5 km), and (b) the resulting Efield produced from this current with the goal of model validation.(c) NBE source current extracted for Rison's NBE3 (h = 6.6 km, l = 500 m, and d = 3.3 km), and (d) the resulting E-field produced from this current with the goal of model validation.