First Map of Coherent Low‐Frequency Continuum Radiation in the Sky

Lightning discharges and radio transmitters emit low‐frequency (∼3–300 kHz) electromagnetic waves with large electric field strengths and stable phases. This phase stability makes it possible to map the source locations of lightning and transmitters in the sky. Electromagnetic waves with smaller electric field strengths generally exhibit a reduced phase stability, caused by numerous simultaneous physical processes that blend into an underlying continuum radiation trapped inside the Earth‐ionosphere cavity. It is therefore currently not known whether the source locations of continuum radiation can be determined. Here we show the first map of coherent continuum radiation in the sky above an array of high‐precision radio receivers. The source locations of the coherent continuum radiation are found at elevation angles ∼30∘−60∘ above the horizon. The identified source locations are attributed to intermittent radio transmitters that emit electromagnetic waves with electric field strengths ∼2 orders of magnitude below the instrumental noise floor. The results demonstrate that it is possible to simultaneously map the signals from coherent continuum radiation, lightning discharges, and radio transmitters in the sky. This work thereby lays the foundation toward the discovery of many more coherent source locations of low‐frequency electromagnetic waves in the sky. It is expected that the identified source locations vary with time as a result of the impact of solar variability on the D‐region ionosphere. Future studies have therefore the potential to contribute to a novel remote sensing and an improved understanding of the D‐region ionosphere, influenced by the near‐Earth space environment.


Introduction
The Stanford ELF/VLF Radio Noise Survey (Fraser-Smith, 1995;Fraser-Smith & Bowen, 1992, and references therein) used instantaneous amplitudes of low-frequency radio waves (Chrissan & Fraser-Smith, 2003; see also Chrissan & Fraser-Smith, 2000) to characterize atmospheric noise, which is a source of interference for radio transmissions and remote sensing (Spaulding, 1995). This atmospheric noise is caused by lightning discharges that occur around the world and therefore reflect the global thunderstorm activity and its diurnal and seasonal variability (Chrissan & Fraser-Smith, 1996). The amplitude variations of low-frequency radio transmissions are also used to monitor the electromagnetic wave propagation in the Earth-ionosphere cavity (e.g., Inan et al., 2010; see also Barr et al., 2000, and references therein). Such kind of measurements serve to remote sense the ionized state of the upper atmosphere, that is, the lower D-region ionosphere, (e.g., Ohya et al., 2018; see also Clilverd et al., 2012;Higginson-Rollins & Cohen, 2017;Macotela et al., 2017, and references therein), which responds to space weather, solar variability, and ionospheric modification by lightning discharges (e.g., Shao et al., 2013; see also Haldoupis et al., 2012;Kotovsky et al., 2016;Smith et al., 2004,, and references therein). The phase of radio transmissions is less often used to determine the ionisation of the lower ionosphere, even though it is an important observation to constrain and enhance the results obtained from amplitude measurements alone (e.g., Thomson et al., 2017; see also Thomson, 2010;Thomson et al., 2011, 2014, andreferences therein). The main reason why phase measurements are less common is that the associated signal processing has a larger degree of complexity (e.g., Gross et al., 2018;Koh et al., 2018, and references therein). More recently, phase measurements of lightning discharges have become of increasing interest, partly to determine the impact of lightning on the lower ionosphere (McCormick et al., 2018) or for their use in lightning detection networks ; see also Dowden et al., 2002;Liu et al., 2016). In particular, three-dimensional lightning detection networks, commonly referred to as lightning interferometers, have made extensive use of the cross-correlation function, which implicitly includes phase information (e.g., Shao et al., 2018; see also Abbasi et al., 2018;Lyu et al., 2016;Rison et al., 2016;Stock et al., 2014;Wu et al., 2014, and references therein). Explicit use of phase information is more commonly used in radar applications, for example, for meteor detection (Younger & Reid, 2017). The phase stability of low-frequency electromagnetic waves from radio transmissions, lightning discharges, and electromagnetic waves generated in near Earth space has been relatively little studied. For example, the phase stability of radio transmissions for submarine communication have been investigated in detail during the absence of solar variability (Thomson et al., 2017). Similarly, the exemplary phase stability of chorus wave packets in near Earth space was determined on board of Van Allen Probe spacecraft by explicit use of the instantaneous frequency (Santolik et al., 2014). More recently, the phase stability has been used to map the source locations of electromagnetic signals in the sky ; see also Füllekrug et al., 2014). This phase stability is an important prerequisite to successfully map well-focused source locations in the sky, in particular for the coherent signals of continuum radiation that are found up to 2 orders of magnitude below the noise floor of the radio receivers used . The aim of this work is to map this coherent continuum radiation in the sky, based on a detailed investigation of the phase stability of electromagnetic waves from lightning discharges and radio transmissions. The organization of this paper is as follows: Section 2 develops the theoretical framework to describe phase stability, section 3 shows exemplary measurements of the phase stability for lightning discharges and radio transmissions, section 4 describes the statistical analysis of the electric field strengths of coherent continuum radiation, and section 5 demonstrates how coherent continuum radiation is mapped in the sky above the array.

