Study of oblique whistlers in the low-latitude ionosphere , jointly with the C / NOFS satellite and the World-Wide Lightning Location Network

We use the C/NOFS satellite’s Vector Electric Field Instrument (VEFI) to study the relationship of impulsive electron whistlers in the low-latitude ionosphere to lightning strokes located by the World-Wide Lightning Location Network (WWLLN). In order to systematize this work, we develop an automated algorithm for recognizing and selecting the signatures of electron whistlers amongst many Very Low Frequency (VLF) recordings provided by VEFI. We demonstrate the application of this whistler-detection algorithm to data mining of a∼ two-year archive of VEFI recordings. It is shown that the relatively simple oblique electron whistler adequately accounts of the great majority of lowlatitude oscillatory VLF waves seen in this study.


Introduction
Transient electric currents in lightning radiate Very-Low-Frequency ("VLF"; 3-30 kHz) pulses known as "sferics" (Volland, 1995).Sferic peak radiated power ranges to >10 11 W, making sferics the most powerful and common impulsive VLF noise source on Earth.Moreover, sferics are effectively ducted in the Earth-ionosphere waveguide (Volland, 1995), so that they can be detected even after propagating over global distances (see, e.g., the review in Said et al., 2010).
The Earth-ionosphere VLF waveguide is not quite perfect, however.There is eddy-current dissipation in the ground, at the lower boundary (Cooray and Ming, 1994).There is also continuous leakage of energy upward through the upper boundary, penetrating the ionosphere, in what has Correspondence to: A. R. Jacobson (abramj@u.washington.edu)been called the "penetrating mode" (for an insightful fullwave treatment, see e.g., Piggott et al., 1965;Pitteway, 1965).The penetrating mode becomes the oblique electron whistler in the ionosphere, where a WKB description is valid.The oblique whistler at ionospheric heights, as opposed to whistlers in the outer magnetosphere, propagates nearly vertically, with horizontal wavevector constrained (by Snell's Law) to match that of the underlying sferic, and with the much stronger vertical wavenumber given by the Appleton-Hartree dispersion relation (see, e.g., Sect.4.9 in Krall and Trivelpiece, 1973).
Whistlers in the ionosphere have been observed both from rockets (see the references in Holzworth et al., 1999;Kelley et al., 1990) and from satellites (recent reviews of two decades of results can be found in Holzworth et al., 2011;Santolik et al., 2009).The present article presents a method to advance and systematize the exploitation of oblique whistlers observed by VLF sensors on low-Earthorbit satellites.These applications include study of ionospheric electron density, and the details of VLF propagation through the lowermost ionosphere, the D-region (Lehtinen and Inan, 2009).order of 10 s.During these burst recordings, the electric field is determined nearly from DC nearly to the Nyquist frequency, 16 kHz.This allows an opportunity for study of lightning-generated waves ("sferics"), although lightning per se is not the satellite's or VEFI's main mission.

Extent of the database used here
The overall period for this study runs from 4 November 2008 through 30 November 2010, roughly two years.We focus here on VEFI electric-field "burst" recordings, digitized at nominally 32 kilosamples s −1 , the fastest rate available, with 16-bits range.We are not exploiting the wave polarization but are only interested in the dispersion of temporal features, and in the dispersion-corrected time-of-arrival of those features.Therefore we limit this study to a single electric-field component "E 34 " on the antenna-boom pair 3-4, or "34".During the two-year duration studied here, there were 15 259 burst-mode recordings of E 34 at the 32-ksample s −1 rate, each recording several seconds duration (variable).

Automated detection of whistler dispersion in VEFI
VLF data

Expected dispersion for oblique whistlers
We shall describe our approach and criteria for recognizing ("detecting") whistler events in the VLF data from VEFI, a posteriori.The search for whistler events is conducted anew in each analysis window handed to the algorithm.Within each analysis window we perform an automated search for dispersed pulses displaying whistler-like spectral dispersion (Krall and Trivelpiece, 1973).This form of dispersion is expected of the upward-propagating "penetrating solution" in the lower ionosphere, or D-region (Piggott et al., 1965;Pitteway, 1965), matching to the oblique whistler wave.The penetrating solution is the simplest way for lightning-generated VLF in the Earth-ionosphere waveguide to leak energy upward toward the satellite on the ionospheric topside.It requires no "duct" structures but rather relies simply upon Snell's law, i.e., the conservation of horizontal wavevector.In the quasi-longitudinal approximation, for frequencies f below both electron gyrofrequency f ce and the electron plasma frequency f pe , we expect the phase to vary with frequency according to where φ is phase (radians), z is the vertical axis, and β is the angle between the magnetic field, B, and the vertical, z.This takes the form that phase is proportional to f 1/2 .For the latter in Hertz, the constant of proportionality (henceforth, "dispersion constant") has units of s 1/2 .For the range of heights of C/NOFS we expect the dispersion constant to lie in the range 10-100 s 1/2 .The physical significance of the constant is that it is a vertical integration of the electron plasma frequency through the ionosphere, from the D-region up to the satellite altitude, with a weighting factor determined by the geomagnetic field vector.This integration is loosely analogous to the more common "Total Electron Content", or TEC, measured e.g., with Very High Frequency (VHF) emissions of lightning (Huang andRoussel-Dupré, 2005, 2006;Jacobson et al., 1999;Roussel-Dupré et al., 2001).The approach to retrieving the "best" dispersion constant will follow the same formalism as developed for retrieving TEC using VHF signals, but with a dispersion relation appropriate to VLF rather than VHF.

Standard analysis windows
To illustrate the methods used, we will focus on one day, 4 January 2010, during which VEFI gathered 29 burst recordings of E 34 .Each recorded burst is identified by the UT trigger time, as yymmdd hhmmss followed by three digits for additional milliseconds.For example, the first burst record on 4 January 2010 was triggered at 00:09:18.971s UT and thus is identified as "100104 000918971".
Figure 1 shows a moving-window spectrum of the entire burst record of E 34 during this first recording of the UT day.There are nearly 4 × 10 5 samples.Numerous features are visible, including an especially intense, apparently dispersed pulse a little after halfway through the record.To provide an automated analysis of all the data in a burst record, we divide it into successive analysis windows, each 16 384 samples long, and overlapping by 50 % (8192 samples).For example, into this record we can place 45 such analysis windows.The small leftover at the end of the record is excluded from the analysis.

Prewhitening
Now to illustrate prewhitening, we consider one analysis window, the 21st window within the burst record of Fig. 1.This gives 16 384 samples of E 34 , and we suppress E 34 in the first half, so that we are left with 8192 zeroes followed by 8192 surviving samples of E 34 .This "zero padding" provides room for the leftward migration of signal energy in the dechirping that follows.We then perform a single, coherent Direct Fast Fourier Transform (DFFT) on these 16 384 samples.Once in the frequency domain, we rank frequencies according to power.The complex Fourier transform for any frequency in the top 10 % of power, that is above the 90th percentile, is then clamped in amplitude at the 90thpercentile amplitude.The phases are not altered.This crude, non-causal "prewhitening" effectively removes most effects of narrow-band anthropogenic carrier signals (Jacobson et al., 1999).We then take the prewhitened Fourier coefficients and return to the time domain.Figure 2c shows the movingwindow spectrum of the prewhitened E 34 , within the 21st window.Most of the energy is on the right side (the latter half) of the time distribution.Some energy has leaked, however, into the first half.That is due to the non-causal manipulation of the Fourier coefficients during prewhitening.

Dechirping against whistler dispersion
Let us once again consider the complex Fourier coefficients, after they are prewhitened in the above manner, but before returning to the time domain.For any choice of the dispersion constant, we can compensate for the spectral phase dispersion according to the f 1/2 variation explained above.Once the phase is compensated for that trial value of the dispersion constant, we can return to the time domain and inspect the result.Figure 2d shows the moving-window spectrogram of prewhitened, phase-shifted E 34 , within the 21st window.This is like Fig. 2c except for the addition of phase-shifting to the manipulation of Fourier coefficients.Obviously the peaks are sharpened and made more intense (vertical), at the expense of the slowly-varying pedestal.The phase-shifting to sharpen the peaks is called "dechirping" (Jacobson et al., 1999).We can optimize this dechirp by maximizing a quantitative figure-of-merit: first, we take the square of the timedomain electric field after it has been dechirped with a trial value of the dispersion constant.Second, we smooth the squared field by a running window of 3 samples.Third, we take the fourth power of the smoothed E 2 .Fourth, we sum this highly nonlinear time series over the entire analysis window.This sum is called the "figure-of-merit".It is a way of favoring highly concentrated peaks at the expense of the pedestal, due to the intense nonlinearity of the eighth power of the electric-field envelope.Figure 2a shows the figure-ofmerit, after normalization by its peak value, as a function of trial value of the dispersion constant.The constant is gridded in steps of 1 s 1/2 , from 0 to 100 s 1/2 .We then iterate the scan, on a finer scale from 80 % to 120 % of the preliminary peak from Fig. 2a.This iterated scan for figure-of-merit is shown in Fig. 2b.Its peak occurs at 36.4 s 1/2 .This is then adopted as the optimum value of dispersion, and is used for the final dechirping shown in Fig. 2d.
The entire division of a record into overlapping analysis windows, then prewhitening, then dechirping, is performed automatically.

Acceptance criteria for automated dechirping
The product of the dechirping for each analysis window is a gridded "figure-of-merit", at two nested resolutions in the dispersion constant.After the first resolution is calculated, two criteria are required to be met, and a decision is made whether to proceed further.The first criterion relates to the width of the peak in the low-resolution figure-of-merit (see Fig. 2a).We sum the self-normalized figure-of-merit over the whole domain (0 to 100 s 1/2 .)This sum is effectively a ponderation width.We require this width to be less than 5 % of the maximum dispersion measure sampled, i.e. to be less than 5 s 1/2 .The second criterion is to require that the normalized figure-of-merit at the left end of the grid, i.e. at zero dispersion, not exceed 0.2.This ensures that the peak in the figure-of-merit not encroach on the zero-dispersion limit, where other forms of electric-field disturbance are likely to dominate.These two criteria are not completely independent.

Automated identification of peaks in dechirped signal
A diagnostic time series is defined to facilitate automated identification of peaks in the dechirped signal.The electric field undergoes a direct Fast Fourier Transform (DFFT), is prewhitened in the frequency domain; Shown here is the spectral density of the prewhitened field, with a moving 256-sample FFT.Most of the energy is missing from the first half of the window, which had been zeroed prior to the first DFFT; there is however some leakage due to the prewhitening.(d) Similar, but with the electric field dechirped for the optimal dispersion measure (36.4 s 1/2 ).The peaks in this window (and all other windows) are determined by the peakfinder algorithm; three peaks are found in this window.The first peak in this window is the 45th peak (out of 98 total) in the overall record.
in two steps.First, the optimally dechirped signal is squared and averaged over 21 samples around each point with a running average window.From this smoothed signal we then subtract a running average of 501 samples around each point, for baseline subtraction.(This vignettes the 250 points at either edge, but since the window contains 8192 samples of non-suppressed data, this edge effect is only a small penalty.)The resulting series is the "indicator" series, which is then sent to a peak-picker algorithm.The peak-picker identifies peaks meeting three criteria: first, the "contrast" of a candidate is calculated as the peak divided by the median with respect to the entire window.We require contrast >3.Second, the width of a candidate peak is defined as a ponderation width, namely the sum over samples within the peak, divided by the maximum peak value.We require width <200 samples.Third, we look at the proximity of neighboring peaks and eliminate neighbors closer than 100 samples to each other.The surviving neighbor is the higher; the lower neighbor is eliminated.In this manner, for example, three peaks are accepted in the window of Fig. 2d.These peaks are #45, 46, and 47 in the entire burst record, which contains a total of 98 peaks automatically determined by these procedures.All three peaks in this particular window are associated with the same optimal dechirp constant, 36.4 s 1/2 .The assigned dechirp constant is applied to all peaks in a given analysis window.

Automated generation of peak files
For each burst recording, we compress the selected peaks and supporting data into a "peakfile".If there are no qualifying peaks, then the file contains no data.If there are peaks, then for each qualifying peak the following is stored: (1) 256 samples of dechirped electric field data (25 % pre-peak, 75 % post-peak), (2) the UT time of the dechirped peak, (3) the index of the analysis window in which the peak is detected, and (4) the optimal dispersion constant for the window in which the peak is detected.The peakfile stores the VLF waveform data of each qualifying event in dechirped form, which is temporally compressed.Thus the peakfiles constitute a small volume of data (e.g., only ∼1.5 Gbyte total for the two-year study period.)Table 1's upper three lines tabulate the characteristics of the peakfile archive by itself.Of the 15 259 eligible burst recordings, 7677 contain at least one identified peak.The total number of peaks in the two-year dataset is 420 265.The example day of 4 January 2010, treated in Figs. 1 and  2, is quite typical of these statistics.The day comprises 29 burst recordings.The highest count of peaks found in any one burst record of this day is 276, and the next highest (98) is found in the burst record of Fig. 1.The other records containing any peaks have 92, 61, 36, 29, and 19 peaks, respectively.Each of the remaining records contains zero qualified peaks.6 Illustrative behavior of qualified peaks

Competing dispersion effects
Consider the 45th qualifying peak within the burst recording "100104 000918971".This peak is the first (of three) within that record's 21st analysis window and is the first of three peaks in Fig. 2c, d.In Fig. 3 we illustrate the efficacy of the dechirping procedure.Figure 3a shows the smoothed square of the dechirped electric field, while Fig. 3b shows the smoothed square after rechirping the electric field, that is, removing the phase correction.(This rechirp is implemented not on the whole analysis window but just on the 256 samples of the peakfile, and thus results in a circular wrap of energy in the 256-sample domain.) Figure 3a, b has the same vertical scale, and illustrate the "gain" of the dechirping procedure.Figure 3c, d repeats this comparison, but for the oscillating electric field, not its smoothed square.Again, the effect of dechirping is to concentrate energy near sample #64.However, in the dechirped electric field (Fig. 3c) there is also an extended coda at lower frequency than the oscillations near sample #64.We explore this further with the moving-window spectrogram of the dechirped electric field, in Fig. 4. The FFT window is 32 samples and is advanced in steps of 2 samples.The color-scale shows the logarithm (base 10) of the spectral density.The extended coda is around 3-4 kHz.This residual dispersion has not been dechirped by the whistler algorithm (Eq.1), which searches for dispersion having a pole at zero frequency, not at finite frequency.Since all of the data in our burst-recording archive is for a nighttime Earthionosphere waveguide, it is not surprising to find "tweek" (Kumar et al., 2008) features such as the coda in Fig. 4.These tweeks are seen to some extent (though infrequently in as dramatic form as in Fig. 4) in much of our dechirped data.The tweeks slightly perturb the whistler dechirp, biasing the optimal dispersion constant to a higher estimate than is true for the whistler alone.This can be seen in Fig. 4, where at higher frequencies well above the coda, say f > 6 kHz, the dechirped pulse is slightly over-corrected for its whistler dispersion.If one found it necessary to remove this slight bias caused by tweeks, then the dispersion constant would need to be corrected on a peak-by-peak basis, not for an entire analysis window.
We have examined power-vs-time within a peak file, for numerous peaks.There is no restriction to simple, monopulse impulses.Rather, some of the dechirped peak curves are quite complex, even though they have been optimally dechirped.This is not unexpected, due to the dispersion of a sferic waveform imposed by lateral propagation in the Earthionosphere waveguide prior to the upward leakage of energy into the oblique whistler.

Variation of peak behavior within a burst record
For the burst record "100104 000918971", we now illustrate the variation of fitted dispersion versus time within the record.The record is only ∼12 s in duration, during which the satellite travels only ∼90 km along its orbit.We see variations within a few seconds, but these probably cannot be interpreted as physically significant in terms of our model.Our consideration of the predicted whistler dispersion (Krall and Trivelpiece, 1973) is based on the WKB approximation, which requires the scale of horizontal variation of the medium to greatly exceed the inverse of the horizontal wavenumber.The horizontal wavenumber is set by Snell's Law governing the coupling of the quasi-guided wave below the ionosphere to the "penetrating solution" propagating upward through the D-region.It is not likely that a WKB interpretation is applicable to explaining the dispersion's finescale variations.
Figure 5 shows the initial normalized figure-of-merit for each of the 45 analysis windows in the burst record "100104 000918971".Each window marked with an "X" fails the test for acceptable figure-of-merit (see Sect. 4.5 above).The dispersion constant for each of the 98 qualifying peaks is shown in Fig. 6, versus peak index.The grouping into several analysis windows is obvious from the plateaux.Figure 6b shows the dispersion variations on a finer scale, and evidently the variations amount to ±5 % of the mean dispersion.

Implications of dispersion variations for peak timing
If the fluctuations in dispersion cannot be interpreted (within a WKB framework) straightforwardly, then the fluctuations must be seen as causing "propagated errors" in estimates of the peaks' occurrence times.We illustrate this in Fig. 7, showing the smoothed square of the dechirped electric field vs. time, for the same peak #45 as is in Figs.3-4.The central peak is for the parent analysis window's optimum dechirp (constant = 36.4s 1/2 .)The two displaced peaks, each lower and wider, are for ±5 % adjustments (relative to 36.4 s 1/2 ) of the dechirp.We note that even a "small" perturbation of the dispersion constant results in a ±1.5-ms "propagated error" in the peak's arrival time.This is quite typical of the arrivaltime propagated errors in our overall two-year dataset.These errors dwarf the uncertainties from ground-truthing systems such as WWLLN, and are the controlling limit of coinci-Fig.7. Effect of dechirp errors on the time-of-arrival of peak, using peak #45 from the burst recording of Fig. 1.The light solid curve is the smoothed square of the dechirped electric field, using the optimal dechirp (36.4 s 1/2 ).Each heavy curve is for dechirping with a different dispersion constant, perturbed from optimal by +5 % (solid) and −5 % (dashed).
dence between VEFI dechirped signatures and those groundtruthing data.

Example of a lower-dispersion record
During the same day, 4 January 2010, the record containing the most (276) qualifying peaks is "100104 141821566".The electric field for the entire record is shown in Fig. 8.This record contains 46 analysis windows. Figure 9 shows the initial normalized figure-of-merit vs. dispersion constant, for each of these 46 windows.Only four windows (marked "X") fail the figure-of-merit criteria.The variation of iterated dispersion for the peaks is shown in Fig. 10 (which is similar to Fig. 6).Now the optimal dechirp is at roughly half the value of dispersion constant seen earlier for record "100104 000918971".The fluctuations in dispersion are several % of the median.The oval marks a series of peaks seemingly following a sustained downward trend in dispersion.
Figure 11 shows the iterated normalized figure-of-merit vs. dispersion constant, for the 42 analysis windows within "100104 141821566" passing the figure-of-merit test (see Fig. 9    oval in Fig. 10) is now marked with asterisks.These highresolution figure-of-merit plots show another perturbing effect: oscillatory variation vs. dispersion constant.This oscillation is due to interference between a limited num- ber of Fourier components contributing to the signal being dechirped.The oscillation-induced errors are still small compared to the random (or at least uninterpretable) variations of several %.

Records yielding no peaks
About half the burst records during the two-year study (Table 1) yielded no peaks.During the example day 4 January 2010, only 8 of the 29 burst records yielded any peaks.
Usually the reason for this is the presence of competing electric-field perturbations that do not exhibit characteristic whistler dispersion.Observing these "competing" perturbations is in fact a principal goal of the C/NOFS mission (de La Beaujardiere, 2004).These perturbations are often related to in situ ionospheric plasma irregularities responsible for scintillations on transionospheric radio-communications links.Figure 12 shows the electric field for the entire record "100104 064015612", a record from 4 January 2010 that yielded no qualifying whistler peaks.Most of the impulsive features in Fig. 12 entail a step change in the baseline, which is of course unlike the start of a dispersed whistler.Amidst these non-whistler events one's eye might detect, perhaps, a couple of whistler-like disturbances.Figure 13 shows the preliminary figure-of-merit plots for all 46 analysis windows in this record.Every analysis window fails the figure-ofmerit test.Most of these failures are quite obvious, as witnessed by the slow slope up to the left (low dispersion) end of most plots.That is symptomatic of features lacking systematic spectral dispersion.The two possible whistler-like   disturbances from Fig. 12 occur during the two windows marked "W" in Fig. 13.Each of these windows' figure-ofmerit shows less baseline slope upward to the left than do their neighboring windows.
In general, the presence of strong non-dispersed electricfield signatures in a burst record tends to prevent identification of dispersed whistlers, at least with the conservative selection criteria imposed in our algorithms.This is certainly borne-out by the burst record "100104 064015612", in which even the two visually-apparent whistler disturbances are rejected.

Propagation model
WWLLN provides lightning-stroke location reports with finer than 30-µs temporal accuracy and finer than 15-km spatial accuracy (Jacobson et al., 2006;Lay et al., 2004Lay et al., , 2007;;Rodger et al., 2004Rodger et al., , 2005Rodger et al., , 2006)).We assume that the radiated VLF "sferic" from the stroke propagates in the Earth-ionosphere waveguide to the vicinity of the subsatellite point, from whence it is upwardly coupled into the penetrating solution (Piggott et al., 1965;Pitteway, 1965) and then propagates vertically to the satellite according to the oblique whistler dispersion relation.This very simple sequence of propagation, first laterally outward in the Earthionosphere waveguide, then vertically through the ionosphere, has been extensively validated elsewhere (see, e.g., Holzworth et al., 1999).It relies on the Quasi-Longitudinal (QL) approximation to the Appleton-Hartree dispersion relation.For typical ionospheric parameters of the equatorial ionosphere, this QL approximation is valid to within a fraction of a degree of the magnetic equator, effectively everywhere sampled by the C/NOFS orbit.
Although the correct effective group speed in the Earthionosphere waveguide is slightly less than c (Dowden et al., 2002), we use exactly c for that leg of the propagation, because the error incurred is still small compared to the much larger errors due to uncertainty in the whistler dispersion constant.As shown earlier in Sects.6.2 to 6.4, the whistlerdispersion uncertainty typically adds a few (to several) milliseconds of blur to our model estimate of the dispersive delay of the whistler.This dwarfs the error incurred by assuming speed exactly c in the Earth-ionosphere waveguide.The error in the whistler dispersive delay similarly dwarfs the errors in the WWLLN stroke location/timing.The dispersionimposed error also exceeds the typical on-board timing accuracy (1 ms) of VEFI timestamping.
The dechirped-peak time in the dechirped VEFI signal is already compensated with respect to the dispersive portion of the propagation delay from the subsatellite point to the satellite.(Notice the leftward shift of the peaks in Fig. 2d relative to Fig. 2c.)To compare the listed WWLLN stroke UT times with the VEFI dechirped-peak UT times, we need to subtract from each VEFI dechirped-peak time two additional terms.The first is the vacuum time-of-flight from the satellite down to the base of the ionosphere (∼85 km at night), and the second is the propagation in the Earth ionosphere waveguide.These two timing corrections, added to the correction already done of the dispersive delay, effectively transport the VEFI time back to the candidate lightning.We call this the "corrected" VEFI time.

