© P. Flisek et al., Published by EDP Sciences 2023

1 Introduction

The propagation of radio waves through the Earth’s ionosphere is affected by the presence of spatial inhomogeneities in the electron density distribution. After propagating through ionospheric plasma inhomogeneities radio waves exhibit fluctuations on their received amplitudes and phases as a result of their wavefront being scattered by the change in refractive index associated with the irregularities. A relative drift between ray path and irregularities leads to the observation of temporal fluctuations on the intensity and phase of the received radio waves: this phenomenon is known as radio wave scintillation.

Whilst satellite radio waves are typically limited to specific and discrete frequencies (e.g. beacons at VHF/UHF such as the historical Wideband, NNSS, and Tsikada, and spread-spectrum signals in the L-band such as GNSS), the Low-Frequency Array (van Haarlem et al., 2013) is capable of detecting and measuring radio wave frequencies closely spaced and over a wide range of frequencies at VHF. Although its primary use is as an interferometer for astronomical imaging of the radio sky, it can also be utilised to detect scintillation from irregularities occurring in the inner-heliosphere and in the mid-latitude ionosphere: in the latter case, irregularities can form, for example, as a consequence large-to-small scale Travelling Ionospheric Disturbances (TIDs) (Hernandez-Pajares et al., 2012) in conjunction with instabilities such as the Perkins’ mechanism (Kelley, 2009; Fallows et al., 2020).

The amount of radio wave scintillation depends upon various factors, including the spatial gradient in the electron density and the radio wave frequency. Generally, the amount of radio wave scintillation decreases with increasing radio wave frequency. For example, the scintillation level was observed to decrease as f⁻ⁿ, with n determined experimentally on the basis of early beacon satellite measurements (Crane, 1977). As different scintillation levels are typically observed at different radio wave frequencies propagating along the same line of sight, a general frequency dependence suggests that there can be various gradients in the ionospheric electron density spatial distribution where the change in electron density can take place over different spatial scales and, thus, originate different levels of scintillation at different radio wave frequencies.

As an example, irregularities in the mid-latitude ionosphere are capable of inducing scintillation on radio wave frequencies in the VHF range: however, the gradients associated with these irregularities are typically too low to induce scintillation on radio wave frequencies in the L-band, typical of Global Navigation Satellite Systems (GNSS). Therefore, as scintillation on different radio wave frequencies is sensitive to different scales (i.e. electron density gradients) in the ionosphere, the observation of scintillation between the VHF and L band enables better detection of irregularities and their spatial and temporal evolution (Fallows et al., 2014; de Gasperin et al., 2018).

The LOFAR (Low-Frequency Array) radio telescope is an interferometer and it is composed currently of 52 stations located in many parts of Europe (van Haarlem et al., 2013). Most of the stations are located in the Netherlands, forming a dense network referred to as the Core and a number of so-called Remote stations. Currently, there are also 14 stations forming an extensive network called ILT (International LOFAR Telescope). Two more ILT stations are under construction.

Each station consists of up to 3264 omnidirectional dipole antennas (in full ILT configuration). The antennas are divided into two separate types: Low Band Antennas (LBA) operating in the frequency range from 10 to 90 MHz and High Band Antennas (HBA) which are able to receive signals at the frequency from 110 to 240 MHz. HBA antennas are grouped in 16 pairs of dipoles within special tiles. Signals from individual dipoles are sampled and digitized by using a 200 MHz clock, generating raw data in the region of approximately up to 10 Gbits/s: raw data are distributed to an analysing system through a dedicated network. In addition, an electronically-controlled beamforming is also utilized.

Due to the frequency range in which LOFAR detectors work, the state and dynamics of the ionosphere can have a significant impact on the result of the observations. During the observation mode, the Stokes I parameter is utilized to estimate the intensity of radio waves at various frequencies, thus creating various channels of time series of radio wave intensities. The display of radio wave intensities as a function of time simultaneously over different radio wave frequencies (or channels) is known as the dynamic spectrum of the observations.

Poland operates three LOFAR stations (in Borowiec – PL610, Łazy – PL611 and Bałdy – PL612) (Fig. 1) together with a dense network of GNSS geodetic permanent stations (also complemented by a GNSS ionospheric and scintillation monitor) – ASG-EUPOS. This particular configuration of instruments allows the presence of multiple ionisation scales in the mid-latitude ionosphere to be investigated by comparing measurements of scintillation from LOFAR with measurements of scintillation and rate of change of TEC from GNSS.

Figure 1 The positions of Polish LOFAR stations: PL612, PL611, PL610.

The analysis presented here intends to address the following questions: (1) what are the ionisation scales which GNSS and LOFAR are sensitive to and (2) how scintillation varies between the VHF and the L-band. Based upon two distinct case studies characterised by disturbed magnetic conditions, a methodology was developed for the analysis and detection of ionospheric structures by estimating the amount of scintillation originating from ionospheric irregularities as observed through LOFAR radio telescopes. The methodology was then validated by means of a comparison between co-located observations from LOFAR and GNSS. This validation provides insights on how traditional GNSS ionospheric observations can be augmented by means of LOFAR VHF scintillation measurements.

2 Data and methodology

During scintillation observation campaigns in the international mode, LOFAR typically observed 3 natural astronomical targets: Cassiopeia A (CasA), Cygnus A (CygA), and Taurus A (TauA). The observations utilised here are part of the program “Monitoring Scintillation above LOFAR”, and are stored on the LOFAR Long Term Archive (LTA) https://lta.lofar.eu/Lofar under project codes LC7_001 and LC8_001. LOFAR collected data with 10 ms time resolution for 100 frequency channels, each with a bandwidth of 195 kHz, sparsely covering the range 21.77–76.07 MHz. This analysis focused on data collected through the Polish LOFAR stations at Baldy (PL612), Borowiec (PL610) and Łazy (PL611) (Krankowski et al., 2014).