Theory of Phase Stability
The time-dependent complex trace (Taner et al., 1979) recorded at the n'th radio receiver in the array is composed of the electromagnetic source field E(t) that propagates as an electromagnetic wave e −i(k(t)r n − t) with the wavenumber k(t) and radian frequency across the geographic receiver locations r n such that (1) The electromagnetic source field is composed of the time-dependent amplitude A(t) and phase (t), which is calculated from an inversion of equation (1) by use of the recordings with N receivers in the array The complex temporal coherency t of the electromagnetic source field is calculated from the individual samples l = 1, 2 … L, which are not necessarily sequential where l is the instantaneous phase of the electromagnetic wave at the sample l. The temporal coherency of electromagnetic waves is often quite large such that it is convenient to use the temporal quality instead

Instantaneous Frequency
The instantaneous phase is a product of the instantaneous, or apparent, radian frequency a (t) and time t where f a is the instantaneous, or apparent, frequency and the first derivative of the apparent radiant frequency with respect to time is neglected in this first-order approximation. The instantaneous frequency can be calculated explicitly from the phase of the electromagnetic source field

Radio Science
10.1029/2018RS006705 such that (Taner et al., 1979). For the numerical implementation of the temporal derivatives for the real and imaginary part of the electromagnetic source field, the symmetric difference quotient is used. For example, the discrete derivative of the real electromagnetic source field is calculated numerically from where the sample time t l = lΔt is the product of the sample number l and the sampling time interval Δt = t l+1 − t l . The main argument for choosing the symmetric difference equation for a numerical implementation of the temporal derivatives is given by the following reasoning: If the first-order difference quotient would be used, the time series of the instantaneous frequency would be shortened by one sample. On the other hand, if the accuracy of the derivative is increased by using higher-order terms with samples that are more distant from the sample for which the derivative is calculated, the accuracy of the derivatives at the beginning and at the end of the time series differs from those in the rest of the time series. As a result, the symmetric difference equation is considered to be a pragmatic and computationally efficient compromise for the numerical calculation of the instantaneous frequency. In this case, the instantaneous frequency f a is finally calculated from It is interesting to note that equation (9) assigns an instantaneous frequency to each amplitude and phase calculated for each sample of the time series. As a result, each sample carries an instantaneous amplitude, phase, and frequency, similar to the spectral coefficients of a Fourier transform. In this sense, one might think about this methodology as some sort of Fourier transform for each individual sample.

Instantaneous Frequency Distribution
The instantaneous frequency is calculated for each individual sample from the measured electromagnetic source field (equation (9)). As a result, an ensemble of many instantaneous frequencies can be characterized by its distribution function. For example, radio transmissions use phase modulation schemes to convey information such that the distribution of instantaneous frequencies results from the particular phase modulation scheme used. For continuum radiation, the phase is random in nature and it is interesting to identify the corresponding distribution of instantaneous frequencies. In first-order approximation, the instantaneous frequency results from a ratio of normally distributed random numbers, which is Cauchy distributed (Johnson et al., 1994, chapter 16, p. 318). In this case, the instantaneous frequency distribution is where f 0 is the center frequency and s is the half width half maximum of the distribution. The maximum of the distribution is located at This maximum drops to half of its value at The center frequency f 0 and width 2s are commonly used to describe the Cauchy distribution, which has no converging statistical moments such as a mean and standard deviation. The width of the distribution encapsulates 50% of all the measured instantaneous frequencies in the interval ranging from f 0 − s to f 0 + s. The approximation of the instantaneous frequency distribution with the Cauchy distribution has been tested by using normally distributed random numbers as an input to the array processing. It is found that the Cauchy distribution approximates the instantaneous frequency distribution almost exactly, even though neither the nominator nor the denominator in equation (9) are normally distributed. This means that the Cauchy distribution varies relatively little when random numbers that are not normally distributed are used in the nominator and denominator, in agreement with theory (Johnson et al., 1994, chapter 16, p. 319).

Understanding of Instantaneous Frequency
The physical meaning of the instantaneous frequency is twofold: (1) The instantaneous frequency can be regarded as the momentary frequency of an electromagnetic wave. If the electromagnetic wave is down converted by multiplication with the time-varying part of the plane wave exp(−i t) in equation (1), then the instantaneous frequency is the frequency offset with respect to the center frequency of the electromagnetic wave.
(2) The instantaneous frequency conveys information on the phase stability of an electromagnetic wave. For example, a constant positive or negative instantaneous frequency offset from 0 corresponds to a constantly changing phase. When the instantaneous frequency fluctuates around 0, the phase can be regarded as relatively stable. In other words, the smaller the offset of the instantaneous frequency from 0 and the smaller its fluctuation around 0, the more stable is the phase. In this sense, one can think about the instantaneous frequency as an indicator for the phase stability of an electromagnetic wave. It is noted that electromagnetic waves are commonly composed of a continuum of numerous frequency components. In this case, the instantaneous frequency reflects the dominating frequency content of the electromagnetic wave, which is time dependent, in the sense of an instantaneous dynamic spectrum. The detailed analysis of instantaneous frequencies and phase stability thereby offers scientific insights into the physics of electromagnetic waves, which are illustrated with some examples in the following section.