VEFI-WWLLN coincidences
We will now compare the corrected VEFI times with the times of candidate lightning strokes provided by WWLLN.During the cumulative time covered by the 7677 burst records with at least one qualified VEFI peak, there are 179 730 WWLLN-located strokes within 10 000 km of the subsatellite point (see Table 1).Figure 14a shows a histogram of the corrected time difference (VEFI -WWLLN), in 1-ms bins, for all stroke locations within 10 000 km of the subsatellite point.Figure 14b is similar, but restricted to stroke locations within 1000 km.The baseline in such a histogram gives a measure of the level of "accidental" coincidence (Jacobson et al., 2000), so that the peak-to-baseline ratio is effectively the ratio of "real + accidental" to "accidental" coincidences, or the coincidence signal-to-noise ratio (SNR).Figure 14 shows that the SNR exceeds 8 for distance <10 000 km, and approaches 20 for distance <1000 km.
There is also a distinct secondary peak in Fig. 14 near −10 ms.This corresponds to a signal occurring before the main stroke by 10 ms, such as a leader step.Although WWLLN detects and locates many return strokes, it normally cannot detect/locate low-amplitude events such as leader pulses.The VEFI burst records in this study, on the other hand, are not triggered by the VEFI signal, and thus are not limited to signals from strong strokes.In fact we have seen many VEFI dispersion-corrected signals whose cadence and complexity resemble leader steps preceding a negative cloud-to-ground initial stroke.Note that the central peak in Fig. 14's histograms has a full-width at half-maximum of around 5-8 ms.This width is primarily a propagated effect from the uncertainty of the dispersion constant.If that uncertainty could be relieved, then the peak could be made narrower, and this would give a higher histogram SNR.
Somewhat arbitrarily, we set the criterion for "coincidence" to be the central portion of the peak in Fig. 14a, from the peak out to, and including, the bins whose count is greater than 5× the median.This translates to VEFI-WWLLN time differences in the range −3.0 to 0.0 ms.For range <10 000 km, there are 12 170 such coincidences, and for range <1000 km, there are 4104 (see Table 1).For all of these coincidences, we attribute the WWLLN-supplied stroke location as the source of the signal observed by VEFI.

