© B. Boyde et al., Published by EDP Sciences 2022

1 Introduction

Radio waves emitted by astronomical sources and observed from the Earth are affected by any plasma through which they propagate. This means that observations of such radio signals contain information on the structure of the medium along the whole line of sight. As a result, these observations can be used to infer the properties of the interplanetary medium, the solar wind and the Earth’s ionosphere (e.g. Cordes et al., 2006; Brisken et al., 2009; Tokumaru et al., 2019; Fallows et al., 2020). While structures at all points along the line of sight can cumulatively contribute to the observed signal variation, the inferred properties of the scattering plasma can allow the approximate location of the structure to be determined (e.g. Fallows et al., 2016).

The effects of these plasma populations on radio waves arise due to the plasma density dependence of the refractive index, meaning that variations in plasma density across the line of sight cause distortions of the phase fronts. Depending on the size of the plasma irregularities, the effects can be refractive or diffractive in character, with diffractive effects associated with smaller irregularities. The characteristic length scale associated with diffractive scattering is called the Fresnel scale D_F (e.g. Basu et al., 1998), given by

$$D_{\mathrm{F}}=\sqrt{2\lambda L},$$(1)

where λ is the wavelength of the radio source and L is the distance from the observer to the scattering medium. The interference of the scattered waves with distance from the scattering region leads to rapid variations of received signal intensity as a function of time and space at a distance from the scattering medium. Due to the frequency dependence of the refractive index and Fresnel scale, the relative importance of refractive and diffractive effects changes when observing at different radio frequencies, and hence broadband observations allow much more detailed characterisation of the structures than is possible with a single isolated frequency.

The LOw-Frequency ARray (LOFAR; van Haarlem et al., 2013) is a low-frequency (10 MHz–250 MHz) radio telescope comprising a network of stations spread across Europe from Ireland to Latvia, with a dense “core” of stations in the northeast of The Netherlands. Although it is primarily intended for radio astronomy, at these frequencies, the ionosphere has a significant impact on the propagation of radio waves from astronomical sources (e.g. Mevius et al., 2016; de Gasperin et al., 2018). The wide bandwidth available to LOFAR means that the frequency dependence of the ionospheric effects can be determined, which is useful for characterising the structure responsible for the observed intensity variations.

Travelling ionospheric disturbances (TIDs) are horizontally propagating waves in the ionosphere, generally accepted to be manifestations of atmospheric gravity waves (AGWs) (Hines, 1960). These AGWs can be generated by various processes, including weather in the lower atmosphere, auroral precipitation and heating from ionospheric current systems (Hunsucker, 1982). Most studies of TIDs focus on either large-scale TIDs (LSTIDs) with wavelengths of order 1000 km or medium-scale TIDs (MSTIDs) with wavelengths of roughly 100–300 km (e.g. Hunsucker, 1982; Habarulema et al., 2013; Oinats et al., 2016; Koval et al., 2017; Terra et al., 2020; Themens et al., 2022). For example, Koval et al. (2017, 2018, 2019) showed a strong focus on solar radio emission due to the passage of MSTIDs. However, much less attention has been paid to small-scale TIDs (SSTIDs) (i.e. those with wavelengths <~100 km).

The scarcity of observations of SSTIDs is largely a result of the relative difficulty of observing such small-scale waves, as they are below the typical grid size of Global Navigation Satellite System (GNSS) Total Electron Content (TEC) maps (e.g. Otsuka et al., 2013, used 80 × 80 km to observe MSTIDs) and the typical time cadence of regular ionosonde measurements (e.g. Amorim et al., 2011; Fedorenko et al., 2011, both used 15 min time cadence to observe MSTIDs). Some studies have focused on small-scale ionospheric disturbances; for example, Baskaradas et al. (2014) observed many small-scale irregularities (time scales ~10 s) in the ionosphere using a single ionosonde operated at a high cadence with a fixed frequency, but they did not repeat periodically as would be expected for a TID. Small-scale TIDs were also reported by Ivanova et al. (2011) using oblique incidence sounding. These SSTIDs were observed in all cases for which LSTIDs were detected, as well as in some cases with no larger-scale disturbances. However, they did not specify any parameters of these observed TIDs. Ionospheric oscillations with periods <10 min were reported by Lan et al. (2018) using a single high-cadence ionosonde. These short-period variations were restricted to the lower altitudes <~200 km, which they ascribed to shorter period AGWs being filtered out at higher altitudes due to the increase in the Brunt-Väisälä frequency. They also confirmed that these shorter period variations were not detected in GNSS TEC maps. However, they were unable to conclusively determine if this was due to their short wavelength or to not reaching the altitude of the F2 peak and hence having relatively small TEC perturbations associated with them. However, none of these studies observed the SSTIDs from multiple locations and therefore could not determine velocities or wavelengths associated with them, only the period.

2 LOFAR observations