Measurements of Phase Stability
This section aims to illustrate the phase stability of lightning discharges and radio transmissions by use of the instantaneous frequency as described in section 2. We make use of recordings with an array of N = 10 radio receivers distributed over an area ∼1 km × 1 km on Charmy Down airfield near Bath in the United Kingdom, which is well documented (Füllekrug et al., 2014;Mezentsev & Füllekrug, 2013). The electric field strengths are measured with a sampling frequency of 1 MHz, an amplitude resolution of ∼25 μV/m, and a timing uncertainty ∼20 ns (Füllekrug, 2010). The measurements from 15:00:00 to 15:00:10 UT on 13 May 2011 are extracted from the original recordings for an exemplary analysis of the phase stability conveyed by the instantaneous frequency. The seemingly short duration of 10 s is chosen here because it enables an analysis without the need for high-performance computing. The measurements are subsequently filtered for lightning pulses from 2-18 kHz and intermittent radio transmissions, that is, the U.K. radio clock from 59 to 61 kHz and LOng-Range Aid to Navigation (LORAN) transmitters from 90 to 110 kHz ; Table S1).

Phase Stability of Lightning Pulses
For long-range lightning detection networks, the phase of the maximum pulse amplitude is often used to determine the lightning occurrence time (e.g., Liu et al., 2018, and references therein) because the rising and falling edge of lightning pulses do not exhibit the same phase stability. For example, the instantaneous frequency of lightning pulses tends to be larger during the rising edge of the pulse when compared to the falling edge of the pulse (Figure 1a). In this case, the instantaneous frequency is always positive, that is, larger than the center frequency of 10 kHz. Other lightning pulses exhibit negative instantaneous frequencies throughout the pulse ( Figure 1a). The polarity of the instantaneous frequency determines the direction of the phase progression, which is either clockwise or anticlockwise ( Figure 2, right, Füllekrug et al., 2016). In these circumstances, the phase stability depends on the total amount of the phase progression during the entire lightning pulse. It is currently not known with certainty whether this phase stability is a property of the causative lightning discharge or a consequence of the electromagnetic wave propagation within the Earth-ionosphere cavity. However, it is speculated that the distance between the lightning discharge and the array contributes significantly to the observed instantaneous frequencies .

Relation Between Phase Stability and Amplitudes of Lightning Pulses
It was recently reported that the instantaneous frequency of lightning pulses implicitly depends on the instantaneous amplitude because the phase is better defined during relative maximum amplitudes when compared to relative minimum amplitudes . We extend this analysis here by comparing the instantaneous frequencies of all the measured relative amplitude maxima and minima of the amplitude envelope. It is found that the majority of instantaneous frequencies for amplitude maxima range within the physically The amplitude envelope of electric field strengths exhibits lightning pulses superimposed on continuum radiation. The instantaneous frequencies of lightning pulses and continuum radiation mainly range from 2 to 18 kHz (colored line). The rising edges of the pulses tend to have larger instantaneous frequencies than the falling edges. The inset shows a different lightning pulse with instantaneous frequencies that are mainly smaller than the center frequency of 10 kHz. (b) The relative maxima of the amplitude envelope have a distribution of instantaneous frequencies that mainly scatter between ∼2 and 18 kHz around the center frequency of 10 kHz (black). The distribution of instantaneous frequencies calculated from relative minima exhibit a much larger scatter (gray), which extends significantly beyond ∼2-18 kHz (dashed lines). The inset shows the distribution functions of relative amplitude maxima (black) and relative amplitude minima (gray). The relative minima extend well into the continuum radiation below the noise floor of the radio receiver.
meaningful interval ∼2-18 kHz defined by the filter used during the signal processing. For relative amplitude minima, the instantaneous frequencies vary more widely, that is, well beyond the boundaries of this interval ( Figure 1b). In other words, the phase stability is indeed larger for relative amplitude maxima when compared to relative amplitude minima. The distribution function of relative amplitude maxima is roughly log normal, whereas the relative amplitude minima are almost uniformly distributed below the instrumental noise floor ∼25 μV/m (Figure 1b, inset). As result, the distribution function of all the measured instantaneous amplitudes can be explained with a Rayleigh distribution for small amplitudes up to the instrumental noise floor ( Figure 2a). To validate this result, normally distributed random numbers are used as an input to the array processing. This simulated amplitude distribution follows almost exactly the entire Rayleigh distribution. The remaining small differences are attributed to the fact that the real and imaginary part of the simulated electromagnetic source field are not normally distributed, as pointed out in the discussion of the Cauchy distribution (section 2.2). Above the instrumental noise floor, the amplitudes of the electromagnetic source field do not follow a Rayleigh distribution any more, which indicates that the amplitudes convey physically meaningful information. It is proposed here that up to ∼92% of the observed amplitudes could possibly result from lightning discharges. Only ∼8% of the measured amplitudes might have a different origin, which is more random in nature.