Distribution of distance to coincident lightning
We expect stronger signals to be more detectable by the dispersed-whistler algorithm, relative to weak signals.Apart from the natural random variability of power of lightning VLF emissions, a more deterministic control is exerted by the distance from the lightning to the subsatellite point.Several effects converge: first, as the distance increases, then the received electric field should scale as (distance) −1/2 , as befits a spherical-shell waveguide such as the Earth-ionosphere waveguide.Second, as the distance increases, lossy attenuation further degrades the signal signal arriving at the receiver.Third, as the distance increases, it becomes more likely for the lightning source to be located in daylight (even though the satellite is in night), so that at least part of the waveguide path involves a relatively lossy, daylit D-region as the waveguide upper boundary.
Figure 15 shows (solid curve) the distribution of distances from the subsatellite point to WWLLN-located lightning contemporaneous with those 7677 burst records containing at least one peak.The distribution for VEFI-coincident WWLLN-located lightning is shown as the dashed curve.Each curve is separately normalized to its own peak, for ease of comparibility.The dashed curve is basically a distribution of source distances for VEFI-detected peaks.It is clearly biased in favor of shorter distances, consistent with our a priori expectations just outlined.

Distribution of integrated energy in VEFI-detected peaks
We have seen in Sect.7.3 clear evidence of the short-range bias for lightning sources of VEFI-detected peaks.Presumably the underlying cause for this bias is the obvious one: that the energy in the VEFI-detected pulses would tend to decrease with distance from the lightning source, assuming similar source-energy spectra for all locations.We now test this explicitly for the WWLLN-coincident VEFI pulses, integrating E 2 to obtain a measure of the received energy for each peak whose source has been reliably located.For this exercise we lift the distance limit (<10 000 km) that has been used so far.The prior limit, 10 000, is roughly the arc distance equal to one-quarter of the Earth's circumference.Beyond this threshold, the path extends into the antipodal hemisphere, where, strictly speaking, the curvature of the waveguide fronts reverses sign.If there is coherent propagation in a spherical-shell Earth-ionosphere waveguide, then we expect this one distance-dependence control (listed first in Sect.7.3) to reverse sign beyond 10 000 km (Dowden et al., 2002).Thus, we expect one (of the three) biases against long-range detection to reverse sign, but the other two biases are unaffected and continue to mitigate against long-range detection.When we lift the distance limit, we find only about 1/3 as many coincidences from opposite-hemisphere lightning (>10 000 km) as from same-hemisphere lightning <10 000 km (see Table 1).Considering both hemispheres together, we show in Fig. 16 a scatter plot of the integrated square of the electric field recorded by VEFI, vs. distance from the subsatellite to the coincident source lightning.
For distance <10 4 km, Fig. 16's trend in energy is downward vs. distance.There is much scatter, but the linear loglog trend might be consistent with a power law: energy ∼ (distance) λ , with λ in the range −3 to −5.Note that the effect of the geometrical spreading within the waveguide can account for only (at most) a −1 contribution to λ at close distances, expected to reverse sign at longer distances >10 4 km.Indeed, the data, despite much spread, seems to change its downward trend at ∼10 4 km, as if we might in fact be seeing the expected behavior of the spherical-shell effect.Future studies of this effect would need to impose better controls over other variables, such as the portion of the propagation www.ann-geophys.net/29/851/2011/Ann.Geophys., 29, 851-863, 2011 path in daylight, etc.One useful step would be to conduct a study of this sort for lightning in which the peak current, or even better, the integrated power emitted, is independently ground-truthed.Then scatter from the random variation of source energy could be reduced.