The methodology developed here was validated by comparing LOFAR VHF scintillation observations with ionospheric measurements of scintillation and TEC fluctuations obtained through GNSS ground receivers. The understanding of the sensitivity of the two instruments to ionospheric irregularities enables us to appreciate how modern LOFAR observations of ionospheric scintillation can augment traditional GNSS ionospheric observations. Standard RINEX 30 s observables available from the IGS network https://cddis.nasa.gov/archive/ (Johnston et al., 2017) as well as from the ASG-EUPOS network (over 100 stations distributed evenly across the area of Poland) (Bosy et al., 2007) were utilised to estimate the spatial and temporal distribution of ionospheric structures in conjunction with LOFAR observations. The ionospheric structures were estimated by means of temporal variations (and geographical maps) of the rate of change of the TEC.

Furthermore, a GNSS scintillation receiver (Septentrio PolaRxS Pro) located 200 m from the LBA section of the LOFAR Baldy station further provided finer scale observations of ionospheric structures occurring during LOFAR observations. The scintillation receiver provided multi-constellation observations with a time resolution of up to 100 Hz.

Two case studies were selected for this analysis on the basis of the simultaneous availability of data from both the GNSS scintillation monitor and the LOFAR Polish stations, as well as the general geomagnetic conditions which were characterised through the Kp and Dst indices. In each of these two case studies, one disturbed day was considered and compared with a quiet reference day.

2.1 Intensity scintillation from LOFAR

The amount of scintillation observed with LOFAR on the intensity of radio waves was quantified by means of the S₄ index, according to Briggs & Parkin (1963):

(1)

where I is the intensity of the received signal; 〈 〉 in general denotes ensemble averaging, but it is in practice approximated with time average.

LOFAR observations of radio wave intensities utilised in this analysis were sporadically affected by RFI. In the very first step of estimating the level of scintillation originated by ionospheric structures on the radio wave frequencies detected by LOFAR, RFI-induced outliers in the data (i.e. spikes in the estimates of the radio wave intensity) need to be mitigated. The method described by Fallows et al. (2020) was incorporated for this purpose. In the RFI mitigation process, the median filter for each frequency band was applied. The threshold for the RFI detection was set at the level of the 10th percentile (5σ) for each channel. Spikes remaining after the filtering and larger than the threshold were finally cut out from the dynamic spectra.

“Clean”, RFI-free intensities in each channel were detrended and normalized in order to obtain zero-mean normalized intensity, which allows us to estimate the temporal fluctuations on the radio wave intensities induced by scintillation. Detrending was done by subtracting a moving average with a 3-minute window. The zero-mean intensity fluctuations I_normalized are given by:

(2)

Detrending removes possible trends from the intensity observations whereas normalization removes the zero-frequency spectral component (Forte et al., 2022). Finally, the S₄ index was calculated by taking the standard deviation of the zero-mean intensity fluctuations in equation (2). The standard deviation was calculated over 3 min in the case of a 3-minute moving average (as utilized in the detrending) but output over a sliding window every one minute: a value of S₄ every minute was output in order to compare LOFAR VHF S₄ values with GNSS measurements (typically output every minute). The subsequent steps utilized in the estimate of the S₄ scintillation index are illustrated in Figure 2.

Figure 2 LOFAR scintillation spectra at subsequent stages of elaboration: (a) the dynamic spectrum (i.e., raw intensity measurements), (b) RFI-free dynamic spectrum, (c) zero-mean normalized intensities and (d) the estimated S4 index.

The time interval of 3 min utilized in the moving average for the estimate of the zero-mean intensity fluctuations was considered as a compromise between the need to appreciate how scintillation varied in time (e.g., due to the variability of the ionospheric irregularities traversed) and the approximation of ergodicity. This compromise is illustrated in Figure 3 which shows the 1D Power Spectral Density (PSD) of the zero-mean normalised intensity fluctuations for DOY271 of 2017 (this was one of the case studies considered, as detailed below), collected during the time interval 16:45–18:30 UT. The 1D PSDs refer to the zero-mean intensity fluctuations from the channel corresponding to a radio wave frequency of 48.92MHz. From Figure 3 it can be noticed that the 1-min PSDs (Fig. 3a) do not exhibit a fully developed low-frequency limit as opposed to the PSDs estimated over 2–5 min (Figs. 3b–3e). The PSDs calculated over 2–5 min contain a better resolution of the low-frequency limit and are characterized by a higher frequency resolution due to a higher number of samples. Differences between the PSDs calculated over different temporal intervals can be noticed from Figure 3: for example, the 2-min PSDs show differences when compared with the PSDs calculated over 3–5 min. These differences are due to the lack of ergodicity in the zero-mean intensity fluctuations. The longer the time interval the larger the spatial distances over which the observations would average through (Forte et al., 2022): for example, in the case of plasma irregularities in the Earth’s ionosphere, a relative drift (i.e., irregularities drifting relative to the ray path) of 100 m/s would imply that irregularities distributed over approximately 18 km (transverse to the ray path direction) contributed to the 3-min PSDs (i.e., 100 m/s · 60 s · 3). The largest spatial scale contributing to the zero-mean intensity fluctuations (or the low-frequency limit in their PSD) corresponds to the outer scale: the turbulent spectrum that characterizes the irregularities extends between an outer scale and an inner scale. However, under the weak scattering approximation irregularities inducing scintillation have a spatial scale smaller than the Fresnel scale and it can be assumed that radio waves traverse plasma density irregularities distributed along a phase-changing screen transverse to the ray path. The Fresnel scale is given by:

(3)