Instantaneous Frequencies of Lightning Pulses
The observation that the vast majority of instantaneous amplitudes could be physically meaningful strongly suggests that the distribution function of the measured instantaneous frequencies conveys information on the phase stability of lightning pulses. To quantify this phase stability, the distribution for instantaneous frequencies of lightning pulses with spatial qualities q s > 1.0 is calculated. This sample encompasses ∼32.5% of all the validated instantaneous frequencies from 2 to 18 kHz. It is found that the resulting instantaneous frequency distribution peaks near ∼11 kHz ( Figure 2b). As a result, the lightning pulses are consistently offset by ∼1 kHz from the center frequency of ∼10 kHz used for the signal processing, in agreement with the instantaneous frequencies of the relative amplitude maxima, which exhibit a similar offset (Figure 1b). This observation means that lightning pulses are more likely to exhibit a positive polarity of the instantaneous frequency after down conversion. To ensure that this result is not a consequence of the array processing, the distribution for instantaneous frequencies of lightning pulses with spatial qualities q s < 0.5 is calculated. This sample encompasses ∼32.7% of all the valid instantaneous frequencies from 2 to 18 kHz. The resulting reference distribution is centered on ∼10 kHz, it does not show any positive offset when compared to the The instantaneous frequencies for the radio clock at 60 kHz (black) are almost perfectly normal distributed (red). The inset shows that the first deviations from the normal distribution occur at the 3 level near ±3 Hz.
distribution of the instantaneous frequencies from lightning pulses, and it follows almost exactly a Cauchy distribution with a width of ∼4.4 kHz as described in section 2.2 ( Figure 2b). The difference between the peaks of the observed instantaneous frequency distributions at 11 and 10 kHz therefore indicates that the observed instantaneous frequencies of lightning pulses are physically meaningful. These instantaneous frequencies of lightning pulses possibly depend on the average distance between the array and the lightning discharges . In summary, the phase stability of lightning pulses is determined by the instantaneous frequencies that exhibit a significant offset from the center frequency chosen for the signal processing. It is expected that an adjustment of the center frequency to ∼11 kHz would remove the observed offset while the width of the distribution remains unchanged.

Phase Stability of Radio Transmissions
Radio transmitters often use phase modulation schemes to convey information such that the phase stability is relatively small unless the phase modulation is removed from the original recordings (e.g., Koh et al., 2018, and references therein). For example, submarine communication transmitters tend to use a phase modulation 10.1029/2018RS006705 scheme that is equivalent to instantaneous frequency changes ±50 Hz (Figure 2c). This variability of the instantaneous frequency corresponds to a relative stability of ∼50 Hz/22.1 kHz = 0.23%. This is relatively large when compared to the U.K. radio clock, which exhibits a relative stability ∼3 Hz/60 kHz = 5 ⋅ 10 −3 %, which is equivalent to ∼50 parts per million at the 3 level (Figure 2d), where the distribution starts to deviate from a normal distribution. It is possible to express this phase stability as a timing uncertainty as described in the following section.

Phase Stability and Timing Uncertainty
The nominal time period T 1 to complete one cycle of the carrier frequency is compared to the time period T 2 which is influenced by a relatively small deviation of the instantaneous frequency f a = f c + f a from the carrier frequency f c by f a , such that The corresponding timing uncertainty ΔT = T 1 − T 2 is then given by after an expansion of the quotient in first-order approximation. The resulting timing uncertainty for the radio clock is then ΔT ≈ 3 Hz/(60 kHz) 2 = 0.8 ns at the 3 level. This timing uncertainty is relatively small when compared to the timing accuracy of the radio receiver ∼25 ns. This is so, because the measurement of the instantaneous frequency is based on the phase derivative, which corresponds to a determination of the timing uncertainty that accumulated between two consecutive samples within 1 μs. This timing uncertainty is much smaller than the accumulated timing uncertainty ∼25 ns of the GPS clock, which was determined between two consecutive GPS pulses separated by 1 s (Füllekrug, 2010). As a result, it is not easily possible to distinguish whether the observed timing uncertainty of the radio clock results from the transmitter, the electromagnetic wave propagation from the transmitter to the receiver, or the GPS clock, which disciplines the radio receiver. But it is possible to use the calculated timing uncertainty as a lower limit for the timing uncertainty of the radio receiver. For example, the timing uncertainty observed for the submarine communication transmitter (Figure 2c) and lightning pulses (Figures 1b and 2b) are clearly larger. For example, the equivalent timing uncertainty for the submarine communication transmitter is ΔT ≈ 10 Hz/(20 kHz) 2 = 25 ns. This timing uncertainty thereby indicates that it has a physical cause because it is much larger than the lower limit of the timing uncertainty of the radio receiver, as determined by the radio clock. It is therefore possible to conclude that the observed phase stability for the submarine communication transmitter is either a consequence of the electromagnetic wave propagation and/or the timing accuracy of the transmitter. Similarly, the phase stability for lightning discharges ΔT ≈ 2 kHz/(10 kHz) 2 = 20 μs is either a consequence of the electromagnetic wave propagation and/or the time-dependent lightning current.