Discussion
The overall success of the correlation between WWLLN and corrected VEFI times tends to confirm the approximate validity of our simplified propagation model responsible for the signal transfer from the source lightning to VEFI.In the low-magnetic-latitude conditions sampled by C/NOFS, the dominant propagation path seems to be the simplest, requiring no ducts: the oblique whistler coupled vertically from the Earth-ionosphere-waveguide's "leaky top" (Piggott et al., 1965;Pitteway, 1965).The simple propagation allows VLF signals from extremely distant lightning (∼10 4 km) to arrive at C/NOFS in a predictable manner.This simple sequence of propagation has already been observed (Holzworth et al., 1999); all we add in the present study is copious statistical support.We caution, however, that more complex propagation modes that have been precisely documented at higher magnetic latitudes by the DEMETER satellite (Chum et al., 2009), such as the "subprotonospheric (SP) whistler", may be present but may simply elude the simplistic search algorithm we have adopted.The SP whistler would undergo internal reflections between the ionosphere and the protonosphere, and would be associated with relatively nearby lightning.The work using DEMETER was a case study, and it is hard to see how to automate a search for SP whistlers in order to treat hundreds-of-thousands of candidate events.Our search algorithm penalizes multi-dispersive propagation and favors simple, single-mode whistler dispersion.More work will be required to determine whether, on the one hand, we are simply missing SP whistlers which are actually present due to lightning near the subsatellite point, or whether, on the other hand, the SP whistlers are not present at the low latitudes of C/NOFS, in keeping with early predictions of SP latitude dependence (Carpenter et al., 1964).We mention that, at least in mid-latitudes, the unducted, oblique whistler has already been shown to account for much of lightning-induced electron precipitation (Bortnik et al., 2006a, b;Lauben et al., 1999;Peter and Inan, 2004), so this seems to add to the original confirmation of the role of the unducted, oblique whistler (Holzworth et al., 1999).
The retrieved dispersion constant for the detected VEFI peaks will in a forthcoming publication be applied to a study of low-latitude ionospheric electron-density structure and climatology.

