Determining Ionospheric Drift and Anisotropy of Irregularities from LOFAR Core Measurements: Testing Hypotheses behind Estimation

Abstract

We try to assess the validity of assumptions taken when deriving drift velocity. We give simple formulas for characteristics of the spatiotemporal correlation function of the observed diffraction pattern for the frozen flow and the more general Briggs model. Using Low-Frequency Array (LOFAR) Cassiopeia intensity observation, we compare the experimental velocity scaling factor with a theoretical one to show that both models do not follow observations. We also give a qualitative comparison of our drift velocity estimates with SuperDARN convection maps. The article is essentially an extended version of the conference paper: “Determining ionospheric drift and anisotropy of irregularities from LOFAR core measurements”, Signal Processing Symposium 2021 (SPSympo 2021).

1. Introduction

LOFAR (Low-Frequency Array) is an international research facility aiming mainly at interferometry in low frequencies to observe the early stages of the universe’s evolution, especially the so-called reionization epoch [1]. Astronomical interferometry is a powerful observation tool that allows reconstructing the radiation intensity angular distribution based on radio observations taken at many different observation points simultaneously forming pairwise many interferometers. It can achieve a fantastic angular resolution determined by the largest baseline and wavelength of the observed radiation by the basic formula: $2\lambda /D$, where: $\lambda$ is the wavelength and D is the distance between a pair of receiving stations forming an interferometer (the baseline) [2]. The distribution of points is crucial—while long baselines give desired resolution, a subset of stations is needed to be closely spaced (core stations), reducing the aliasing. The basic idea described above coins the name for the technique as VLBI (Very Large Baseline Interferometry). The reconstruction of the sky from such measurements becomes quite challenging but is possible. Nevertheless, the reconstructed image often contains some distortions caused by the irregular cosmic plasma through which the radiation propagates [2]—the phenomenon is often described as the flickering of observed radio sources. The angular deflection of a source image is due to an ionization gradient across the direction to the source [2,3].

Scintillation is the name of the described phenomenon and, in general, is defined as fluctuations in signal parameters after passing through a medium with heterogeneous distribution of the refractive index [4]. In the case of space plasma, the factors that determine scintillation are electron density distribution and varying magnetic. There are two scintillation regimes: diffractive and refractive scintillation. In fact, both contribute to the so-called diffraction pattern in the observation plane. The former is responsible for the focusing/defocusing effect whilst the latter denotes linear superposition of contributions from the scattering medium [4].

Let us briefly present a formal description of scintillation theory in the so-called diffusion approximation for electromagnetic wave propagation. If the electromagnetic wave’s frequency is far from any characteristic plasma mode frequency (plasma or gyro frequency), and the size of plasma irregularities is much larger than the wave’s wavelength one can apply the scalar equation for the complex wave amplitude u [4]:

(1)$-2i{k}_0\frac{\partial u}{\partial z}+{\nabla }_{\perp }^{2}u+k_0^2ϵ\left(\mathbf{r}\right)u=0$

where u is (possibly complex) slowly varying wave amplitude, $k_0$ is the wavenumber of the incident wave, coordinate z is along wave propagation and perpendicular to irregularities layer, ${\nabla }_{\perp }^2$ is the Laplace operator in the plane perpendicular to ${\mathbf{k}}_0$ and $ϵ\left(\mathbf{r}\right)$ is the fluctuating part of plasma dielectric permittivity.

If the thickness L of the irregularity layer (Figure 1) satisfies the condition $\sqrt{\lambda L}\ll r_0$, where $\lambda$ is the wavelength of incident wave and $r_0$ is the outer scale of irregularities (the largest size of irregularities), the layer distorts mainly the phase of the incident wave. It is the so-called phase screen approximation. The wave subsequently propagates in the free space and generates a complicated interference pattern at a distance from the irregularity layer [4]. Wave propagation in the phase screen approximation can be expressed using the Fresnel diffraction formula [4]:

(2)$u({\mathbf{r}}_{\perp },z)=u_0\frac{i{k}_0}{2\pi z}\int d{\mathbf{r}}_{\perp }^{\prime }exp\left[-i\left(\frac{k_0{({\mathbf{r}}_{\perp }-{\mathbf{r}}_{\perp }^{\prime })}^2}{2z}+\delta \varphi \left({\mathbf{r}}_{\perp }^{\prime }\right)\right)\right].$

where $u({\mathbf{r}}_{\perp },z)$ is the complex wave field observed at a distance z from the irregularity layer, $u_0$ is an amplitude of incident plane wave characterized by the wavenumber $k_0$ and ${\mathbf{r}}_{\perp }$ is vector perpendicular to ${\mathbf{k}}_0$. The initial phase fluctuations $\delta \varphi$ on a screen are related to fluctuations of the total electron content $\Delta N_e$ by the formula:

(3)$\delta \varphi \left({\mathbf{r}}_{\perp }\right)=-r_e\lambda \int dl\frac{\Delta N_e({\mathbf{r}}_{\perp },l)}{1-{\left(\frac{f_p}{f}\right)}^2},$