Phase Stability and Temporal Quality
It was shown in the previous section that the phase stability of the measurements is relatively large and the influence of the timing uncertainty of the GPS clock is relatively small. It is therefore possible to investigate the phase stability of intermittent radio transmissions using the temporal quality (equation (4)), which is a measure of the complex temporal coherency of electromagnetic waves (equation (3)). Here we make use of LORAN transmitters for marine navigation at 100(±10) kHz, which transmit repetitive pulse sequences Füllekrug et al., 2009; see also Mezentsev & Füllekrug, 2013, and references therein). For example, the pulse sequences emitted by the four strongest LORAN transmitters exhibit various instantaneous frequencies (Figure 3a). Most of the pulses exhibit instantaneous frequencies that are ∼0-3 Hz below 100 kHz and the instantaneous frequencies of one pulse sequence vary between ∼0 and 1 Hz above 100 kHz. These small deviations from the carrier frequency are comparable to those observed for the U.K. radio clock and indicate a large phase stability. It is noted that the rising and falling edges do not exhibit larger and smaller instantaneous frequencies as observed for the lightning pulses (compare to Figure 1a). The phase stability of the radio transmissions is evaluated by calculating the instantaneous frequency from the complex temporal coherency for 876 repetitions of the same pulse sequence and a subsequent calculation of the temporal quality. This temporal quality varies from ∼3 down to continuum radiation with a temporal quality of ∼0.5 (Figure 3b). The decrease of the temporal quality with each pulse sequence mainly results from the distance between the transmitters and the receivers, but it also depends on the quality of the transmission. For example, the third pulse sequence starts with two pulses with large temporal qualities, followed by six pulses with smaller temporal qualities. The same effect is observed in a less pronounced way for the fourth pulse  Figure 3a from ∼33 to 42 ms. (d) The temporal quality of coherent LORAN transmissions is constant with increasing integration steps, while the temporal quality of incoherent radiation decreases with each integration step. The LORAN pulses from ∼19 to 26 ms become apparent after ∼150-200 integration steps. The secondary pulses from ∼0 to 10 ms with significant temporal qualities result from sky waves.
sequence. The reason for this difference in temporal qualities is that the two leading pulses are fixed in time, while the six consecutive pulses are time shifted by ±1 μs to convey some information. This imposed time shift reduces the temporal quality of the transmissions.

Phase Stability of Coherent Continuum Radiation
The phase stability described in the previous section can be used to detect intermittent radio transmissions, which are part of the continuum radiation with amplitudes that are normally too small to be detected because they fall below the instrumental noise floor. In this case, the temporal coherency is used for the detection of the radio transmissions in the first instance. For example, three LORAN transmitters are detected by calculating the temporal quality from 654 repetitions of three pulse sequences (labeled 5-7 in Figure 3c). The instantaneous frequency is calculated from the complex temporal coherency, which has the same phase as the electromagnetic source field (equation (3)). It is found that all three transmissions exhibit a surprisingly large phase stability, even though the corresponding temporal qualities are very small. It is interesting to determine the minimum number of repetitions of the pulse sequence needed to enable a detection of the weakest LORAN transmission (labeled 6 in Figure 3c). This aim can be achieved by successively increasing the number of integration steps L in equation (3). It is found that the first and third pulse sequences (labeled 5 and 7 in Figure 3c) can be identified after ∼5-15 and ∼20-30 repetitions, respectively (Figure 3d). For the

10.1029/2018RS006705
identification of the second transmission with the smallest temporal quality (labeled 6 in Figure 3c), ∼150-200 integration steps are required for a reliable detection. It is noted that the first pulse sequence exhibits secondary pulses, which result from sky waves that arrive at a later time. Overall, the results demonstrate the extraordinary phase stability of coherent continuum radiation. The corresponding amplitudes of the electric field strengths below the noise floor of the receiver are determined by use of a thorough statistical analysis described in the following section.

Statistical Analysis of Coherent Continuum Radiation
The main challenge for the determination of electric field strengths as part of continuum radiation is that the amplitudes are smaller than the noise floor of the radio receiver. It is common practice to estimate the reduction of the noise floor of the instrument 0 by using the standard deviation of the mean, which is determined from the number of repetitions L used for calculating the average of a measurement of interest. In this case, the reduced noise floor of the instrument is L = 0 ∕ √ L. For example, the standard deviation of the mean for L = 654 repetitions with an instrumental noise floor of 0 = 25 μV/m is L ≈ 0.98 μV/m. In other words, an impulse with an electric field strength of 9.8 μV/m would not exceed the instrumental noise floor, but it would exceed the standard deviation of the mean by a factor of 10. This reasoning can be used to determine how many samples are needed to determine a mean value , which exceeds the standard deviation of the mean by the factor For example, to detect an average electric field strength of = 0.6 μV/m in the presence of an instrumental noise floor of 0 = 25 μV/m with a significance of = 5 standard deviations, more than L = 43, 403 repetitions of the electric field strengths would be required. This number is relatively large when compared to the ∼150-200 repetitions required for the detection of coherent continuum radiation described in section 3.2.3 (Figure 3d). In addition, this first-order approach only works if the average converges to the true mean value which is not necessarily guaranteed in the presence of strong interference. Therefore, a more robust estimation of the mean and standard deviation is required as described in the following section.

