Multipath propagation of low‐frequency radio waves inferred from high‐resolution array analysis

The low‐frequency radio sky shows the locations of electromagnetic radio sources with a characteristic dilution of precision. Here we report a thorough high‐resolution analysis of radio waves from low‐frequency (∼20–150 kHz) radio communication transmitters which are recorded with a small aperture array of radio receivers during the day. It is found that the observed dilution of precision results from the array geometry of the radio receivers, a birefringent wave propagation, and the correlated multipath propagation of low‐frequency radio waves. The influence of the array geometry on the dilution of precision is reduced by taking into account the impulse response of the array. This procedure reveals for the very first time the splitting of one single radio source into two distinct source locations separated by ∼0.2°–1.9° which result from a birefringent wave propagation. The two locations are yet more clearly identified by using the polarity of the modulated wave number vectors of the radio waves. This polarity is also used to quantify the dilution of precision arising from correlated multipath propagation which is discriminated against wave number fluctuations arising from the timing accuracy of the radio receivers. It is found that ∼69% of the wave number variability is of natural origin and ∼31% originates from the timing accuracy of the receivers. The wave number variability from correlated multipath propagation results in a standard deviation ∼2–8% relative to the source location. This compact measurement of correlated multipath propagation is used to characterize the uncertainty of source locations in the radio sky. The identification of correlated multipath propagation strongly suggests the existence of very fast processes acting on time scales <1 ms in the D region ionosphere with physically meaningful effects on low‐frequency radio wave propagation. This important result has implications for practical applications in that the observed multipath propagation enables the determination of natural limits for the accuracy of navigation and lightning location methods using low‐frequency radio waves.


Introduction
The low-frequency radio sky from ∼20 to 250 kHz exhibits anthropogenic electromagnetic noise sources with a finite angular resolution of unknown physical origin. It was speculated that this angular resolution is determined by a birefringent wave propagation and random fluctuations of the wave number vectors in the sky [Füllekrug et al., 2015]. In other geoscience applications, such as lightning detection networks [e.g., Lyu et al., 2014;Wu et al., 2014;Stock et al., 2014;Rakov, 2013;Bitzer et al., 2013;Mezentsev and Füllekrug, 2013, and references therein], it is common practice to estimate the lightning location accuracy by a forward simulation of the network performance taking into account known instrumental measurement uncertainties, [e.g., Stock et al., 2014;Bitzer et al., 2013]. This general strategy can be applied to small aperture arrays which determine the arrival angle and length of the wave normal vector to locate radio sources in the sky. For example, a small aperture array was operated on Charmy Down airfield near Bath in the UK [Füllekrug et al., 2014;Mezentsev and Füllekrug, 2013]. The observed wave number vectors exhibit a significant scatter in the wave number domain such that only the median wave number vectors were used to determine the locations of radio sources in the sky. This scatter manifests itself in bivariate distribution functions with a characteristic ellipsoidal shape, which is independent of the radio source [Füllekrug et al., 2014, Figure 2, left]. These observations strongly indicate that the specific array geometry of the radio receivers is responsible for the general shape of the wave number scatter, while the scatter itself is of physical origin. This effect is known as the dilution of precision in navigational service provision, e.g., Global Navigation Satellite Systems (GNSS), where the dilution of precision is attributed to the momentary geometric configuration of the satellite constellation [Langley, 1999].

10.1002/2015RS005781
As a result, this contribution describes the influence of the array geometry on the observed bivariate distributions of the scattered wave numbers and aims to quantify the physical processes leading to the scatter toward an improvement of the location accuracy of small aperture arrays.