where λ is the wavelength and z is the distance to a hypothetical phase-changing screen in the ionosphere that approximates the scattering experienced by radio waves during their propagation through ionospheric irregularities. Figure 4 shows the Fresnel scale for LOFAR VHF radio wave frequencies in Figure 2. In this case, it was assumed that scintillation originated in the F region (i.e., z = 350 km). Under the weak scattering approximation for irregularities in the F region, a time interval of 3 min was assumed to provide a good compromise between the need to appreciate the temporal variability of scintillation, the spatial scales inducing scintillation, and the approximation of ergodicity. A shorter time interval (e.g., 2 min) would indeed have a low-frequency limit closer to the Fresnel scale, whereas larger time intervals (e.g., 4–5 min) would imply averaging over spatial scales larger than the Fresnel scale with the implication that the temporal variability would be smoother. Hence, the moving average in equation (2) and the standard deviation in equation (1) were estimated over a time interval of 3 min for the LOFAR VHF scintillation observations considered here.

Figure 3 Power Spectral Density (PSD) estimated for the radio wave frequency 48.92 Mhz on zero-mean intensity fluctuations observed on DOY271 in 2017 during the time interval 16:45–18:30 UT. The PSDs were calculated over different time intervals: (a) 1 min, (b) 2 min, (c) 3 min, (d) 4 min and (e) 5 min.

Figure 4 Fresnel scale variation with frequency.

In the case of GNSS, the estimate of the S₄ scintillation index follows a similar method: S₄ is estimated by means of the standard deviation of the normalized intensity fluctuations. The normalization and standard deviation are typically estimated over a time interval of 1 min, which removes any trend related to satellite motion. The Fresnel scale in the case of the GNSS L1 radio wave frequency is approximately of the order of 365 m (assuming again a phase screen at a distance of 350 km from the receiver) and, in the presence of weak scattering, spatial scales smaller than 365 m contribute to scintillation. In addition, the contribution from thermal noise present in the receiver is automatically removed (van Dierendonck et al., 1993) to provide a corrected estimate of the S₄ scintillation index.

In the case of GNSS observations, the corrected S₄ scintillation indices are output every minute. Both under the weak scattering approximation and in the presence of multiple (strong) scattering, the difference between the spatial scales contributing to scintillation at VHF and L band make the two instruments sensitive to irregularities forming over different spatial scales in the ionosphere because the PSD of intensity fluctuations covers different ranges of spatial frequencies (Forte, 2008, 2012b; Forte et al., 2022). Therefore, the comparison between co-located scintillation observations from GNSS and LOFAR allows us to detect irregularities forming over a wider range of spatial scales in the ionosphere, than typically detected by using GNSS alone.

On the other hand, irregularities forming over spatial scales larger than the outer scale contribute to phase fluctuations, as measured for example on GNSS signals. Therefore, the comparison between LOFAR and GNSS was extended to include not only scintillation indices but also a measure of phase fluctuations. Phase fluctuations can be quantified by utilising dual-frequency phase observations from GNSS: the geometry-free combination of dual-frequency carrier phase observations provides an estimate of the TEC. The difference of TEC over consecutive epochs (or rate of change of TEC) provides a measure of temporal fluctuations in phase observations induced by irregularities in the ionosphere.

2.2 Intensity scintillation and TEC fluctuations from GNSS

Co-located with the LOFAR station PL612 at Baldy was a GNSS ionospheric monitor which outputs scintillation indices and uncalibrated Slant TEC every minute, together with standard RINEX observables and 50-Hz raw estimates of slant TEC, signal intensity, and signal phase. In the analysis described here the presence of ionospheric irregularities was detected by comparing enhancements in the S₄ scintillation index observed through LOFAR with any enhancements in both scintillation and the rate of change of TEC as observed through GNSS. In the case of PL612, LOFAR scintillation indices at VHF, GNSS scintillation indices at L-band, and GNSS Rate of Change of TEC (ROT) were compared to provide insights on the presence of ionospheric irregularities as well as on the variation of the S₄ index over a wide interval of radio wave frequencies.

GNSS ROT for a specific satellite in view was calculated as (Pi et al., 1997):

(4)

where TEC_k and TEC_i+k are the (uncalibrated) Slant TEC at epochs k and k + 1, respectively. GNSS ROT was estimated over 1-minute intervals. In the case of PL612, GNSS ROT was estimated over two different temporal intervals (i.e., 1 min and 1 s) to provide more information on the scale size of the ionospheric structures detected.

The comparison of LOFAR 1-minute S₄ (i.e., output every 1 min and based on a 3-minute time interval for moving average and standard deviation) across all radio wave frequencies with co-located GNSS TEC fluctuations estimated over 1 min and 1 s enabled an estimate of the scale size of the ionospheric structures detected, thus providing insights on the sensitivity of the two instruments to ionisation gradients present in the ionosphere and on the possibility to combine their measurements for space weather monitoring purposes.

2.3 General considerations about signal levels

In the validation of the methodology through comparison of scintillation observations collected through LOFAR and GNSS it is worth considering the following aspects in relation to the overall signal levels involved. The two instruments collect radio signals in two different ways: whilst LOFAR observations of particular radio objects are based on the process of beamforming (van Haarlem et al., 2013; Błaszkiewicz et al., 2016), GNSS extracts and estimates information from radio signals transmitted from satellites through a demodulation process.

In the case of LOFAR, the width of the beam created during the beamforming process has an angular size of 4° for the LBA antennas (Błaszkiewicz et al., 2021). The signal-to-noise of a particular object drives the detection: typically, the signal from a given object is considered as detected when the signal level exceeds the 3σ threshold (where σ indicates the average noise level). The signal-to-noise ratio can be increased by using long exposures, however, in the case of scintillation observations (such as those illustrated here) only very bright sources are used. In particular, standard LOFAR scintillation observations consider three bright radio sources from the A-team: i.e., supernova remnants Cassiopeia A (CasA) and Taurus A (TauA) as well as radio galaxy Cygnus A (CygA). The absolute flux density of these sources is well established: at the radio wave frequency of 50 MHz, the flux density measures 27,104 Jy for CasA, Cyg A 22,146 Jy for CygA, and 2008 Jy for Tau A (where 1 Jy = 10⁻²⁶W · m⁻²Hz⁻¹) (de Gasperin et al., 2020). As a comparison, the flux density of the Sun at 50 MHz is approximately of the order of tens of kJy during quiet conditions and approximately of the order of several MJy during disturbed conditions (Ho et al., 2008). Out of these three sources, in this study, only observations from CasA and CygA were presented because TauA was at very low elevation angles during the time intervals considered. In general, as the noise level varies from a few Jy up to tens of Jy as a function of the zenith angle (Błaszkiewicz et al., 2018), the signal-to-noise ratio for CasA and CygA is approximately of the order of 2710 and 2214, respectively: assuming a noise level of few tens of Jy, these values amount approximately to 34 dB and 33 dB, respectively.