Resampling Coherent Continuum Radiation
It is suggested here to resample the electric field strength distribution, which offers a more robust estimate of its statistical moments, such as the mean and the standard deviation. For example, the jackknife method uses randomly chosen subsets from the original data to calculate statistical moments (e.g., Tukey, 1958, and references therein). In addition, the bootstrap method allows to use the original data multiple times during the resampling process (e.g., Efron, 1979, and references therein). Here we resample the electric field strength distribution with the bootstrap method after mitigating the influence of exceptionally large values arising from interference. The smallest 68.27% absolute electric field strengths are extracted from the original data. The chosen percentage roughly corresponds to the electric field strengths that lie within ±1 of the original distribution. In this way, outliers caused by strong interference are excluded from the analysis at the outset. Subsequently, the resampled mean and standard deviation are calculated by use of the discrete Fourier transform. The discrete Fourier transform is used here because it offers a computationally efficient way to calculate numerous standard deviations from the spectral coefficients in a single transform. More specifically, the spectral coefficients c m at the harmonic numbers m are calculated from a discrete Fourier transform of the physically meaningful real part of the electromagnetic source field E r where n is the sample number, N is the total number of samples, and √ 2∕N is a normalization factor, which is chosen such that the real and imaginary part of each spectral coefficient is an estimate of the standard deviation at the harmonic number m.

Resampled Mean
The resampled mean is obtained from the spectral coefficient c 0 at the harmonic number m = 0 by using the normalization factor 1∕N

10.1029/2018RS006705
The resulting estimate of the resampled mean̂= c 0 is normally distributed around the true mean of the original distribution function with the standard deviation of the mean = 0 ∕ √ N.

Resampled Variance
Similarly, the resampled variancê2 is calculated from the real and imaginary parts of the spectral coefficientŝ where M = M 1 + M 2 is the number of the randomly chosen real and imaginary parts c m 1 ,r and c m 2 ,i of the spectral coefficients at the harmonic numbers The resulting estimate of the resampled variance where Γ is the gamma function which is calculated here in its numerical form (Cody, 1976).

Resampled Distribution Functions
The resampling of the distribution function with the bootstrap method makes it possible to repeat the resampling process as many times as required to obtain a smooth distribution function f (̂, , ) to determine the resampled mean ( Figure 4a) and a smooth distribution function f (X, ) for the resampled variance ( Figure 4b).
The true mean and variance are given by the most likely values, or modes, of the respective resampled distribution functions. The numerical accuracy of the resampled mean and variance is determined by the chosen bin size used to calculate the resampled distributions. The smoothness of the resampled distributions is determined by the number of realisations used during the resampling. It is important to keep in mind that the resampled mean and variance are not direct observations. The resampled mean is effectively a measure for the small offset of the instrumental noise floor away from 0. Similarly, the resampled variance is effectively a measure of the uncertainty for the resampled mean. The main advantage of resampling the mean and variance is that the resulting distribution functions can be described almost exactly by a normal distribution for the mean and a 2 distribution for the variance. This is so because the instrumental noise floor of the radio receivers is almost exactly normally distributed (Füllekrug, 2010, Figure 5). Any significant deviations from the normal distribution and the 2 distribution would indicate that the resampled mean and variance do not describe electric field strengths as part of continuum radiation below the instrumental noise floor. For example, the normal distribution of the resampled mean could be evaluated at ±1 , ±2 , and ±3 , which encompass ∼68.27%, ∼95.45%, and ∼99.73% of the data. Similarly, the 2 distribution of the resampled variance could be evaluated at the mode, mean, and variance, that is, at ( −2) 2 , 2 , and 2 4 , which encompass ∼37%, ∼58%, and ∼95% of the data, in agreement with the theory developed for a unit variance (e.g., Kendall & Stuart, 1977, p. 398). Significant experimental deviations from these theoretical relationships could be used to identify electric field strengths, which are not part of continuum radiation. A more traditional approach would be to fit a normal distribution and a 2 distribution to the resampled distributions. In this case, a merit function, for example, the mismatch, could be used for the identification of electric field strengths, which are not part of continuum radiation. However, this approach favors fitting the peak values of the distributions and relatively large deviations from the smaller values in the wings of the distributions would have little, if any, influence on the quality of the fit. However, the main result of this section is that it is possible to quantify the mean and variance of coherent continuum radiation below the instrumental floor by calculating resampled distribution functions with the bootstrap method. This result means that it will be possible to calculate the average energy flux, or Poynting vector, of continuum radiation in future studies. The critical first step toward this long-term objective is to map the wave number vectors of coherent continuum radiation in the sky, which is the subject of the following section.