Fig. 1 .
Fig. 1.Spectral density of electric field on antenna pair 3/4 versus time, during the first ∼12-s-duration burst recording of 4 January 2010, with sampling rate ∼32 ksamples s −1 .The FFT window is 1024 samples wide, advanced by 128 samples.This record will be dissected into 45 windows, each containing 16 384 samples, and overlapping by 8192 samples, for dispersion and peak-finding analysis.

Fig. 2 .
Fig. 2. Example of dispersion and peak-finding analysis, using the 21st window of the record in Fig. 1 above.The window contains 16 384 samples, of which the first 8192 are set to zero before processing.The data in this modified window is then Fourier-transformed, pre-whitened, and dechirped against an oblique-whistler dispersion relation.(a) Figure-of-merit for dechirp vs. preliminary grid of 101 dispersion constants, from 0 to 100 s 1/2 in steps of 1 s 1/2 .(b) Iterated figure-of-merit for a grid of 101 values of the dispersion constant, from 80 % to 120 % of the optimum found in the first scan.(c)The electric field undergoes a direct Fast Fourier Transform (DFFT), is prewhitened in the frequency domain; Shown here is the spectral density of the prewhitened field, with a moving 256-sample FFT.Most of the energy is missing from the first half of the window, which had been zeroed prior to the first DFFT; there is however some leakage due to the prewhitening.(d) Similar, but with the electric field dechirped for the optimal dispersion measure (36.4 s 1/2 ).The peaks in this window (and all other windows) are determined by the peakfinder algorithm; three peaks are found in this window.The first peak in this window is the 45th peak (out of 98 total) in the overall record.