On the other hand, GNSS receivers estimate the signal level from the receiver’s tracking stage by combining in-phase and quadrature samples in time (Van Dierendonck, 1995). Estimates of GNSS signal nominal levels are found to be approximately in the region of 40–50 dB (from arbitrary units) as estimated from tracking stages, with drops observed in the presence of scintillation (Forte, 2012a, b). Typical signal-to-noise ratios of GNSS signals are found to be approximately in the region of 40 dB-Hz in absence of scintillation-induced fading (Parkinson & Spilker, 1996).

Despite differences in the way LOFAR and GNSS process radio signals and estimate signal levels, the above considerations on the overall signal-to-noise ratio allows to assume that the signals detected by both instruments are high enough above the background noise level to give a meaningful ionospheric measurement. Given possible differences in the absolute signal levels, the comparison between scintillation observations from LOFAR and GNSS illustrated here was attempted not by considering the raw observations (i.e., the value of I in Eq. (1)) rather by considering the zero-mean normalized intensity fluctuations (i.e., I_normalized in Eq. (2)). The normalization allows to compare signatures on signals even if the average signals’ levels may be different (e.g., between different GNSS satellites or between GNSS and LOFAR).

3 Results

Two case studies were utilised to illustrate the methodology of detecting ionospheric structures by means of LOFAR VH scintillation measurements: the results were validated by comparing LOFAR scintillation measurements and ground GNSS ionospheric observations (i.e. scintillation and TEC fluctuations). Table 2 describes the cases selected for this analysis, including the time interval of the LOFAR observations as well as Kp and Dst indices which describe the geomagnetic conditions during the observation. Each case study was characterised by measurements from a magnetically disturbed day as well as from a magnetically quiet day: the latter was utilised as a reference which the measurements from the disturbed day can be compared to.

3.1 Case study 1: DOY087 and DOY090, 2017

The results corresponding to these particular case studies are illustrated in Figure 3 (PL612), Figure 5 (PL610), and Figure 6 (PL611). Stronger scintillation was detected on the radio emission from CygA compared to CasA, something that appears to be in common amongst all case studies presented here. TauA was also observed but remained below 20° elevation angle throughout the observation period.

Figure 5 LOFAR S4 of Cassiopeia A (a, b) and Cygnus A (c, d) calculated for days 87 (quiet) and 90 (disturbed) in 2017, recorded on LOFAR station PL612. GNSS ROT values were calculated for 1 s (e, f) and 60 s (g, h) for days 87 and 90 of 2017. GNSS S4 index (i, j) observed on L1 (blue dots), L2 (red dots) and L5 black dots) frequencies recorded with the ionospheric monitor co-located with the PL612; LOFAR scintillation indices for the 48.92 MHz channel (from Figs. 4a to 4d) are also shown as pink area (Cas A) and blue area (Cyg A). Gaps in the GNSS scintillation indices were due to issues related to data downloading.

Figure 6 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated for days 87 and 90 of 2017, recorded on LOFAR station PL610. ROT values (e, f) calculated for 60 s for days 87 and 90 of 2017 from observations recorded by receiver near PL610 station.

PL612 on DOY090 2017 detected strong scintillation in the morning (Fig. 5). Higher S₄ values appeared between 05:30 UT and 06:05 UT, reaching a maximum value of up to approximately 0.35, overall. Localised higher S₄ values appeared at 05:50 UT, in conjunction with a solar flare with maximum flux B3.3: the white stripe in Figure 5 indicated the S₄ value corresponding to the flare that has been removed. DOY087 2017 was the closest quiet reference (Fig. 5) and it shows very low S₄ values, approximately about 0.05 without any variation through the whole observation period. Figures 6 and 7 show observations collected by the other Polish LOFAR stations (PL610 and PL611 respectively). At both PL610 and PL611, S₄ values appeared to be very similar to those observed at PL612.

Figure 7 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated for days 87 and 90 of 2017, recorded on LOFAR station PL611.

Figures 5 and 6 show GNSS ROT for PRNs: G03, G07, G09, G23, G28, G30, R05, R07, R08, R09, R10, R15, R16 – these are the PRNs visible from the co-located GNSS receivers with line of sight closer to the line of sights of CygA and CasA. In the case of the PL611 LOFAR station, no GNSS data were available for the case studies considered.

Although no significant enhancement in 60 s and 1 s ROT across all colocated stations could be appreciated, an increase in the occurrence of cycle slips can be noticed after 6:45 UT with cycle slips visible in both 60 s and 1 s ROT. In comparison, no cycle slips were visible on the quiet reference day, which suggests that cycle slips were of ionospheric origin. In the case of PL610, GNSS observations from the nearby ASG-EUPOS receiver BOR1 (located 600 m from the PL610 station) showed the occurrence of cycle slips in the same time interval as observed from Baldy (PL612) scintillation receiver (co-located with PL612).

Given the proximity of the lines of sight from the PRNs considered for the GNSS measurements and those from CasA/CygA considered for LOFAR measurements, it is plausible to assume that LOFAR stations and co-located GNSS receivers were observing similar ionospheric irregularities. It is interesting to observe that LOFAR seemed to detect enhancements in S₄ scintillation indices somehow earlier than GNSS detected an increase in the occurrence of cycle slips.