where $r_e$ is the classical electron radius, f—the wave’s frequency, and $f_p$—the plasma frequency of the medium.

In the situation where the thickness of the irregularity layer does not satisfy the requirements for the phase screen approximation, it is still possible to use phase screen approximation. However, to do so, certain conditions must be met. The medium needs to be divided into multiple layers, each treated as the phase screen, where radiation emerging from the preceding layer is incident upon the next one [5].

Scintillation of distant cosmic radio emissions can provide interesting information on the cosmic medium itself, its internal spatial structure (spatial power spectra), basic evolution characteristics (drift velocity and its dispersion), etc. [6,7,8] In this paper, we focus on the drift velocity estimation of a diffraction pattern observed on the ground using scintillation mode observation of LOFAR radio telescope [1,9]. We review the most important work on the correlation method of drift estimation. A large number of consistent measurements of LOFAR signal amplitude allows the validation of basic model assumptions when estimating the drift velocity. We discuss possible explanations of our findings concerning existing work and future research perspectives.

2. Materials and Methods

The evolution model of a diffraction pattern evolution can be composed of two parts: evolution of the scatter (the medium on which the radiation is scattered) and a model for radio wave propagation to relate quantities estimated in the experimental reference frame to interesting quantities describing a medium through which waves propagate. Let us consider at first the simplest scenario. The so-called frozen-in motion of a scalar field $\psi$ (i.e., rigid motion of a medium) can be described by:

(4)$\begin{array}{c}\hfill \psi (\mathbf{r},t)=\psi (\mathbf{r}-\mathbf{v}t)\end{array}$

and the phase screen propagation model that results in diffraction pattern drifting with the same velocity as the scatter [10]. In general, the problem of determining medium evolution based on scintillation measurements is an example of inverse problems theory and, to our knowledge, has no general solution. Derivation of medium evolution involves multipoint measurements analysis [11]. When the stationarity (homogeneity) assumption is reasonable, the two methods are used to derive evolution characteristics: the cross-correlation analysis [12], and the dispersive (cross-spectral) analysis [7]. First relates to the characteristic features of auto- and cross-correlation function, second makes use of the cross-spectrum phase, giving information about dispersion. By Wiener-Khinchin theorem, both can be related to one another.

Let us assume the field that is observed undergoes simple rigid motion, which means its temporal evolution follows Equation (4) with constant velocity ${\mathbf{v}}_F$, where subscript ${}_F$ stands for “frozen flow” estimate. It means, under the assumption of homogeneity of the field, that cross-correlation of sampled time series of the field at positions ${\mathbf{r}}_{\mathbf{1}}$, ${\mathbf{r}}_{\mathbf{2}}$ at time instants $t_1$, $t_2$ will be:

(5)$<\psi ({\mathbf{r}}_1-{\mathbf{v}}_F{t}_1)\psi ({\mathbf{r}}_2-{\mathbf{v}}_F{t}_2)>=C(\mathsf{\zeta }-{\mathbf{v}}_F\tau ),$

where C is the auto-correlation of the random field:

(6)$C\left(\mathsf{\zeta }\right)=<\psi \left({\mathbf{r}}_1\right)\psi \left({\mathbf{r}}_2\right)>$

and $\mathsf{\zeta }={\mathbf{r}}_1-{\mathbf{r}}_2$, $\tau =t_1-t_2$ which we call the displacement coordinates.

From (3) and (4), one can see that if irregularities are in a frozen flow, it is so for any quantity formed from u since the formulas do not involve temporal changes. In our analysis, we will use the observed intensity of a radio source, which in our case means: $\psi =u{u}^{*}$.

Next simplification considers modelling of anisotropy of the field. We assume that anisotropy can be described by taking the argument of the correlation function as the quadratic form $q\left(\zeta \right)={\zeta }^T{Q}_s\zeta$, whose matrix $Q_s$ is symmetric positive-definite [6]. Then the correlation function can be written as:

(7)$C\left(\mathsf{\zeta }\right)=\rho \left({\mathsf{\zeta }}^{{}^T}{Q}_s\mathsf{\zeta }\right),$

where $\rho$ describes desired radial dependence of the correlation function. Having established relation (7), we are able to assess the drift velocity of the pattern by relating features of temporal cross-correlations to the features of the field correlation function. One of the most frequently used is maximum of temporal cross-correlation taken at two positions separated by vector $\mathsf{\zeta }$:

(8)$\begin{array}{c}\hfill \frac{\partial }{\partial \tau }C(\mathsf{\zeta }-{\mathbf{v}}_F\tau )=\dot{\rho }\frac{\partial }{\partial \tau }\left[{(\mathsf{\zeta }-{\mathbf{v}}_F\tau )}^T{Q}_s(\mathsf{\zeta }-{\mathbf{v}}_F\tau )\right]=\\ \hfill 2\dot{\rho }\left[{\mathbf{v}}_F^T{Q}_s{\mathbf{v}}_F\tau -{\mathsf{\zeta }}^T{Q}_s{\mathbf{v}}_F\right]=0.\end{array}$