Wave Number Vector Determination
This contribution uses observations with a small aperture array which consists of 10 wideband radio receivers distributed over an area of ∼1 km × 1 km near Bath in the UK. Exemplary measurements from 15:00:00 to 15:00:02 on 13 May 2011 are extracted from the original recordings which are taken with a sampling frequency of 1 MHz for half an hour. The seemingly short duration of ∼2 s is chosen to enable a determination of the wave number vectors without the need for high-performance computing facilities. The wave number vectors are inferred by making the ansatz that the electromagnetic signal from a single source sweeps as a plane wave E(t) exp(−i(kr n − t)) across the network, where E(t) is the complex amplitude of the electromagnetic source field which varies with time t, k is the wave number vector of the electromagnetic wave, r n are the geographic coordinates of the radio receivers n = 1, 2, … N = 10, and is the frequency of the radio wave. The recordings at all individual receivers are multiplied with the reference plane wave exp(i(k 0 r n − t)), where k 0 is the reference wave number vector which assumes a propagation along the great circle path from the transmitter to the receiver at the speed of light. The multiplication of the recorded time series with the reference plane wave results in a complex time series with a phase which is referenced to the start time of the recordings and k 0 . Subsequently, a low-pass filter is applied to the time series to avoid interference with the signals from neighboring transmitters in the spectrum which typically exhibits a signal-to-noise ratio ∼1-2 orders of magnitude measured with an accuracy ∼20 μV/m [Füllekrug, 2010]. The resulting complex time series y can then be explained by the referenced plane wave from a single source where k = k − k 0 . Carrier wave transmitters produce an electromagnetic wave field where all propagation paths, i.e., ground and multihop sky waves, are recorded simultaneously at the receiver with one apparent source location in the sky. These multiple signals at the same frequency might be discriminated with an analysis of the covariance matrix of the complex time series in equation (1), e.g., by use of the MUltiple SIgnal Characterization (MUSIC) algorithm [Stoica and Nehorai, 1989;Bienvenu and Kopp, 1980;Schmidt, 1979]. However, this contribution focuses on an analysis of the entire electromagnetic wave field to study the more delicate multipath propagation which results in a subtle spatial spreading of the apparent wave number vector arising from small variabilities of the propagation path. The wave number vectors are inferred from the phase of equation (1) where x n,1 and x n,2 are the n ′ th receiver location coordinates in the east and north directions, Φ is the phase of the electromagnetic source field, and n is the observed phase of the referenced complex time series. Equation (2) is solved with the Gaussian method of least squares to minimize the phase residuals n such that the resulting wave number vectors are finally calculated from k = k + k 0 . The phase residuals n ≈ 20 ns are in excellent agreement with the instrumental accuracy [Füllekrug et al., 2014;Füllekrug, 2010] which is a remarkable result, given that the array only samples ∼1 km of an electromagnetic wave with a wavelength of ∼15 km (∼7%) at frequencies ∼20 kHz. The inferred wave number vectors exhibit a comparatively large variability k around their median value such that it was speculated that these wave number fluctuations contain a physically meaningful signal [Füllekrug et al., 2015] which is the subject of this contribution.

Theory and Analysis
The observed scatter of the wave number vectors is assessed by use of a perturbation theory with small deviations from the solutions of the array analysis. It is sensible to base this assessment on the array aperture because it is expected that the array geometry defines the ellipsoidal shape of the bivariate scatter of the wave numbers. The array aperture T(k, r n ) is defined by the wave number vector k and the n = 1, 2, … N radio receiver locations r n