Figure 5 also illustrates GNSS S₄ scintillation indices at L1, L2, and L5 (Figs. 5i–5j) observed on the same PRN links considered for the GNSS ROT (Figs. 5e–5h) through the ionospheric monitor co-located with the LOFAR PL612 station. Scintillation at the L band was very low which suggests that electron density gradients were forming mainly over larger spatial scales (i.e., larger than the Fresnel scale at L band).

3.2 Case study 2: DOY271 and DOY275, 2017

The results corresponding to these particular case studies are illustrated in Figure 8 (PL612) and Figure 9 (PL611). DOY271 2017 was the most disturbed day among the cases considered here. The PL612 station detected an enhancement in scintillation on CygA after 17:45 UT with S₄ exceeding 0.5. Similarly to case study 1, CasA showed lower scintillation than CygA (Figure 8). TauA was again below the elevation angle of 20° and too low to provide meaningful comparison with GNSS observations. Interestingly, the scintillation index S₄ showed higher values for higher frequencies on both targets. On the quiet reference day DOY275 2017 (Fig. 8), LOFAR S₄ appeared overall at a lower level. The measurements at PL611 (Fig. 9) appeared very similar to those recorded at PL612. LOFAR S₄ index was high throughout the observation period with some bursts to even higher values. At PL611 higher S₄ indices occurred in the mid-frequency interval (i.e., 40–60 MHz) in contrast to PL612 where S₄ indices were higher at higher frequencies (i.e., 50–75 MHz). On the other hand, Cassiopeia A observations recorded at PL611 showed lower S₄ indices than those recorded at PL612: the highest value of S₄ for PL611 observations was 0.25 (Fig. 9), meanwhile S₄ for PL612 reached 0.45 (Fig. 8) in the case of Cassiopeia A.

Figure 8 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated respectively for days 271 and 275 of 2017, recorded on LOFAR station PL612. ROT values calculated for 1 s (e, f) and 60 s (g, h) for days 271 and 275 of 2017. GNSS S4 index observed on L1, L2 and L5 frequencies recorded with the ionospheric monitor collocated with the PL612 LOFAR station superimposed on CasA and CygA S4 index calculated for the 48.92 MHz channel (i, j).

Figure 9 LOFAR S4 of Cassiopeia A (a, b) and Cygnus A (c, d) calculated respectively for days 271 and 275 of 2017, recorded on LOFAR station PL611.

GNSS 60 s ROT showed higher values at PL612 (Fig. 8) throughout the whole observation time on DOY271. After 18:00UT 60 s ROT increased approximately to 0.5 TECu/min. No cycle slips were observed in this case as opposed to case study 1. GNSS 60s ROT on DOY275 (a quiet reference day) showed lower values between −0.2 and 0.2 TEC/min. Higher 1 s ROT values occurred on DOY271 but not on DOY275.

Figure 8 also illustrates GNSS S₄ scintillation indices (Figs. 8i–8j) observed on the same PRN links considered for the GNSS ROT (Figs. 8e–8h) through the ionospheric monitor co-located with the LOFAR PL612 station. Similarly to case study 1 (Fig. 5), scintillation at L band was very low in case study 2. The very low values of GNSS S₄ scintillation indices suggests again that electron density gradients were forming mainly over spatial scales larger than the Fresnel scale at L band.

Overall, in every selected case study some enhancements in GNSS ROT and in LOFAR S₄ appeared during disturbed conditions. In the case of PL612 enhancements in GNSS ROT could be observed only over a temporal interval of one minute; no noticeable ROT enhancements were observed over a shorter temporal interval of 1 s, suggesting that the ionospheric irregularities should have a spatial scale of the order (and larger than) of few kilometres in the direction across the ray path. This can be seen by assuming a value of approximately 100 m/s for the relative velocity between ray path and the irregularities at F region heights: the distance covered by the ray path in 1 min is 100 m/s · 60 s = 6000 m: however, as irregularities can have higher drift velocities this value represents a lower limit to this estimate.

4 Discussion

The observations presented a methodology capable of detecting ionospheric structures by observing and quantifying the radio wave scintillation that they induce and that it is observed through LOFAR radio telescopes. The methodology was validated by comparing scintillation measurements from three LOFAR stations in Poland with (nearly) co-located GNSS scintillation and ROT measurements. Overall, some enhancements in scintillation detected through the LOFAR stations on radio wave frequencies received from CygA and CasA tended to occur during magnetically active conditions, with S₄ indices generally lower in the case of CasA than in the case of CygA. Similarly, some enhancements in 60 s GNSS ROT for those PRNs with a line of sight closer to CygA and CasA tended to occur during more active conditions (in some cases together with cycle slips), whereas no enhancement was observed on GNSS 1 s ROT. On the other hand, the GNSS S₄ scintillation index (as estimated through a GNSS scintillation monitor co-located with the LOFAR PL612 station) remained at very low values throughout the cases considered.

The comparison between LOFAR VHF scintillation indices, GNSS L band scintillation indices, and GNSS ROT indicates the type of spatial scales over which electron density gradients were forming in the middle-latitude ionosphere during the case studies considered here. The overall ionospheric conditions during the two case studies considered can be appreciated by means of meridional plots of the ROT Index (ROTI) (Cherniak et al., 2014). ROTI was calculated as the standard deviation of ROT values over a time interval of 5 min (Pi et al., 1997):

(5)

where the ROTI at a given epoch was given by the standard deviation taken over the 10 preceding epochs (5 min, given 30 s RINEX observations). The ROTI meridional plots (similar to keograms) were calculated for a selected meridian with the latitudinal step of 0.2° in order to assess the evolution of ionospheric structures both in latitude and in time.

In order to appreciate the longitudinal evolution of ionospheric structures as well, Figures 10–12 show meridional plots for three selected meridians: 20° (longitude referring to the LOFAR station PL612) together with 17° and 24° as east- and westward references. Figures 10 and 11 show meridional plots for DOY087 (quiet reference day) and DOY090 (disturbed day) for case study 1, respectively. Similarly, Figures 12 and 13 show results for DOY275 (quiet reference day) and for DOY271 (disturbed day) for case study 2, respectively.