Provided $\dot{\rho }=\partial \rho /\partial q$ is non zero, we obtain the expression for ${\tau }_m$, which is the time delay for maximum of the correlation function, and instantly for its gradient with respect to the displacement coordinates $\mathsf{\zeta }$:

(9)${\tau }_m=\frac{{\mathsf{\zeta }}^T{Q}_s{\mathbf{v}}_F}{{\mathbf{v}}_F^T{Q}_s{\mathbf{v}}_F}\to {\nabla }_{\mathsf{\zeta }}{\tau }_m=\frac{Q_s{\mathbf{v}}_F}{{\mathbf{v}}_F^T{Q}_s{\mathbf{v}}_F}.$

To derive the explicit formula for ${\mathbf{v}}_F$ let us observe that:

(10)${\mathbf{v}}_F^T{\nabla }_{\mathsf{\zeta }}{\tau }_m=1$

and:

(11)${\left({\nabla }_{\mathsf{\zeta }}{\tau }_m\right)}^T{Q}_s^{-1}{\nabla }_{\mathsf{\zeta }}{\tau }_m=\frac{1}{{\mathbf{v}}_F^T{Q}_s{\mathbf{v}}_F}.$

Combining these we arrive at:

(12)${\mathbf{v}}_F=\frac{Q_s^{-1}{\nabla }_{\zeta }{\tau }_m}{{\left({\nabla }_{\zeta }{\tau }_m\right)}^T{Q}_s^{-1}{\nabla }_{\zeta }{\tau }_m}.$

It means, that knowing the matrix $Q_s$, and ${\nabla }_{\zeta }{\tau }_m$, which are to be estimated from data, one can compute the drift velocity $\mathbf{v}$. The presented formulas align with those given by Briggs [6] when dropped purely temporal decorrelation. A large number of available measurement positions allows a direct approach, while the small number of positions of observations motivated numerous works to overcome this limitation by comparing correlations at many different time lags, see for example [13,14].

We present a simple numerical example illustrating described modelling. A set of virtual receiving stations (black circles in Figure 2) samples a random field that is drifting past. An example of received signals is depicted in Figure 3. Some mutual relations between each other can be observed. In Figure 4a cross-correlations for each pair of signals for zero time lag have been plotted as a function of the separation vector. To assess matrix Q we took least square fit of the function $q({\zeta }_x,{\zeta }_y)=q_{11}{\zeta }_x^2+2q_{12}{\zeta }_x{\zeta }_y+q_{22}{\zeta }_y^2$ to estimated correlations near $\mathsf{\zeta }=(0,0)$. Resulting ellipsis has been plotted at the same figure. Clearly, the correlation captures well basic spatial properties of the sampled field.

Figure 4b,c show time lag for maximum cross-correlation as a function of separation (it is clearly linear as formula (9) predicts). Now, combining the information on structure geometry with time delay gradient, we were able to estimate drift velocity (red strip in Figure 4c). The true drift velocity (black on Figure 4b) is very close to the estimated (red on Figure 4c). For comparison, we also give an isotropic estimate represented as the blue strip (the estimate of drift velocity obtained using the assumption of isotropy of random field—the quadratic form matrix is a scalar multiple of identity matrix in this case).

We will take an important note at this point. The above example shows that we reconstruct the spatio-temporal correlation for a given random field as:

(13)$C(\mathsf{\zeta },\tau )=\frac{1}{2T}{\int }_{-T}^{T}dt\psi (\mathbf{r},t)\psi (\mathbf{r}+\mathsf{\zeta },t+\tau )$

As pointed out by Briggs in the exposition of his full correlation analysis [6], the basic frozen flow drift estimation presented above can lead to serious overestimation of velocity magnitude in the presence of temporal decorrelation. When additional temporal decorrelation is present the evolution of a field can be described by a spatio-temporal quadratic form that takes the following form in the rest frame:

(14)$Q=\left[\begin{array}{cc}{Q}_s& 0\\ 0& \beta \end{array}\right].$

It also describes an elliptic shape of the isosurfaces of spatio-temporal correlation. The $Q_s$ is the spatial part of the full matrix, while $\beta$ describes temporal fading. When the field is moving with the velocity v the matrix of the quadratic form transforms into: $Q^{\prime }=T^{T}QT$, where T is the Galilean transformation to a new coordinate system moving with the velocity $\mathbf{v}$:

(15)$T=\left[\begin{array}{ccc}1& 0& {\mathit{v}}_{\mathit{x}}\\ 0& 1& v_y\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{cc}I& \mathbf{v}\\ {\mathbf{0}}^T& 1\end{array}\right].$

Writing Q’ explicitly, it reads:

(16)$Q^{\prime }=\left[\begin{array}{cc}{Q}_s& Q_s\mathbf{v}\\ {\mathbf{v}}^T{Q}_s& \beta +{\mathbf{v}}^T{Q}_s\mathbf{v}\end{array}\right].$

The important property of the Briggs model is that the correlation function depends on $\mathsf{\zeta }$ and $\tau$ through the same functional form $\rho$:

(17)$<\psi ({\mathbf{r}}_1,t_1)\psi ({\mathbf{r}}_2,t_2)>=C(\mathsf{\zeta },\tau )=\rho \left({[\mathsf{\zeta },\tau ]}^T{Q}^{\prime }[\mathsf{\zeta },\tau ]\right).$

For this model we obtain the gradient of time lag to the maximum of cross-correlation:

(18)${\nabla }_{\zeta }{\tau }_m=\frac{Q_s\mathbf{v}}{\beta +{\mathbf{v}}^T{Q}_s\mathbf{v}}$

For frozen flow and Briggs models, the gradient is linear, so this property cannot be used to discriminate them. The experimentally assessed quantities: ${\nabla }_{\zeta }{\tau }_m$ and $Q_s$ are the same, thus between the frozen flow velocity estimate ${\mathbf{v}}_F$ and $\mathbf{v}$ the scalar relation should hold:

(19)${\mathbf{v}}_F={\gamma }_F\mathbf{v}\to {\gamma }_F=1+\beta /{\mathbf{v}}^T{Q}_s\mathbf{v},$

where the scaling factor $\gamma$ gives the measure of overestimation of the drift velocity, when a pure frozen flow estimator is taken.

The relation (17) provides another way of assessing the drift velocity:

(20)$\mathbf{v}=Q_s^{-1}{Q}_{[1:2;3]}^{\prime }.$

Briggs model enables differences in functional dependence between spatial and temporal correlations [6]; however, the simple form of dependence allows for the construction of relatively simple and practical algorithms for estimating the drift velocity [12,14]. In these algorithms, the arguments of the temporal cross-correlation functions for different positions are compared when the function values are equal. This leads to methods that are essentially similar to the Formula (20). The same form of dependence for spatial and temporal changes facilitates the analysis and construction of practical algorithms, but lacks a solid physical foundation. This problem was discussed in the works of Little and Ekers [15], where the authors considered the dispersion of structures motion, while Wernik et al. [10] studied the importance of ionospheric propagation effects in the case of the vertical drift gradient in the scattering layer.