Fig. 3 .
Fig. 3. Time series for the 45th peak (the first above in Fig. 2d).There are 256 samples for this peak in the distilled peak archive.(a) Smoothed square of electric field vs. time for the dechirped peak.(b) Smoothed square of electric field after rechirping via DFFT on these 256 samples of electric field.(c) Dechirped electric field.(d) Rechirped electric field.

Fig. 5 .
Fig. 5. Preliminary figure-of-merit vs. dispersion constant, for all 45 analysis windows within the burst recording of Fig. 1.Curves for successive windows are displaced downward in each column.Any window marked "X" fails the figure-of-merit test.

Fig. 6 .
Fig. 6.Optimal dispersion constant for each of the 98 peaks selected within the burst recording of Fig. 1, vs. index of peak.Multiple peaks within the same analysis window forcedly have the same dispersion constant.(a) Dispersion constant, full scale on vertical.(b) Dispersion constant, zoomed to 34-42 (s 1/2 ) on vertical, showing variations during burst recording.

Fig. 8 .
Fig. 8. Similar to Fig. 1, but for the burst recording of 4 January 2010 yielding the most peaks for that day.

Fig. 9 .
Fig. 9. Similar to Fig. 5, but for the burst recording of Fig. 8.Only four (marked "X") of the 46 windows fail the figure-of-merit test.