10.1002/2015RS005781
where N = 10 is the total number of radio receivers and each summand on the right-hand side describes a linear wave model, i.e., the basis function of a sparse spatial Fourier transform. The array aperture describes the response of the array to an electromagnetic wave with the reference wave number vector k 0 = (k 0,1 , k 0,2 ) T with k 0 = √ k 2 0,1 + k 2 0,2 = 0 ∕c, where k 0,1 and k 0,2 are the reference wave number components in the geographic east/west and north/south directions, 0 is the center frequency of a radio source, and c is the speed of light. More generally, the array aperture describes how well the array focuses a collimated electromagnetic wave onto the wave number domain. The array aperture is maximum when the exponent of each linear wave in equation (3) vanishes In this idealized case, the array aperture is exactly T(k, r n ) = 1 for k − k 0 = 0. Yet, the observations suggest that real measurements result in small deviations from the ideal case such that residual wave number vectors k occur k − k 0 = k.
Three possible cases need to be distinguished.
1. The residual wave number vectors can result from coherent recordings made by all radio receivers in the array. In this case, the observed wave number vectors have a well-defined physical origin as a result of correlated multipath propagation, where the wave number vectors spread over a certain spatial area. 2. The residual wave number vectors can result from incoherent recordings made by all radio receivers in the array. In this case, the observed wave number vectors result from random fluctuations of the residual phase measurements arising from uncorrelated multipath propagation where the incoherent scatter of the radio wave produces phase residuals which are larger than the residuals measured under laboratory conditions. This uncorrelated multipath propagation is not considered here because the observed phase residuals are the same as under laboratory conditions. 3. The residual wave number vectors can result from random fluctuations of the residual phase measurements arising from the timing accuracy at each individual radio receiver. In this case, the observed phase residuals are as small as those measured under laboratory conditions. The phase residuals are specifically thought to originate from the random clock jitter n = 0 t n where t n ≈ 20 ns [Füllekrug, 2010].
In all the three cases, the residual wave number components in the local east/west direction k 1 and in the north/south direction k 2 are related to the phase residuals n and the relative geographic coordinates x n,1 and x n,2 for each individual receiver by combining equations (4) and (5) such that Note that equation (6) treats the phase residuals n as observations to determine the unknown wave number fluctuations k which is an extension of, and in full agreement with, equation (2), because the phase of the electromagnetic source field Φ is the same at all receivers. Equation (6) is formulated for each receiver and summarized in matrix form, where G is composed of the x n,1 and x n,2 coordinates in each individual row n where k are the wave number vector residuals with zero mean. This linear system of equations is solved by the generalized inverse G −g = (G T G) −1 G T to infer the residual wave number vector k = ( k 1 , k 2 ) T from the observed phase residuals . The wave number residuals are characterized by their covariance matrix where the brackets describe an average, or expectation, value of many wave number residuals with zero mean. This covariance matrix is entirely determined by experimental measurements. For example, a radio transmitter for submarine communication in Rhauderfehn, Germany, operates at a center frequency of

