Array analysis of electromagnetic radiation from radio transmitters for submarine communication

The array analyses used for seismic and infrasound research are adapted and applied here to the electromagnetic radiation from radio transmitters for submarine communication. It is found that the array analysis enables a determination of the slowness and the arrival azimuth of the wave number vectors associated with the electromagnetic radiation. The array analysis is applied to measurements of ∼20–24 kHz radio waves from transmitters for submarine communication with an array of 10 radio receivers distributed over an area of ∼1 km ×1 km. The observed slowness of the observed wave number vectors range from ∼2.7 ns/m to ∼4.1 ns/m, and the deviations between the expected arrival azimuths and the observed arrival azimuths range from ∼−9.7° to ∼14.5°. The experimental results suggest that it is possible to determine the locations of radio sources from transient luminous events above thunderclouds with an array of radio receivers toward detailed investigations of the electromagnetic radiation from sprites.


Introduction
Recent theoretical model calculations predicted that the exponential growth and splitting of sprite streamers near ∼40 km height cause electromagnetic radiation in the frequency range ∼100 kHz [Qin et al., 2012]. This prediction was promptly confirmed by experimental measurements of intense consecutive, or dancing, sprites [Füllekrug et al., 2013b]. Similarly, sprite streamers near ∼75 km height are the source of electromagnetic radiation in the frequency range <1 kHz [e.g., Qin et al., 2012;Cummer and Füllekrug, 2001;Pasko et al., 1998;Cummer et al., 1998, and references therein]. It is hence somewhat surprising that sprites do not radiate significantly in the frequency range near 10 kHz [Franz et al., 1990], which is possibly attributed to the relatively low streamer density of sprites near ∼50-60 km height. On the other hand, recent electromagnetic radiation measurements during sprite occurrences, which are located only ∼80 km away from a radio receiver, show a rapid succession of numerous broadband pulses with pausing time intervals ∼50-100 μs [Füllekrug et al., 2013a], similar to the pulse sequences observed during the initiation of lightning discharges [e.g., Kolmasova et al., 2014;Gurevich et al., 2013;Stolzenburg et al., 2013, and references therein]. In fact, many sprites are associated with such pulse sequences, which blend into electromagnetic noise at large distances from the sprite [Füllekrug et al., 2013b]. Yet the underlying wave number vectors required for a more detailed array analysis to locate sprites are not very well known.
Wave normal, or Poynting, vectors of electromagnetic waves are measured in space with three axis electric and magnetic field sensors [e.g., Santolik et al., 2014Santolik et al., , 2003. At the surface of the Earth, it is not possible to use three axis electric or magnetic field measurements because the boundary conditions for electromagnetic fields near a conductor result in a screening of the vertical magnetic and horizontal electric fields and in a superposition with the electromagnetic fields induced inside the conductive Earth. As a result, this contribution aims to determine the wave number vectors by use of an array of radio receivers. Specifically, the array analyses developed for seismic and infrasound arrays [e.g., Füllekrug et al., 2013c;Rost and Thomas, 2002;Gaffet et al., 1998;Cansi, 1995] are adapted and subsequently exemplarily applied to ∼20-24 kHz radio waves from radio transmitters for submarine communication to validate the array analysis.

Ansatz
An array of 10 radio receivers was deployed over an area of ∼1 km × 1 km on an old airfield at Charmy Down near Bath in Southwest England [Mezentsev and Füllekrug, 2013]. The array is characterized by the sweet FÜLLEKRUG ET AL.
©2014. The Authors. Figure 1. The array of radio receivers. (left) Ten radio receivers are distributed over an area of ∼1 km × 1 km. The sweet spot of the array (red dot) minimizes the average distance <d> to all receivers. The radio waves of four transmitters sweep across the array from different directions (lines). (right) The array response (|T|) is most suitable for radio waves propagating along the north-west to south-east and the south-east to north-west direction as a result of the geometric distribution of the receivers. The Nyquist wave number (dotted circle) and the fundamental wave number (dashed circle) characterize the array. The wave numbers within the quasi-static part of the wave number spectrum are of interest (solid circle). The radio waves from the four transmitters Tavolara (ICV), Rhauderfehn (DHO), Skelton (GQD), and Cutler (NAA) arrive from different azimuths (lines). spot (51.4277 • N, 2.3429 • W), i.e., the location with the smallest average distance ∼419 m to all the radio receivers ( Figure 1, left). The aperture of the array is defined by the geometry of the radio receiver locations that are constrained by logistic access to Charmy Down airfield. The measurements at the receivers y c (r, t) are explained by an electromagnetic source field E(t), which propagates as a plane wave e −i(kr− t) across the array of radio receivers (1) In this description, y c is the complex trace that is composed of a real part and imaginary part, i.e., the originally recorded time series and the time series shifted by 90 • [Bracewell, 1990, and references therein]. This ansatz enables a measurement of the envelope |y c | and the instantaneous phase c when formulated in polar form y c = |y c |e −i c [Taner et al., 1979]. The complex trace depends on time t and on the receiver locations r n = (x 1,n , x 2,n ) T in the east-west (x 1 ) and north-south (x 2 ) directions with respect to the sweet spot of the network, where the receiver locations are numbered from n = 1, 2, … N = 10. The electromagnetic source field E(t) is composed of an amplitude A(t) and phase (t) such that E(t) = A(t)e −i (t) is independent of any receiver location. The wave number vector k = (k 1 , k 2 ) T lies in the horizontal x 1 × x 2 plane because the array has no significant vertical gradient. The product t of the radio wave frequency and time t describes a continuous phase progression that is taken out of equation (1) by multiplication with e −i t such that This step ensures that the phase of y n (t) is referenced to the phase 0 at the start time of the recordings after the application of a low pass filter with a cutoff frequency of 150 Hz.

Array Aperture
The time independent array aperture or transfer function, describes the impulse response of the array with respect to the electromagnetic source field (Figure 1, right). This array transfer function is characterized by the Nyquist wave number k N and the fundamental wave number k F . The Nyquist wave number k N = 2 ∕(2d min ) is determined by the smallest distance d min ≈ 220.2 m between two receivers in the array and the fundamental wave number k F = 2 ∕(2d max ) is determined by the largest distance d max ≈ 1261.4 m between two receivers in the array. It is convenient to FÜLLEKRUG ET AL.
©2014. The Authors. normalize the wave numbers k with the free space wave number k 0 = ∕c, where is the frequency of the radio wave and c is the speed of light in vacuum such that k N ∕k 0 ≈ 30.9 and k F ∕k 0 ≈ 5.4 for a radio wave with a frequency of 22 kHz. However, this study is focused on the determination of relative wave numbers k∕k 0 ranging from 0 to 2 in the quasi-static part of the wave number spectrum (Figure 1, right) to ensure that there is no ambiguity in the interpretation of the results.

Wave Number Vectors
The wave number vectors of the radio waves are determined by writing equation (2) in polar form and by arranging the phases into one equation for each of the receivers (x 1,n , x 2,n , 1) where the phase residuals n need to be minimized. Equation (5) is an overdetermined linear system of equations that can be solved with the Gaussian method of least squares to obtain the best fitting solutions for k 1 , k 2 , and .
Four terrestrial radio transmitters for submarine communication are recorded with the array of radio receivers from 15:00:00 to 15:00:02 UTC on 13 May 2011 (Table 1). The wave number vectors are inferred by the use of equation (5). Subsequently, the wave number vectors are binned in a two-dimensional k 1 ∕k 0 × k 2 ∕k 0 grid with a resolution Δ(k∕k 0 ) = 0.01. This distribution function exhibits four distinct maxima, which correspond to the wave number vectors of the four radio communication transmitters (Figure 2, left).

Horizontal Slowness
The length of the wave number vector k = |k| determines the apparent velocity v a = ∕k, which describes the radio wave propagation velocity as it appears to the array in a projection of the three-dimensional wave number vector onto the two dimensions of the array. The apparent velocity of the radio waves is often expressed in the corresponding horizontal slowness s = 1∕v a , which enables a convenient linear assessment of the travel time t d of the radio wave across the array of receivers t d = sd, where d is the distance traveled across the network. For example, an electromagnetic wave in vacuum has the slowness s 0 = 1∕c ≈ 3.34 ns/m. A larger slowness s > s 0 indicates an apparent velocity which is less than the speed of light, and a smaller slowness s < s 0 indicates that the radio wave arrives at the array from an elevation angle Θ. This elevation angle is determined by the ratio between the radio wave propagation speed in the medium near the array v m and the apparent velocity v a such that cos Θ = v m ∕v a , as inferred from a geometric analysis of the plane wave arriving at the array [Rost and Thomas, 2002, Figure 1]. The horizontal slowness of the radio waves from the four radio communication transmitters range from ∼2.65 ns/m to ∼4.13 ns/m (Table 1). These results show that the radio waves of the transmitters arrive from varying elevation angles that result from the distance-dependent superposition of the transverse electric and transverse magnetic radio wave propagation modes (TE n /TM n ), possibly assisted by an influence of the geomagnetic field on the radio wave propagation.

Arrival Azimuth
The orientation of the wave number vector reflects the arrival azimuth Φ of the radio wave, which is determined by the ratio between the east-west (k 1 ) and the north-south (k 2 ) component of the wave number vector such that Φ = arctan (k 2 ∕k 1 ). These arrival azimuths are measured with respect to the geographic east and compared to the expected arrival azimuths that are calculated assuming a propagation from the transmitters to the array along the great circle paths (Table 1). The expected arrival azimuths exhibit deviations from the measured arrival azimuths ΔΦ, which range from ∼−9.7 • to ∼14.5 • ( Table 1). The observed deviations are indicative of an anisotropic wave propagation [e.g., Füllekrug and Sukhorukov, 1999], possibly as a result of the distance traveled at varying angles with the geomagnetic field [Barr et al., 2000, and references therein].

Electromagnetic Source Field
The wave numbers can be used to determine the electromagnetic source field E(t) from the recordings by rewriting equation (2) e −ik r n E = y n + y n ∧ n = 1, 2, … N = 10 (6) where the residuals of the measurements y n need to be minimized. Equation (6) is an overdetermined linear system of equations that can be solved with the complex Gaussian method of least squares, which results in a sparse spatial Fourier transform to determine the best fitting solution for the electromagnetic source field. It is convenient to normalize the electromagnetic source field E(t) with the temporal average of its amplitude < |E(t)| > such that E(t)∕< |E(t)| > describes relative changes of the amplitude and phase of the source field. The resulting normalized source field is binned on the complex plane with a resolution Δ(E∕ < |E| >) = 0.01. This distribution of the electromagnetic source field shows the constellation diagram of the phase modulated transmissions, e.g., for the transmissions from Rhauderfehn (Figure 2, right). The rotation of the constellation diagram with respect to the real and imaginary axes of the complex reference frame results from the fact that the phase is measured relative to the phase 0 at the start time of the recordings. The normalized source field is closely related to the coherency of the electromagnetic source field y n e ik r n |y n e ik r n | FÜLLEKRUG ET AL.
©2014. The Authors. which ranges between 0 ≤ E c ≤ 1 and measures the significance of the electromagnetic source field. The coherency can be used, for example, as a weight for phase coherent stacks to increase the signal-to-noise ratio of average electromagnetic source fields [Schimmel and Paulssen, 1997].
The phase modulation of the radio transmissions results in a deterministic phase progression which is conveniently expressed in the apparent frequency [Taner et al., 1979] where E r and E i are the real and imaginary part of E (Figure 2, right) and Δt=1 μs is the sampling time interval of the recordings. The apparent frequency is measured relative to the center frequency of the transmission (Table 1) because it results from the temporal derivative of the phase in equation (9). This relative apparent frequency covers a total range of ∼±50 Hz, and it exhibits distinct oscillations with a repetition period of ∼5 ms (Figure 3, left), which is observed in the spectrum of |E| as a spectral peak at 200 Hz with seven harmonics (Figure 3, right). Note that the spectral amplitudes of |E| drop over many orders of magnitude down to the numerical limit of the spectral analysis.

Residuals
The residual phases n are composed of two components. The first component is independent of time and frequency, and it is attributed to quasi-static clock offsets. The second component is not correlated between the receivers of the array, and it is attributed to the jitter of the clocks used at the individual receivers. The residual phases are expressed as temporal uncertainties t = ∕ and follow a normal distribution function with a standard deviation ∼15.3 ns (Figure 3, left, inset). The residuals of the electromagnetic source field y n are also composed of two components. The first component is independent of time and frequency, and it is attributed to the calibration of the receivers in the array. The second component is not correlated between the receivers of the array, and it is attributed to the quantization noise of the analogue-to-digital converters used at the individual receivers. The residuals follow a normal distribution function with a standard deviation ∼3.3 μV/m (Figure 3, right, inset). The residual distributions for the time jitter and the quantization noise are in excellent agreement with the technical specifications of the radio receiver [Füllekrug, 2010, Figures 5 and 6]. This important result means that each sample of the measurements is explained by the wave number vector (k 1 , k 2 ) T , and the electromagnetic source field E with an accuracy down to the instrumental limits.

Discussion
In seismic and infrasonic array analysis, the elevation angles are determined from the wave propagation velocity in the medium v m near the array that is not yet very well known for radio waves. This knowledge FÜLLEKRUG ET AL.
©2014. The Authors. could be gained by designing a volumetric, or three-dimensional, array to provide experimental evidence for the elevation angles, for example, by placing radio receivers on a self-configuring array of unmanned aerial vehicles, e.g., quadcopters, or by placing the array into a volcanic crater of suitable size. Alternatively, it might be possible to assess all the radio transmitters within the bandwidth of the radio receiver to determine a most likely, or minimum, wave propagation velocity in the medium near the array toward a determination of the elevation angles. Such a comprehensive assessment of all the elevation angles of electromagnetic radiation from known radio transmitters enables the opportunity to establish a standard reference map of the low-frequency radio sky, similar to known stellar constellations in radio astronomy. The main application of this reference map is to determine the relative location of electromagnetic radiation from radio sources above thunderclouds, such as sprites, gigantic jets, and possibly electron beams. For example, sprites are thought to cause impulsive radio noise associated with the growth and branching of individual streamers. The associated impulsive radio signals have a broad spectrum such that the array analysis to locate radio noise from sprites would use that part of the electromagnetic spectrum which is not occupied by known radio transmissions. The reference map will thereby enable a determination of the location of electromagnetic radiation from sprites relative to the known locations from radio transmitters. Another possibility to assess the radio noise from sprites is to slightly modify the array analysis to apply beam forming. Beam forming steers the antenna aperture toward the expected location of the electromagnetic radiation from sprites above thunderclouds to maximize the sensitivity of the array transfer function in that direction. This strategy requires a prediction of the propagation path of electromagnetic radiation from sprites which is well understood.

Summary
The data analysis concepts developed for seismic and infrasonic arrays are successfully adapted and applied to an array of radio receivers to infer the electromagnetic source fields and the wave number vectors of ∼20-24 kHz radio waves from transmitters for submarine communication. The wave number vectors provide novel information on the horizontal slowness of the radio waves, which range from ∼2.7 ns/m to ∼4.1 ns/m depending on the arrival azimuths of the radio waves. It is expected that more extended array analyses will provide a wealth of novel insights into the physics of low-frequency radio waves, and possibly radio noise from sprites and yet unknown radio sources, with experiments that remain to be conducted in the future.