Fig. 10 .
Fig. 10.Similar to Fig. 6, but for the burst recording of Fig. 8.An oval marks a sequence of five consecutive analysis windows apparently showing a sustained trend in dispersion.

Fig. 11 .
Fig. 11.Iterated figure-of-merit vs. dispersion constant for the burst recording of Fig. 8.Each trace is for a successive analysis window, and only the 42 windows containing at least one selected peak are shown.The five traces marked with an asterisk correspond to the sustained trend in dispersion seen in Fig. 10 above.

Fig. 12 .
Fig. 12.Similar to Fig. 1, but for a burst recording of 4 January 2010 yielding no peaks.

Fig. 13 .
Fig. 13.Similar to Fig. 5, but for the burst recording of Fig. 12.Each of the 46 successive analysis windows fails the figure-of-merit test.Some windows (e.g., the lowest in the left column) show evidence of whistler dispersion, but the competing non-dispersed activity still causes the dispersion figure-of-merit test to fail.

Fig. 14 .
Fig. 14.Histograms of VEFI peak time minus WWLLN stroke time, after correction for WWLLN stroke's propagation to the satellite in two path segments: First, from the stroke location to the subsatellite point in the Earth-ionosphere waveguide, at speed c, and second, vertically upward to the satellite as an oblique electron whistler.(a) Including WWLLN stroke locations out to 10 000 km from the subsatellite point.(b) Including WWLLN stroke locations out to 1000 km from the subsatellite point.The histogram bins are 1-ms wide.

Fig. 15 .
Fig.15.Histograms of arc distance from stroke location to subsatellite point.Each histogram is normalized to unity peak for ease of comparison.Solid curve: all WWLLN strokes.Dashed curve: only WWLLN strokes whose propagation-corrected arrival time is coincident (−3 to 0 ms in Fig.14above) with the VEFI peak.

Fig. 16 .
Fig. 16.Scatter plot of the integral of E 2 at the satellite vs. arcdistance to the coincident WWLLN stroke location.