Radio Science
10.1002/2015RS005781 23.4 kHz and the observed scatter of the wave number residuals can be summarized in a bivariate distribution function (Figure 1a). This bivariate distribution exhibits an ellipsoidal shape and the scatter can be characterized by the covariance matrix in equation (8).
The aim of the forward simulation is to explain the observed covariance matrix of the wave number residuals C k from fundamental principles. Toward this end, the covariance matrix can be estimated, or simulated, by use of suitable phase residuals. In this case, the estimated covariance matrixĈ k results from a multiplication of the right-hand side of equation (7) with its transpose and taking the expectation valuê where the simulated covariance matrix of the wave number residuals is calculated from the covariance matrix of the observed phase residuals C =< T >. The covariance matrix of the phase residuals can be expanded into the identity matrix I multiplied by the variance of the phase residuals 2 assuming the phase residuals are uncorrelated and have the same variance at each radio receiver in the array such that This important result means that the covariance matrix of the wave number residuals can be estimated from the variance of the phase residuals and the matrix (G T G) −1 which is characterized by the array geometry. In other words, the observed scatter of the wave numbers is entirely determined by the array geometry and the statistical properties of the phase residuals. The covariance matrix of the phase residuals is simulated here by using the observed phase residuals where the spectral amplitudes are preserved and the phases of the spectrum are scrambled by a random permutation. The application of this method ensures that the statistical properties of the phase residuals are as similar as possible to the observed phase residuals and satisfy the assumptions leading to equation (10), i.e., that the phase residuals are uncorrelated and have the same variance at each radio receiver. The random permutation is necessary because using the observed phase residuals would result in the observed wave number residuals and hence not provide any new information.
The estimated bivariate distribution functionf is calculated from the simulated covariance matrix where the |Ĉ k | is the determinant of the covariance matrix estimated from the phase residuals in equation (10). The resulting bivariate distribution function of this forward simulation with the estimated covariance matrix of the wave number residuals (Figure 1b) is much smaller than the observed distribution function ( Figure 1a). It is interesting to note that the simulated distribution function is similar to the array aperture in equation (3) when a suitable scaling factor is applied [cf. Füllekrug et al., 2015, Figure 1]. As a result, the simulated distribution function is understood as the response of the array to a radio wave arriving from one specific location with uncertainties determined by the clock jitter t n of the radio receivers. The comparison between the observed and simulated wave number distribution functions (Figures 1a and 1b) therefore suggests that ∼31% of the standard deviation from the observed wave number scatter can be explained by measurement uncertainties and ∼69% of the scatter has a physical origin.
The large scatter of the observed wave number distribution (Figure 1a) can be reduced significantly by use of the simulated array response (Figure 1b). This sharpening of the image is achieved by filtering the observed image with the inverse impulse response of the array which needs to be determined. The simulated array responsef is related to the Dirac function d by the two-dimensional convolution f = h * d [Gonzalez and Woods, 2008], where h is the impulse response, or transfer function, in the wave number domain. This convolution is conveniently linearized by making use of a two-dimensional discrete Fourier transformF where H∘D is the element-wise, or Hadamard, product of the Fourier transformed impulse response H and Dirac function D and .∕ is the element-wise matrix division. The calculated impulse response can then be applied to the Fourier transform of the observed bivariate distribution function F D = F.∕H (14) which results in an estimated location of the radio source in the wave number domain after an inverse two-dimensional discrete Fourier transform ofD. In ideal circumstances, the result would be one single point d in the wave number domain. Yet the application of the impulse response to the observations reveals two different locations for the same radio source (Figure 1c). This surprising result strongly suggests that the observed scatter of the wave numbers ( Figure 1a) partly results from a birefringent radio wave propagation.
In searching for a physical explanation of the observed splitting of the wave number scattering, the apparent frequency f a was investigated in more detail [Taner et al., 1979], where E r and E i are the real and imaginary part of the electromagnetic source field. The apparent frequency exhibits very small deviations ∼0.2% = ±50 Hz/23400 Hz for the radio transmitter in Rhauderfehn (Figure 1d). These small frequency variations result from the modulation scheme of the radio transmissions known as Minimum Shift Keying (MSK) [Thomson, 2010, and references therein]. Wave number vectors with a positive apparent frequency are associated with a clockwise phase progression of the electromagnetic source field in the constellation diagram and negative apparent frequencies are associated with a counterclockwise phase progression. Another way of expressing this phase progression is to state that the modulating wave number vector exhibits an alternating polarity. As a result of this modulation, it is necessary to segment the observed wave numbers with respect to the polarity of the apparent frequency ( Figure 1d, dashed lines). This segmentation results in two different wave number distributions after correction for the impulse response ( Figure 1e). The two wave number distributions with different apparent frequency polarity can clearly be distinguished. The wave numbers associated with positive apparent frequencies are ∼2% smaller than the wave numbers associated with negative apparent frequencies. The segmentation is also applied to the electric field strength of the electromagnetic source field which exhibits a distinct splitting into two different distributions (Figure 1f ). Electric field strengths with positive apparent frequencies exhibit ∼7% smaller amplitudes than electric field strengths with negative apparent frequencies. In summary, radio waves with positive apparent frequencies during a clockwise phase progression in the constellation diagram result in smaller wave numbers and smaller electric field strengths than radio waves with negative apparent frequencies.
The segmentation with respect to the polarity of the apparent frequency results in a reduction of the original variance of the wave number scatter by ∼10-20%. Yet given that the reduced wave number scatter is still much larger than the difference between the wave numbers with positive and negative apparent frequencies, it is concluded that the remaining wave number fluctuations are correlated across the entire array. This correlated multipath propagation can be quantified by use of an eigenvalue analysis of the covariance matrices for the observed and estimated wave numbers where R andR are the rotation matrices used to determine the variances along the principal axes of the positive semidefinite covariance matrices. The principal axes of these two matrices are almost identical as a The transmitters are characterized by their call sign, center frequency f , geographic location latitude (lat), longitude (lng), and distance d to the array. The relative wave numbers for negative apparent frequencies k − ∕k 0 are larger when compared to the wave numbers for positive apparent frequencies k + ∕k 0 . The relative standard deviations of the correlated multipath propagation for negative apparent frequencies − m ∕k − and positive apparent frequencies + m ∕k + range from ∼2 to 8%. confirmed by the calculated angles between the principal axes of the scatter ellipse and the geographic reference frame. The variances for the correlated multipath propagation can then be calculated from the differences between the observed and simulated variances of the wave numbers for the major and minor axis of the bivariate distributions when assuming that the simulated variances are statistically independent from the variances of the correlated multipath. In this way, the observations of the wave numbers are corrected by removing the contribution of measurement uncertainties to the observed variances.
It is possible to determine a compact measurement of the correlated multipath propagation by calculating the average standard deviation directly from the observed and simulated covariance matrices where tr is the trace of the covariance matrix. This description of correlated multipath propagation is well suited for practical applications when the correlated multipath propagation along the principal axes is of minor interest. For example, the multipath propagation of five low-frequency radio communication transmitters was investigated in detail (Table 1). It is found that all the normalized wave numbers are larger for negative apparent frequencies than for positive apparent frequencies. The normalized standard deviation of the wave number scatter resulting from the correlated multipath propagation m ∕k ranges between ∼2 and 8% for negative and positive apparent frequencies. This combined standard deviation of the correlated multipath propagation is hence a sensible measure of the uncertainty for source locations in the radio sky.

