Lightning Interferometry Uncertainty, Beam Steering Interferometry, and Evidence of Lightning Being Ignited by a Cosmic Ray Shower

We present an uncertainty analysis on correlation‐based time delay estimate, the basis for broadband lightning interferometry. A normal interferometry would yield much higher uncertainties than the theoretically predicted lower bound due to finite correlation window lengths. However, if the signals are aligned in time, the uncertainties approach the lower bound and can be used to indicate the interferometric uncertainties. Based on this, we introduce a beam steering interferometry technique. It first estimates a direction centroid for the lightning sources with a normal interferometry and computes the time delays among the signals. It then shifts the raw data with the time delays to align the signals roughly in time and reprocess the aligned signals. It finally pushes the reprocessed results to their correct positions, based on the time delays estimated in the first step. We apply this technique on a fast positive breakdown (FPB) process that started a normal intracloud lightning. The FPB process is shown to have a much more complex structure and development than a normal interferometry would provide. More importantly, from both interferometric and polarization analyses, the FPB appears to be ignited by a cosmic ray shower (CRS). We estimate the radio frequency strength and frequency content related to a presumed CRS in a thunderstorm electric field and find that they are in agreement with the observations. We examine the electric field effect of the CRS front and find that it could raise the field above the threshold for positive breakdown and is capable of igniting the FPB (and the lightning).


