Detection of Low‐Frequency Continuum Radiation

The electromagnetic spectrum at low frequencies from ∼3 to 300 kHz is dominated by impulses from lightning discharges and anthropogenic radio transmissions used for communication. Electromagnetic waves generated in near‐Earth space exhibit generally smaller amplitudes that are attenuated when travelling through the ionosphere before they can be observed at high and midlatitudes. Electromagnetic waves with yet smaller amplitudes contribute to the overall electromagnetic energy trapped within the Earth‐ionosphere cavity. At this point, the electromagnetic waves from all possible sources blend into an unstructured continuum radiation near the instrumental noise floor, which is often considered to be a fundamental limit to scientific discovery. As a result, the sources of continuum radiation have been little studied and are essentially unknown. Here we show how low‐frequency continuum radiation is detected and discriminated against known radio sources and instrumental noise by use of rigorous criteria inferred from novel precision measurements with an array of radio receivers. In particular, it is found that coherent continuum radiation from intermittent radio transmitters exhibits electric field strengths 0.5–0.7 μV/m, which are almost 2 orders of magnitude smaller when compared to the noise floor of the radio receivers ∼25 μV/m. Another part of the continuum radiation is found at local zenith above the array when it is caused by random instrumental noise. The results exemplify the possibility to extract sources from continuum radiation to study their origin and physical properties, which can contribute to an improved understanding of the impact of space weather and solar variability on the Earth's upper atmosphere.


Introduction
The recordings of low-frequency radio waves are commonly described as signals, if they are wanted, and noise, if they are unwanted. For example, the Stanford ELF/VLF Radio Noise Survey (i.e., Fraser -Smith, 1995;Fraser-Smith & Bowen, 1992, and references therein) mainly reflects the average spectral amplitudes of intermittent lightning discharges (Chrissan & Fraser-Smith, 1996) measured around the world to characterize global thunderstorm activity (Füllekrug & Fraser-Smith, 1997). The wording noise survey is adopted from larger radio frequencies where lightning discharges are described as atmospheric noise that interferes with radio transmissions used for communication or remote sensing (Spaulding, 1995). On the other hand, lightning researchers often describe radio transmissions as noise (e.g., Füllekrug et al., 2015;Jacobson et al., 1999;Smith et al., 2004, and references therein) because anthropogenic radio transmissions limit the detection threshold for natural weak lightning processes such that the influence of disturbing radio transmitters needs to be mitigated. For example, in long-range lightning detection networks spectral filter techniques are applied to increase the detection efficiency (Dowden et al., 2002;Liu et al., 2016, 2018, andreferences therein). When electromagnetic signatures from intermittent lightning discharges and radio transmissions in specific frequency ranges are absent, or removed, small radio signals originating from wave particle interactions in near-Earth space can be observed at high latitudes and from satellites in space (e.g., LaBelle & Anderson, 2011;Manninen et al., 2016;Parrot et al., 2015;Titova et al., 2017, and references therein). In this context, the wording noise is used to describe the unstructured behavior of radio waves observed in dynamic spectra close to the noise floor of the radio receiver. Similarly, the subtle radio signals from sprite streamers in the middle atmosphere appear as irregularly spaced pulses in radio recordings described as low-frequency radio noise Qin et al., 2012). In general, the classification of low-frequency radio measurements in two categories named signal and noise is commonly based on recordings with a single radio receiver. Recent recordings with an array of 10 radio receivers have shown that small electromagnetic signatures are well correlated over an area of 1 km × 1 km   Figure 1) and it was suggested to use the spatial quality, 10.1029 or flatness, of low-frequency radio waves to discriminate lightning pulses from smaller electric field strengths with unknown physical origin .