Figure 10 ROTI meridional plot for DOY 87 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).

Figure 11 ROTI meridional plot for DOY 90 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).

Figure 12 ROTI meridional plot for DOY 275 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).

Figure 13 ROTI meridional plot for DOY 271 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).

Tenuous enhancements in ROTI (up to approximately 0.1 TECU/min) tended to occur in the disturbed days of the two case studies considered. Diagonally-shaped structures on the meridian plots (e.g., DOY 090 2017 03:00–06:00 UT, Figs. 11a and 11b) seem to suggest that ionospheric structures were generally moving southward with time, whereas widespread enhancements could indicate structures forming locally. Whilst in some cases enhancements in LOFAR VHF scintillation appeared to be consistent with the tenuous enhancements in ROTI (e.g., PL612 in DOY271, 18:00-18:30 UT, Fig. 8), in other cases, LOFAR scintillation enhancements appeared during the absence of noticeable enhancements in ROTI (e.g., PL612 in DOY090, 05:30-06:00 UT, Fig. 5).

The case studies discussed here indicate that enhancements in scintillation on LOFAR VHF radio wave frequencies did not always correspond to enhancements in scintillation on GNSS L-band radio wave frequencies. This implies that electron density gradients developed over spatial scales smaller than the LOFAR Fresnel scale (under the assumption of weak scattering), but they were not prominent over spatial scales smaller than the GNSS Fresnel scale: i.e., ionisation gradients over spatial scales smaller than few hundred metres were not intense enough to induce enhancements in scintillation on GNSS L-band signals.

In order to illustrate the comparison between the S₄ computed from LOFAR and GNSS scintillation receiver the correlation coefficients have been calculated. The coefficients were made with the Pearson method, where values vary between −1 and 1 (−1 is linear anticorrelation, 1 is linear correlation and 0 means no correlation). As the PRNs passes were short, which made the coefficients irrelevant, the average values of GNSS S₄ have to be taken into the correlation. The averages were made for each system (GPS, GLONASS and GALILEO) individually. Each S₄ value of the PRN has been averaged in time. In the case of S₄ produced by the LOFAR, the middle frequency has been chosen.

Results of the correlations are presented in Table 2. The coefficients were estimated for each case and each source separately. The strongest correlation (around 0.4) between GNSS and LOFAR S₄ is visible in case 271/2017 on CygA, however, it still should be considered low. The rest of the cases show no correlation for both of the quiet days as well as the disturbed once. It should be noted that the level of S₄ for both the LOFAR and GNSS is different, from what is visible in Figure 8.

As the inertial subrange for LOFAR VHF (i.e., the spatial scales developing in the presence of turbulence between an outer and an inner scale) is different from the inertial subrange for GNSS L-band, scintillation detected by the two instruments is sensitive to irregularities forming over different spatial scales. Tenuous enhancements in GNSS 60 s ROT can be attributed to ionospheric gradients forming over spatial scales larger than approximately 6 km in the horizontal direction (i.e., assuming a relative drift between ray path and irregularities of 100 m/s, the spatial scale would be 100 m/s × 60 s = 6000 m) and extending over a wider range of altitudes because irregularities with spatial scales larger than the outer scale can originate phase fluctuations (Forte et al., 2017; John et al., 2021). Given that GNSS ROT does not show any enhancement over 1 s intervals, then these ionospheric gradients (inducing tenuous enhancements on GNSS 60 s ROT) are likely to have a spatial scale, transverse to the ray path direction, that is of the order of or larger than approximately 6 km (by accounting for the GNSS ray path scan velocity relative to the ionospheric drift); their spatial scale along the ray path direction is likely to be of tens of kilometres (Forte et al., 2017; John et al., 2021).

However, LOFAR VHF scintillation was induced by irregularities with electron density fluctuations distributed over smaller spatial scales, in the VHF inertial sub-range. These gradients seemed to be not intense enough over smaller spatial to induce scintillation at the L band.

Therefore, in relation to question (1), this aspect suggests that LOFAR VHF scintillation measurements can detect ionospheric irregularities with spatial scales of approximately up to 3 km and distributed over at least 6 km horizontally and over several tens of kilometres vertically. These irregularities are not necessarily detected by means of the GNSS ROT, which implies that LOFAR VHF scintillation measurements offer a higher sensitivity to weaker electron density gradients occurring in the ionosphere than GNSS ROT (or scintillation) measurements, where these gradients originate a rather small signature. That is, whilst the enhancement in LOFAR VHF scintillation tends to be distinct (and higher than the noise level), the enhancement in GNSS scintillation and ROT tends to remain low (and closer to the noise level).

This aspect is connected with the evidence on question (2), where LOFAR VHF scintillation can show enhancements on the S₄ scintillation index (induced by ionospheric irregularities with spatial scales in the VHF inertial sub-range of the electron density spatial fluctuations), whereas GNSS L-band scintillation shows very low S₄ scintillation index values (indicating that electron density fluctuations with spatial scales in the L-band inertial sub-range tend to be not intense enough to induce scintillation at L band in the mid-latitude ionosphere). These aspects need to be considered if observations from LOFAR and GNSS were to be combined for wider ionospheric monitoring.

The considerations above apply whenever LOFAR VHF scintillation originated by ionospheric structures, which introduces a further point in the discussion: i.e., whether the enhancements in LOFAR VHF scintillation that were not consistent with GNSS ROT observations were of ionospheric or other (e.g., interplanetary) origin. In the presence of weak scattering, the spatial scales originating scintillation are those smaller than the Fresnel scale. The scattering occurring on the wavefront of radio waves propagating through plasma density irregularities translates into temporal fluctuations when there is a relative drift between the ray path and the irregularities. In this case, the Fresnel temporal frequency is given by:

(6)

The Fresnel frequency depends upon the distance z to the hypothetical phase screen (which approximates weak scattering) and the relative drift V^REL. Various values of f_F can be determined by different combinations of the parameters V^REL and z (Forte et al., 2022).

Figure 14 illustrates the PSDs from CasA and CygA scintillation observations considered within the two case studies. The Fresnel frequency (i.e., the frequency at which the PSDs start to roll off according to a power law in double logarithmic scale) appears to remain consistent throughout the observations, with a value approximately between 10⁻² and 10⁻¹ Hz. Although these values of the Fresnel frequency can be originated by plasma density irregularities both in the ionosphere (with moderate-to-high relative drift) and in the inner heliosphere (with moderate-to-high drift), the fact that differences can be observed in the observations from different LOFAR stations seems to suggest that the observations considered throughout case studies 1 and 2 are more likely to be of ionospheric origin (Forte et al., 2022).

Figure 14 Power Spectra density plots for presented cases. Every PSD plot is made individually for every LOFAR station (PL610, PL611 and PL612) and for each target (CasA and CygA). All PSDs are calculated for middle channel.

In order to ascertain the origin of the LOFAR scintillation observations presented here, the crosscorrelation functions between pairs of available LOFAR stations for each source were estimated. From the peak of the cross-correlation function it is possible to estimate the drift of the scintillation pattern, which coincides with the drift of the irregularities originating scintillation. The cross-correlation function was estimated by means of the cross-power spectral density: before taking its inverse Fast Fourier Transform, the cross-power spectral density was band-pass filtered (between 0.02 Hz and 1 Hz) in order to remove noise, following the method described in Fallows et al. (2016, 2020).

The cross-correlation functions corresponding to different combinations of all available stations and sources for the two case studies presented here are shown in Figures 15 and 16: these figures only show the cross-correlation function for the LOFAR radio wave frequency that exhibited the highest peak (48.92 MHz). The lags corresponding to the cross-correlation peak are approximately between 372 s and 652 s and correspond to the differences in the observed LOFAR VHF S₄ (Figs. 5–9). Considering the different baselines distances, the drift velocities corresponding to the cross-correlation peaks are approximately between 538 m/s and 994 m/s, which are values typical of ionospheric drifts (Tsugawa et al., 2004; Borries et al., 2009; Panasenko et al., 2019). ACE and Wind observations of the solar wind speed at the time of the observations reveal a flow of the order of approximately 600 km/s (Fig. 17): typical solar wind speeds are of the order of a few hundred km/s (Asai et al., 1998; Ondoh & Marubashi, 2001).

Figure 15 Cross-correlation functions (CCF) values for each baseline and both sources, DOY 090/2017.

Figure 16 Cross-correlation functions (CCF) values for each baseline and both sources, DOY 271/2017.

Figure 17 Proton density (upper panels), proton density ratio (middle panels) and solar wind speed (lower panels) obtained from ACE and Wind satellites for DOY 87 of 2017 (a), DOY 90 of 2017 (b), DOY 275 of 2017 (c) and DOY 271 of 2017 (d). The green bands indicate the exact times during which the LOFAR observations were collected.

The size of the ionospheric structures, their shape, and their direction of motion with respect to the size and orientation of the baselines considered here can account for differences in the cross-correlation peak values and their occurrence at positive/negative lags. The estimate of the components of the ionospheric drift based on the baselines considered here (such as discussed in Fallows et al. (2020)) was not attempted here: the main aspect to consider is indeed the order of magnitude of the drift which is plausible to be of ionospheric origin for the case studies considered here.

5 Conclusions

LOFAR radio telescopes constitute a cutting-edge instrument in modern radio astronomy, operating at several tens of sites and providing a pathfinder to the Square Kilometre Array Observatory (SKAO). The latest upgrade allows for systematic measurements aimed at space weather monitoring. This study established a novel methodology that allows LOFAR to detect and characterize ionospheric irregularities by measuring the VHF radio wave scintillation that they induce. This novel methodology is capable of estimating the S₄ scintillation index attributable to ionospheric irregularities by accounting for non-ergodicity in the measurements in conjunction with the typical VHF inertial sub-range where electron density irregularities can induce scintillation.

Measurements from co-located ground GNSS receivers and LOFAR stations in Poland were compared in the presence of ionospheric irregularities to validate the detection of ionospheric irregularities by means of LOFAR VHF scintillation observations. GNSS L-band scintillation indices and GNSS ROT were compared with scintillation indices measured through LOFAR over a wide range of VHF radio wave frequencies received from the radio objects CasA and CygA. Some enhancements in LOFAR VHF S₄ indices and in GNSS 60 s ROT tended to occur during moderately disturbed magnetic conditions (not necessarily in a consistent way), however the electron density gradients associated with these ionospheric irregularities were too weak to enhance GNSS L-band scintillation.

Measurements of LOFAR VHF scintillation, GNSS L-band scintillation, GNSS 60 s ROT, and GNSS 1 s ROT evaluated in two case studies seem to suggest that the corresponding ionospheric irregularities appeared to form over spatial scales of the order of at least a few kilometres across the ray path and extending over a wider range of altitudes: some of these structures can be detected through LOFAR better than through GNSS. Measurements of LOFAR VHF scintillation can indeed be utilised for the detection of ionospheric irregularities characterised by spatial scales of approximately up to 3 km and distributed over at least 6 km horizontally and over several tens of kilometres vertically. The gradient in electron density associated with these structures may be enough to induce scintillation at VHF and enhancements in GNSS ROT or it may be enough to induce scintillation at VHF but not enough to induce enhancements in GNSS ROT. This aspect suggests that LOFAR VHF scintillation measurements have a higher sensitivity to ionospheric gradients than GNSS.