Mapping Coherent Continuum Radiation in the Sky
The previous sections have clearly shown that coherent continuum radiation has a remarkable phase stability (section 3), which can be used to calculate the electric field strengths of coherent continuum radiation below the instrumental noise floor (section 4). This promising result strongly suggests that it is possible to map the wave number vectors of coherent continuum radiation in the sky. These wave number vectors can be calculated from equation (1) by use of the complex temporal coherency n (t) at the n'th radio receiver in the array where n (t) is the phase of the complex temporal coherency measured at the n'th receiver in the array. It is possible to replace the complex trace by the complex temporal coherency because it is assumed here that the phase and wave number vector of the electromagnetic source field remain constant during the integration time, at least in first-order approximation. The transfer function between the complex temporal coherency n (t) and m (t) at the m'th radio receiver in the array is then The phase of this transfer function is used to formulate a linear system of equations to determine the wave number vector where x n,1 and x n,2 are the geographic coordinates of the radio receivers in the east and north direction, measured in m with respect to the sweet spot of the array, and n and m are the phase uncertainties of the complex temporal coherency.

Arrival Azimuth and Elevation Angle of Wave Number Vectors
The wave number vector measured in two horizontal dimensions (equation (26)) is determined here by a projection of the wave number vector in three dimensions onto the horizontal plane in two dimensions. This projection can be expressed in Cartesian or spherical coordinates where k 0 = ∕c is the wave number of an electromagnetic wave that sweeps across the array at the speed of light c, Φ is the arrival azimuth of the wave number vector with respect to the geographic east direction, and Θ is the elevation angle of the wave number vector measured from the horizon upward. Equation (27) can be used to determine the arrival azimuth and elevation angle of the wave number vector. For example, the arrival azimuth can be calculated from

10.1029/2018RS006705
Similarly, the elevation angle can be calculated from Note that the reference wave number k 0 = ∕c used here is chosen based on a best fit to the observed dispersion relation (Füllekrug, Mezentsev, et al., 2015, Figure 2, right). A more thorough assessment of the reference wave number can be achieved by using a three-dimensional array as discussed in the following section.

Wave Number Vector in Three Dimensions
For a three-dimensional array with significant geographic height differences x n,3 − x m,3 between individual radio receivers, it is possible to expand equation (26) into three dimensions such that all three components of the wave number vector are measured directly In this case, the elevation angle can also be calculated from the vertical component of the wave number vector The main advantage of a three-dimensional array is that the length of the wave number vector becomes an observable quantity which enables a direct measurement of the wave propagation velocity This is an important measurement because it was previously pointed out that the elevation angle determination in two dimensions (equation (29)) depends on the wave number k 0 of the electromagnetic wave. For example, in free space, k 0 = ∕c, where c is the speed of light. In the presence of the Earth's conductivity, it is more realistic to use k 0 = k m , where k m is the limiting wave number for the local medium k m = ∕v m , which is determined by the local wave propagation velocity v m . If confirmed by corresponding measurements, a significant difference between the elevation angles determined by a two and three-dimensional array (equations (29) and (32)) could indicate that the local wave number number vector deviates from the Poynting vector, for example, as a result of an anisotropy introduced by local conductivity anomalies (e.g., Bor et al., 2016, and references therein). As a result, it is interesting to consider the construction of a three-dimensional low-frequency array, for example, by placing radio receivers on prominent topographic features or unmanned aerial vehicles. Given that the wave number vector can be calculated for each sample from equation (30), a three-dimensional array would be able to measure the instantaneous wave propagation velocity for each individual sample by use of equation (34).

Mapping Wave Number Vectors
The arrival azimuth and elevation angle of the wave number vector calculated from equations (28) and (29) are mapped on the celestial sphere above the array to visualize the arrival direction of the wave number vector. Four areas with enhanced numbers of source locations can be distinguished (Figure 5a). These areas are related to the pulse sequences of the LORAN transmissions as part of the coherent continuum radiation discussed in section 3 (Figure 3c). Each pulse sequence that belongs to one individual transmitter is extracted

10.1029/2018RS006705
and mapped individually in the sky above the array (Figures 5b-5d). The first transmitter with the largest phase stability exhibits the most focused source locations that are located at the largest elevation angles ∼60 ∘ (Figure 5b). A secondary maximum occurs at lower elevation angles and exhibits arrival azimuths that are slightly shifted to the west. These two areas of source locations are attributed to sky waves, as the ground wave tends to get quickly absorbed when observed beyond the horizon during day time (e.g., Figure 9c in Mezentsev & Füllekrug, 2013, and references therein). The second transmitter with the smallest phase stability exhibits the largest scatter in the sky around elevation angles ∼40 ∘ (Figure 5c). The source locations of this transmitter were originally concealed by the sky waves of the first transmitter when all source locations were mapped simultaneously (compare Figure 5c to Figure 5a). The source locations of the third transmitter are found at elevation angles ∼35 ∘ -40 ∘ toward the east of the other two transmitters. The results clearly indicate that the size of the area covered by the source locations is inversely proportional to the phase stability of the electromagnetic waves. In other words, for the mapping of coherent continuum radiation with well-focused source locations in the sky, large phase stabilities are a desirable prerequisite. The overall phase stability is determined by the measured instantaneous amplitude, frequency, and phase and their corresponding uncertainties, which are discussed in the next section.