Introduction
Since its introduction to lightning observations , broadband radio frequency (RF) interferometry has become a widely used technique in the research community for detailed lightning studies Dong et al., 2001;Lyu et al., 2019;Rison et al., 2016;Shao et al., 2018;Stock et al., 2014Stock et al., , 2017Sun et al., 2013;Tilles et al., 2019). In the last two decades, this technique has been continuously developed and improved at a number of research institutes in the areas of broadband data acquisition and interferometric data processing. Broadband interferometers digitize and record raw broadband signals from spatially separated receiving antennas and postprocess the digital signals to form lightning maps. Although a large amount of data has to be stored for each lightning flash, one has the flexibility to process the raw data in a variety of ways to achieve the best possible results. As a result, this technique has been contributing to new and greater understandings of lightning breakdown processes, noticeably the observations of lightning-initiating fast positive breakdown (FPB) processes (e.g., Lyu et al., 2019;Rison et al., 2016;Shao et al., 2018). FPBs propagate at a high speed of several times 10 7 m/s and effectively carry positive charge in the direction of propagation. Similar to FPBs but opposite in polarity and propagation direction, Tilles et al. (2019) reported with broadband interferometer observations that some intracloud (IC) flashes started with fast negative breakdown (FNB).
For increasingly more detailed spatial and temporal lightning analysis with broadband interferometer observations, it is critical to understand the mapping uncertainties of such a technique. Stock et al. (2014Stock et al. ( , 2017 presented an initial uncertainty analysis for their interferometer observations based on Carter's (1987) lower bound uncertainty formula on coherent time delay estimates. However, Carter's theoretical prediction is based on assumptions of stationary random signals and on sufficiently long correlation time windows that are often not satisfied in lightning interferometric processing. In this paper, we conduct a series of numerical simulations to show that a normal interferometric process that uses cross correlation to determine the time delay between a pair of signals will yield much larger uncertainties than the theoretical lower bound. This is especially true for lightning analysis in which shorter correlation windows are always desired for high time resolutions. Nevertheless, when the pair of raw signals is aligned in time, the uncertainties approach the lower bound and the lower bound can be used to indicate the interferometric mapping uncertainties.
Based on the uncertainty analysis, we introduce a beam steering interferometric technique to improve the broadband interferometric mapping capability and to more reliably estimate the mapping uncertainties/ errors for each lightning source. Such a technique processes the same raw data with two iterations: (1) finding the rough direction of the lightning source with a normal interferometric process for a selected time interval and (2) aligning in time the raw data from different antennas based on the rough source direction and reprocessing the data to produce a more accurate and detailed interferometric image. It is demonstrated in this paper that such a two-step process can significantly improve the map quality with the same raw data.
We then apply the beam steering interferometry technique to a FPB process that initiated a normal polarity IC flash. FPB is observed to occur in the initial microseconds of both IC and cloud-to-ground lightning (CG) discharges (Rison et al., 2016;Shao et al., 2018) and are considered to initiate the lightning flashes. With the beam steering interferometric analysis, it is found that such positive breakdown can be much more complex than they originally appear. We found that the FPB presented in this study shows two apparent groups of breakdown processes that propagate downward slightly differently in time and follow two apparent paths. More importantly, the two groups appear to start in front of an arc-shaped source front that resembles the front of a cosmic ray shower (CRS), as will be discussed in detail in this paper.
In addition to the interferometric analysis, we analyzed the polarization state of the corresponding RF signals. The polarization observations further support the development of the interferometric sources and indicate that the sources for the arc-shaped front are likely related to a CRS front.
To examine the CRS hypothesis, we calculate, to order of magnitude, the RF signal strength and frequency content related to a presumed CRS process that enters a region of electric field in the storm. With an assumed field strength of 110 kV/m at 8-km altitude, the field strength for relativistic runaway electron avalanches (RREA), we find that the estimated RF magnitude and frequencies are in agreement with that observed for the arc-shaped feature. We further estimate the electric field effect due to the presumed CRS front and find that the field can be enhanced by three times or more of the ambient field. Such an enhanced field is high enough to ignite conventional positive breakdown that leads to a normal IC flash.

Uncertainties in Broadband Lightning Interferometry
For a broadband interferometer, the direction of the RF signal can be obtained by evaluating either time delays or phase differences between pairs of receiving signals. With a baseline d formed by two separated antennas, the 1-dimensional (1-D) angular direction can be determined by cos α ¼ cτ d Here, α is the angle between the source direction and the baseline; c is the speed of light; f is signal frequency; τ is the time delay, and ΔΦ is the phase difference measured at f. In a normal broadband interferometric process, a common and synchronized time window is chosen for all pairs of digitized signals from which a corresponding τ or ΔΦ is computed. For the chosen time window, the resultant cosα for each baseline is computed based on τ or ΔΦ and the source's 2-D direction is determined by combining the different baseline measurements (e.g., Shao et al., 2018).

Uncertainties for Time Delay Estimates
The Los Alamos broadband interferometric mapping and polarization observation system (BIMAP) detects the source direction by computing the time delay τ between pairs of receiving antennas. The time delay is determined by examining the cross correlation of the received signals for each chosen, finite time window as the window slides forward along the data series, as described in Shao et al. (2018). To interpret the detailed lightning processes obtained by BIMAP, it is necessary to understand the mapping uncertainties in the observations. Fortunately, studies for similar direction-finding problems had been conducted decades ago in the acoustic community (e.g., Carter, 1987;Hahn, 1975) and has recently been used in Stock et al. (2014) for broadband lightning interferometric observations. A pair of time series data from a broadband interferometer can be assumed to be where n 1 (t) and n 2 (t) are noises and s(t) is the signal. If the noises and the signal are stationary random processes and are uncorrelated with each other, the lower bound for the uncertainty σ τ on time delay estimate τ can be expressed as (e.g., Hahn (1975)), based on the Cramer-Rao Lower Bond (CRLB) theory, where T is observation period; f is frequency; N 1 ( f ), N 2 ( f ), and S( f ) are power spectra of the noises and the signal; and f 1 and f 2 are the lower and higher cutoff frequencies of the receivers. In a simplified case where the signal and noises power spectra are constant across the frequencies and the two noise sequences have the same power, Equation 2 becomes where SNR = S/N is the signal-to-noise ratio. This simplified formula has been commonly used in the research communities to estimate the best possible performance for the time delay analysis (e.g., Carter, 1987;Stock et al., 2014).
Nevertheless, Equations 2 and 3 are only accurate when T is large, for example, hundreds times of 1/( f 2 f 1 ), and with an idealized square-shaped filter between f 1 and f 2 . In practice, the latter can never be achieved with realistic filters, and for the impulsive nature of lightning sources and the desire of high time resolution analysis, the former is often not satisfied. To obtain a more realistic lower bound σ τ with BIMAP's specific receiver response and data processing, numerical simulations have been conducted.
In the simulation, two independent Gaussian random series of σ 2 n ¼ 1 are generated to represent the noises, and one Gaussian random series with selected σ 2 s values is generated for the signal. The sample rate of the signal and the noises are assumed to be 190 MHz to represent BIMAP's digitizing rate. The signal is then added to the two noises respectively, and the second series is shifted with a chosen number of data points to represent the time delay. A steep 80-MHz low-pass and a soft 20-MHz high-pass digital filter are applied on the data series to simulate BIMAP's receiver and antenna responses. Time delay between the two filtered series is computed with cross-correlation analysis, and a fine, subsample time delay is interpolated with a parabolic fit of 3 points around the cross-correlation peak. It is noted that linear upsampling of the original time series data, as was done by BIMAP previously, does not significantly and monotonically improve the time delay estimate, especially when the data are upsampled by more than a factor of 2. This is understandable because the linear upsampling does not add any more information to the original data. In this study and future BIMAP processing, we will not upsample the original data by more than factor of 2.
Each run of the simulation described above gives a single estimate of the time delay b τ. To investigate the statistical uncertainty of the time delay estimate, we repeat such a simulation 1,000 times with independently generated noises and signals and compute the standard deviation b σ τ of the 1,000 simulatedb τ s. The solid lines in Figure 1a show the simulated b σ τ as function of input SNR σ 2 s =σ 2 n À Á for four selected time window sizes: 64, 128, 256, and 2,048 BIMAP-equivalent data points (T: 0.337, 0.674, 1.347, and 10.779 μs), from top to bottom, respectively. In these simulations, the input time delay was set to be zero. Figure 1a represent the corresponding theoretical lower bound σ τ for the four window sizes. The lower bound in each simulation run is computed with Equation 2, and a final σ τ for each given window size and SNR are obtained with the mean of the corresponding σ τ s of the 1,000 runs. Equation 2 instead of 3 is used for the lower bound calculation because the time window under consideration is not long enough to ensure a flat power spectrum for the noises and the signal or to ensure the same spectra for the two noises in each individual runs. Since the applied filter of 20-80 MHz do not cutoff all the signals and noises beyond the 3-dB passband, the integrations in Equation 2 are carried out from 0 to the Nyquist frequency of 95 MHz to include all the signal and noise power. Figure 1a shows that above a certain SNR, the simulated uncertainties b σ τ (solid lines) of the time delay estimate for each selected window size are in good agreement with the predicted lower bound σ τ (dashed lines), and the theoretical lower bound σ τ can be used to indicate the errors in the time delay estimate. On the other hand, below certain SNRs the simulated b σ τ increases sharply far above the predicted lower bound for each of the selected window sizes. This transition feature at lower SNR is in agreement with a more rigorous theoretical prediction of Ianniello (1982). For broadband lightning interferometry analysis, time delay estimate in this lower SNR region is apparently not trustworthy and the interferometric results should be ignored.

The dashed lines in
Nevertheless, it needs to be emphasized that the results in Figure 1a are based on zero time delay in the two data series. In a realistic observation with a normal interferometric processing, the signal in the two series will have certain time delays. For instance, with BIMAP's 19-m baselines and 190 MHz sampling rate, one can have a maximum of 12 data points delay between two received signals. Statistically, if the signal is stationary and the source is not moving, one can estimate the time delay with a long enough time window to achieve the lower bound uncertainty. Unfortunately, neither of those conditions is typically applicable to lightning observations. Figure 1b shows the effect of the signal time delays on the uncertainty estimate for a selected window size of 64 data points (T: 337 ns). In the figure, the line at the bottom represents that of the zero delay, the same as the black solid line in Figure 1a. The lines from lower to upper represent the uncertainties for time delays of 2, 4, 8, and 12 BIMAP-equivalent data points. As evident in the plot, the uncertainties increase with increasing signal delays across all the SNR values. This is due to the shortening of the common portion of the signal in the selected time window between the two data series. At larger SNR values, each of the lines becomes flat and shows little improvement with increasing SNR. This is because that with larger SNR the cross correlation is increasingly dominated by the signal itself and the uncertainty is dominated by the misalignment of the signal in the two series. Figure 1c is the same as Figure 1b but for a longer window size of 256 data points. As expected, it shows similar features as that in Figure 1b. However, the corresponding uncertainties are lowered proportionally due to the increased time window. Figures 1b and 1c show that for normal interferometric lightning analysis, the uncertainties for the time delay estimate could be much greater than the corresponding theoretical lower bound and they do not necessarily follow the approximate 1/SNR trend if the time window is not long enough.
So far, we have not used the simplified Equation 3 for lower bound uncertainty analysis, partly because it is not clear how to choose the proper values for f 1 and f 2 with the realistic filters. However, we can empirically determine the effective f 1 and f 2 by using the results based on the more general Equation 2. Figure 1d shows that by choosing f 1 = 10 MHz and f 2 = 82 MHz for the BIMAP antenna and receiver response of 20-to 80-MHz 3-dB passband, the lower bound uncertainties from Equation 3 are in good agreement with that of Equation 2. The upper lines (black for Equation 2, red for Equation 3), where the red line overlaps the black line, are for time window of 64 data points and the lower lines, which are also overlapped, for time window of 256 points. In practice, since the signal is mixed with the noise in a random manner, it is difficult to extract precisely the signal from the data series, and it is difficult to apply Equation 2 for the lower bound analysis since it requires explicit knowledge of the signal and the noise. Nevertheless, with a proper determination of the effective f 1 and f 2 , Equation 3 can be used in practical data analysis to estimate the theoretical lower bound. The signal power for a selected time window is represented by the peak value of the cross correlation between the two data series. The noise power can be determined with a data section that contains no lightning signals. Therefore, SNR can be reliably evaluated for Equation 3.

Angular Uncertainties
BIMAP consists of two orthogonal baselines of d = 19 m that are aligned in East-West (X) and North-South (Y) directions. Based on the fundamental relation of cos α ¼ cτ d for each baseline, the uncertainties of the source position projected on the horizontal XY plane for BIMAP can be expressed as where σ τx and σ τy are the uncertainties for the time delay estimate between the two pairs of antennas in the EW and NS baselines. After some simple geometrical analyses, the uncertainties for the azimuth and elevation angles can be expressed as where τ x and τ y are measured time delays for the EW and NS baseline measurements. In general, σ τ depends on τ for each baseline, as shown in Figures 1b and 1c and may not be assumed the same between the two baseline measurements. However, if they are assumed the same, Equation 5 can be simplified to Equations 11 and 12 in .

Beam Steering Interferometry
The analysis in last section shows that if the time delay between signals from a pair of antennas approaches zero, the uncertainty of the time delay estimate approaches the theoretical lower bound for the time window sizes that are often used in interferometric lightning analysis. However, lightning sources rarely occur directly above an interferometer and there is almost always a time delay between a pair of signals. The delay is greater for sources at lower elevation angles and for interferometer with longer baselines.
Figures 2a-2c show realistic 1-μs lightning RF data measured by the three BIMPA antennas. BIMAP consists of three antennas that form approximately a right triangle with the antennas separated by 19-m along North-South (NS) and East-West (EW) legs to form three baselines. The antenna at the right angle is considered the reference antenna, and the other two antennas that are placed toward the north and the west are considered outer antennas. Signals from the three antennas are time synchronized and digitized at a rate of 190 megasamples per second. Figures 2a-2c compare the signals for the three pairs of the antennas, north (red) and reference (black), west (blue) and reference, and north and west antennas, respectively. BIMAP evaluates the time delay for each of the three pairs by cross-correlating the corresponding signals over a chosen time window (e.g., the shaded boxes) and determines the 2-D source direction using the three time delays. The time window is then shifted consecutively forward to determine source directions at later times. If the time window encloses the RF signals for a temporally isolated and intense source, accurate and consistent time delays can be obtained among the three pairs, and a reliable source direction can be determined. However, the sliding window often encloses slightly different signals among the three pairs due to their different arrival times at the three antennas, as illustrated with "window n." Near the left edge of "window n," a distinct pulse (in the oval) is enclosed for the reference (black) and west (blue) antennas but not for the north antenna (red). In this case, the cross correlations and the resultant time delays for the three pairs are not applied to the same signals, and the correlation coefficients will be lower for the first and the last pairs and the subsequent source direction will be less accurate or erroneous. For each individual pair, the greater the time delay between the two signals, the less the common portion of the signals is enclosed in the time window, resulting in a greater uncertainty for the time delay estimate, as previously shown in Figures 1b and 1c.
For a given baseline, the direction error due to this "window edge effect" is the smallest when the source is orthogonal to the baseline (α = π/2, when the signals are aligned in time) and is the largest when the source is in the direction of the baseline (α = 0, when the signals are separated by the maximum amount of time). On the other hand, for a source in a certain direction, a longer time window will reduce the error in the time delay estimate, as comparing Figure 1c against Figure 1b, since the RF signals at the window edges have a relatively smaller contribution to the overall correlation estimate. However, a longer correlation window has a greater chance of enclosing multiple RF sources in the window, and the resultant source direction is more likely to be an estimate of the power-weighted direction among the multiple sources. For instance, during the 1-μs time interval in Figure 2, there appear four or more distinct, temporally separated pulses overriding the more continuous lower level signals. If the correlation window is chosen to be 1 μs or longer, it will result in an averaged source direction for all the pulses and the lower level signals in this window.
In lightning interferometric studies, it is always desired to resolve the sources in the finest time and spatial resolution possible with the available raw data. In most reported studies (e.g., Shao et al., 2018;Stock et al., 2014), the time window is chosen to be around 1 μs to provide a good time resolution and high-quality map for a variety of lightning processes. With a shorter time window, the uncertainties from the window edge effect become more significant and the source positions become less accurate or erroneous.
The left column in Figure 3 shows a normal interferometric analysis with a time window of 337 ns (64 BIAMP data points) for the first 500 μs of a small IC. To compute the time delay for each time window, we upsample the raw data by a factor of 2, apply forward and inverse Fourier transformations, and fit the correlation peak with a parabolic function, as that in Stock et al. (2014Stock et al. ( , 2017 and Shao et al. (2018). The time window was shifted 21 ns (4 data points) forward at each step (see Appendix A for such a selection). Figure 3a shows the time waveforms of the broadband RF signal (black) and the electric field changes (red). Figures 3b and 3c show the time development of the source directions in azimuth (AZ) and elevation (EL) angles, and Figure 3d plots the sources in the AZ-EL format. This IC occurred in the cloud nearly due north to BIMAP, and the time delay between the signals from the north and reference antennas is the largest among the three baselines, as shown in Figure 2. As a result, the window edge effect has a bigger impact in the NS direction.
To reduce the location uncertainty with the same raw data and at the same time resolution, we introduce here a two-stage beam steering interferometry that will minimize the erroneous window edge effect. We first process the data with the normal interferometric analysis and find the centroid position for the sources. With this first-stage source position, we compute the corresponding time delays for the outer antennas as compared to the reference antenna and shift backward the signals with corresponding integer numbers of raw points to best represent the respective time delays. Since the sources in the 500 μs in this case are relatively localized as compared to the sample interval of the discrete raw data (5.3 ns, equivalent to 0.084 of the [−1,1] horizontal projection range), a single centroid position is used for the time shifts. The time delay was estimated to be about −7 sample points (−37 ns) for the north antenna and about −3 sample points (−16 ns) for the west antenna. Figures 2d-2f show the signal pairs after shifting the signals 7 and 3 data points for the outer antennas. We then reprocess the temporally aligned signal pairs with the same window size (337 ns, 64 data points) and the same progression steps (21 ns, 4 data points) as that for Figures 3a-3d, as if the sources were pulled overhead. After this second-stage interferometric processing, we move the sources back to their original positions based on the time shifts used for the signal alignment. The right column in Figure 3 shows the results with such a two-stage interferometric processing. As compared to the left column, it reduces substantially the scattering of the source positions. Some detailed time-evolving structures that appear to be noise in the left column can now be better resolved after the two-stage processing, as evident in the initial 100 μs of the discharge.
The beam steering process ensures that all the baselines observe the same portion of the signal, and for each individual baseline the signal pairs are aligned in time that will minimize the uncertainty for its time delay estimate. Because of the latter, the time delay uncertainties approach closely to the theoretical lower bound when the SNR is above a certain level, as shown in Figure 1a. On the other hand, if one uses the lower bound uncertainty in a normal interferometry analysis, it would falsely underestimate the error, as discussed in last section (Figures 1b and 1c).
With the raw data aligned closely in time (Figures 2d-2f) after the beam steering processing, we can now use the lower bound time delay uncertainties to approximate the angular uncertainties in our analysis. In this analysis, Equation 3 is used to compute the time delay uncertainties for the EW and NS baselines (σ τx , σ τy ) with the empirically determined f 1 and f 2 and with data-based SNR estimate for each correlation time windows, as discussed in last section. With the time delays (τ x , τ y ) for the EW and NS baselines and their associated uncertainties (σ τx , σ τy ), the angular uncertainties are computed with Equation 5. Based on the simulation result in Figure 1a for the time window of 64 data points, interferometric results with corresponding SNR greater than 6 are considered trustworthy and are presented in this analysis. Figures 4a-4c show the time waveforms of the RF signal and the estimated angular uncertainties for AZ and EL for the 500 μs of the IC process and show that the uncertainties are anticorrelated with the signals strength (SNR), as expected. Figure 4d shows the source positions in the AZ-EL format with the corresponding error bars for each of the mapped source points.
It should be noted that after the signals are aligned in time, the uncertainties of σ τx and σ τy are approximately the same, and the angular uncertainties can follow the same relation as that reported in Stock et al. (2014), for example, The two-step interferometry described here is new to lightning interferometric observations since the broadband interferometry was introduced over two decades ago . It improves the mapping quality from processing the original out-of-phase signals to processing the time-shifted, approximately in-phase signals. Nevertheless, time (or phase) shift technique itself has long been used in radio astronomy (Ellingson et al., 2013; Chapter 2 in Thompson et al., 2017) and in space-based terrestrial imaging applications (Shao et al., 2004). Similar techniques have also been used in nowadays commonly available phased-array radars, in which phase shift (equivalent to time shift in time domain) to an array of transmitting/receiving elements is adjusted to actively steer the radar beam to desired directions. For lightning observations, the RF signals could arrival from any directions over the whole sky at a given time instant. The technique presented here needs to first find the general direction of the sources with a normal interferometric process, shift the signals to roughly align them in time (equivalent to steering the beam of the antenna array to the general direction of the sources), and then reprocess the data with a second-stage interferometry. Although a simple and straightforward practice, it is to our surprise that it has not been implemented earlier by lightning researchers.

Interferometric Observations
We now apply the beam steering interferometry on a downward FPB process that lasted about 14 μs and started a normal IC flash over Los Alamos on 6 June 2017. The overall BIMAP result for this flash has been reported in Shao et al. (2018). In this study, we process the initial 100 μs of the flash with the normal interferometry, determine the centroid of the source positions in the 100-μs interval, align the raw data by shifting 3 and −6 raw data points for the EW and NS baselines, and reprocess the aligned data to produce the final source positions. BIMAP's sampling resolution is 5.3 ns, so a single data point corresponds to about 8°AZ and 6°EL intervals for sources at 54°EL, at which the lightning started ( Figure 5). Since the lightning sources in this 100-μs interval only extended a few degrees near the initiation region of 280°AZ and 54°EL ( Figure 5), the same single centroid is enough for proper alignment of the raw data records for the time interval.
In the analysis of this FPB process, the correlation time window was chosen to be 674 ns (128 raw data points) and the correlation window advances 21 ns (4 data points) consecutively. The longer 128 points (as compared to 64 data points in the last section) is chosen to reduce the time delay (and angular) uncertainties for weaker signals, especially for the initial 1-2 μs of the lightning process when the RF signal just emerges from the background noise ( Figure 5a).
Figures 5b and 5c show the time-dependent AZ and EL results with a normal interferometry for the 14-μs FPB process. Figure 5h plots the corresponding source positions in the AZ-EL format. Although noisy, it can be seen that the sources in general extend downward. BIMAP, as other single-station lightning interferometers, does not provide 3-D positions for lightning RF sources but only a 2-D directional map. Therefore, our results are limited to 2-D observations. Nevertheless, it is well known that normal ICs in summer New Mexico storms start between 7 and 10 km above mean sea level (e.g., see the supporting information in Rison et al., 2016) or 5-8 km above BIMAP at Los Alamos. For sources around 54°EL in this case, they are about 6-10 km away from BIMAP. The~1°EL general downward extension of the sources (as indicated in Figure

Journal of Geophysical Research: Atmospheres
Figures 5d, 5e and 5i show the corresponding results after the two-stage interferometric processing. It is immediately clear that such an approach provides a much more detailed image than that of the normal interferometry, and some finer spatial structures and more complicated temporal developments start to emerge. For the beam steering interferometry, as well as for the normal interferometry, we only plot the source positions when their SNR is above 3. Below that, the uncertainties are expected to be far above the theoretical lower bound and cannot be trusted, as indicated by the green lines in Figure 1a.
With the time alignment of the raw data in the beam steering interferometry, it is possible to use the lower bound uncertainties for the time delay estimate to deduce the corresponding angular uncertainties (Equation 5). Figures 5f and 5g show the angular error bars on top of the time-dependent AZ and EL results. Comparing with Figure 5a, larger errors are associated with lower signal strength (lower SNR) in the first 1-to 2-μs interval when the lightning signal barely emerged from the background noise. In later time when the SNR are larger, the error bars are shortened and are not clearly visible in the plots. Figure 5j plots the results in the AZ-EL format, same as that in Figure 5i, but with the AZ-EL errors indicated by the white crosshairs.

Journal of Geophysical Research: Atmospheres
It should be noted that the actual errors for each individual source point would be greater than the lower bound, but the beam steering interferometry is the best process to allow the errors to approach the lower bound, especially for lightning analysis in that the RF sources are impulsive and dynamical and short correlation time windows are required. It should also be noted that with any selected window size it is possible that multiple, spatially separated lightning sources could be occurring simultaneously with respect to the window size. The resultant source position is then a centroid location among the possible simultaneous sources, as illustrated in Appendix B. This suggests that the lightning map obtained with a broadband interferometer is likely a spatially smoothed and more concentrated map as compared to the spatial scatterings of the actual sources. This is why shorter window size is always desired to reduce this possible effect, as long as the associated uncertainties are small enough such that the mapped positions are meaningful and trustworthy.
To examine in more detail the spatial and temporal development of this FPB process, Figure 6 shows the time sequence of the source positions in every 1-μs interval. The overall sources as that in Figure 5i are shown as gray dots as the background, and the sources in a specific 1-μs interval are shown by the red dots with their corresponding 1-σ error bars.
During the first 0.4 μs (Figure 6a), the sources appear to increase in EL as though the sources move upward (from the beginning to −7 μs, Figure 5g). However, at this stage the angular errors (1-σ lower bound errors are shown, actual errors can be larger as discussed earlier) are comparable to the source extension and it cannot be certain if the sources actually move upward. In fact, during the first 1.4 μs (Figures 6a and 6b, up to −6 μs), the sources are widely scattered and the position errors are large due to lower SNR, and we have low confidence on the precise positions for each individual points. Nevertheless, it is conceivable that the sources as a group propagates in general from upper left toward lower right from Figures 6a and 6b during the 1.4 μs interval. Starting from the third time interval of −6 to −5 μs (Figure 6c), the errors become small enough and the mapped positions are considered more trustworthy.
It is interesting to see that in Figure 6c, the sources appear to occur along an arc below and in front of the earlier sources in Figures 6a and 6b. Figure 6c shows that in this time interval, the sources extend about 1°A Z from left to right in a 0.5-μs time period (between the first two vertical lines in Figure 5d in which AZ changes from lower to higher values). A 1°AZ interval at 54°EL corresponds to about 60-100 m in 1-D distance, resulting in about one half to two thirds c of 1-D propagation speed. Figure 6d shows that in the next 1 μs, more sources are mapped along the right half of the arc. In this time interval the sources are mapped to propagate backward from the right to the left at about the same high speed, as indicated by the fast AZ decreasing between the second and the third vertical lines in Figure 5d. The apparent fast, near the speed of light back and forth source movement in Figures 6c and 6d along the arc feature is unlikely due to actual propagation of the lightning sources. Given that interferometer maps the centroid of sources within its correlation time window, it is more likely that in Figures 6c and 6d simultaneous multiple sources occur either at or beyond the two extremes of the arc, and possibly with additional sources along the arc between the two extremes, as illustrated in Appendix B. In this case, the apparent rightward propagation in Figure 6c (with increasing AZ) is likely due to relative intensity increase of the sources toward the right, and the apparent leftward propagation in Figure 6d is likely due to intensity increase at the left. With these considerations and simultaneous polarization observations discussed later, Figures 6c and 6d suggest that a simultaneous, spatially spread source front (similar to a wave front), starting from Figure 6a from the upper left, is the reason for the observed arc feature. In the rest of Figure 6, this apparent source front is depicted with a blue dashed line.
In Figure 6e, 3.4 μs into the beginning of the mapped lightning process, the lightning sources start to propagate predominately downward. For the convenience of description, we call the region below the arc the downward region. We also divide the downward region to the left and right sides with the dashed black line. Figure 6e shows that the sources start at the midpoint below the arc feature and propagate downward along the right side of the downward region. Figure 6f shows that new sources appear back to the top of the downward region, with some sources occur at the left side. In the next 2 μs (Figures 6g and 6h), the sources are mapped mostly to the right side and the sources propagate further downward. Figures 6i and 6j show that more sources occur in the upper region at the right side, while in Figure 6j the left side appears to extend further downward as continues from Figure 6g. In Figures 6k and 6l Figure 6m, the sources reach the bottom of the FPB at both sides. Figure 6n shows that some sources occur above the bottom at both the right and left sides during the time period between 5 and 6 μs, apparently due to the start of simultaneous negative breakdown at the top of the downward region. This is illustrated in Figure 7 that shows the observations in an extended time period of about 50 μs, including 35 μs of upper negative breakdown following the initial FPB. As shown in Figure 7c, starting from the time of Figure 6n (5-6 μs, shaded) to 9 μs the sources start to transition from the lowest EL (bottom of FPB) to the upper EL above the arc feature (Figure 7d), with a couple of up and down source oscillations. After the 9-μs mark, sources are all mapped above the arc feature (Figure 7c and 7d) and the sources extend slowly (on the order of 10 5 m/s) upward as normal negative breakdown, as shown in Figure 7c with the slow upward trend after the 9-μs mark. Normal negative breakdown is associated with breakdown in virgin air at an average speed of several times 10 5 m/s that effectively carries negative charge in its propagation direction, as reported extensively in all lightning interferometric observations, either with narrowband (Rhodes et al., 1994;Shao et al., 1995;Shao & Krehbiel, 1996) or broadband systems (e.g., Rison et al., 2016;Shao et al., 2018;Stock et al., 2014) and in time-of-arrival observations (e.g., Rison et al., 2000).
Although complex, the FPB presented in this paper nevertheless extended downward in general. The downward propagations of the left and the right sides appear to be asynchronous in time, suggesting that two different groups of FPB develop separately from the bottom edge of the arc feature. At the right side, the FPB reaches near the bottom in Figure 6h, while at the left side it almost reaches the bottom in Figure 6i, delayed about 4 μs. If we consider the right side group from Figures 6e-6h, it extends about 1°EL in about 4 μs, giving Comparing to the complex spatial structure and temporal development presented here, we notice that FPBs and FNBs in ICs reported in Rison et al. (2016) and Tilles et al. (2019) are mostly simply linear in structure and monotonic in propagation. Their observations show that even if there were multiple sources in a chosen time window (0.7-1.4 μs for the former and 0.7 μs for the latter, with some extent of averaging among the windows), the dominant breakdown sources propagate monotonically downward for FPBs and upward for FNBs in their ICs. Their observed source propagation together with consistent electric field change analysis shows that the inferred breakdown developments are physically real and are not a result of false Decompose of (f) in three time intervals. Time development in each of the decomposed plots is indicated by colors from back to red.

Journal of Geophysical Research: Atmospheres
interferometric imaging of possible multiple sources. Nevertheless, event NBE2 in Rison et al. (2016) does show a more complex development that depicted "several attempted downward events" in their supporting information Figure S5, similar to the asynchronous developments of the breakdown processes between the right and left paths in our Figures 6e and 6f.
In Figure 6o, we attempt to depict and summarize the development of this FPB process. In the first 1-2 μs the lightning sources propagate from upper left toward lower right. The sources are widely scattered due to their lower signal intensity. At about 2-3 μs, a source front of an arc-shaped feature is observed that extended mostly horizontally. At about 3-4 μs, downward positive breakdown started at the bottom of the arc and developed more vertically downward. In this process, there appear to be two groups of FPB that both started at the bottom of the arc-shaped source front and propagated downward slightly asynchronously in time.

Polarization Observations
In addition to the broadband interferometric observation, BIMAP simultaneously detects the polarization state of the mapped RF signals. A detailed description of BIMAP's full capability, including its polarization capability, has been reported in Shao et al. (2018). More details on effect of antenna pattern and effect of data processing on polarization analysis, and the first space-based polarization observation of lightning RF signals, have been reported in Jacobson (2001, 2002). Figure 8 shows the corresponding polarization observation for the FPB process. Figure 8a shows the time waveform of the RF signal, same as that in the earlier figures. Figure 8b shows the power spectrogram of the RF signal in dB scale. Figure 8c shows the spectrogram of the degree of polarization, indicating the fraction of the total RF power that is polarized (coherent). The degree of polarization ranges between 0 and 1, corresponding to totally unpolarized and totally polarized signals. In the plot, the degree of polarization is indicated by the colors and is shown in the range of 0.5 (black) to 1 (red). Figure 8d shows the ellipticity of the polarization with a measure of open angle formed by the ratio of the minor and major axes of the polarization ellipse. A linear polarization corresponds to a 0°open angle, and a circular polarization corresponds to a 45°open angle. An elliptically polarized signal will have an open angle between 0°and 45°. The sign of the open angle indicates the sense of the polarization rotation, either clockwise or anticlockwise. Figure 8d shows that most of the RF signals for this event are nearly linearly polarized, with the open angles close to 0°(green). Figure 8e shows the orientation of the polarization. The orientation is measured with the major axis of the polarization ellipse to BIMAP's EW baseline, and it has an unambiguous range of 180°, that is, from −90°to +90°in Figure 8e. Detailed description of the polarization analysis and presentation is reported in Shao et al. (2018) and Jacobson (2001, 2002).
Since the RF signals are mostly linearly polarized, we will focus on the orientation of the polarization for the mapped lightning sources. For a polarized signal, its orientation is expected to be parallel to the direction of the discharge current or the drift direction of the underlying electrons. Normally, the electron drift is driven by the local electric field, and the polarization orientation will provide a measure of the direction of the local electric field. However, for MeV and higher-energy electrons, local electric field has less an effect on the motion of the electrons, and the RF radiation is mostly produced by electron collisions in the air and electron deflection by the Earth's magnetic field, as will be discussed in a future report.
To obtain the polarization orientation for the mapped lightning sources, we require the degree of polarization ( Figure 8c) at a chosen time and frequency to be greater than 0.65 and at the same time require the signal strength (Figure 8b) to be greater than the background noise. Figures 8d and 8e show only the selected results above the two thresholds. At each time step related to the interferometric process, we find the median of the polarization orientation between 25 and 75 MHz (Figure 8e) and use it to represent the orientation for the specific mapped source. Figure 8f overlaps the polarization orientation over the mapped sources for the entire FPB process, and Figures 8g-8i decompose the results in three time sections, −7.6 to −5, −5 to 0, and 0 to 6 μs, respectively. Polarization orientations shown in these plots are geometrically transformed from the antenna plane to the AZ-EL plane. As discussed early, during the beginning 1-2 μs (Figures 6a and 6b), the lightning sources appear to propagate form upper left toward lower right although the angular errors are large and the sources are widely spread. During the next 1 μs (Figure 6c), an apparent arc-shaped source front is observed. The polarization analysis supports these interferometric observations. Figure 8g shows that the polarization 10.1029/2019JD032273

Journal of Geophysical Research: Atmospheres
orientations during the initial time interval (−7.6 to −6 μs, Figures 6a and  6b) are in parallel to the presumed propagation direction, suggesting that the underlying electrons move or drift in the direction of propagation. For the arc-shaped feature, the polarization orientations are orthogonal to the arc front, agreeing with the occurrence of a simultaneous source front. If the sources had traveled laterally back and forth along the arc, the polarization would be expected to be horizontally orientated, parallel to the arc. Figures 8h and 8i show that in the downward region the polarization is more vertically orientated, parallel to the direction of the downward source propagation. This indicates that the local electric field in this region is approximately directed uniformly downward.
It is interesting to note that the negative breakdown processes that followed the FPB and occurred atop the arc feature show randomly oriented polarization ( Figure 7e). As discussed in Shao et al. (2018), such observations indicate that the corresponding breakdown processes are random in direction, apparently due to a nonuniform and complex in situ electric field structure in that region. Given that the preceding virgin air FPB was along rather uniform downward field (Figures 8h and 8i), the field structure on top is either affected by the FPB or affected locally by the negative breakdown processes themselves.

Possible Explanations of the FPB Process
The polarization observation supports and further confirms the overall picture outlined in Figure 6o. The downward processes below the arc-shaped source front are similar to that reported in Rison et al. (2016), although there appear to be two groups of FPB in this event instead of a single FPB process. The main question is that what physics processes are responsible for the earlier RF sources that apparently propagate more horizontally and form an apparent arc-shaped source front before the downward FPB started?
Both the propagation and the polarization show an apparent change of direction before and after the occurrence of the arc-shaped feature. Given that discharges follow the direction of the ambient electric field, it is difficult to understand how the electric field would change its direction in such a small scale (100-200 m, Figure 6o) in the electrically undisturbed, virgin air in the cloud. On the other hand, the arc-shaped source front closely resembles a pancake front of a cosmic ray air shower (CRS; e.g., Chapter 16 in Gaisser et al., 2016) and the thunderstorm electric field will have minimal effect on the path of a CRS due to the ultrahigh energy (greater than tens to hundreds of megaelectron volts) of the CRS particles near the core. If the sources prior to the arc and for the arc itself are related to a CRS, the change of propagation and polarization directions can be understood. In this case, the CRS is likely entering the electrified storm region from the upper left direction (Figure 6o), and the following FPB will appear to start near the core and in front of the CRS (at the midpoint of the arc, Figure 6e) and propagate more vertically downward in the direction of the local thunderstorm field, as conceptually illustrated in Figure 9.
However, in a fair weather condition, only high-energy CRSs (>10 16 eV primary energy) can produce a strong enough RF signal detectable on the ground, if the RF antenna is in the forward direction of the shower. Such a RF signal is produced by deflections of high-energy electrons and positrons in the shower by the Earth's magnetic field (Schröder, 2017). The detected RF signal will be a brief, 10s ns coherent pulse due to the speed of light propagation of the CRS front (Shao et al., 2018). Nevertheless, the high-energy electrons in a CRS will ionize the air to produce a swarm of a large number of low-energy, eV-level electrons as they travel forward. When this process occurs in the thunderstorm's electric field, the swarm of low-energy electrons will collectively drift in the direction parallel to the electric field and can produce much stronger RF signals. Let us now examine the possibility of such a scenario.
When a 10 16 -eV CRS enters deep into the atmosphere, it produces about 10 7 electrons at any instant near the shower maximum (above~4-5 km asl, depending on its incident angle) with an average energy of about 20 MeV (e.g., Chapter 15-16 in Gaisser et al., 2016). Each CRS electron will produce more lower energy electrons in a cascade manner and eventually result in an average of 20 MeV/34 eV = 5.8 × 10 5 eV-level, low-energy electrons before the seeding CRS electron loses all its energy, where 34 eV is the ionization energy for the air. The low-energy free electrons are then lost quickly in 10s ns by attaching to the air atoms (Dwyer & Babich, 2011). Since the CRS continuously supplies about the same amount of the seeding electrons as it traverses an~1-km range near its shower maximum, it is safe to assume that there is a steady amount of 5.8 × 10 5 low-energy free electrons for each CRS electron at any given time instant. As the CRS traveling at the speed of light forward, a swarm of low-energy electrons is continuously generated behind the CRS front and the swarm follows the front and travels at an apparent speed of the speed of light.
In a fair weather condition without electric field, the low-energy electrons drift randomly in the air and will produce little RF signals. However, if the CRS enters an RREA field region in a storm, for instance, E = 110 kV/m at 8-km altitude, the low-energy electrons will drift in parallel to the field and result in a current of i = N le eμE, where N le is the total number of the low-energy electrons, e is the electron charge, μ ≅ 0.3 (m 2 /Vs) is the electron mobility at 8 km. In this case, i = 10 7 · 5.8 × 10 5 · 1.6 × 10 −19 · 0.3 · 1.1 × 10 5 = 3.1 × 10 −2 A. When such a current emerges after the CRS enters the RREA field from outside of the storm, the radiated RF signal can be estimated with Equation 6 in Shao et al. (2005) where r is the distance from the CRS to the probe, v the propagation speed of the current (the swarm of low-energy electrons in this case), θ is the angle between the directions of v and the probe, and i(0,t') is the emerging current in the retarded time frame. In a scenario that the CRS enters the presumed RREA field region from a zero field region, the magnitude of i(0,t') can be approximated as that estimated above, 3.1 × 10 −2 A. For r = 6-10 km as for the case of the FPB, v ≅ c for the swarm of the low-energy electrons, and a representative viewing angle of θ = 45°, we have dE = 0.37-0.23 mV/m. This is in general agreement with the observed RF magnitude at the time of the arc-shaped sources (e.g., at −5 μs in Figure 8a), in which 100 digit units are calibrated to correspond to 0.64 mV/m for BIMAP observations in Los Alamos. Here, the estimate is done for 8-km altitude, but the order of magnitude result will apply to the expected range of IC initiation heights (7-10 km) for New Mexico storms.
In addition, it is useful to estimate the thickness of the swarm of the low-energy electrons behind the CRS front for examination of possible radiating RF frequencies. The thickness for a CRS front itself is typically a few meters (e.g., Chapter 16 in Gaisser et al., 2016). If the low-energy electrons are produced and attached to the air atoms instantaneously, the thickness of the swarm would be also about a few meters. However, the thickness of the low-energy electrons is likely greater than the CRS thickness since the original CRS seeding electrons would be slowed down behind the CRS front during the ionization process and since it takes a finite time period for the free electrons to attach to the air atoms. For a rough estimate, we can use the attachment time to constrain the upper limit for the thickness. At 8 km, the attachment time is about 6 × 10 −8 s for the low-energy electrons (Gurevich et al., 2004). This suggests that the low-energy electrons produced at the CRS front will be captured by air atoms 18 m behind the speed of light front, giving a rough thickness of 18 m. Therefore, the likely thickness for the swarm of the low-energy electrons is in the range of few meters to about 18 m. If we choose 10 m as a representative thickness, the time period for the low-energy electron current will be 33 ns. Such a current pulse will produce RF signals around 30 MHz (1/33 ns), as observed by BIMAP in the initial microseconds of the discharge process (Figure 8b).
These order of magnitude analyses for a 10 16 eV CRS entering an RREA field indicate that both the RF magnitude and the RF frequencies are in general agreement with BIMAP's observations at the occurrence of the arc-shaped feature. Therefore, it is conceivable that the sources for the arc feature and that in prior are associated with the CRS in an RREA field.
In the above discussion, we have not considered electron multiplication effect due to possible RREA process seeded by the CRS electrons. First, it is not known whether or not the local ambient field is above the RREA threshold. Past balloon-borne observations showed that in-cloud electric fields are almost always bounded by the breakeven field in thunderstorms (Marshall et al., 1995), which is about 0.7 of the RREA field at any altitude. Second, if the local field is indeed above the RREA threshold but less than the field strength for conventional breakdown, the electron multiplication e-folding length would be many hundreds of meters and the e-folding time would be several microseconds at the IC height (Dwyer, 2010). Given the spatial scale observed in Figures 5 and 6, such a large e-folding length may not exist here, possibly due to a limited spatial 10.1029/2019JD032273 Journal of Geophysical Research: Atmospheres extent of the RREA field in the cloud. In this case, RREA will still increase the number of high-energy electrons from the CRS seeding electrons, but the increase will not be significant to change the rough estimate of the low-energy electrons and the result of the RF signal magnitude.
The observation and the analysis discussed above appear to suggest that the downward FPB below the arc-shaped feature (Figures 6 and 9) was ignited by a CRS. However, it has been convincingly argued by Dwyer (2010) and Dwyer and Babich (2011) that a CRS will not be able to initiate lightning, even with the maximum possible RREA electron multiplication in place. They showed that the density of low-energy electrons produced by RREA can never reach a level to significantly increase the conductivity to start a hot conventional lightning breakdown. For instance, Dwyer (2010) reported that with a large 10 avalanche length in a RREA field (sea level-equivalent of 284 kV/m), the density of the low-energy electrons at 8 km is as low as 4 × 10 10 m −3 for a high-energy CRS of 10 17 eV. With such an electron density, the air cannot be heated enough to transform from a diffused electron cloud to a hot lightning channel, as originally proposed by Gurevich et al. (1992Gurevich et al. ( , 1999. Instead of focusing on the possibility of a direct transition from RREA to a hot lightning channel, we examine here how the CRS-related low-energy electrons would affect the electric field ahead of the CRS front and examine if the altered field would be sufficient to initiate positive breakdown. As discussed earlier, a continuous and steady 5.8 × 10 12 low-energy electrons would be produced by the forward moving front of a 10 16 eV CRS, without the consideration of RREA multiplication. Since the majority of the high-energy CRS electrons are produced near the core (along the axis of shower propagation), within about r = 10 m from the core (90% of CRS electrons), the resultant low-energy electrons can be assumed to be within the same radius. With a thickness of dz ≅ 10 m the average density of the low-energy electrons can be estimated to be ne = 1.8 × 10 9 m −3 at 8 km, agreeing on order of magnitude with a more thorough calculation in Rutjes et al. (2019). This gives an average conductivity of σ = n e eμ = 8.6 × 10 −11 S/m, which is extremely small by any standards. However, it is still several orders of magnitude greater than the air conductivity of 10 −15 -10 −13 S/m. Since the low-energy electrons are expected to have energies in the range of eVs to <34 eV, the Debye length for the electron swarm is estimated to be a fraction of a meter, less than one tenth of the swarm depth. Therefore, even with the extremely low electrical conductivity, the electron swarm will behave as a finite conductive volume (shaded area in Figure 9) in the storm's electric field. In a simplified case that the CRS moves in parallel to the thunderstorm's ambient electric field and the front profile of the low-energy electron density is assumed spherical (the real profile near the core should be much sharper than a sphere), the electric field E at the midpoint in front of the of the CRS front can be approximated as (Page 284, Ward & Hohmann, 1988) where E 0 is the ambient field and σ 0 and σ e are the conductivities for the ambient air and the region of the low-energy electrons. Since σ 0 is much less than σ e , we have E = 3E 0 , a factor of 3 enhancement of the ambient field. Equation 8 is for a stationary finite conductive sphere. In the case of a CRS, the ionized region moves near the speed of light, and the electric field at the front is compressed and will be greater than the corresponding stationary field (Appendix C). In any case, the electric field in front of the CRS is expected to be enhanced due to the swarm of the low-energy electrons, likely by more than a factor of 3. If a RREA field exists in the cloud, a CRS of about 10 16 eV that sweeps through this region would likely enhance the field strength beyond the positive breakdown field, which is about 1.4 times of the RREA field in the air. Therefore, it is conceivable that the FPB below the arc-shaped feature in Figures 6 was ignited by such a CRS-related mechanism, as illustrated in Figure 9.
After the FPB started, the RF signals will be dominated by the FPB process which are apparently consisted of a great number of positive breakdown, as discussed in Rison et al. (2016) and as shown in this paper. In the meantime the CRS continues to propagate downward at the speed of light, but its RF emission is much weaker than that of the positive breakdown and cannot be simultaneously mapped by the lightning interferometer. However, the downward CRS will weakly ionize the region in front of the FPB and will help to increase the propagation speed of the FPB. As reported here and in Rison et al. (2016), FPB travels at an apparent speed greater than 10 7 m/s, about 2 orders of magnitude faster than any known experimental or 10.1029/2019JD032273

Journal of Geophysical Research: Atmospheres
theoretically predicted speed for virgin air positive breakdown. The preionization by the CRS may provide an explanation for the seemingly unphysical high propagation speed. Figure 9 summarizes our hypothesized concept for the initiating discharge process presented in this paper. A CRS enters the electrified storm region from upper left. Due to the electric field, the swarm of eV-level electrons generated by CRS' high-energy electrons in the CRS front produces detectable RF signals off of the core direction. This swarm of eV-level electrons at the same time enhances the electric field in front of the CRS front near the core to above the positive breakdown threshold and ignites the downward FPB process. Unlike the trajectory of the CRS that is little affected by the storm field, FPB is expected to propagate in the direction of the electric field, mostly downward in this case.
If the above analysis is valid, it can be inferred that CRSs as low as 10 14 eV could enhance the field to the positive breakdown threshold at the IC initiation height if the ambient field reaches the level of RREA or roughly the level of the breakeven field (which is 0.7 RREA, on the same order of the RREA field).

Summary and Discussion
In this paper, we analyzed the uncertainties on broadband lightning interferometry, based on previously reported statistic theories of coherent time delay estimates. Comparing the theoretical prediction with our numerical simulation, we found that when the raw signals are aligned in time (zero time delay), the uncertainties reach the theoretical lower bound (Figure 1a). However, if the raw signals are not aligned in time, as typical for lightning interferometric observations, the uncertainties can be much higher than the theoretical lower bound. The uncertainties increase as the size of the correlation window decreases, as demonstrated in Figures 1b and 1c. This is expected since the theoretical prediction is based on the assumption of large correlation window. For lightning observations, we desire to have as short a window as possible due to the impulsive, dynamical, and possible multiple-source nature of the lightning processes. In this case a normal interferometric analysis cannot use the lower bound as an indicator for the mapping uncertainties.
We then introduced a beam steering interferometry technique that process the raw data in two steps. It first processes the data with the normal interferometry, finds the general centroid direction of the sources, and estimates the corresponding signal time delays at the antennas. It then shifts the raw data in time as if to pull the sources overhead, applies interferometry to the shifted data, and finally pushes the reprocessed results back to the correct positions. With the alignment of the raw data in such a process, the uncertainties of the time delay estimate for each baseline approach the theoretical lower bound and the corresponding angular uncertainties can be reasonably estimated. In the two-stage interferometry analysis, a centroid direction is determined for a time interval during which all the lightning sources are located in a small angular region. If the sources travel a large angular region, the centroid direction cannot be used for the raw data shift, since it would align the signals for a certain section of the sources but would possibly increase the time delay for the other sections. Therefore, the two-stage technique is specifically useful for detailed analysis of spatially compact events. Or, if desired, one can apply such a technique consecutively through an entire lightning flash with consecutively determined centroid directions.
The beam steering interferometry improves the location accuracy for sources at lower elevation angles at which the time delays among the signals are greater. This should not be confused with the uncertainties due to dα j j ¼ c d sin α dτ . The technique improves the τ estimate for sources at smaller α by pulling the sources close to α ¼ π 2 , such that τ is equally accurate for sources from α = 0 to π 2 . With the same τ accuracy across the sky, the angular uncertainty still increases with 1 sin α due to the geometry effect.
We apply the beam steering interferometry on a FPB process that initiated a normal IC. It is immediately clear that the lightning-initiating FPB is much more complex than we have previously observed. As described in a great detail in last section, the FPB analyzed in this paper shows a more horizontally orientated arc-shaped feature 1-2 μs into the beginning of the mapped sources. Then two apparent groups of FPB started from the midpoint of the arc propagated mostly vertically downward at the speed of several times 10 7 m/s. The speed estimated for each apparent group is similar to that reported in Rison et al. (2016).
In addition to the detailed source map, we also analyzed the polarization of the corresponding RF signals. Together with the source structure and development, the polarization result supports that the sources of the arc feature and that in prior are consistent with CRS-related RF sources when the CRS entered a high field region from the left (Figures 6 and 8). The two groups of the downward FPB are then likely ignited at the front of the CRS.
To investigate the CRS hypothesis, we examined the possible RF signal strength and frequency content for a 10 16 eV CRS entering a region of RREA field in the storm. With rough order-of-magnitude analyses, we found that such a CRS is capable of producing the same level of RF signal at about the same frequencies as that observed for the arc-shaped feature. We further examined the enhancement effect of the CRS front on the ambient electric field and found that it is possible to increase the field strength by a factor of 3 or more at the midpoint and in front of the CRS front. Such an enhancement could raise the field strength above the conventional positive breakdown field threshold and would be able to ignite the downward FPB. In our analysis, we have not considered the electron multiplication effect due to possible RREA process seeded by the CRS electrons, but only the low-energy eV-level electrons produced by the CRS electrons. Figure 9 summarizes our hypothesized concept based on the observations presented in this paper. To have a significant RREA electron multiplication, there has to exist a region of many hundreds of meters to a couple of kilometers of RREA field in the cloud (e.g., Dwyer & Babich, 2011), which appears not to be the case given the small spatial scale of the mapped RF sources. Furthermore, we hypothesize that even with the lower breakeven field (0.7 of RREA field), a CRS with a primary energy as low as 10 14 eV might be able to enhance the field enough to ignite positive breakdown at 8-km height.
In our analysis, we did not attempt to examine the possible transition from a RREA process to a hot lightning channel, as has been studied exhaustedly in a large number of previous papers (e.g., Gurevich et al., 1992Gurevich et al., , 1999. It is conceivable such a process may never occur, as discussed in Dwyer (2010) and Dwyer and Babich (2011). If our analysis is valid, it will not need a RREA process to initiate a lightning, and the electric field in the cloud will not need to be unphysically high. Apparently, further thorough theoretical and observational studies are needed to validate our preliminary analyses, especially observations with sensitive, accurate, higher time resolution and polarized broadband lightning interferometry.

Appendix A: Overlap of Consecutive Correlation Windows in Interferometric Analysis
In broadband lightning interferometry analysis, the correlation window slides consecutively forward along the time series data. In the lightning community, the step size of the forward sliding is chosen arbitrarily without any fundamental physics or mathematics basis. For a stationary random process that is generated by a single slow moving source, the step size is commonly chosen to be half of the correlation window length. In this ideal case, a half window slide will provide a complete measure for the source statistics. However, lightning sources are impulsive, dynamical, and cannot be treated strictly as a stationary process, especially when the correlation window length is long enough to include a sequence of different discharge processes. Therefore, for any correlation window that is different than the other, even by a single data point, the signals (discharge processes) covered by the two windows can be different, as shown in Figure 2. If they do cover the same signals/processes, one would expect the sources to be mapped at the same position. In this case the source position is over determined but not falsely determined, and there is no extra information generated. In our analysis, we choose to slide 4 data points each time to skip over a full signal cycle at the middle of BIMAP's frequency band.
To demonstrate the step size effect, Figure A1 shows the same FPB map with step size of 4, 8, 16, 32, 64, and the correlation window size of 128, respectively from a to f. The overall skeleton structure across the figures is the same, but with a larger step size, it misses the apparent source positions that would be covered by an intermittent window. For the lightning process, it is safe to assume the sources vary from one raw data point to the next, and it is better to overdetermine the positions than to miss the possible positions.

Appendix B: Interferometric Centroid of Simultaneous Sources
With the time delay-based interferometric analysis presented in this paper, if two or more lightning sources occur within the same time correlation window and are within the array's angular resolution, a 10.1029/2019JD032273

Journal of Geophysical Research: Atmospheres
power-weighed centroid position will be obtained. Simultaneous sources that are separated by a larger angle would not pass either our correlation-coefficient or close-loop time delay criterion and will not yield a viable source location. Figure B1 illustrates a scenario with two sources that separate by 1°in a presumed 1-D direction of 54°from the baseline. Applying the first-step in the beam steering interferometry, we steer the phase center to the middle between the two sources, the two sources are −0.5°and 0.5°from the phase center, as indicated be the two circular dots (S1 and S2) in the plot. In this case, the effective baseline length is d sin 54°= 15.4 m, where d = 19 m for BIMAP.
We generate two independent noise series and two independent signal series, all in Gaussian. The standard deviations for the noises are set to be 1. The standard deviation for S1 is set to be 10, giving a SNR of 100 for S1. We select a range of standard deviations for S2, between 1 and 100 (SNR: 1 to 10 4 ), to test effect on the centroid position due to different power ratios between the two sources.
For Antenna 1, we added the noise and the two signals directly. For Antenna 2, we compute the time delay Δt corresponding to 0.5°with the effective baseline length, shift S1 and S2 −Δt and Δt respectively, and add the shifted signals to the other noise series. Since a small angular deviation of 0.5°yields a time delay much less than BIMAP's sampling resolution (5.3 ns, 190 MHz), all the noises and signals are initially generated with a sampling rate of 100 times of the BIMAP sampling rate. We then apply 20-to 80-MHz filtering (see the main text) to the simulated data for the two antennas, subsample the data to BIMAP's 190 MHz, and process the data with 128 data points (as for Figures 5-7) to estimate the resultant angular position.
For each selected S2 to S1 power ratio, we repeat the above process 10,000 times and determine the mean and standard deviation of the corresponding centroid positions. For instance, for power ratio of 0.01, the position was determined to be about −0.4°, near the stronger S1 position, with a standard deviation of about 0.06°, as indicated by the leftmost data point Figure A1.
(a-f) Source maps with 4, 8, 16, 32, 64, and 128 slide points for a correlation window of 128 points. Figure B1. Interferometric centroid position for two simultaneous sources, as a function of relative power ratio between the two sources.
in Figure B1. When the power ratio equals 1, the centroid is determined to be about 0°, at the midpoint between the two sources. As the power for S2 increases, the centroid position moves closer to S2, as would be expected. The decreasing angular error (vertical bar) toward the right is due to the increasing SNR.
It is interesting to note that the centroid position does not vary linearly across the large range of the power ratio between the two sources. As indicated in Figure B1, when the power ratio changes from 1 to about 10 (or equivalently to 0.1), the centroid position moves relatively quickly from midpoint toward 0.4°(−0.4°), near the position of the stronger source. This feature suggests that the centroid position between two sources can be readily dominated by the stronger source when its power is 10 times (3 times in RF amplitude) or greater than the other source.