When scintillation observed through LOFAR radio telescopes is of ionospheric origin, LOFAR VHF scintillation observations can be utilized for the identification of the presence of ionospheric irregularities. The methodology for the calculation and comparison of LOFAR S₄ presented here forms the basis for automated and rapid monitoring of ionospheric irregularities, which can be applied to all LOFAR radio telescopes and which can augment traditional GNSS ionospheric observations. A disadvantage of LOFAR VHF scintillation observations is the need to ascertain whether the origin of scintillation is due to irregularities in the ionosphere or elsewhere. Given the propagation geometry, there also is the possibility that intensity fluctuations originating in the inner heliosphere could overlap with those originating in the ionosphere (Forte et al., 2022). However, a clear advantage in using LOFAR VHF scintillation for ionospheric studies is that (once the ionospheric origin is verified) these observations have a higher sensitivity to weaker electron density gradients than GNSS and the potential to detect ionospheric structures typically not detectable by only using traditional ionospheric GNSS measurements. This aspect allows us to take into account a wider variety of ionisation scales occurring in the ionosphere, which is essential for modelling purposes.

Acknowledgments

This paper is based on data obtained with the International LOFAR Telescope (ILT) under project codes LC7_001 and LC8_001 available through LOFAR Long Term Archive (LTA): https://lta.lofar.eu/Lofar. LOFAR (van Haarlem et al., 2013) is the Low-Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland.

UWM would like to thank the Ministry of Education and Science of Poland for granting funds for the Polish contribution to the International LOFAR Telescope (decision number 2021/WK/02) and for maintenance of the LOFAR PL-612 Baldy (MSHE decision no. 28/530020/SPUB/SP/2022). The UWM contribution is also supported by the National Centre for Research and Development, Poland, through grant ARTEMIS (decision numbers DWM/PL-CHN/97/2019 and WPC1/ARTEMIS/2019) and the National Science Centre, Poland, through grant 2017/27/B/ST10/02190 and the National Science Centre, Poland for granting “LOFAR observations of the solar corona during Parker Solar Probe perihelion passages” in the Beethoven Classic 3 funding initiative under project number 2018/31/G/ST9/01341.

The work carried out by BF at the University of Bath was supported by the UK Natural Environment Research Council [Grant number NE/R009082/1, Grant number NE/V002597/1, and Grant number NE/W003074/1]. The editor thanks four anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Figure 1 The positions of Polish LOFAR stations: PL612, PL611, PL610.
In the text

	Figure 2 LOFAR scintillation spectra at subsequent stages of elaboration: (a) the dynamic spectrum (i.e., raw intensity measurements), (b) RFI-free dynamic spectrum, (c) zero-mean normalized intensities and (d) the estimated S4 index.
In the text

	Figure 3 Power Spectral Density (PSD) estimated for the radio wave frequency 48.92 Mhz on zero-mean intensity fluctuations observed on DOY271 in 2017 during the time interval 16:45–18:30 UT. The PSDs were calculated over different time intervals: (a) 1 min, (b) 2 min, (c) 3 min, (d) 4 min and (e) 5 min.
In the text

	Figure 4 Fresnel scale variation with frequency.
In the text

Figure 5 LOFAR S4 of Cassiopeia A (a, b) and Cygnus A (c, d) calculated for days 87 (quiet) and 90 (disturbed) in 2017, recorded on LOFAR station PL612. GNSS ROT values were calculated for 1 s (e, f) and 60 s (g, h) for days 87 and 90 of 2017. GNSS S4 index (i, j) observed on L1 (blue dots), L2 (red dots) and L5 black dots) frequencies recorded with the ionospheric monitor co-located with the PL612; LOFAR scintillation indices for the 48.92 MHz channel (from Figs. 4a to 4d) are also shown as pink area (Cas A) and blue area (Cyg A). Gaps in the GNSS scintillation indices were due to issues related to data downloading.

In the text

	Figure 6 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated for days 87 and 90 of 2017, recorded on LOFAR station PL610. ROT values (e, f) calculated for 60 s for days 87 and 90 of 2017 from observations recorded by receiver near PL610 station.
In the text

	Figure 7 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated for days 87 and 90 of 2017, recorded on LOFAR station PL611.
In the text

Figure 8 LOFAR S4 of CasA (a, b) and CygA (c, d) calculated respectively for days 271 and 275 of 2017, recorded on LOFAR station PL612. ROT values calculated for 1 s (e, f) and 60 s (g, h) for days 271 and 275 of 2017. GNSS S4 index observed on L1, L2 and L5 frequencies recorded with the ionospheric monitor collocated with the PL612 LOFAR station superimposed on CasA and CygA S4 index calculated for the 48.92 MHz channel (i, j).

In the text

	Figure 9 LOFAR S4 of Cassiopeia A (a, b) and Cygnus A (c, d) calculated respectively for days 271 and 275 of 2017, recorded on LOFAR station PL611.
In the text

	Figure 10 ROTI meridional plot for DOY 87 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).
In the text

	Figure 11 ROTI meridional plot for DOY 90 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).
In the text

	Figure 12 ROTI meridional plot for DOY 275 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).
In the text

	Figure 13 ROTI meridional plot for DOY 271 for meridians: 17° E (a), 20.5° E (b) and 24° E (c).
In the text

	Figure 14 Power Spectra density plots for presented cases. Every PSD plot is made individually for every LOFAR station (PL610, PL611 and PL612) and for each target (CasA and CygA). All PSDs are calculated for middle channel.
In the text

	Figure 15 Cross-correlation functions (CCF) values for each baseline and both sources, DOY 090/2017.
In the text

	Figure 16 Cross-correlation functions (CCF) values for each baseline and both sources, DOY 271/2017.
In the text

	Figure 17 Proton density (upper panels), proton density ratio (middle panels) and solar wind speed (lower panels) obtained from ACE and Wind satellites for DOY 87 of 2017 (a), DOY 90 of 2017 (b), DOY 275 of 2017 (c) and DOY 271 of 2017 (d). The green bands indicate the exact times during which the LOFAR observations were collected.
In the text