Assessment of Uncertainties 5.4.1. Uncertainty of the Instantaneous Amplitude
The instantaneous peak amplitudes of the LORAN transmissions vary from ±27−29 μV/m, ±0.5−0.7 μV/m, and ±2-4 μV/m (Füllekrug et al., 2018, Figures 4c and 4d) and correspond to temporal qualities ∼0.5, ∼0.1, and ∼0.2 (Figure 3c). The uncertainty of the instantaneous peak amplitude is determined by the mode of the resampled variance ±4.8 μV/m (Figure 4b), which is almost one order of magnitude larger than the resampled mean ±0.6 μV/m (Figure 4a). Alternatively, the uncertainty of the instantaneous peak amplitude could be characterized by the standard deviation of the resampled mean ±0.44 μV/m (Figure 4a). In either case, the uncertainty is determined by the time-independent instrumental noise floor. As a result, it is concluded that the temporal qualities of the LORAN transmissions, and hence their phase stabilities depend on their instantaneous peak amplitudes, in agreement with the results in section 3.1.1.

Uncertainty of the Instantaneous Frequency
The instantaneous frequencies of the LORAN pulses vary from 0-3 Hz below 100 kHz to 0-1 Hz above 100 kHz, which is similar to the phase stability of the radio clock (Figure 3c and section 3.2.2). The uncertainty of the instantaneous frequency can be characterized by the width of the corresponding Cauchy distribution (equation (12) and Figure 2b). However, the uncertainty of the instantaneous frequency of the LORAN transmissions has not been quantified in more detail because it is comparable to the accuracy of the radio clock, and it is relatively small when compared to the phase stability that results from the instantaneous amplitude (section 5.4.1).

Uncertainty of the Instantaneous Phase
The phase uncertainties that result from the instantaneous amplitude and the instantaneous frequency both contribute to the uncertainty of the wave number vectors that are calculated from equation (26). A large phase stability results in a small scatter of the wave number vector locations in the sky and small phase stabilities result in a large scatter of the wave number vectors locations in the sky. A more rigorous treatment of the phase stability shows that the covariance matrix of the wave number vector scatter can be calculated from the covariance matrix of the phase residuals (Füllekrug, Mezentsev, et al., 2015, equation (9)). This wave number scatter can be reduced significantly by use of the array response, which results in a sharpening of the sky map (Füllekrug, Mezentsev, et al., 2015, equation 14). This image sharpening is not applied here because the phase stability for coherent continuum radiation can also be improved by calculating averages over longer time intervals in the first instance (Figure 3d).

Increasing the Phase Stability
The phase stability can be increased by use of a longer integration time to calculate the temporal coherency. However, this strategy implicitly assumes that the wave propagation conditions within the Earth-ionosphere cavity do not change significantly during the integration time. The Earth-ionosphere cavity is bound by a fixed ground conductivity and a highly variable conductivity of the lower D-region ionosphere. In other words, it is expected that a most suitable integration time does exist which is (1) long enough to ensure a sufficient temporal coherency for mapping well-focused source locations in the sky and (2) yet short enough to avoid blurring of the sky caused by time-varying ionospheric conditions.

Summary and Conclusions
The phase stability of low-frequency electromagnetic waves is assessed for the first time in unprecedented detail with a rigorous analysis based on novel recordings with an array of radio receivers. It is found that the coherent continuum radiation with electric field strengths below the instrumental noise floor exhibits a remarkable phase stability, as measured by the temporal coherency, or temporal quality, of the electromagnetic waves. The phase stability of coherent continuum radiation increases in proportion to the integration time used to calculate the temporal coherency. The source locations of intermittent radio transmitters with electric field strengths up to 2 orders of magnitude below the instrumental noise floor are successfully mapped in the sky above the array. The size of the area in the sky covered by the identified source locations is found to be inversely proportional to the phase stability of the electromagnetic wave. These results clearly demonstrate that it is possible to simultaneously map the source locations of coherent continuum radiation, lightning discharges, and radio transmissions in the sky. This work thereby opens a novel avenue to quantify the energy fluxes, or Poynting vectors, for many more coherent sources of continuum radiation in the celestial sphere above the array. It is expected that these energy fluxes and their apparent source locations in the sky are time dependent. Such variability is imposed by the constantly changing wave propagation conditions controlled by the ionized upper atmosphere from ∼60-90 km height. This D-region ionosphere is influenced by upwelling neutral density waves and solar variability intruding from the near-Earth space environment. We therefore anticipate that this work has the potential to lead to novel remote sensing of the D-region ionosphere in the years to come and a corresponding improved understanding of the ionized upper atmosphere.