On the 15th of September 2018, between 10:15 and 11:48 UT, LOFAR was used to observe the radio sources Cygnus A (3C 405: RA 19h59m28s, Dec. 40.73°) and Cassiopeia A (3C 461: RA 23h23m24s, Dec. 58.82°) (observation ID L667596 under project code LT10_001, data available at https://lta.lofar.eu/). Each LOFAR station recorded independent intensity spectra of each source at a time resolution of 10.49 ms and a frequency resolution of 195.3 kHz between 24.99 MHz and 63.86 MHz. This provided 200 frequency channels for each source, and the observed spectra were combined into a dynamic spectrum of intensity as a function of frequency and time for all analyses.

In order to isolate the effects of the ionosphere, the dynamic spectra were filtered to remove radio frequency interference (RFI) and detrended to remove variation in intensity due to the antenna bandpass and changes in source elevation. The approach to this is described in Fallows et al. (2020) and involves flattening the data with a median filter with a window of (1.95 MHz × 0.5 s). RFI was then identified as any point which exceeded the median by more than 5 standard deviations after flattening. Any channels which contained more than 20% RFI were entirely removed, and isolated RFI points were replaced by linear interpolation. Each channel was then normalised by dividing by a fitted third-order polynomial, which removed two effects. Firstly it removed the difference in magnitude across channels due to the source spectrum and bandpass of the antenna, and secondly, it removed the temporal variation in intensity that arises from the difference in antenna sensitivity at different source elevations. For the purposes of plotting, the data were then down-sampled by a factor of 100 in time, giving a time resolution of approximately 1.05 s for the plots.

Figure 1 shows the dynamic spectrum of Cygnus A observed from the UK LOFAR station (51.1°N, 1.4°W, 176 m above sea level; note that lower frequencies are towards the top of the plot). Cygnus A was at a very low elevation, rising from 4.8° to 11.7° over the course of the observation and moving from 19.9° to 36.3° in azimuth clockwise from north. The features are pseudo-periodic, with one appearing roughly every 10 min over approximately an hour. These repeated intensity enhancements suggest focusing on a succession of electron density minima, such as would be associated with a TID and has been reported in solar radio observations by Koval et al. (2019). The features consistently display a broadening towards lower frequencies with clear diffraction fringes, as well as a time asymmetry with lower frequencies leading the higher frequencies. The broadening towards lower frequencies supports the idea that these are a result of focusing from an ionospheric lens, as the lower frequencies will reach their focus before the LOFAR antenna and therefore are increasingly spread out as observed by LOFAR.

Fig. 1 The filtered dynamic spectrum of Cygnus A from the LOFAR station in the UK (station UK608: 51.1°N, 1.4°W, 176 m above sea level). The frequency scale is inverted, with lower frequencies at the top of the plot. The normalised intensity is the intensity after the filtering process, with a normalised intensity of 1 being the expected value for an undisturbed ionosphere. Individual features are numbered to distinguish them in further discussion. The horizontal white lines are the channels removed due to RFI contamination.

Some individual features are shown in Figure 2, illustrating the fine structure that is present. Although the fringes at lower frequencies are consistent in all features, there are significant qualitative differences. The type of feature appears to alternate, as features 2 and 4 both show less dispersive intensity enhancements at higher frequencies, with a sharp transition to more dispersive behaviour below ~40 MHz, similar to a transition between weak and strong scattering behaviour (e.g. Rino et al., 1981), and the presence of diffraction fringes in the more dispersive frequency range. Features 1, 3 and 5 do not share this sharp transition and display a much more symmetrical form in the low frequencies. Features 1 and 5 also show a clear splitting into a “doublet” at frequencies above the apparent focal frequency, but it is unclear if this is also the case for feature 3 due to the proximity of the focal frequency to the upper-frequency limit of the observation.

Fig. 2 The dynamic spectra of the five features highlighted in Figure 1. Note that each panel uses a different colour scale for intensity to provide maximum contrast in each panel and that each has a different time scale.

Besides the lower frequencies leading the higher frequencies in time, another more subtle asymmetry in the features is the fringe spacing on the leading and trailing edges of the envelope in features 1, 3 and 5. In features 3 and 5, the fringes are more closely spaced in frequency along the trailing edge of the feature (later in time) than the leading edge (earlier in time), whereas for feature 1 the opposite is true. This is likely the result of some asymmetry in the ionospheric structures responsible for the features.

To the best of our knowledge, features such as these in dynamic spectra have not been reported elsewhere. The “spectral caustics” described by Koval et al. (2017) in solar radio observations have a similar overall form but do not display the fine structure of diffraction fringes present here. In their classification, feature 3 would be described as “inverted V like”, features 1 and 5 as “X like” and features 2 and 4 appear similar to their “fiber like” caustics. Many of their observations have been shown to correlate with the minima of passing MSTIDs Koval et al. (2019), and they were able to explain at least the “inverted V like” caustics as arising from focusing on the minima of such TIDs (Koval et al., 2018). The absence of fringes in their observations is likely due to the larger angular size of the Sun compared to Cygnus A.

The observations of Kuiack et al. (2021) of source magnification between 57.6 and 62.5 MHz using LOFAR may reflect the same phenomenon as observed here, although it is over a greatly reduced frequency range. They reported one case with a main peak of intensity with secondary peaks to either side, which they showed was described well by a first-order Bessel function of the first kind. Most of their observations do not display the secondary peaks, more in line with the observations of Koval et al. (2017), suggesting that the features shown in Figure 2 are a special case of a more general phenomenon. Kuiack et al. (2021) argue that these magnifications were most likely a result of focusing on electron density minima in the ionosphere, such as TIDs.

3 Modelling

3.1 Analytic modelling

Many of the features observed here using LOFAR are strongly similar to those produced by an analytic single-phase screen model in Meyer-Vernet (1980). This model replicates the envelope broadening towards lower frequencies and the pattern of overlapping fringes by assuming an infinite one-dimensional thin phase screen containing a single sinusoidal variation. With some minor simplifying approximations, Meyer-Vernet (1980) showed that an analytic solution could be found for the observed intensity, meaning that the parameter space can be quickly and easily sampled to compare the model to the observed dynamic spectrum. The parameter space that is sampled can be constrained by other observations, such as ionosondes and GNSS TEC measurements, as will be explored in Section 3.2.

The model uses co-ordinates x and z, where x is aligned with the phase screen and z along the line of sight (which is normal to the screen), with the source at z = −∞ and the screen at z = 0 as shown in Figure 3. The phase change Φ imparted by a given point x on the screen at time t is given by

$$\mathrm{\Phi }\left(x,t\right)=\frac{{\mathrm{\Phi }}_0}{f}\cos \left(\frac{x-vt}{d}\right)+\frac{\alpha x}{f}=\frac{{\mathrm{\Phi }}_0}{f}\cos \left(X\left(x,t\right)\right)+\frac{\alpha x}{f},$$(2)

where f is the frequency of the radio wave, and the other variables are defined in Table 1. The inverse frequency dependence of the phase screen amplitude assumes that the fluctuations in density are small and that the radio frequency remains well above the plasma frequency. This is necessary to ensure that the collisional and magnetic field terms are negligible and so the refractive index n can be approximated by

$$n\approx \sqrt{1-\frac{f_p^2}{f^2}}\approx 1-\frac{f_p^2}{2f^2}$$(3)

where f_p is the plasma frequency. The variation in the background ionosphere (i.e. the ionosphere before being perturbed by the TID) is assumed to be linear as this is the most variation that can be included while retaining the analytic solution. This corresponds to a linear variation in TEC with position along the propagation direction of the TID.

Fig. 3 A schematic diagram of the model, showing the co-ordinate system and several parameters. Dashed arrows represent raypaths for a given frequency.

Rather than directly inputting a value for z, it is more physically meaningful to define the altitude of the phase screen h and use the known elevation of the source θ. The value of z can then be calculated from these and the radius of the Earth R_E as

$$z=-R_E\sin \left(\theta \right)+\sqrt{\left(R_E\sin \right(\theta ){)}^2+(R_E+h{)}^2-R_E^2},$$(4)

assuming that the Earth is a perfect sphere and the observer is at a negligible height above sea level. The source is assumed to have a Gaussian brightness distribution (intensity as a function of angle) of

$$B\left({\theta }_k\right)\propto \exp \left(\frac{{\theta }_k^2}{\Delta {\theta }^2}\right),$$(5)

where θ_k is the angle between the radio wavevector and the z-axis, and for Cygnus A, the angular size Δθ is assumed to be 0.5 arcmin (Carilli et al., 1991; Skrutskie et al., 2006). While the actual structure of Cygnus A is not Gaussian, actually consisting of two distinct lobes, this simple model is adequate to explain the observed features.

The normalised intensity I that would be observed at the point (x, z) and time t is calculated in the Fresnel approximation. By assuming that the dominant contributions come from the regions of the screen close to the line of sight, the analytic expression for I can be shown to be (Meyer-Vernet, 1980)

$$I=1+2\sum _{p=1}^{\mathrm{\infty }}\mathrm{}{J}_{\mathrm{p}}\left(\frac{2{\mathrm{\Phi }}_0}{f}\sin \left(\frac{\mathrm{pZ}}{2}\right)\right)\cos \left(\mathrm{pX}\left(x,t^{\mathrm{\prime }}\right)\right)\exp \left(-(\mathrm{pu}{)}^2\right),$$(6)

where J_p is a Bessel function of the first kind, $Z=\frac{\mathrm{zc}}{2\pi f{d}^2}$ is the normalised distance to the screen, $t^{\text{'}}=t+\frac{\mathrm{zc\alpha }}{2\pi f^{2}v}$ is the time adjusted for the refraction in the background ionosphere (Meyer-Vernet et al., 1981), and $u=\frac{\Delta \theta z}{2d}$ is the normalised angular size of the source. We note that the value of x chosen for the observer does nothing except shift the observed feature in time, and so we can assume x = 0 without loss of generality. The infinite sum is, in practice, possible to compute due to the rapid decay of both the Bessel function term and the exponential for high p.

The two terms that both decays with high p provide the explanation for fringes appearing in certain cases and being absent in others. If the sum is truncated by the decay of the exponential term (i.e. by the finite size of the source), then the fringes are smoothed out. If, on the other hand, the Bessel function is primarily responsible for the truncation, then the fringes are observed. The characteristic p values at which the terms start to decay are $p\approx \frac{1}{u}$ and $p\approx \frac{\mathrm{e}{\mathrm{\Phi }}_0}{f}$ for the exponential and Bessel function terms, respectively, where e is Euler’s number (Meyer-Vernet, 1980). This means that the approximate condition for the presence of fringes can be expressed as

$$\frac{\mathrm{e}{\mathrm{\Phi }}_0}{f}<\frac{2d}{z\Delta \theta },$$(7)

meaning that for a given phase screen the presence of fringes requires the source size to be below a certain threshold.

In practice, we do not have a given phase screen. Instead, the source size will be taken from previous observations of the source, and the focal frequency is the simplest physical observable. Considering the phase screen as acting like an idealised thin parabolic lens, we find using geometric optics that $\frac{{\mathrm{\Phi }}_{0}z}{d^2}=C$ where C is some constant determined by the focal frequency. This represents the relationship between lens curvature $\frac{{\mathrm{\Phi }}_0}{d^2}$, which is equivalent to focusing power, and focal length at the given frequency z; as focusing power increases, focal length decreases and vice versa. Combined with equation (7), this means that for a given observation, the presence of fringes places an upper limit on the wavelength of the TID responsible.

In order to directly compare the model to observations, it is important to account for the background intensity recorded by LOFAR. This is a combination of diffuse background emission and instrumental effects, and in the LBA corresponds to a total system effective flux density (SEFD) of ~30 kJy (van Haarlem et al., 2013), which is unaffected by ionospheric irregularities. The flux density of Cygnus A at the frequencies considered here is ~20 kJy (de Gasperin et al., 2020). However, the received flux density from Cygnus A will be elevation dependent due to the difference in the projected area of the antenna array, whereas the SEFD will be approximately elevation independent as the increase in the field of view at low elevations counterbalances the reduced sensitivity to background emission. As a result, the observed intensity is given by

$$I=1+2\frac{I_{\mathrm{CygA}}\sin \left(\theta \right)}{I_{\mathrm{CygA}}\sin \left(\theta \right)+I_{\mathrm{sys}}}\sum _{p=1}^{\mathrm{\infty }}\mathrm{}{\mathrm{J}}_{\mathrm{p}}\left(\frac{2{\mathrm{\Phi }}_0}{f}\sin \left(\frac{\mathrm{pZ}}{2}\right)\right)\cos \left(\mathrm{pX}\left(x,t^{\mathrm{\prime }}\right)\right)\exp \left(-(\mathrm{pu}{)}^2\right),$$(8)

where I_CygA is the flux density of Cygnus A, I_sys is the SEFD and θ is the elevation.

3.2 Applying the model

In order to compare the observed features to this model, some parameters were fixed based on complementary observations discussed below, and others were adjusted to match the observed intensity distribution. The screen altitude h, velocity v and source size Δθ were fixed to 200 km, 50 ms⁻¹ and 0.5 arcmin, respectively, using these complementary observations. For the remaining parameters: d and Φ₀ were adjusted to match the focal frequency and fringe spacing by eye (lower values of d lead to wider fringe spacing), and α was set to match the time asymmetry that was observed. In order to give a quantitative comparison between the model and observations, several characteristic values were defined: the focal frequency, the frequency of the 5th fringe along each edge of the envelope (leading and trailing) and the time between the 5th fringe on the leading and trailing edges. The fringe location was defined as the maximum intensity within the fringe, and the values for each feature and their modelled replicas are given in Table 2.

The strength of this model in replicating the types of features observed by LOFAR is clearly illustrated in Figure 4. The overall envelope shape, the pattern of overlapping fringes and time asymmetry are all present, and the required wavelength suggests that the features are caused by an SSTID (Λ = 20 km). The modelled intensity variations also match very closely in magnitude to the observation. If the background is not accounted for, the modelled focal intensity is ~15 times the mean, demonstrating the importance of correcting for this effect. The periodicity of the model, in this case, is shorter than the observations, meaning that multiple features are shown, but the individual modelled features match well. The width of the modelled feature given in Table 2 is significantly lower than the observed feature, but this is due to the edges of the observed feature being approximately straight lines in the time-frequency plane, whereas the modelled feature curves outwards significantly, meaning that it is not possible to match the width across a range of frequencies. There is a small discrepancy in fringe spacing between observation and model, as given in Table 2, indicating that the parameters used in this case are not optimal, but nevertheless it illustrates that features of this type can be represented by this model.

Fig. 4 Feature 3 along with its modelled replica. The time axes cover the same duration, with t = 0 in the model corresponding to the minimum electron density (i.e. the maximum of phase change) lying on the line of sight. The parameters used here are: Φ0 = 8 × 108 rad Hz, Λ = 20 km, θ = 7° and α = 1.89× 105 rad Hz m−1.

Figure 5 shows several data sources that were used to constrain various parameters of the model. The ionosonde data were used to constrain h, as TID amplitudes typically maximise around the altitude of maximum electron density (Fedorenko et al., 2011). This gave an approximate value of h = (199 ± 7) km based on the mean and variance of the manually scaled hmF2 values between 10:15 and 11:45 UT (the duration of the LOFAR observation), and so a value of h = 200 km was used for all modelling. The TEC data from MADRIGAL were used to attempt to constrain α, as it provides, in principle, a maximum value if the TID propagation aligns perfectly with the TEC gradient. The gradients were calculated from the difference in TEC between two cells separated by 4° of latitude, selected to minimise the noise while ensuring the values are still sufficiently local (latitudinal gradients are neglected as they are approximately zero). All gradients calculated between 10:15 and 11:45 UT centred on latitudes between 58° and 61°N and longitudes between 5° and 7°E are considered, as this was roughly the region of the ionosphere corresponding to the observation of Cygnus A from UK608 assuming a thin shell ionosphere at an altitude of 200 km. The distribution of these gradients is shown in panel 2 of Figure 5 and could then be mapped to α by converting from vertical TEC to the slant path observed by LOFAR (assuming an elevation of 7° corresponding to the third feature and a thick shell ionosphere between 150 and 500 km). The conversion from TEC to phase change Δϕ normalised to 1 Hz is given by

$$\Delta \varphi =\frac{-\mathrm{TEC}{e}^2}{4\pi c{m}_e{ϵ}_0},$$(9)

where e is the charge of an electron, c is the speed of light, m_e is the mass of an electron, and ϵ₀ is the permittivity of free space. As with equation (2), this assumes that the fluctuations in electron density are small and that the plasma frequency remains well below the radio frequency. The mean and variance of the TEC gradients suggested that α = (1.20 ± 0.39) × 10⁴ rad Hz m⁻¹ assuming a northward velocity.

Fig. 5 The additional data used to constrain the model parameters. The leftmost panel shows the altitude of peak ionospheric electron density (hmF2) measured from the Juliusruh ionosonde (54.6°N, 13.6°E). The middle panel shows the distribution of north-south vTEC (vertical TEC) gradients (positive implies increasing to the north) calculated from MADRIGAL GNSS TEC based on the difference over 4° of latitude centred on points between 58°–61°N, 5°–7°E, and 10:15–11:45UT. The third panel shows vTEC measured from the Southampton (SOTN) GNSS receiver (50.9°N, 1.5°W) for all azimuths between 0° and 60° clockwise from the north.

The GNSS TEC data from Southampton were used to constrain the amplitude of variations in the phase screen, as the TEC perturbation amplitude ΔTEC along the line of sight can be related to Φ₀ using equation (9). By inspection of Figure 5, the waves in the TEC data have ΔTEC ≈ 0.05 TECu, but the corresponding TEC observed along the LOFAR line of sight cannot be precisely calculated as it is highly sensitive to the relative orientation of the line of sight and TID phase front, especially for short wavelength TIDs. However, as an order of magnitude estimate, this is still useful, suggesting Φ₀ = 5 × 10⁸ rad Hz.

As well as the data shown in Figure 5, the LOFAR observation itself can be used to constrain the velocity v of a structure if it is observed by several stations. The dense network of stations in The Netherlands provided observations of a single feature from all 24 stations in the LOFAR core, with an example shown in Figure 6 (below). Following the approach described in Fallows et al. (2020), this allowed the propagation velocity of the structure to be estimated by considering the cross-correlations between stations to estimate the propagation time along each baseline (remote stations were not included as the feature evolved in frequency over these larger scales). Combining this gave an estimate for the velocity of the structure relative to the motion of the LOFAR lines of sight of (57 ± 5) ms⁻¹ at an azimuth of 104° clockwise from north. For all modelling, the velocity was assumed to be v = 50 ms⁻¹ as a rough estimate, although the velocity of the structures observed from the UK station is unlikely to be exactly the same as those observed from The Netherlands. However, an error in the assumed velocity will only rescale the time axis and will not distort the overall intensity distribution otherwise.

Fig. 6 The feature observed in the dynamic spectrum of Cygnus A by one of the stations in The Netherlands (station CS001: 52.9°N, 6.9°E), alongside the modelled replica (note that different colour scales are used in each panel). As in Figure 4, the time scale of both panels is the same, with t = 0 corresponding to the minimum of electron density (i.e. the maximum of phase change) lying on the line of sight. The parameters used here are: Φ0 = 8 × 108 rad Hz, Λ = 25 km, θ = 8.8° and α = 7.58 × 104 rad Hz m−1.

Correcting the apparent source motion from each station allowed us to estimate the physical velocity of the ionospheric structure. Assuming an altitude of 200 km suggests a velocity of (75 ± 5) ms⁻¹ relative to the Earth. Knowing the direction of propagation allowed us to estimate the resulting value of α projected onto this direction, which is the value significant to the model. In this case, this suggested a small negative value of α = (−2.9 ± 1.0) × 10³ rad Hz m⁻¹ (neglecting any uncertainty contribution from the velocity estimate).

The amplitude of the phase perturbation used in Figure 4 is consistent with that implied by the GNSS TEC measurements shown in Figure 5, but the value of the gradient term α is greater in magnitude and positive despite the expected small negative value. This suggests that although the chosen value of α can approximately replicate the observed time asymmetries, the actual physical mechanism for the asymmetry is not a large-scale horizontal electron density gradient. One candidate to explain the asymmetry besides the horizontal TEC gradient is the low elevation of the source (7°), meaning that the incident radio waves are significantly refracted as they enter the ionosphere by the vertical gradients in electron density. Because this refraction is frequency-dependent, the different frequencies pass through the peak of ionospheric density at different locations and so “see” the TID passing at different times, and would cause the lower frequencies to lead to higher frequencies if the structures propagate towards the observer. However, assuming that the propagation direction estimated from the stations in The Netherlands is consistent with the features observed from the UK, the waves are propagating slightly away from the observer as viewed from the UK station. Another possible factor is the alignment of the line of sight with the phase fronts of the TID, which was shown by Koval et al. (2018) to create significant asymmetries in their modelled dynamic spectra.

As well as the features observed from the UK station shown in Figure 2, one similar feature was observed in the spectrum of Cygnus A by many of the LOFAR stations in The Netherlands (which was used to estimate v), with an example shown in Figure 6. This feature is the one which most closely matches the model in its form, with no doublet above the focal frequency and its asymmetry very closely replicated by the effect of the chosen value of α. However, this value of α is again inconsistent with the observed TEC gradients and velocity, which imply a negative α, so one or more other physical mechanisms must be responsible. As in Figure 4, the feature can only be replicated in the model by an SSTID (Λ = 25 km). The absence of a doublet in both the model and this observation suggests that the doublet in features 1 and 5 may arise from some perturbation to the simple sine wave; this will be discussed in the following section. Unlike in Figure 4, the modelled intensity is noticeably higher than is observed, suggesting that the model may be overly idealised in many cases. Also, despite the observed feature showing distinct curvature in the edges of its envelope, the modelled replica is still significantly narrower, as shown in Table 2. This could be corrected by increasing the TID length scale d, but this would cause the fringe spacing to decrease, meaning that the model cannot precisely match the observation in this case (unless v is an overestimate).

3.3 Numerical modelling

The analytic model derived by Meyer-Vernet (1980) provides a very clear replication of several of the observed features, but the most prominent aspect that is absent is the doublet above the focal frequency shown in Figure 2. This intensity distribution cannot be replicated by the single sine wave model, clearly suggesting that there is scope to extend the model. In order to investigate this, it was necessary to move from the analytic model to a numerical phase screen approach.

The numerical phase screen approach was outlined in detail in Sokolovskiy (2001) and allows for an arbitrary number of parallel screens to be considered. However, for the purposes of this work, we continued to consider a single phase screen so as to minimise the number of parameters required to define the model. While multiple phase screens may provide a more detailed representation of the ionosphere, it will be shown that a single phase screen is nevertheless capable of replicating the features observed by LOFAR in this case.

The use of numerical phase screen models in the ionosphere has previously focused on anthropogenic radio sources, such as ground- or space-based radar systems and GNSS signals (e.g. Hocke & Igarashi, 2003; Wang et al., 2014; Ludwig-Barbosa et al., 2019; Carrano et al., 2020; Ding et al., 2021), meaning that care must be taken in defining the input spectrum of a natural source as considered here. As the source is incoherent, each angular component must be treated separately to avoid the spurious phase relationships between the components that occur if they are combined into a single spectrum. The intensity contributions from each component are then summed to give the total observed intensity distribution. As in the analytic model, the source is assumed to be Gaussian for simplicity, although any intensity distribution is possible in the numerical model. The phase screen is sampled with 2 m resolution to ensure that the full angular spectrum is available throughout the observing band of LOFAR (the upper-frequency limit of ~65 MHz used in this observation corresponds to a wavelength of ~4.70 m). The screen contains 2¹⁸ samples, giving an extent of 524.288 km. These values are chosen to ensure that the wavelengths of the TIDs considered are far lower than the size of the screen and to provide a high angular resolution when the Fourier transform is applied. The numerical phase screen calculation gives intensity as a function of x, which is then converted to a function of t by assuming a constant velocity v.

The similarity of the low-frequency behaviour of features 1 and 5 to the analytic model suggests that the simple sinusoid is a good starting point. The doublet at high frequencies is indicative of focusing from two discrete points, corresponding to two minima of phase curvature, rather than the single minimum per period present with a single sine wave. The simplest way to modify the simple sinusoid to achieve this is to add a second harmonic of lower amplitude to equation (2), meaning that the screen’s effect becomes

$$\mathrm{\Phi }\left(x,t\right)=\frac{{\mathrm{\Phi }}_0}{f}\left[\cos \left(X\left(x,t\right)\right)+A\cos \left(2X\left(x,t\right)+\mathrm{\varphi }\right)\right]+\frac{\alpha x}{f},$$(10)

where A and ϕ are the relative amplitude and phase offset of the harmonic, respectively. As the observed features are approximate time-symmetric (except for the distortion represented by α), the phase screen should retain its symmetry, corresponding to ϕ = 0. This means that this modification to the phase screen introduces a single additional free parameter to the model, A. Although the comparison to the observed features in Section 3.2 suggests that the physical cause of the asymmetry was not the horizontal TEC gradients represented by α, it was retained as a simple and effective term to replicate the observations.

In order for this phase screen to contain two minima of phase curvature, the constraint on the value of A is A < −0.0625. A further constraint is suggested by the fact that the intensity between the two peaks of the doublet is not suppressed (i.e. is ≥1), indicating that there is no de-focusing in this region. This lack of de-focusing indicates that the curvature is not positive at any point between its minima, corresponding to the constraint A > −0.25. An example of a perturbation satisfying these constraints is shown in Figure 7 alongside the unperturbed phase variation.

Fig. 7 The phase perturbations associated with the simple sine wave perturbation given by equation (2) (red dashed line) and the perturbation given by equation (10) (black solid line). The minima of curvature for both curves are also marked. Both are normalised to an amplitude of 1 and wavelength of 1 and have α = 0. The harmonic has a relative amplitude of A = −0.12.

The presence of significant harmonics in TID waves has been reported elsewhere, for example, van de Kamp et al. (2014) identified a TID using incoherent scatter radar and GNSS TEC maps, with the observed periods differing by a factor of two. They explained this as the two methods identifying different harmonics of the same underlying TID structure. Harmonic generation in the ionosphere has also been reported in numerical modelling work by Kirchengast (1997), where it occurred even for a monochromatic driver due to non-linear processes such as frictional heating. While these works both focused on TIDs with periods on the order of an hour, theoretical work has shown harmonic generation, from both spectra of AGWs and non-linear ionospheric responses to monochromatic AGWs, without any assumptions about the period of the waves considered (Chao-Song & Jun, 1991; Nekrasov et al., 1995; Nekrasov & Shalimov, 2002).

The addition of this harmonic to the modelled phase screen is effective in replicating the doublet, as shown in Figure 8. The qualitative agreement between the model and observation remains strong, but as in Figure 6 the intensity is overestimated (the observed focal intensity is ~1.75, model is ~2.8). Another discrepancy visible in Figure 8 is that the model predicts slightly higher intensities for the fringes along the centre of the envelope than on the boundaries, which is not apparent in the observation. The width of the observed feature matches very closely with the model in this case, as shown in Table 2.

Fig. 8 Feature 1 and its modelled replica. Both plots use the same time scale and t = 0 corresponds to the TID minimum lying on the line of sight (note that different colour scales are used in each panel). The parameters used here are: Φ0 = 1.1 × 109 rad Hz, Λ = 30 km, θ = 5°, α = 1.45× 105 rad Hz m−1 and A = −0.12.

The addition of the harmonic is also able to replicate feature 5, as shown in Figure 9. The TID amplitude and wavelength are both lower than those used to replicate feature 1 (wavelength reduced by factor 2, amplitude by factor 3), which is required to increase the fringe spacing while retaining the approximately same focal frequency. In this case, the qualitative agreement is not as strong as for feature 1 (Fig. 8), as the modelled feature spreads out far more towards the low frequencies than towards the high frequencies, which is not the case for the observed feature. The overestimate in intensity for this feature is much smaller than in Figures 6 and 8. While there is clearly notable variation in the parameters required to replicate the different features of the observation, all are consistent with an SSTID.

Fig. 9 Feature 5 and its modelled replica. Again the time scale of the model and observation are the same, but the intensity scales used are different. The parameters used here are: Φ0 = 3.5 × 108 rad Hz, Λ = 15 km, θ = 8.4°, α = 7.37× 105 rad Hz m−1 and A = −0.12.

4 Discussion

The modelling described above is able to represent many of the features observed by LOFAR in this observation, despite being an extremely simplified representation of the actual ionosphere. Whereas the actual ionosphere varies in three spatial dimensions and in time, the phase screen model of Meyer-Vernet (1980) is able to replicate these observations by considering only variation in a single spatial dimension and restricting time variation to propagation with a single constant velocity. The phase screen model allows approximate TID parameters such as wavelength and TEC perturbation amplitude to be estimated and indicates that LOFAR is picking up signatures of TIDs too small to be resolved by the most common techniques of ionospheric monitoring.

The extension to numerical modelling allows us to build upon the initial model of Meyer-Vernet (1980) and suggests that LOFAR can infer the presence of non-linearities in the TID perturbation. This means that as well as identifying TIDs with shorter wavelengths than those observed by other techniques, these observations can provide information on the fine structure of such TIDs. The numerical modelling approach also enables non-wave-like perturbations to be considered, meaning that other seemingly unrelated features in dynamic spectra can potentially be replicated with this approach. The main limitation of this is that the model, as currently defined, could not replicate a situation where multiple structures are present and propagate in different directions or with different velocities, as observed by Fallows et al. (2020).

The observations from the stations in The Netherlands allow us to confirm that these small-scale disturbances are propagating horizontally, supporting the arguments of both Ivanova et al. (2011) and Lan et al. (2018) that their observations of short-period ionospheric variations at a single point were manifestations of SSTIDs. Although only one feature was observed from The Netherlands, and hence no period is directly observable in this case, the phase screen modelling and observed period of the similar features from the UK both strongly suggest that this feature is a manifestation of an SSTID.

The introduction of non-linearities to the modelling also builds upon the explanation provided by Koval et al. (2018) for their “X-shaped” caustics. They observed some intensity enhancement above the focal frequency when modelling focusing from a perfectly sinusoidal TID and proposed that this focusing could therefore explain both their “inverted V like” and “X like” caustics. However, their model did not replicate the two sharp intensity enhancements but rather a more diffuse intensity enhancement above the focal frequency, and also did not extend as far above the focal frequency as the observed “X like” features do. Explaining the “X like” caustics as arising from non-linearities, such as the harmonic described in Section 3.3, also provides a clear explanation for why they occur in some cases and not others. However, there are likely other non-linear perturbations that would lead to similar features in dynamic spectra.

The effectiveness of the phase screen model, despite clearly being an oversimplification, suggests that the TID perturbation may be highly localised in altitude. This would be consistent with the observations of Lan et al. (2018) that short-period perturbations are restricted to lower altitudes, up to around the observed F2 layer peak in this case, and if the AGW driving the observed TID can only just reach the altitude of the F2 layer peak, the resulting plasma density perturbations would be expected to be extremely localised in altitude (i.e. not extending significantly into the topside). The low altitude of the F2 layer peak (~200 km) may also explain why these features were observed in this case but are not seen more regularly. If the F2 peak was at a higher altitude, then any TID-driven variations in electron density and hence TEC would be dominated by AGWs capable of reaching this higher altitude, thereby reducing the significance of the shorter wavelength waves observed here. This low altitude also increases the importance of observing at low elevation in order to provide a sufficient path length for focusing to occur well within the LOFAR observing band.

However, there are still significant discrepancies between the observed and modelled dynamic spectra, most notably the inadequacy of the horizontal TEC gradient to explain the asymmetry. Although the parameter α does not provide a physically reasonable explanation for the observed time asymmetry, it is possible to closely replicate the observations with an arbitrary value. This suggests that although the physical mechanism is not the large-scale horizontal gradients in the ionosphere identified from MADRIGAL GNSS TEC data, whatever mechanism is responsible is likely to cause similar frequency dependence. The asymmetries arising from the low elevation and the alignment between the line of sight and TID phase fronts demonstrated by Koval et al. (2018) are both candidates. However, this was not explored here as neither effect can be quantified in a way that allows it to be directly added to the simple single-phase screen model.

As well as the discrepancies between observed and modelled spectra, there are two features in the UK observation of Cygnus A that the model does not replicate. Features 2 and 4 are distinct from the others largely due to the nature of their time asymmetry. Whereas the others remain approximately symmetric around a curve described by some value of α, these features show little dispersion in the higher frequencies and strong dispersion in the lower frequencies. This suggests that any perturbation capable of explaining these features must either be spatially asymmetric or change significantly in time as it passes through the line of sight. Therefore, in order to attempt to replicate these features, the model perturbation would need to be extended beyond what has been considered here to allow for strongly asymmetric perturbations.

5 Conclusions

We have identified features in dynamic spectra observed by LOFAR on the 15th of September 2018 that indicate the passage of an SSTID across the line of sight. They are characterised by a broadening towards the lower frequencies overlaid with diffraction fringes and a focal frequency typically around 40–60 MHz. Phase screen modelling can replicate the observed intensity distribution in time and frequency, allowing the wavelength and amplitude of the SSTID to be estimated given its altitude. The structures that explain the observed features are too small to show up in GNSS TEC maps, but TEC time series from individual GNSS receivers show waves with the expected amplitude and period present at the time. Favourable conditions for observing the features require a compact source to resolve the fringes and a low elevation to allow for a long focal length.

The simplicity of the model required to replicate the observations is striking because it implies that the ionospheric variations were dominated by a single structure with a well-defined altitude and wavelength. There is no evidence of scattering from a spectrum of irregularities, as has been reported by Fallows et al. (2020), suggesting that either condition was not favourable for the significant growth of irregularities or that such irregularities had not had time to develop. Work is in progress to determine whether this single-phase screen approach can be applied to other features observed in LOFAR dynamic spectra.

These observations show that LOFAR can resolve ionospheric structure on a scale other approaches tend to miss. GNSS TEC maps and ionosondes are typically too coarse, and GNSS scintillations arise from much smaller structures (Basu et al., 1998). This suggests LOFAR observations may be able to complement others to provide a complete description of TIDs, both by covering TID wavelengths that are otherwise not commonly observed and by revealing the presence of certain non-linearities which may not be apparent from other observing methods. However, a more systematic study would be required to determine whether these features are a regular occurrence or if the observation considered here happened under specific favourable conditions that are rarely repeated.

Acknowledgments

This work is supported by the Leverhulme Trust under Research Project Grant RPG-2020-140. Ben Boyde acknowledges receipt of a PhD studentship from the same grant. This paper is based on data obtained with the International LOFAR Telescope (ILT) under project code LT10_001. LOFAR (van Haarlem et al., 2013) is the LOFAR 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 benefited 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. LOFAR data are available at https://lta.lofar.eu/. This paper uses data from the Juliusruh Ionosonde, which is owned by the Leibniz Institute of Atmospheric Physics, Kuehlungsborn. The responsible Operations Manager is Jens Mielich. The ionosonde data were accessed through the Global Ionospheric Radio Observatory (GIRO; Reinisch & Galkin, 2011), accessible at http://spase.info/SMWG/Observatory/GIRO. The services of the Natural Environment Research Council (NERC) British Isles continuous GNSS Facility (BIGF), https://www.bigf.ac.uk, in providing archived GNSS data (and/or products) to this study are gratefully acknowledged. GPS TEC data products and access through the Madrigal distributed data system are provided to the community by the Massachusetts Institute of Technology under support from US National Science Foundation grant AGS-1952737. Data for the TEC processing is provided from the following organisations: UNAVCO, Scripps Orbit and Permanent Array Center, Institut Geographique National, France, International GNSS Service, The Crustal Dynamics Data Information System (CDDIS), National Geodetic Survey, Instituto Brasileiro de Geografia e Estatística, RAMSAC CORS of Instituto Geográfico Nacional de la República Argentina, Arecibo Observatory, Low-Latitude Ionospheric Sensor Network (LISN), Topcon Positioning Systems, Inc., Canadian High Arctic Ionospheric Network, Institute of Geology and Geophysics, Chinese Academy of Sciences, China Meteorology Administration, Centro di Ricerche Sismologiche, Système d’Observation du Niveau des Eaux Littorales (SONEL), RENAG: REseau NAtional GPS permanent, GeoNet – the official source of geological hazard information for New Zealand, GNSS Reference Networks, Finnish Meteorological Institute, SWEPOS – Sweden, Hartebeesthoek Radio Astronomy Observatory, TrigNet Web Application, South Africa, Australian Space Weather Services, RETE INTEGRATA NAZIONALE GPS, Estonian Land Board, Virginia Tech Center for Space Science and Engineering Research, and Korea Astronomy and Space Science Institute. B. Dabrowski and A. Krankowski thank 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 UWM authors also thank the Ministry of Education and Science (MES), Poland for granting funds for the Polish contribution to the International LOFAR Telescope (agreement no. 2021/WK/02). The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 The filtered dynamic spectrum of Cygnus A from the LOFAR station in the UK (station UK608: 51.1°N, 1.4°W, 176 m above sea level). The frequency scale is inverted, with lower frequencies at the top of the plot. The normalised intensity is the intensity after the filtering process, with a normalised intensity of 1 being the expected value for an undisturbed ionosphere. Individual features are numbered to distinguish them in further discussion. The horizontal white lines are the channels removed due to RFI contamination.

In the text

	Fig. 2 The dynamic spectra of the five features highlighted in Figure 1. Note that each panel uses a different colour scale for intensity to provide maximum contrast in each panel and that each has a different time scale.
In the text

	Fig. 3 A schematic diagram of the model, showing the co-ordinate system and several parameters. Dashed arrows represent raypaths for a given frequency.
In the text

	Fig. 4 Feature 3 along with its modelled replica. The time axes cover the same duration, with t = 0 in the model corresponding to the minimum electron density (i.e. the maximum of phase change) lying on the line of sight. The parameters used here are: Φ0 = 8 × 108 rad Hz, Λ = 20 km, θ = 7° and α = 1.89× 105 rad Hz m−1.
In the text

Fig. 5 The additional data used to constrain the model parameters. The leftmost panel shows the altitude of peak ionospheric electron density (hmF2) measured from the Juliusruh ionosonde (54.6°N, 13.6°E). The middle panel shows the distribution of north-south vTEC (vertical TEC) gradients (positive implies increasing to the north) calculated from MADRIGAL GNSS TEC based on the difference over 4° of latitude centred on points between 58°–61°N, 5°–7°E, and 10:15–11:45UT. The third panel shows vTEC measured from the Southampton (SOTN) GNSS receiver (50.9°N, 1.5°W) for all azimuths between 0° and 60° clockwise from the north.

In the text

Fig. 6 The feature observed in the dynamic spectrum of Cygnus A by one of the stations in The Netherlands (station CS001: 52.9°N, 6.9°E), alongside the modelled replica (note that different colour scales are used in each panel). As in Figure 4, the time scale of both panels is the same, with t = 0 corresponding to the minimum of electron density (i.e. the maximum of phase change) lying on the line of sight. The parameters used here are: Φ0 = 8 × 108 rad Hz, Λ = 25 km, θ = 8.8° and α = 7.58 × 104 rad Hz m−1.

In the text

	Fig. 7 The phase perturbations associated with the simple sine wave perturbation given by equation (2) (red dashed line) and the perturbation given by equation (10) (black solid line). The minima of curvature for both curves are also marked. Both are normalised to an amplitude of 1 and wavelength of 1 and have α = 0. The harmonic has a relative amplitude of A = −0.12.
In the text

	Fig. 8 Feature 1 and its modelled replica. Both plots use the same time scale and t = 0 corresponds to the TID minimum lying on the line of sight (note that different colour scales are used in each panel). The parameters used here are: Φ0 = 1.1 × 109 rad Hz, Λ = 30 km, θ = 5°, α = 1.45× 105 rad Hz m−1 and A = −0.12.
In the text

	Fig. 9 Feature 5 and its modelled replica. Again the time scale of the model and observation are the same, but the intensity scales used are different. The parameters used here are: Φ0 = 3.5 × 108 rad Hz, Λ = 15 km, θ = 8.4°, α = 7.37× 105 rad Hz m−1 and A = −0.12.
In the text