Rationale
The main aim of this work is to advance the theory of low-frequency array analyses to enable a scientific treatment of small electric field strengths and to demonstrate with exemplary array measurements that this novel theory can be used to identify and quantify meaningful physical processes. The irregular behavior of small electric field strengths is specifically referred to continuum radiation throughout this contribution. A detailed investigation of the scientific merit of continuum radiation is thought to be important for three major reasons: (1) Continuum radiation is more common than specific events. An event normally describes only a brief excursion from continuum radiation. Therefore, it is important to identify the physical causes of continuum radiation which is more abundant than individual events. (2) The exclusion of continuum radiation from an analysis of signals of interest increases the quality of the results. A rigorous scientific discrimination between continuum radiation and signals of interest by use of the spatial quality, or flatness, of low-frequency radio waves serves to either detect continuum radiation for further analysis or to exclude continuum radiation from further analysis. (3) The wave number vectors of continuum radiation reflect the near-Earth space environment. The main scientific merit of array analyses is the determination of wave number vectors. For example, it was recently demonstrated that the D-layer ionosphere is a birefringent plasma for low-frequency electromagnetic waves (Füllekrug et al., 2015). This birefringence results in the splitting of wave number vectors, which provides detailed information on the reflection heights in the D-layer ionosphere.

Long-Term Objectives
Based on these promising results, it is expected that the arrival angles of wave number vectors will change over time and thereby reflect solar variability and space weather induced changes in the near-Earth space environment, which affect the D-layer ionosphere. The work presented here is intended to contribute toward laying the theoretical and experimental foundation for using low-frequency array analyses to remote sense the D-layer ionosphere for the benefit of near-Earth space research. The ionospheric remote sensing itself is not the subject of this contribution and remains one of the long-term objectives for the application of array analyses in the future.

Structure
The organization of this paper is as follows: Section 2 introduces the novel array analysis to treat small electric field strengths, section 3 describes the physical properties of electromagnetic source fields with small electric field strengths for comparison with intermittent lightning pulses and radio transmissions, and section 4 demonstrates how the statistical properties of small electric field strengths can be used to extract meaningful physical processes from the array measurements.

Array Analysis
The electromagnetic recordings Y n (t) at each receiver n = 1, 2 … N in the array are used to form the complex trace y n = Y n + i(Y n ), where the Hilbert transform (Y n ) shifts the recorded waveforms by ∕2 (Taner et al., 1979). The complex trace is explained by an electromagnetic source field E(t) = A(t)e i (t) with amplitude A(t) and phase (t) that propagates as a plane wave e −i(k(t)r n − t) with the real wave number vector k(t) to the geographic locations of the receivers r n in the east and north direction at the radian center frequency of interest such that This ansatz works well for large observed electric field strengths y n = |y n |e i( n + n ) when the observed phase n + n has a small uncertainty n . Smaller electric field strengths tend to exhibit larger phase uncertainties that can exceed the maximum possible phase difference Δ = | n − m |, between two receivers n and m in the array. This maximum phase difference is determined by the baseline b, i.e., the largest possible distance between two receivers in the array, such that Δ = Δt = b∕c, where Δt is the corresponding maximum time delay and c is the speed of light. In extreme cases, the observed phase difference can exceed the maximum possible phase difference |( n + n )−( m + m )| > Δ , or the phases n and m are found across the + ∕ − boundary of their principal value such that the observed phase difference can appear to be larger than it actually is. In either case, large phase uncertainties can occur and hinder an accurate determination of the wave number vector by solving the system of equation (1). As a result, it is suggested here to infer the wave number vector from the transfer functions that are calculated for mutual pairs of receivers.

Analysis of Small Electric Field Strengths
The complex trace at the receiver location m can be described by The transfer function between the complex traces in equations (1) and (2) is then whereŷ n (t) andŷ * m (t) are the complex traces normalised to unit amplitude and * denotes complex conjugation. The phase of the transfer function is finally used to formulate a system of equations that determines the wave number vector where x n,1 and x n,2 are the geographic coordinates in the east and north direction, respectively. The main merit of formulating this overdetermined system of equations is that the phase differences n − m are centered around zero within the physically possible range ±Δ . This constraint is used here to define a weight for the observed phase difference between two receivers with the raised cosine probability Large phase differences | n − m + n − m | that stray beyond the physically possible range Δ are thereby weighted less when equation (4) is solved for k.