Interpretation and Discussion
The analysis of the wave number scattering from a radio wave transmitter distinguishes four contributions: (1) the array geometry of the radio receivers which determines the ellipsoidal shape of the scatter, (2) birefringent wave propagation which results from the polarity of the modulating wave number vector of the electromagnetic source field, (3) correlated multipath propagation of the radio waves across the entire array, and (4) random fluctuations which result from measurement uncertainties arising from the clock jitter at each individual receiver. Each of these contributions is identified with a distinct methodology, but the interpretation of the correlated multipath propagation is not necessarily unique. An alternative explanation could be, for example, that the GNSS timing signals vary coherently across the entire array which has a size of only ∼1 km × 1 km. In this case, similarly sized variabilities in the F layer ionosphere could influence the transionospheric propagation of the radio waves from the GNSS satellites and thereby result in the observed correlated multipath. While this explanation is plausible, it is considered to be unlikely, because the correlated multipath reported here is observed on the microsecond time scale which seems too fast for ionization variability in the F layer ionosphere because the scintillation index is typically given on the minute time scale. In addition, scintillation is less likely to occur in middle latitudes when compared to high latitudes.
The birefringent wave propagation results in two distinct source locations in the wave number domain. The source location associated with a positive apparent frequency occurs at a smaller wave number such that it would appear at a larger elevation angle in the radio sky when compared to the source location associated with a negative apparent frequency. This angular resolution ranges from ∼0.2 ∘ to 1.9 ∘ . The larger elevation angle for positive apparent frequencies indicates that the sky wave contribution to the radio wave might have penetrated deeper into the lower D layer ionosphere when compared to the source location associated with a negative apparent frequency. The smaller electric field strengths support this interpretation because the longer propagation path would result in a larger absorption of the radio wave. This observation is consistent with the fact that the reflection height of a radio wave penetrating into a plasma is controlled by the plasma frequency which increases with altitude as a result of an increasing electron density. A radio wave with a larger frequency would then penetrate deeper into the lower ionosphere associated with a larger absorption and occurs at a larger elevation angle. Another possibility to explain the birefringence might be that ordinary and extraordinary mode radio waves penetrate to varying depths into the D region ionosphere, where the ordinary mode is unaffected by the geomagnetic field and the properties of the extraordinary mode are influenced by the geomagnetic field [Budden, 1988, chapter 13]. Considering for simplicity the case of waves vertically incident on the horizontally stratified lower ionosphere, the ordinary wave is reflected at the height where the magnetoionic variable X = 1, which is where the wave frequency f is equal to the electron plasma frequency f N [Rishbeth and Garriott, 1969, p. 49], i.e., f N ≈ 20 kHz. The extraordinary wave is reflected at the height where MHz for southern England, reflection takes place at a height where f N ≈ 160 kHz which is appreciably higher than the ordinary wave reflection height.
It is interesting to note that the normalized standard deviation of the correlated multipath propagation supports a possible influence of the geomagnetic field (Table 1). The smallest multipath propagation (∼1.8-2.3%) is observed for the two transmitters toward the east of the array in Rhauderfehn (DHO) and Pinneberg (DDH47) while the largest multipath propagation (∼8.0-8.2%) is observed for the transmitter toward the west of the array in Cutler (NAA). The correlated multipath propagation observed for the transmitters toward the north in Skelton (GQD) and the south in Tavolara (ICV) is found between the two extremes (∼3.5-5.5%). This result might be explained with observations of anisotropic wave propagation as a result of the ground conductivity [e.g., Hutchins et al., 2013;Taylor, 1960, and references therein] and/or theoretical model calculations of low-frequency radio waves propagating in the east/west and west/east directions at different angles with the geomagnetic field [Galejs, 1972, chapter 6 and references therein]. It seems therefore possible to establish a long-wave scintillation index for amplitude and phase measurements of low-frequency radio waves to describe the natural variability of the D region ionosphere. However, it is unlikely that this scintillation index would be able to distinguish between kinetic reaction times and transport phenomena, because this scintillation index will require a certain integration time to reach a stationarity limit. Radar observations are more likely to provide this level of detail for limited spatial areas [e.g., Barabash et al., 2012;Osepian et al., 2009;Friedrich and Rapp, 2009, and references therein].
The multipath propagation observed here has an immediate application toward a physics-based assessment of the natural limits for the accuracy of navigation and lightning location methods at low frequencies. The angular resolution of radio sources in the sky spans ∼0.2 ∘ -1.9 ∘ which corresponds to a relative location accuracy ranging from tan(0.2 ∕180) ≈ 0.4% to tan(1.9 ∕180) ≈ 3.3%. For navigation, the bivariate normal distribution functions enable the calculation of an average over time such that the accuracy will decrease proportional to ∼ 1∕ √ N, where N is the number of microseconds used for calculating the average. For example, if the stationarity limit of the scintillation index is reached after ∼1 ms, the maximum navigation accuracy that can be achieved ranges from ∼0.4/ √ 1000% to ∼3.3/ √ 1000% such that a transmitter at a distance of 1000 km could be located with an accuracy ∼0.1-1.0 km. The current of lightning discharges varies on the time scales of microseconds such that the prospects of temporal averaging are very limited. On the other hand, the spectrum of lightning discharges is very broad such that an average over many frequency bands can be calculated. For example, if 33 spectral bands with statistically independent multipath propagation are used, the maximum lightning location accuracy that can be achieved ranges from ∼0.4/ √ 33% to ∼3.3/ √ 33% such that a lightning discharge at a distance of 1000 km could be located with an accuracy ∼0.7-5.7 km. As a result, navigation and lightning location methods will work best when mainly ground waves are used for the analysis to minimize ionospheric influences.

Summary
The source locations of low-frequency radio communication transmitters are investigated in detail in the wave number domain with a high-resolution analysis of measurements with a small aperture array of distributed radio receivers. It is found that the scatter of the wave numbers results in a physically meaningful dilution of precision of the location as a result of four main factors: (1) the array geometry of the radio receivers, (2) birefringent wave propagation arising from the polarity of the modulating wave number vector of the radio wave, (3) correlated multipath propagation of radio waves across the array, and (4) random fluctuations arising from the clock jitter at each individual receiver. In particular, it is found that radio sources with fractionally larger frequencies occur at larger elevation angles and exhibit more absorption as a result of the deeper penetration into the D region ionosphere when compared to lower frequency radio waves. The correlated multipath propagation occurs on the microsecond to millisecond time scale which strongly suggests that the observed wave number scatter originates from the variability of the lower D region ionosphere with implications for the limiting accuracy of navigation and lightning location methods using low-frequency radio waves.