Data used for this analysis come from LOFAR core and remote station observations [1]. The Low-Frequency Array (LOFAR) is an excellent astronomical instrument, and a handy tool for studying irregularities in the ionosphere [1]. Due to its operational frequency range (10–270 MHz), LOFAR is very sensitive even to tiny changes in ionospheric electron density. The instrument’s interferometric nature allows for multi-point observations, thus giving the possibility for ionospheric scintillation measurements over distances ranging from tens of meters to hundreds of kilometres. The project for scintillation monitoring over the LOFAR stations has been carried out for several years, and a large amount of data has been collected and stored in the Long Term Archive (LTA, https://lta.lofar.eu (accessed on 13 July 2022)). Available data contains signal amplitude for a few strongest radio sources measured at all core and remote stations in beam forming mode (not as an interferometer). Based on the LTA data, correlation analysis between stations can be done in order to obtain information about the characteristics of ionospheric structures.

The Cassiopeia signal amplitude was observed under three various geomagnetic conditions (the key parameter describing geomagnetic conditions—Dst index, is shown in Figure 5: a for L547449; b for L547785; c for L552177 set respectively). L547449 observations took place during the calm period, the L547785 is a phase of relaxation of a small geomagnetic disturbance, and L552177 was measured during the main phase of the magnetic storm (Table 1 gathers the details concerning the data sets). The amplitude dynamic spectra for these periods are shown on Figure 6—one can see an increasing density of the structures observed which may be related to the increasing speed of ionospheric structures. Positions of ionospheric pierce points (IPPs) are shown as a red line pieces in Figure 7 (IPPs were calculated at an altitude of 350 km).

3. Results

The Equation (7) implies that the correlation function should satisfy the basic equation:

(21)$\frac{\partial }{\partial \tau }C-\mathbf{v}·{\nabla }_{\mathsf{\zeta }}C=0,$

which describes the homogenous drift of the correlation function in the displacement coordinates. Figure 8 illustrates such behaviour of the correlation function for the L547785 data set. One can see that the correlation function seems to be in movement (the plots show correlations for subsequent time instants).

The observed signals were analysed as described earlier in the numerical example. Figure 9a,b,c show correlations for zero time shift as a function of stations position difference (the displacement coordinates) with the estimated quadratic form ellipse overlaid for different datasets, respectively. Geometry for all three cases seems to be well-defined, and their sizes do not depend on geomagnetic conditions. Figure 9a1,b1,c1 show the time shift of the maximum cross-correlation as the function of the separation, the determined velocity vectors (red stripe) and the isotropic estimate, where Q is a scalar multiple of the identity matrix, (blue strip). Black strip means vector speed with magnitude successively: 100 m/s, 1000 m/s, 1000 m/s for the reference. The exact results are gathered in Table 2. One can observe the linear dependence of the time shift as a function of the separation vector.

The frozen drift values shown in the table seem to be overestimated, especially for disturbed conditions. It may result from non-stationarity. The conditions are changing during the time period or an incorrect model used for the estimation. To check this, we divided the fragments into shorter segments of 3.4 min (2048 points), overlapping by half their length. We also estimated $Q^{\prime }$ matrix to compute drift velocity according to the Formula (20). In order to do so, for large correlation values of $C(\mathsf{\zeta },\tau )$ around 0, we fit the polynomial of the second degree in three variables (${\zeta }_x,{\zeta }_y,\tau$):

(22)$C({\zeta }_x,{\zeta }_y,\tau )\approx q_{11}{\zeta }_x^2+q_{22}{\zeta }_y^2+q_{33}{\tau }^2+2q_{12}{\zeta }_x{\zeta }_y+2q_{13}{\zeta }_x\tau +2q_{23}{\zeta }_y\tau .$

This allows calculating Q’ and the velocity ${\mathbf{v}}_T$. The Figure 10, Figure 11 and Figure 12 show the results for both estimators and the source elevations for considered periods. In green, we marked the periods where we made the comparison with SuperDARN radar convection pattern assessments. We limited the plots to values less than 2000 m/s.

Having many pairs with a significant correlation at hand, it is possible to check the validity of the assumption of frozen drift. According to the Formula (5), cross-correlation in the frozen flow case should be the intersection of the spatial correlation that goes through the separation vector for a given pair in the direction of the velocity vector. We constructed such an intersection by taking computed spatial correlations in a strip distant from the mentioned straight line no more than d. We introduced spatial coordinates $\mathsf{\zeta }={\mathsf{\zeta }}_0+{\mathbf{v}}_F\tau$ for the temporal cross-correlation to compare the two quantities. Figure 13 shows such a comparison for the selected data segment. One can see that the red points (spatial correlations) do not lie on the blue line showing temporal cross-correlation between two stations. As we pointed out, the frozen drift assumption can lead to overestimating the drift velocity. We tried to determine the factor by which ${\mathbf{v}}_F$ should be multiplied so that the curves overlap. For this purpose we minimise the following quantity:

(23)$\sum _k{\left[C_s\left({\mathsf{\zeta }}_k\right)-C_t\left(({\mathbf{v}}_T/\gamma ){\tau }_k\right)\right]}^2$

with respect to $\gamma$, where $\gamma$ gives the experimental velocity scaling factor. We also have given the velocity scaling factor for Briggs model (20) ${\gamma }_F$, which we can get from quadratic form $Q^{\prime }$:

(24)${\gamma }_F=\frac{q_{33}}{{\mathbf{v}}^T{Q}_s\mathbf{v}}$

At this point we used the estimate (20) of the drift velocity that - in the Briggs model - should be consistent with (24). Figure 11a shows experimental velocity scaling factor $\gamma$ vs. ${\gamma }_F$ (20) given by Briggs model. The different markers refer to different data sets: black diamonds—L552177; blue circles—L547785; red triangles—L547449. One can see a quite large discrepancy between the theoretical scaling and experimental one, moreover, the markers representing different data sets are distributed nonuniformly. For a comparison Figure 11a gives experimental velocity scaling factor $\gamma$ vs. the ratio of velocity estimators magnitude $v_F/v_T$. In the latter case, one can see that the velocity modulus quotient is much closer to the experimental scaling factor than the theoretical value, and the points corresponding to the different data sets are evenly distributed. In Figure 12 the relation between azimuths of ${\mathbf{v}}_F$ and ${\mathbf{v}}_T$ is shown. They roughly agree, however a minority has large deviations, which may suggest that the relation between estimates is not scalar. The points grouping suggests that estimates have mainly zonal component, which is in agreement with the SuperDARN observations (Figure 14, Figure 15, Figure 16 and Figure 17).

We also present a comparison of the obtained results with the velocity estimates that use other experimental principles. For comparison, we will use the available plots showing the convection of the northern polar ionosphere modelled on the basis of measurements by the SuperDARN system [16]. Figure 14, Figure 15, Figure 16 and Figure 17 show the velocities overlayed on the ionospheric convection maps given by SuperDARN system around the moments of time marked with green lines in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 respectively. The green colour strip gives the frozen flow estimate ${\mathbf{v}}_F$ while the blue one the ${\mathbf{v}}_T$. We keep exact mutual relation between estimates while abandoning accurate comparison to the SuperDARN drift velocity magnitudes. The reason is twofold: we would like to highlight our results that otherwise would be obscured if given in scale. The second reason is that the location of our estimates is outside of the validity of SuperDARN modelling yet still close enough to bear some features of large-scale polar ionosphere convection. Inevitably, the comparison is qualitative for the reasons we mentioned, but they can be used to trace possible drift at mid-geomagnetic latitudes. Figure 14 shows that for chosen time instants of the L547449 data set our estimates agree with distant polar convection at least around 18:22 UT. Similarly, Figure 15 shows an agreement of our estimates with the polar ionosphere convection for the L547785 data set and Figure 16 for L552177. Figure 17 shows signatures of reaching a flow stagnation point around 16:45 UT. Here, additionally, we display SuperDARN convection assessment before and after the event. A tongue of eastward circulation reaches far west.

Although the estimates of the ionospheric drift velocity follow the convection pattern assessed by SuperDARN system for mid-latitude ionosphere may seem exaggerated for sets L547785 and L552177. To illustrate, that such values are still possible we present the drift estimations made by the ionosphere soundings at Pruhonice (49°59′N 14°32′E), Juliusruh (54°37′N 13°22′E) and Fairford (51°41′N 1°46′W) during the storm of 13.10.2016 which are respectively: 976 +/− 130 m/s (at 16:49 UT), 405 +/− 18 m/s (at 17:15 UT), 833 +/−146 m/s (at 17:53 UT) (https://giro.uml.edu/driftbase/ (accessed on 13 July 2022)). This illustrates possible ionospheric variability in space and time during disturbed conditions at mid-latitudes.

4. Conclusions

This article presents an experimental test of the assumptions of the models used to determine the ionospheric drift velocity. In our particular case of determining the drift velocity of the diffraction pattern created at the Earth’s surface by Cassiopeia radiation and measured by the LOFAR system. We briefly presented the specifics of the observed field. The basic models and methods for estimation of drift velocity have been described. We also designed our original algorithm version that uses many measurement points and a direct correlation function approximation in the vicinity of the origin of the reference frame. This approach may be a source of some approximation errors, but if only their effect can be described by the multiplication of the quadratic form by a constant scalar, it should not affect the drift velocity estimate. We did not perform any detailed analysis of errors and their sources, nevertheless, the consistency of the results and their similarity to the results of other experimental methods serve as a vindication of some sort. The algorithm itself is independent of the hypotheses taken by other methods, and this property allowed us to test them.

In many cases, we have realised that the assumption about the frozen-in evolution of ionospheric irregularities is violated. We adopted a model with the diffraction pattern’s temporal decorrelation, which gave more reliable drift velocity estimates. This was accomplished by fitting a three-dimensional polynomial to the spatio-temporal correlations obtained from LOFAR’s scintillation amplitude measurements. We compared our estimates with SuperDARN convection maps. The scintillation-based velocity assessment reproduces quite well ionospheric convection features captured by SuperDARN. Using a comparison of theoretical scaling relation between frozen and decorrelation estimators to actual scaling, we have also shown that the Briggs model doesn’t describe spatio-temporal behaviour of the correlation function correctly. Nevertheless, it gives a method of estimating the drift velocity that agrees reasonably well with other independent estimates (SuperDARN, ionosondes) and has been used often for this purpose and still is [17].

Author Contributions

Conceptualization, M.G. and H.R.; Formal analysis, M.G. and K.B.; Investigation, M.P.; Methodology, M.P.; Supervision, H.R.; Writing—review & editing, B.M., D.P., L.T. and H.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

LOFAR data used in this work is available through the LOFAR Long Term Archive (LTA): https://lta.lofar.eu/ (accessed on 13 July 2022). SuperDARN convection maps are available through http://vt.superdarn.org/tiki-index.php?page=Conv+map+overview (accessed on 13 July 2022). Ionosonde drift data are accessible at https://giro.uml.edu/driftbase/ (accessed on 13 July 2022).

Acknowledgments

The authors acknowledge the use of SuperDARN data. SuperDARN is a collection of radars funded by national scientific funding agencies of Australia, Canada, China, France, Italy, Japan, Norway, South Africa, the United Kingdom, and the United States of America. This paper uses ionospheric data from the USAF NEXION Digisonde network, the NEXION Program Manager is Annette Parsons.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Illustration of the scattering by a screen.

Figure 2. Numerically generated random field and the receiving stations marked as black circles.

Figure 3. The data sets obtained from virtual receiving stations from Figure 2.

Figure 4. Colour-coded estimate of the correlation function as a function of separation—(a), and time lag to a maximum of temporal correlation in separation coordinates, black, red and blue lines give accordingly: true velocity, velocity estimate that takes into account anisotropy of structures and isotropic velocity estimate—(b,c).

Figure 5. Dst index for considered datasets: (a)—L547449; (b)—L547785; (c)—L552177.

Figure 6. Amplitude dynamic spectra: (a)—L547449 set; (b)—L547785 set; (c)—L552177 set.

Figure 7. Daily positions of IPPs: (a)—L547449 set; (b)— L547785 set; (c)—L552177 set. The red line segment gives IPPs for observation intervals.

Figure 7. Daily positions of IPPs: (a)—L547449 set; (b)— L547785 set; (c)—L552177 set. The red line segment gives IPPs for observation intervals.

Figure 8. Correlation for subsequent time delays: (a)—0; (b)—10 s; (c)—20 s.

Figure 9. Correlations for zero time lag as a function of separation vector for datasets: (a)—547449, (b)—L547785, (c)—L552177, with quadratic form ellipsis overlaid; time lag for maximum of cross-correlation as a function of separation vector for datasets: (a1)—547449, (b1)—L547785, (c1)—L552177, with velocity estimates: red - taking into account anisotropy of structures; blue - isotropic. Black vectors indicate velocities with magnitude: 100 m/s, 1000 m/s, 1000 m/s for sets 547449, L547785, and L552177, respectively.

Figure 10. The drift velocity magnitude and azimuth as a function of time for the L547449 dataset for two estimators: (12)—(a) and (20)—(b); panel (c) shows the elevation of the source.

Figure 11. Collective plot for all considered datasets (blue circles—L547449, red triangles—L547785, black diamonds—L552177): (a)—experimental scaling factor [Forumla omitted. See PDF.] vs the theoretical derived from Briggs model, (b)—experimental scaling factor [Forumla omitted. See PDF.] vs. ratio of estimators magnitude [Forumla omitted. See PDF.].

Figure 12. Azimuth of the frozen flow drift velocity estimate [Forumla omitted. See PDF.] vs. azimuth of [Forumla omitted. See PDF.].

Figure 13. An example of the procedure for obtaining the experimental velocity scaling factor: (a)—a bar containing spatial correlations taken for comparison, (b)—spatial correlations superimposed on temporal cross-correlation between a pair of stations (they do not overlap), (c)—spatial correlations superimposed on temporal cross-correlation between a pair of stations after the optimal speed scaling—minimising the expression (22).

Figure 14. Qualitative caption. comparison of drift estimation by two methods, according to formula (12)—green—and formula (20)—blue—with ionospheric convection obtained by the SuperDARN system for moments: marked in Figure 10a,b with green vertical lines, for the quiet time (set L547449); panel (a) for 17:42 UT, and panel (b) for 18:22 UT of 2016.09.20.

Figure 15. Qualitative caption. comparison of drift estimation by two methods, according to formula (12)—green, and formula (20)—blue, with ionospheric convection obtained by the SuperDARN system for moments: marked in Figure 18a,b with green vertical lines, for the quiet time (set L547785); panel (a) for 3:54 UT, panel (b) for 5:30 UT, and panel (c) for 6:54 UT of 2016.09.25.

Figure 15. Qualitative caption. comparison of drift estimation by two methods, according to formula (12)—green, and formula (20)—blue, with ionospheric convection obtained by the SuperDARN system for moments: marked in Figure 18a,b with green vertical lines, for the quiet time (set L547785); panel (a) for 3:54 UT, panel (b) for 5:30 UT, and panel (c) for 6:54 UT of 2016.09.25.

Figure 16. Qualitative caption. comparison of drift estimation by two methods, according to formula (12)—green—and formula (20)—blue—with ionospheric convection obtained by the SuperDARN system for moments: marked in Figure 19a,b with green vertical lines, for the quiet time (set L552177); panel (a) for 16:24 UT, and panel (b) for 16:48 UT of 2016.10.13.

Figure 17. Details of the convection pattern change around 16:48 UT for the set L552177; panel (a) for 16:42 UT just before estimation time 16:48 UT (panel b), and panel (c) for 16:50 (just after estimation time).

Figure 18. The drift velocity magnitude and azimuth as a function of time for the L547785 dataset for two estimators: (12)—(a) and (20)—(b); panel (c) shows the elevation of the source.

Figure 19. The drift velocity magnitude and azimuth as a function of time for the L552177 dataset for two estimators: (12)—(a) and (20)—(b); panel (c) shows the elevation of the source.

Table 1

Datasets selected for the analysis and description of the geomagnetic conditions over the time of observation.

Data Set	Date & Time of Start	Date & Time of Stop	Short Characteristics	Frequency
L547449	2016-09-20 17:30:00	2016-09-20 18:40:00	quiet	61.8 MHz
L547785	2016-09-25 03:30:00	2016-09-25 06:55:00	small disturbance	27.7 MHz
L552177	2016-10-13 15:55:00	2016-10-13 17:30:00	storm —main phase	61.8 MHz

Table 2

Drift velocity estimates obtained from selected datasets. Table results correspond to the graphical presentation in Figure 4.

Data Set	${\mathit{v}}_{\mathit{x}}$	${\mathit{v}}_{\mathit{y}}$	Magnitude	${\mathit{v}}_{\mathit{x}}$ Isotropic	${\mathit{v}}_{\mathit{y}}$ Isotropic	Magnitude Isotropic
	[m/s]	[m/s]	[m/s]	[m/s]	[m/s]	[m/s]
L547449	142	−35	147	86	−90	124
L547785	694	408	804	549	552	778
L552177	−2160	795	2370	−1526	1586	2200

References

1. Haarlem, M.; Wise, M.; Gunst, A.; Heald, G.; McKean, J.; Hessels, J.; de Bruyn, A.; Nijboer, R.; Swinbank, J.; Fallows, R. et al. LOFAR: The low-frequency array. Astron. Astrophys.; 2013; 556, 53. [DOI: https://dx.doi.org/10.1051/0004-6361/201220873]

2. Thompson, A.; Moran, J.; Swenson, G., Jr. Interferometry and Synthesis in Radio Astronomy; Springer: Berlin/Heidelberg, Germany, 1991; 1, pp. 11-36. [DOI: https://dx.doi.org/10.1007/978-3-319-44431-4]

3. Mevius, M.; van der Tol, S.; Pandey, V.N.; Vedantham, H.K.; Brentjens, M.A.; de Bruyn, A.G.; Abdalla, F.B.; Asad, K.M.; Bregman, J.D.; Brouw, W.N. et al. Probing ionospheric structures using the LOFAR radio telescope. Radio Sci.; 2016; 51, pp. 927-941. [DOI: https://dx.doi.org/10.1002/2016RS006028]

4. Yeh, K.C.; Liu, S.J. Radio wave scintillation in the ionosphere. Proc. IEEE; 1982; 70, pp. 324-360.

5. Spivack, M.; Uscinski, B. The split-step solution in random wave propagation. J. Comput. Appl. Math.; 1989; 27, pp. 349-361. [DOI: https://dx.doi.org/10.1016/0377-0427(89)90022-8]

6. Briggs, B.H. On the analysis of moving patterns in geophysics. I: Correlation analysis. J. Atmos. Terr. Phys.; 1968; 30, pp. 1777-1788. [DOI: https://dx.doi.org/10.1016/0021-9169(68)90097-4]

7. Briggs, B.H.; Vincent, R.A. Spaced-antenna analysis in the frequency domain. Radio Sci.; 1992; 27, pp. 117-129. [DOI: https://dx.doi.org/10.1029/91RS03051]

8. Coles, W.A.; Kaufman, J.J. Solar wind velocity estimation from multi-station IPS. Radio Sci.; 1978; 13, pp. 591-597. [DOI: https://dx.doi.org/10.1029/RS013i003p00591]

9. Fallows, R.A.; Bisi, M.M.; Forte, B.; Ulich, T.; Konovalenko, A.A.; Mann, G.; Vocks, C. Separating nightside interplanetary and ionospheric scintillation with LOFAR. Astrophys. J.; 2016; 828, L7. [DOI: https://dx.doi.org/10.3847/2041-8205/828/1/L7]

10. Wernik, A.W.; Liu, C.H.; Yeh, K.C. Modeling of spaced-receiver scintillation measurements. Radio Sci.; 1983; 18, pp. 743-764. [DOI: https://dx.doi.org/10.1029/RS018i005p00743]

11. Grzesiak, M.; Wernik, A.W. Dispersion analysis of spaced antenna scintillation measurement. Ann. Geophys.; 2009; 27, pp. 2843-2849. [DOI: https://dx.doi.org/10.5194/angeo-27-2843-2009]

12. Costa, E.; Fougere, P.F. Cross-spectral analysis of spaced-receiver measurement. Radio Sci.; 1988; 23, pp. 129-139. [DOI: https://dx.doi.org/10.1029/RS023i002p00129]

13. Costa, E.; Fougere, P.F.; Basu, S. Cross-correlation analysis and interpretation of spaced-receiver measurements. Radio Sci.; 1988; 23, pp. 141-162. [DOI: https://dx.doi.org/10.1029/RS023i002p00141]

14. Rino, C.L.; Livingston, R.C. On the analysis and interpretation of spaced-receiver measurements of transionospheric radio waves. Radio Sci.; 1982; 17, pp. 845-854. Available online: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/RS017i004p00845 (accessed on 13 July 2022). [DOI: https://dx.doi.org/10.1029/RS017i004p00845]

15. Ekers, R.D.; Little, L.T. A method for analysing drifting random patterns in astronomy and geophysics. Astron. Astrophys.; 1971; 10, pp. 306-309.

16. Ruohoniemi, J.M.; Baker, K.B. Large-scale imaging of high-latitude convection with Super Dual Auroral Radar Network HF radar observations. J. Geophys. Res. Space Phys.; 1998; 103, pp. 20797-20811. Available online: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/98JA01288 (accessed on 13 July 2022). [DOI: https://dx.doi.org/10.1029/98JA01288]

17. Wang, J.; Morton, Y.T. A comparative study of ionospheric irregularity drift velocity derived from a GNSS receiver array and Poker Flat Incoherent Scatter Radar measurements during high-latitude ionospheric scintillation. J. Geophys. Res. Space Phys.; 2017; 122, pp. 6858-6881. Available online: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017JA024015 (accessed on 13 July 2022). [DOI: https://dx.doi.org/10.1002/2017JA024015]