Weighted Inversion
The system of equation (4) is solved with the Gaussian method of weighted least squares to minimize the phase uncertainties n,m = n − m , which cannot be explained by plane waves propagating across the array where W is the matrix with the weights from equation (5) on its diagonal, G is the matrix with the relative geographic coordinates, k is the wave number vector, and p is the vector with the measured phase differences. The weighting makes the inversion computationally expensive, because the weighted model matrix A = G T WG and the solution vector b = G T Wp need to be computed for each individual sample. As a result, it is more efficient to calculate the inverse matrix explicitly from where the matrix and vector elements are given by

Radio Science
such that the wave number vector is finally calculated from ( k 1 k 2 ) = 1 a 1,1 a 2,2 − a 1,2 a 2,1 ( a 2,2 −a 1,2 −a 2,1 a 1,1 Once the wave number vector is determined for each sample, equation (1) can be reformulated to solve for the electromagnetic source field by use of the complex Gaussian method of least squares The multiplication of the complex trace at each receiver y n (t) with the conjugate plane wave e i(k(t)r n − t) on the right-hand side of equation (16) means in practice that the phase of the complex trace n is shifted by the phase of the plane wave k(t)r n − t to the phase of the electromagnetic source field (t).

Spatial and Temporal Quality
It is now possible to define a spatial coherency of the plane wave across the array This spatial coherency of the plane wave is s = 1 when the phase uncertainties n do not vary across the array and the coherency becomes smaller when the phase uncertainties vary across the array. The spatial quality, or flatness, of the plane wave is then given by (Füllekrug et al., 2016, equation (4)). The spatial quality is large when the observed electric field strengths have a large spatial coherency and the spatial quality decreases when the spatial coherency becomes smaller. Similarly, a temporal quality of the electromagnetic source field can be defined by use of the temporal coherency where the index k denotes selected sample numbers, K is the total number of selected samples to calculate the coherency, and k is the phase of the electromagnetic source field. If the selected sample numbers are sequential E k = E(t k ) = E(kΔt), where the sampling time t k = kΔt is calculated from the sampling time interval Δt. It is noted, however, that the selected sample numbers are not necessarily sequential. In either case, the temporal quality is

10.1029/2018RS006654
This temporal quality of the electromagnetic source field is large when the phase of the electromagnetic source field k does not vary for the selected samples and the temporal quality decreases when the phase varies for the selected samples.

Significance of Spatial and Temporal Quality
The significance of the spatial and temporal quality of low-frequency radio waves can be assessed by comparison with the expected qualities when normally distributed random numbers are used as an input to the array analysis. The complex trace of normally distributed numbers has uniformly distributed phases Φ l in the interval [− , + ]. For the definition of the spatial and temporal coherency in equations (17) and (19), the cosine and the sine of the phases are averaged to form new statistical variables where l is the number of a receiver in the array for the spatial coherency or a sample number for the temporal coherency. Similarly, L is the number of receivers in the array for the spatial coherency or the total number of samples used to calculate the temporal coherency. For many realizations of X c and X s with L = 1, the distribution function for X = X c and X = X s is the same in the open interval X ∈ (−1, 1) (Papoulis, 1985, p. 100). The standard deviation of this distribution function √ 2 for a sinusoid with unit amplitude A = 1. It is noted that the distribution function f (X) and its standard deviation X do not change when the sinusoid is randomly sampled, i.e., calculated from uniformly distributed random phases.
When L ⪰ 5 the influence of the central limit theorem increases such that the distribution function for X starts approaching the normal distribution function (Papoulis, 1985, p. 194-195) with zero mean μ = 0 and a standard deviation that reduces the standard deviation X = 1∕ √ 2 by the standard deviation of the mean 1∕ √ 2L for many realizations of X. The approximation of the normal distribution function is better for increasing L, i.e., the number of receivers in the array or the number of selected samples used. For example, for an array with N = 10 receivers, the standard deviation of X is N = 1∕ √ 2N = 0.22. This reasoning applies to X c and X s in equation (21) such that Γ = X c + iX s has the same standard deviation in the real and imaginary part. As a result, (Siddiqui, 1962(Siddiqui, , 1964, and references therein) with a scale parameter which is the standard deviation N of X, i.e., = N = 1∕ √ 2N. The most likely value, or mode, of the Rayleigh distribution for the coherency is given by mod f ( , ) = , the mean is μ f ( , ) = √ ∕2 and the standard deviation is f ( , ) = √ (4 − )∕2 (Papoulis, 1985, p. 112). For example, for an array with N = 10 receivers, the most likely value of the coherency is mod f ( , ) = 1∕ √ 20 = 0.22, the mean of the coherency is μ f ( , ) = 0.22 √ ∕2 = 0.28 and the standard deviation of the coherency is f ( , ) = 0.22 √ (4 − )∕2 = 0.15. Deviations from these characteristic parameters, or moments, to describe the Rayleigh distribution indicate that the input to the array analysis was not a set of uniformly distributed random numbers.
The array analysis described in section 2 is specifically designed to infer robust estimates for the wave number vector, electromagnetic source field, and the spatial and temporal quality of low-frequency radio waves with small electric field strengths observed with an array. Here 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 UK (Füllekrug et al., 2014;. 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 array analysis. 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 discharges from 5 to 15 kHz and intermittent radio transmissions, i.e., the UK radio clock from 59 to 61 kHz and four LORAN transmitters from 90 to 110 kHz (Füllekrug et al., 2015, Table S1). The intermittent radio transmissions are deliberately chosen to enable a direct comparison between the results with an array analysis for small and large electric field strengths. For example, in-between individual impulses from lightning discharges with spatial qualities q s > 1, irregular waveforms with qualities q s < 1 can be observed (Figure 1a). These irregular waveforms are up to 2 times larger than the sensitivity of the radio receiver ∼ 25μV/m, for example, during a sudden increase of irregular waveform ripples, which last for ∼3 ms (Figure 1a, inset). As a result, this irregular behavior is considered to be physically meaningful even though its cause is currently unknown. We suggest here to classify such irregular behavior of low-frequency radio waves as continuum radiation. The wording continuum serves to indicate that the irregular behavior is observed across the entire bandwidth of the radio receivers and can occur at any time of the recordings. The wording radiation is used to indicate that the irregular behavior is caused by physically meaningful radio waves.

Detection of Continuum Radiation
A more rigorous criterion for the detection of continuum radiation can be inferred from an array analysis of lightning discharges in comparison to an array analysis with normally distributed random numbers on input. For example, the distribution function of wavefront qualities measured in the frequency range from 5 to 15 kHz (Figure 1b, black line) differs significantly from the distribution function of wavefront qualities inferred from a simulation with normally distributed random numbers on input to the array analysis (Figure 1b, gray line). To be more precise, ∼ 60% of the observed wavefront qualities are considered to be significant because they clearly exceed the largest simulated spatial quality q s ≈ 1 inferred from random numbers. From a comparison with the recordings in Figure 1a, it is concluded that waveforms with spatial qualities q s > 1 are associated with impulses from lightning discharges. However, ∼ 40% of the observed wavefronts with qualities q s < 1 could be associated with distant lightning discharges, yet unknown sources, or they could be even random in nature. Whatever the physical cause of the continuum radiation is, for this array with 10 radio receivers, a significance threshold for the spatial quality of q s = 1 appears to be useful to identify continuum radiation and to discriminate continuum radiation against signals of known origin. In other words, continuum radiation can be detected by use of the criterion q s < 1 for the spatial quality of the observed radio waves.
It is interesting to note that the simulated distribution of random qualities peaks at q s = 0.12 (Figure 1b, gray line), while the expected spatial quality for normally distributed random numbers is q s = − log 10 (1 − 1∕ √ 2N) = 0.11, where N = 10 is the number of radio receivers in the array (equation (18) and section 2.4). The small difference between the expected and measured spatial quality can be explained as follows: When normally distributed random numbers are used for the real and imaginary part of a random complex measurement without any array analysis, the spatial coherency s follows exactly the expected Rayleigh distribution with the mode mod s = 0.22 (Figure 1b, inset, gray line). When the normally distributed random numbers are used as an input for the array analysis, the resulting electromagnetic source field is also Rayleigh distributed, but with a mode mod s = 0.26 (Figure 1b, inset, black line). This means that the array processing described in section 2 enhances the spatial coherency by ∼18% as a result of the weighted inversion (section 2.2). The net effect of weighting individual samples therefore results in the observed small increase of the spatial quality.

Electromagnetic Source Field of Continuum Radiation
The detection of continuum radiation is based on a threshold for the spatial quality of the observed radio waves where the electromagnetic source field is not considered. Yet the electromagnetic source field incorporates all the information about the physical processes that cause the observed radio waves. The electromagnetic source field can be displayed in a constellation diagram, which is a bivariate distribution function for the Figure 1. Low-frequency continuum radiation. (a) The envelope of electric field strengths exhibits lightning pulses superimposed on continuum radiation. The wavefront qualities for lightning pulses range from 2 to 3 (red) and 1 to 2 (blue), and for continuum radiation from 0 to 1 (black). The inset shows the ripples of continuum radiation for comparison with the instrumental noise floor (dashed line). (b) The distribution function of wavefront qualities (black) is composed of ∼40% continuum radiation (light blue) and ∼60% lightning pulses (dark shades of blue). The spatial qualities inferred from normally distributed random numbers are always smaller than one (gray) and demonstrate the significance of the lightning pulses. The inset shows the spatial coherencies inferred from normally distributed random numbers after the array processing (black) and without array processing (gray) for comparison with their corresponding Rayleigh distribution functions (red lines). The small difference between the modes of the Rayleigh distributions indicates that the array processing has only a minor effect on the inferred coherencies. See text for more details. normalized real and imaginary part of the electromagnetic source field (Füllekrug et al., 2014, Figure 2, right). Here we make use of the pulsed radio transmissions from the UK radio clock at 60 (±1) kHz, LORAN transmitters for marine navigation at 100 (±10) kHz Füllekrug et al., 2009, and references therein), and lightning discharges at 10 (±5) kHz toward a more detailed assessment of continuum radiation. It is found that continuum radiation occurs at the center of the constellation diagram for all frequencies and it can therefore be considered to be frequency independent (Figures 2a, 2c, and 2e). For example, the constellation diagram at 60 kHz clearly exhibits a bimodal pattern, either the observed electric field strengths are very large, or very small (Figure 2a). This observation reflects that the radio clock either transmits a carrier wave or it is switched off. The transition period between these two states is relatively short and appears as a line in the constellation diagram, which connects the continuum radiation at the center of the constellation diagram with large electric field strengths. It is important to note that these large electric field strengths always exhibit the same phase because the radio clock transmissions are controlled by a precise atomic clock. The constellation diagram at 100 kHz is more complex. Each of the four LORAN transmitters exhibits pulses with positive and negative polarity, which follow ellipsoidal shapes (Figure 2c). One specific transmitter exhibits six straight lines in the constellation diagram, which reflect time shifts of the pulses by ±1μs to convey information with a constant phase throughout each pulse. The varying amplitudes in the constellation diagram are indicative of the distances between the LORAN transmitters and the array. The continuum radiation is found in the center of the constellation diagram associated with small electric field strengths. The same applies to the constellation diagram for lightning pulses where the phases vary randomly such that the constellation diagram exhibits a rotational, or radial, symmetry with a maximum in the center (Figure 2e).
Another way to assess continuum radiation is by using bivariate distribution functions inferred from the absolute electric field strengths and spatial qualities of the observed radio waves (Figures 2b, 2d, and 2f ). These q s ∕E diagrams enhance the visualization of continuum radiation and include individual distribution functions for the spatial qualities and electric field strengths of the electromagnetic source field. For example, the q s ∕E diagram for the radio clock is clearly bimodal where the continuum radiation is characterized by spatial qualities q s < 1 associated with electric field strengths <10 μV/m. The q s ∕E diagram for the four LORAN transmitters clearly exhibits varying peak amplitudes for each of the four different LORAN pulses (Figure 2d). The rising and falling edge of the LORAN pulses blend into the continuum radiation, similar to the pulses from lightning discharges (Figure 2f ). Overall, it is concluded from this comparison of the radio clock, LORAN and lightning discharges that the existence of continuum radiation is independent of frequency. The same results have been obtained during an entire minute of recordings with the array, which confirms this conclusion.

Mapping Continuum Radiation and Signals in the Radio Sky
The main aim of array analyses is to map the apparent locations of arriving electromagnetic waves in the radio sky above the array. It has been shown that this is possible for lightning discharges and radio transmitters (Füllekrug et al., 2015, but it is currently unknown whether continuum radiation can also be mapped in the radio sky. As one example, the detected continuum radiation at 60 kHz with spatial qualities q s < 1 (Figures 2a and 2b) is extracted from the recordings to investigate its appearance in the radio sky. The wave number vectors of this continuum radiation are summarised in a bivariate distribution function in the k 1 ∕k 2 plane. This map exhibits a scatter that is reduced by use of the array aperture scaled with the covariance matrix of the wave number vectors to sharpen the image (Füllekrug et al., 2016, equations (6) and (7)). The resulting map of the radio sky shows that the continuum radiation is found at local zenith above the array (Figure 3a). This particular location indicates that on average the observed phase differences, or time delays, between pairs of radio receivers are zero as a result of using transfer functions and weights for the formulation of the Figure 3. Maps of low-frequency continuum radiation and signals. (a) Scatter of wave numbers from continuum radiation with wavefront qualities q s < 1 at 60 kHz after correction with the array aperture scaled by the covariance matrix of the wave numbers. The continuum radiation is found at local zenith above the array. The horizon is shown for comparison at k∕k 0 = 1.2 (solid circle) and k∕k 0 = 1.0 (dashed circle), where k 0 = ∕c is the reference wave number and c is the speed of light. (b) Scatter of wave numbers for the radio clock at 60 kHz and LORAN transmitters at 100 kHz after removing continuum radiation with wavefront qualities q s < 1. The four LORAN transmitters are found close to the horizon (solid and dashed circles) while the radio clock signal is found below the LORAN transmitter in the north direction (white circle). It is concluded that a spatial quality of q s = 1 is a useful threshold for discrimination, either to detect or to exclude continuum radiation. inversion in equation (6). In other words, the weighted inversion focusses continuum radiation onto the local zenith above the array instead of producing randomly scattered wave number vectors across the entire k 1 ∕k 2 plane. This observation is confirmed by using normally distributed random numbers on input for the array analysis which show the same pattern, i.e., the focussing of source locations at local zenith above the array. This pattern can therefore be used fortuitously to identify continuum radiation that is random in nature.
The spatial quality can also be used to remove continuum radiation from the recordings to enhance maps of the radio sky that are supposed to show signals of known origin. For example, the detected continuum radiation at 60 and 100 kHz with spatial qualities q s < 1 (Figures 2a, 2b, 2c, and 2d) is removed from the recordings and the remaining wave numbers vectors are summarized in a bivariate distribution function in the k 1 ∕k 2 plane (Figure 3b). The resulting map shows the locations of four LORAN transmitters, one toward the north, one toward the east, and two toward the south. The radio clock appears as a well confined location below the LORAN transmitter toward the north. In addition, the radio sky is covered by seemingly random source locations that are partly attributed to the rising and falling edges of individual LORAN pulses and the interference between the pulses from independently operating LORAN transmitters. Overall, it is concluded that a threshold of q s = 1 for the spatial quality of radio waves is a suitable criterion to either detect continuum radiation for further analysis or to exclude continuum radiation from further analysis.

Detection of Coherent Continuum Radiation
The detection of continuum radiation with spatial qualities q s < 1 and small electric field strengths ∼10 μV/m near the instrumental noise floor seems to imply at first glance that continuum radiation is of little, if any, scientific interest. On the other hand, a weak enhancement of the number of mapped radio sources is observed above the northern LORAN transmitter which could be physically meaningful (Figure 3b). A forensic visual inspection of the amplitude envelope of the electromagnetic source field reveals the presence of eight consecutive pulses from a LORAN transmitter and a trailing ninth pulse, which indicates that this pulse sequence was emitted by a master transmitter. This pulse sequence is repeated with a period of 90.07 ms and exhibits small electric field strengths ∼50 μV/m (Figure 4a). These small electric field strengths indicate that this LORAN transmitter is not associated with the other four detected LORAN transmitters that exhibit electric field strengths >500 μV/m (Figure 2d). The detected signal can be averaged over many repetition periods to increase the signal to noise ratio. Whilst this method works well for the detected transmitter (Figure 4a), there are no signs of other transmitters that would normally contribute to a LORAN pulse sequence. Several statistical procedures have been implemented with the aim to reveal the missing transmitters from the repetition intervals, i.e., calculating the average, median, and the logarithmic average from K = 654 repetitions of the pulse sequence. Yet neither of these statistical procedures can reveal the missing LORAN transmitters. Finally, the temporal coherency is calculated from the 654 repetition periods, which now clearly reveals two The envelope of electric field strengths at 100 kHz exhibits nine consecutive pulses from a LORAN transmitter with spatial qualities from 2 to 3 (red) and from 1 to 2 (blue), reaching down into continuum radiation from 0 to 1 (black). (b) The temporal qualities inferred from a repetition interval of 90.07 ms exhibit three LORAN transmitters (labels 1, 2, and 3). The theoretical significance limit of the qualities for normal distributed random numbers (dashed line) and five time its value (dotted line) is shown for reference. The inset shows two double pulses that are attributed to ground and sky waves. (c) The distribution functions of the electric field strengths taken at the pulse maxima of the two strongest LORAN transmitters (labels 1 and 3) clearly show negative and positive polarity at ∼ ±28 μV/m (black) and at ±2-4 μV/m (red). (d) The distribution function of the electric field strengths taken from the smallest LORAN pulse maxima (black and label 2) exhibits two relative maxima at ±0.5-0.7 μV/m (dashed lines) while the distribution function of electric field strengths in-between consecutive pulse maxima is single peaked (red). The lower panel shows the difference between the two distribution functions with a relative maximum at ∼0.0 μV/m that indicates a meaningful splitting of the distribution function as part of continuum radiation well below the noise floor of the radio receiver. more LORAN transmitters (Figure 4b). These two LORAN transmitters appear to emit eight or seven consecutive pulses while the first LORAN transmitter emits nine pulses. For comparison, the temporal quality arising from normally distributed random numbers is q t = − log 10 (1 − 1∕ √ 2K) = 0.012, which acts as a detection limit for the temporal qualities. All the detected LORAN transmitters exhibit temporal qualities that are at least 5 times larger than the detection limit such that there is no doubt about their significance. In addition, the pulses from the first LORAN transmitter exhibit a secondary maximum that is attributed to sky waves arriving after the pulse from the ground wave (Figure 4b, inset).

Quantifying Small Electric Field Strengths
After the detection of small LORAN pulses by use of the temporal quality, it becomes now possible to quantify the associated small electric field strengths within the continuum radiation. The maxima of the temporal qualities for all pulses in the repetition period are used as a reference to extract ±30 μs, i.e., ±30 samples, of the electric field strengths around individual pulse maxima. The distribution function of the extracted electric field strengths is subsequently calculated with a resolution of 10 nV/m for each transmitter. These distribution functions exhibit two separate maxima, one for negative pulses and one for positive pulses (Figure 4c). This polarity of the LORAN pulses is embedded in the originally recorded electric field strengths, i.e., the real part of the complex trace, while the envelope of the LORAN pulses reflects the absolute value of the complex trace (Figures 4a and 4b). The LORAN transmitter with the largest temporal qualities exhibits the largest electric field strengths ∼ ±28 μV/m. The LORAN transmitter with the second largest temporal qualities exhibits a distribution of electric field strengths with relative maxima at ∼ ±3 μV/m. This double peaked distribution results from the overlap of the distribution for pulses with negative and positive polarity. It is remarkable that these electric field strengths can be detected and quantified even though the electric field strengths are smaller than the instrumental noise floor ∼ 25μV/m. For the LORAN transmitter with the smallest temporal qualities, the two relative maxima for the distribution of electric field strengths are found at ∼ ±0.6 μV/m and are difficult to distinguish (Figure 4d, upper panel). As a result, the distribution function of electric field strengths measured in-between two consecutive pulse maxima is used as a reference. This distribution is single peaked with a maximum at ∼0.0 μV/m as it would be expected for continuum radiation that is random in nature. The difference between the distribution functions of electric field strengths at the pulse maxima and in-between pulse maxima therefore exhibits a clear maximum at ∼0.0 μV/m (Figure 4d, lower panel). This result indicates a meaningful splitting of the distribution function of the electric field strengths as part of coherent continuum radiation almost 2 orders of magnitude below the noise floor of the radio receiver.

Discussion
A novel array analysis is developed to detect continuum radiation with small electric field strengths and spatial qualities of the radio waves q s < 1 across the array. It is shown that continuum radiation is composed of impulses from lightning discharges (Figures 1a and 2f ), intermittent radio transmissions (Figures 2d and 4), and that it can be random in nature (Figures 2b and 3a). In particular, the temporal coherency of radio waves is used to detect LORAN transmissions that otherwise would have remained hidden below the instrumental noise floor of the recordings. The detection of this coherent continuum radiation is assisted by the analysis of distribution functions (Figures 1b, 2b, 2d, and 2f ). Amplitude distribution functions and interarrival times of lightning pulses have been characterized before in the framework of the Stanford ELF/VLF Radio Noise Survey (Chrissan & Fraser-Smith, 2000, 2003. The aim of that work was mainly to characterize the parameters of the observed distribution functions for lightning, where the largest analyzed frequency was 32 kHz (Chrissan & Fraser-Smith, 1996). The work presented here is focussed on the discrimination between continuum radiation, lightning, and transmitters as measured with an array of radio receivers. It is expected that lightning could also contribute to the continuum radiation observed at 60 kHz, but this contribution appears to be relatively small (Figure 3a). The detection of two LORAN transmitters with small electric field strengths is entirely based on a stable phase of the continuum radiation over time where the actual amplitudes are not relevant (Van Vleck & Middleton, 1966). In fact, the closest four LORAN transmitters have such large amplitudes and frequent pulse repetitions such that the average, median and logarithmic average of the pulse repetitions do not converge. As a result, a stable phase, or large temporal coherency, is considered to be of critical importance for the detection of physically meaningful signals within continuum radiation. This phase stability depends on three main factors, the phase stability of the source, the stability of the Earth-ionosphere cavity in which the radio waves propagate, and the phase stability of the recordings determined by the internal oscillator disciplined by GPS clocks . The LORAN signals with small electric field strengths are detected below the instrumental noise floor, but these small electric field strengths are not recorded directly. In fact, the small electric field strengths are quantified from their most likely values, or modes, of their distribution functions. As a result, more sophisticated statistical methods to determine the phase stability of the recorded electric field strengths and an improved characterisation of their distribution functions are clearly warranted to achieve the long-term objectives stated in section 1.2.

Summary and Conclusions
A novel array analysis for low-frequency radio waves with the capability to treat small electric field strengths has been developed. The rigorous array analysis is applied for the first time to lightning pulses and intermittent radio transmissions. During the absence of lighting pulses and radio signals, an underlying continuum radiation has been detected that is present in all frequency ranges. The continuum radiation is characterized by the spatial quality, or flatness, of radio waves across the array. The coherent continuum radiation is used to identify intermittent radio transmissions with electric field strengths below the noise floor of the radio receivers. Another part of the detected continuum radiation is found at local zenith above the array when it is random in nature. The presented work is thought to be only the starting point toward a determination of many more sources and properties of continuum radiation. For example, the phase stability of coherent continuum radiation could be used to assess the state of the lower ionosphere, which is controlled by space