© ESO 2018

1 Introduction

The ultra-low frequencies (10−100 MHz) are the last poorly explored window available for ground-based astronomical observations. Some attempts have been made in the past to cover this frequency range; notably, the 38 MHz 8th Cambridge survey (8C; Rees 1990) and the 74 MHz Very Large Array Low-frequency Sky Survey (VLSS; Cohen et al. 2007; Lane et al. 2014) pioneered this exploration. More recently, the GaLactic and Extragalactic All-sky MWA Survey (GLEAM; Hurley-Walker et al. 2017) produced images down to 72 MHz. A major limitation when observing at these frequencies is the presence of the ionosphere, a layer of partially ionised plasma, surrounding our planet.

The ionisation of the ionosphere is driven by the UV and X-ray radiation generated by the Sun during the day and is balanced by recombination at night. A lower level of ionisation is maintained during the night by the action of cosmic rays. The peak of the free electron density lies at a height of ~300 km but the ionosphere extends, approximately, from 75 to 1000 km. The free electron column density along a line of sight (LoS) through the ionosphere is generally referred to as the total electron content (TEC). The TEC unit (TECU) is 10¹⁶ m⁻², which is the order of magnitude typically observed at zenith during the night; in normal conditions the TEC during daytime hours is ten times higher. When observing radio emission at long wavelengths, the ionosphere introduces systematic effects such as reflection, refraction, and propagation delay of the radio waves (Mangum & Wallace 2015). For interferometric observation, propagation delay is the main concern (Intema et al. 2009). The effect is caused by a varying refractive index n of the ionospheric plasma along the wave trajectories. The total propagation delay, integrated along the LoS at frequency ν, results in a phase rotation given by $${\Phi }_{\text{ion}}=-\frac{2\pi \nu }{c}{\displaystyle {\int }_{\text{LoS}}\left(n-1\right)}\text{d}l.$$(1)

An n constant in time and space would impose a coherent phase error that would result in a spatial shift of the observed image compared to the true sky. The problem becomes more complicated as n depends strongly on time and position.

Neglecting the frictional force and assuming a cold, collisionless, magnetised plasma (such as the ionosphere), the refractive index n can be calculated exactly (Davies 1990). For signals with frequencies ν ≫ ν_p (the plasma frequency, that for the ionosphere is around 1−10 MHz), it can be expanded (see e.g. Datta-Barua et al. 2008) into a third-order Taylor approximation retaining only terms up to ν⁻⁴ : $$\begin{array}{lll}\hfill & n\approx \hfill & 1-\frac{q^2}{8{\pi }^2{m}_e{ϵ}_0}\cdot \frac{n_e}{{\nu }^2}\pm \frac{q^3}{16{\pi }^3{m}_e^2{ϵ}_0}\cdot \frac{n_{e}B\cos \theta }{{\nu }^3}\hfill \\ \hfill & \hfill & -\frac{q^4}{128{\pi }^4{m}_e^2{ϵ}_0^2}\cdot \frac{n_e^2}{{\nu }^4}-\frac{q^4}{64{\pi }^4{m}_e^3{ϵ}_0}\cdot \frac{n_e{B}^2(1+{\cos }^2\theta )}{{\nu }^4},\hfill \end{array}$$(2)

where n_e is number density of free electrons, B is the magnetic field strength, θ is the angle between the magnetic field B and the electromagnetic wave propagation direction, q is electron charge, m_e is electron mass, and ϵ₀ is electric permittivity in vacuum. In red we show the parameters related to ionospheric conditions and the Earth’s magnetic field. The first term is associated with a dispersive delay proportional to the TEC along the LoS. This is the dominant term; for most radio-astronomical applications at frequencies higher than a few hundred megahertz, higher-order terms can be ignored. The second term is related to Faraday rotation, the positive sign is associated with left-hand polarised signals and the negative sign with right-hand polarised signals. This term depends on TEC and the Earth’s magnetic field. The last two terms are usually ignored but can become relevant for observations at frequencies below 40 MHz. Of these last two terms, the first is dominant and depends on the spatial distribution of the electrons in the ionosphere (Hoque & Jakowski 2008); it becomes larger if electrons are concentrated in thin layers and not uniformly distributed.

Using Eq. (2) we can give order of magnitude estimates of the expected effects at first and second order (see also Petit & Luzum 2010, Ch. 9): $$\begin{array}{lll}\hfill & \delta {\Phi }_1\hfill & =-8067{\left(\frac{\nu }{60\text{MHz}}\right)}^{-1}\left(\frac{\text{dTEC}}{1\text{TECU}}\right)[\text{deg}];\hfill \\ \hfill & \delta {\Phi }_2\hfill & =\pm 105{\left(\frac{\nu }{60\text{MHz}}\right)}^{-2}[\left(\frac{\text{dTEC}}{1\text{TECU}}\right)\hfill \\ \hfill & \hfill & +\left(\frac{\text{TEC}}{1\text{TECU}}\right)\cdot \left(\frac{\text{d}B}{40\mu \text{T}}\right)][\text{deg}];\hfill \end{array}$$(3)

where we adopt a magnetic field B = 40 μT with θ = 45°. The total TEC in quiet geomagnetic conditions can vary from ~1 to ~ 20 TECU from night to day, respectively, and influences the second-order term. Considering a differential TEC (dTEC) of ≃ 0.3 TECU, which is a plausible number for baselines of ~50 km, and observing at 60 MHz, the first-order term produces phase variations of several times 2π. The Faraday rotation instead produces an effect of around ± 30°∕50° (with different signs for the two circular polarisations) at night/day assuming d B = 1%. This effect is not negligible and needs to be corrected (see also Table 1). The effect quickly becomes relatively more severe at lower frequencies because of the 1∕ν² dependency. Higher-order effects can be mostly ignored at 60 MHz. However, at a frequency of ~ 20 MHz, the third-order effect can produce large phase errors.

The paper is outlined as follow: in Sect. 2 we introduce the Low Frequency Array and in Sect. 3the data that we used. In Sect. 4, we present the effect of the ionosphere on radio-interferometric observations at low-frequency. In Sect. 5, we show how our observations can be used to infer ionospheric properties. Conclusions and consideration for future low-frequency experiments are outlined in Sect. 6.

2 The Low Frequency Array

LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) is a radio interferometer that operates at low (including ultra-low) frequencies: 10− 240 MHz. It has 38 stations (aperture arrays capable of multi-beam forming) in the Netherlands, divided into 24 “core stations”, concentrated within 4 km, and 14 “remote stations”, providing baselines up to ~ 120 km¹ . The six innermost stations are packed within 1 km² and are collectively called the “superterp”. Thirteen “international stations” are spread across Europe, but they are not considered in this paper. LOFAR uses two antenna types: High Band Antenna (HBA, used to observe in the frequency range 110− 240 MHz) and Low Band Antenna (LBA, used to observe in the frequency range 10− 90 MHz). In this paper we consider only data from the LBA system, that is the most strongly affected by measurement corruptions induced by the ionosphere. However, most of the results can be extended to higher frequencies.

The LOFAR core is located at 52°54′32^″ N, 6°52′08^″ E. The telescope is marginally affected by ionospheric gradients generated at low latitude by Rayleigh–Taylor instabilities and ionospheric irregularities typical of the auroral regions. At high latitudes, strong refractive index gradients due to field-aligned ionisation structures can cause severe scintillation conditions. The ionosphere conditions in the polar regions are regularly monitored by the Kilpisjärvi Atmospheric Imaging Receiver Array (KAIRA; McKay-Bukowski et al. 2015), a station built using LOFAR hardware in arctic Finland (see e.g. Fallows et al. 2016). This said, the ionospheric impact on a specific LOFAR observation is strongly dependent on the global ionospheric conditions at the time of observation combined with the location of the target in the local sky.

3 Data

For this paper we analysed a set of LOFAR LBA observations pointed at three calibrators: 3C196, 3C295, and 3C380 (see Table 2). For the detailed analysis of ionospheric systematic effects we used a 5.5 hr observation pointed at 3C196 and obtained on May 3, 2013 (18:00 → 23:30 UTC). The observation was taken covering a large continuous bandwidth (22− 70 MHz) so as to evaluate the effect of the ionosphere down to the lowest frequencies. Furthermore, the large bandwidth is an essential tool for separating the various effects based on their different frequency dependencies. The dataset has been calibrated using procedures that will be described in de Gasperin et al. (in prep.). Here we focus on the outcome of the calibration to describe the influence of the ionosphere on the signal measured by the antennas. For each station pair, radio interferometers record streams of data, called visibilities. While ionospheric effects are clearly present in the visibilities, it is easier to analyse them by looking at the solutions. Solutions give a station-based representation of the systematic effects in the form of complex gain factors derived when observing a point source with a known position and flux density (i.e. a calibrator).

We also analysed seven 8-hour-long observations obtained between November 17 and 28, 2017. These observations were taken during the day, pointing the LOFAR beam towards the calibrator sources 3C295 (5 hreach run) and 3C380 (3 hr each run). One observing run of 3C380 failed after two hours and was combined with an additional hour taken 3 days later. These observations have a narrower frequency coverage (42− 66 MHz), centred around the bandpass peak response of the LBA dipoles (~58 MHz). Both the frequency and time resolution of these observations are slightly higher than the one described above; however, this does not affect our results.

Becauseof the low antenna sensitivity and the high sky temperature, the LOFAR LBA system is often in a relatively low signal-to-noise ratio (S/N) regime when observing celestial radio sources. A common way to compensate for this is to average data into larger time or frequency intervals when finding solutions. At low frequencies, the optimal solution intervals are a trade-offbetween S/N and decorrelation. The ionosphere tends to vary very quickly and averaging over more than ~ 5 s in time often does not allow for these changes to be tracked. Furthermore, combining too many frequency channels is also not advisable. As shown in Fig. 1, between the edges of a single LOFAR subband (1SB ≃ 0.195 MHz) centred at 30 MHz, there is a differential phase of 100° (assuming two stations with a differential TEC value of 0.5 TECU). For compact arrays (baselines shorter than a few kilometer, with dTEC ≲ 0.1 TECU), this constraint can be relaxed. To facilitate accurate calibration for all baselines and ionospheric conditions, our calibration has been performed at a relatively high time and frequency resolution, as described in Table 2.

All ionospheric-related systematic effects described in this paper induce direction dependent errors (DDE). This means that the effects vary appreciably as a function of viewing direction, even across the field of view (FoV) of the telescope. LOFAR LBA is characterised by a full width half maximum (FWHM) of ~ 4°. Errors across the FoV are particularly problematic to correct as they require either simultaneous estimation in multiple directions across the FoV (Tasse et al. 2018), or an iterative approach like “peeling” (Noordam 2004; van Weeren et al. 2016). For this paper, we observed only fields whose total flux density is strongly dominated by a central compact source. As a consequence, also solutions are dominated by the systematic effects from that source direction. Therefore, DDE can be treated as direction independenterrors (DIE) whose values are an approximation of the DDE value in the direction of the dominant source.

4 Systematic effects

In Fig. 2 we show phase solutions for each station for 5.5 h of observation of 3C196, in the frequency range 22− 70 MHz. All stations labelled CS are “core stations” – these stations share the same clock and are close together. Stations labelled RS are “remote stations” and each of them has an independent clock. Clocks are not perfectly synchronised and may drift in time with respect to one another. This imprints a time-variable systematic error on the phases of remote stations which is linear in frequency ( ∝ ν). This effect dominates the phase error in the closest remote stations such as RS 106 that in fact shows a rather uniform phase wrapping across the sampled bandwidth. In the most distant remote stations (such as RS 508 or RS 509) the ionospheric dispersive delay ( ∝ 1∕ν) dominates the error budget, and as a consequence phases wrap differently at different frequencies. We can use the different frequency dependency of the two effects (clock errors and ionospheric dispersive delay) to disentangle their contribution to the phase solutions. This process is known as clock/TEC separation and is performed by the LOFAR Solution Tool (LoSoTo² ). The outcome is time-streams of values representing the various systematic effects that we want to disentangle. It is important to note that the clock/TEC separation procedure does not impose any time coherency in the expected systematic effects; each time stamp is treated separately. The fact that the outcome looks well-correlated in time is an indicator that the procedure works as expected. As a final verification we also subtracted the derived phase effect of all the systematic effects combined from the original solutions, this produced the second panel of Fig. 3 which shows rather uniform residuals, meaning that the vast majority of the LOFAR phase systematic effect can be described with a simple model of clock delays and ionospheric effects.

Fig. 1 Ionospheric-induced phase variations between the beginning and the end of a band of 1/2, 1 and 10 LOFAR sub bands (1SB ≃ 0.195 MHz). The dTEC is assumed to be 0.1, 0.3, and 0.5 TECU. These are typical values for distances of a few tens of kilometres. A phase variation larger than ~ 100° creates strong decorrelation. This plot can be used to estimate the maximum amount of averaging (or the maximum band usable to find a single solution) before decorrelating the signal. Coloured bands are the frequency ranges observed by LOFAR.

4.1 Ionospheric (dispersive) and instrumental delays

In the top two panels of Fig. 3 the clock delay and the ionospheric delay are plotted for each station after separation. These values are differential with respect to CS 002, and as a consequence core stations (not plotted) have very small dTEC and a constant differential “clock” due to other instrumental delays. Conversely, remote stations show dTEC values larger than 0.3 TECU for stations that are around 50 km away from CS 002. Furthermore, a strong correlation between the dTEC as measured by neighbouring stations is also visible.

The initial part of the observation is affected by large ionospheric traveling waves. We cross-correlated the dTEC values for RS 407, RS 508, and RS 509 to estimate the wave direction and speed. We obtained a time lag of 137 s (RS 509–RS 508), 177 s (RS 509–RS 407), and 47 s (RS 508–RS 407). These values are compatible with waves traveling at ~ 200 m s⁻¹ from south-west to north-east, values compatible with previous measurements (e.g. Fallows et al. 2016). The speed of these waves might be compatible with the propagation of traveling ionospheric disturbances (TIDs) even if the direction of the propagation is usually equatorward (Cesaroni et al. 2017). In fact, TIDs are associated with traveling atmospheric disturbances (TAD) originated at auroral latitude by Joule heating (Hunsucker 1982). Traveling ionospheric disturbances can have periods of minutes to hours (Hocke & Schlegel 1996). In the second half of the observation the ionosphere becomes less ordered and likely more dominated by turbulent motions. Large variations in the dTEC are not necessarily an indication of a dataset that is difficult to calibrate. An ionosphere that is very active, but relatively coherent across the FoV and with relatively slow variation in time (the first half of the test dataset) is easier to calibrate than an ionosphere that is highly direction dependent with fast scintillations (the second half of the test dataset).

To demonstrate the spatial coherence of the dTEC solutions, we construct spatial screens of the dTEC variations across the array. One screen is made for each time slot solved for in the calibration, fitting to the projected locations (pierce points) of 3C196 onto a thin plane located at a height of 200 km above the ground. Generally, there is one pierce point per station, per direction, per time slot. In the case of the fits discussed here, we use a single direction (that of 3C196); therefore we can have up to 38 pierce points per screen (some of which are flagged; see Fig. 2). We adopt the same method as Intema et al. (2009) which uses Karhunen–Loève (KL) base functions to model the spatial variations of the dTEC values. During the fitting, we adopt a power-law dependence for the phase structure function (see Sect. 5.1), with a power-law index, β = 5∕3. Lastly, we fix the number of KL base vectors fit per screen to five, resulting in five free model parameters per screen.

We plotthe screen fits at two times in Fig. 4. The spatial coherence of the dTEC values at each time is clearly visible in these plots. Therefore, it should be possible to interpolate spatially between pierce points usingthe dTEC screens to obtain dTEC values in directions other than that of 3C196. The screens can also be usedduring fitting as constraints that enforce spatial coherence of the dTEC solutions, thus reducing the number of free parameters solved for during the calibration. We are currently investigating the use of screens in this way.

Fig. 2 Gain phase solutions (from + π: blue to − π: red) for the XX polarisation referenced to station CS 002 (located in the superterp centre). Each panel is one station. White pixels represents bad data that were removed. Stations CS 013, CS 026, CS 031, CS 101, and RS 310 were fully removed due to hardware issues or strong radio frequency interference (RFI) contamination. In the bottom image we show the residuals after subtracting all the effects shown in Fig. 3.

Fig. 3 From top-left to bottom-right: instrumental clock delay (in s). Total electron content variation along the observation (in TECU). Faraday rotation (in rad m−2 ). Ionospheric third-order effect (in rad m−3). All values are differential between CS 001 (assumed constant at 0) and all remote stations (from blue to red in alphabeticalorder).

Fig. 4 Two examples of TEC-screen fits from the observation shown in Fig. 3. Values are in differential TECU withrespect to station CS 002. Each circle represents a LOFAR station; colour coded within the circle is the dTEC value for that station (we note that the markers are filled with the colour corresponding to the measured dTEC). The background colour is the TEC screen fitted across the array as described in the text. Stations crossed with an X are the same excluded from Fig. 2. The full movie is available in the online material.

Fig. 5 As in Fig. 2, but here we show the differential phase solutions in circular polarisation (RR – LL). The plot isolates the effect of Faraday rotation. Since it depends on dTEC and dB, the effect is stronger on stations further away from the reference (CS 002). By fitting a 1∕ν2 dispersive delay we can recover the rotation measure displayed in the third panel of Fig. 3.

4.2 Higher-order terms

With this dataset we were able to measure the second-order ionospheric effect due to Faraday rotation and the third-order dispersive delay effect. Faraday rotation is not easily obtainable from the XX and YY phase solutions. Instead, we took advantage of its different sign in the right and left polarisations (see Eq. (2)). This required a conversion of the dataset from linear to circular polarisation and the extraction of the RR–LL phase solutions (see Fig. 5). Since clock delay and ionospheric first and third terms are scalar, their effect is the same for the RR and LL polarisation and they cancel out if we subtract one from the other. At this point, one can now easily fit the second-order term as a differential delay between the two polarisations with a 1∕ν² dependency.In the third panel of Fig. 3 we show the result in terms of rotation measure (RM). $$\Phi =\text{RM}{\lambda }^2=\frac{2\pi \nu }{c}\cdot \left(\frac{q^3}{16{\pi }^3{m}_e^2{ϵ}_0}\cdot \frac{1}{{\nu }^3}{\displaystyle {\int }_{\text{LoS}}{n}_e}B\cos \theta \text{d}l\right).$$(4)

The term in parentheses comes from the second term of Eq. (2). Again we observe that nearby stations have a similar behaviour and that there is a correlation between the dTEC and the dRM. This is a consequence of the presence of n_e in both the first and second terms of Eq. (2). The measured RM of 0.03 rad produces ~ 45° phase error, which is compatible with expectations given a measured dTEC of 0.3 TECU.

In certain observations, such as the one used in this example, when the ionospheric variation is strong enough, the third-order term can also be extracted by fitting a 1∕ν³ term together with the clock delay and the ionospheric first-order dispersive delay at the clock/TEC separation time. The third-order term becomes relevant only below ~40 MHz but it can rarely be ignored at frequencies close to the plasma frequency (< 20 MHz). In the last panel of Fig. 3 we show the estimated third-order effect using: $$\Phi ={\text{TEC}}_{\text{3}}{\lambda }^3\approx \frac{2\pi \nu }{c}\cdot \left(\frac{q^4}{128{\pi }^4{m}_e^2{ϵ}_0^2}\cdot \frac{1}{{\nu }^4}{\displaystyle {\int }_{\text{LoS}}{n}_e^2}\text{d}l\right).$$(5)

As the first- and third-order terms both depend on the number of free electrons, n_e , not surprisingly they trace each other. In the initial part of the observation, the effect is stronger due to the large dTEC and possibly other ionospheric characteristics (e.g. a thin dominant layer). In the second part of the observation, the third-order term becomes less prominent and harder to fit. As a consequence, imperfect separation in the parameter determination generates small noise-like perturbations in the other estimated parameters. The second-order term is not affected by this problem as its estimation is done with different methods as described above.

4.3 Amplitude scintillations

Scintillations are caused by electromagnetic waves scattered in a non-uniform medium with small changes in the refractive index, such as the ionosphere. A plane wave that enters such a medium with a spatially uniform phase exits the medium with a spatially irregular phase. After propagation to a station, the irregular phases may combine either constructively or destructively (see Kintner et al. 2007). As a consequence, the wave amplitude is increased or decreased and the gain amplitude solutions of the station compensate for the effect by producing an exact opposite trend. To some degree, at frequencies below 100 MHz, we observe that scintillations are always present in LOFAR amplitude solutions.

An example of this is visible in Fig. 6, where an “amplitude wave” crosses the core area (4 km) in 1 min, therefore travelling at a velocity of ~ 240 km h⁻¹. Projecting the linear size of these coherent structures to the height of 300 km, we can estimate an angular size over which amplitude corrections can be considered fairly uniform, which is ~ 20′. Scintillations therefore create a strong direction-dependent amplitude error across the FoV.

Interestingly, wave-like structures are a recurrent pattern across the core stations in the second half of our test observation (see the online material). On the other hand, in the first part of the observation that is dominated by large waves, amplitude scintillations seem to be quasi-simultaneous across the entire array (including remote stations). The amplitude of the scintillations is around ± 10% during both halves of the observation; therefore, it does not seem to depend strongly on the size of the wave phenomena. Although scintillations with timescales larger than a few seconds are highly unlikely to be due to the interplanetary medium, the possibility cannot be entirely ruled out: Kaplan et al. (2015) demonstrated that interplanetary scintillation (IPS) was still visible in observations taken with the Murchison Widefield Array (MWA; Tingay et al. 2013) with a time resolution of 2 s, and Fallows et al. (2016) demonstrated that IPS appears near-simultaneous across the Dutch stations of LOFAR. Furthermore, dedicated observationsof IPS taken with LOFAR show good IPS signal for the radio sources observed here, with both 3C295 and 3C380 used in an IPS observing campaign run during October 2016 and 3C196 used in spring and summertime observations.

Fig. 6 Left to right: time evolution of amplitude solutions for LOFAR core stations for the observation shown in Fig. 3. Each plot is separated from the next by 10 s. Each point represents a core station and each is positioned in the plot according to its geographical location. Colour (blue to red) and size (small to large) are related to the value of the amplitude correction (from − 10% to + 10%). A wave travelling south–north is visible. Similar structures were also reported with single-station observations. The full movie is available in the online material.

5 The ionosphere as sensed by LOFAR

5.1 Structure function

We divided the observations into five time chunks of ~1 h each and calculated the ionosphere phase structure function as described in Mevius et al. (2016). For Kolmogorov turbulence, the phase structure function is a power-law of the form: $$D(r)={\left(\frac{r}{r_{\text{diff}}}\right)}^{\beta }.$$(6)

We fit this function to the data to obtain an estimation of β (expected to be 5∕3 = 1.67 for pure Kolmogorov turbulence) and of r_diff, the spatial scale over which the phase variance is 1 rad², referred to as the diffractive scale (Narayan 1992). Since the phase errors of the core stations are dominated by ionospheric errors, we did not use those phases directly, instead first converting to dTEC. This gives a fast rough estimate of the diffractive scales and thus of the ionospheric quality during an observation. The air mass factor was not corrected for, resulting in a slightly larger variance depending on the elevation angle. We obtained β =1.88, 1.88, 1.92, 1.94, 1.88 and r_{diff@150 MHz} = 7, 8, 16, 16, 22 km, for time chunks 1 to 5, respectively (see Fig. 7). These values for β are larger than expected for Kolmogorov turbulence (β = 5∕3), and are probably due to the effect of large-scale waves, which were not filtered. The baselines that were used for Fig. 7 are much smaller than the typical wavelength of such waves (~ 100−200 km), therefore the effect closely resembles that of a linear gradient over these baselines, which would correspond to β = 2.0 in the phase structure function. The results are similar to those reported in Mevius et al. (2016). The smaller diffractive scales correspond to the first half of the observation, where the magnitude of the dTEC variations is larger.

After removing the large ionospheric gradient (e.g. with a direction-independent calibration), higher orders of the refractive index expansion in Eq. (2) are small, and are assumed to be uniform across the LOFAR beam. To test this assumption we can use the structure function. The LOFAR primary beam size is ~ 4°, that corresponds to 20 km at a 300 km-high ionospheric layer. With a diffractive scale of 10 km and assuming β = 1.7, Eq. (6) gives a phase variance of ~3 rad² at 150 MHz. This corresponds to a dTEC of about 0.03 TECU. Consequently, even at 60 MHz, only the first-order term varies substantially across the FoV, while other terms can be considered constant. This finding can be used to simplify the calibration strategies on the target fields.

Fig. 7 Phasestructure function divided into five time chunks for the observation of Fig. 3. The initial part of the observations is visibly more affected by ionospheric disturbances also in Fig. 3. The phase variance is converted to the expected value at 150 MHz to compare it with other experiments such as LOFAR HBA and MWA.

5.2 Multi-epoch observations

In November 2017, we performed a series of eight observations pointed at the calibrators 3C295 and 3C380. The observations were carried out during daytime hours. We extracted the dTEC and dFR measurements for all of them (see Figs. 8 and 9). A basic expectation of ionospheric models is the presence of an increasing north–south TEC gradient; in all observations, the dTEC values show its presence (red lines are stations in the North, blue lines are stations in the South, all referenced to station CS 002 in the superterp). The wave-like behaviour of the ionospheric systematic effects are clear, although during some time intervals, waves appear more structured and on larger scales than at other times (as seen with different techniques with MWA; Loi et al. 2015a,b).

Ionospheric conditions are driven by a number of factors (e.g. season, solar activity, and latitude). The dominant factor is the day–night cycle. During the night, the ionospheric TEC is reduced; this has an impact on the magnitude of the second-order term of Eq. (2), which is expected to be smaller. However, single-station observations taken under an ionospheric scintillation monitoring project demonstrate that scintillation is often much stronger at night, particularly prior to midnight. We report a large fraction of data loss due to intense scintillation events during night-time observations. During these events we were not able to disentangle the various phase effects, and amplitude solutions show clear signs of decorrelation. Along a seven-night campaign that we carried out in February and March 2015, around 30% of our observations had to be discarded for this reason. Conversely, in 14 days of daytime observations taken in 2017 (7 of which are presented in this paper), only a small percentage of our data were affected by strong scintillations. Since the observations were taken far apart in time, it is difficult to derive strong conclusions. A downside of daytime observations is the (unavoidable) presence of the Sun, combined with the relatively poor side-lobe suppression of a phased array like LOFAR compared to dish-based radio interferometers. However, the low-frequency radio spectrum of the (quiet) Sun is strongly inverted, with S_ν ∝ ν³, and as a consequence, at 50 MHz, the Sun flux density is expected to be a few thousand janskies, which is comparable to the flux density of some other sources present at night. Moreover, solar emission is extended on scales of several arc-minutes; therefore, it is resolved out on moderately long baselines.

5.3 Absolute TEC and satellite comparison

An interesting application is the possibility to estimate the absolute TEC from the differential TEC and differential Faraday rotation. We make the following approximation: $${\text{dRM}}_{12}\approx c\cdot B|{|}_1\cdot {\text{TEC}}_1-c\cdot B|{|}_2\cdot {\text{TEC}}_2,$$(7)

with dRM₁₂ differential Faraday rotation between stations 1 and 2, B||_1,2, the parallel magnetic field, and TEC_1,2, the integrated electron content along the LoS for stations 1 and 2, respectively. We approximate the integral in Eq. (1) with the thin-layer approach. This can be rewritten in terms of absolute TEC at station 1: $${\text{TEC}}_1\cdot \text{d}B|{|}_{12}=C\cdot {\text{dRM}}_{12}-B|{|}_2\cdot {\text{dTEC}}_{12},$$(8)

with dB||₁₂ being the differential parallel magnetic field and dTEC₁₂ the differential integrated electron content between stations 1 and 2. We used the World Magnetic Model (Chulliat et al. 2014) and the measurements of dRM and dTEC to estimate the absolute vertical TEC for the first observation in Fig. 9 in this way. The thin layer model places a single ionospheric layer at an altitude of 450 km. Slant TEC to vertical TEC conversions are done assuming a spherical ionosphere at the same altitude. The resulting vertical TEC at the position of CS 001 as a function of time is shown in Fig. 10. The error bars reflect the spread of the measurements using all different station combinations. In the same plot, we show the interpolated vertical TEC values from three different TEC maps created using the Global Navigation Satellite System (GNSS) and provided by the International GNSS service (IGS). In particular, CODE, UPC, and combined final products have been used (Schaer 1999; Orús et al. 2005; Hernandez-Pajares et al. 2009). The B field parameters and satellite-based TEC values were determined using the RMextract package³ . The measured absolute TEC values are within 10% and following the trend of the GNSS based products from IGS. The time resolution is 10 s, two orders of magnitude better than that of the IGS models which is 2 hr for codg and igsg and 15 min for uqsg. The main uncertainty in the LOFAR measurement comes from the value of d B||, which can be very small. A small unmodelled perturbation of dB|| could artificially enhance the amplitude of the waves as observed in Fig. 10. Here we merely illustrate the main concept. A further discussion of systematic uncertainties of this method will be the subject of a subsequent publication (Mevius et al. in prep.).

6 Conclusions

The ionosphere is the main limiting factor of the quality of low-frequency radio-interferometric observations. The time/space variable refractive index of the ionospheric plasma generates highly direction-dependent dispersive delays that affect phases recorded by the interferometer. At ultra-low frequencies, differential Faraday rotation and higher-order terms (in frequency) also become prominent. The most important results of this paper are as follows.

• We show that LOFAR station-based gain phase can be decomposed into a small number of systematic effects: clock delays, ionospheric effects of first, second (Faraday rotation), and third order. The third-order effect is only important for observations below ~40 MHz.

• We showthat the ionospheric parameters we derived with our decomposition are consistent with expectation (e.g. they are time and spatially coherent) and with independent measurements from satellites. This procedure also demonstrates that LOFAR can be used to obtain independent measurements of the absolute TEC.

• We show that visibility amplitudes are affected by scintillations and that at ultra-low frequencies (<100 MHz) they are always present. We show that amplitude scintillations also follow patterns both in time and space and in certain periods behave like travelling waves across the array.

• We showed that scintillations at night can result in a data loss of ~30% and thatthis phenomenon is reduced during the day. Based on these results we would recommend observation during daytime hours. However, the origin of the scintillations (interplanetary or ionospheric) needs to be clarified and more statistics need to be gathered for a definitive answer.

We also gather here some considerations for future low-frequency experiments:

• Low-frequency telescopes operate in a regime of low S/N due to the high sky temperature. The ionosphere itself contributes with one or more parameters to be estimated every second. It is important that no further degrees of freedom are added to the solving process (e.g. by using different clocks in different stations).

• Again, to maximise the S/N, modern solvers need to be conceived in order to exploit all known coherent structures of the ionosphere, including spatial, time, and frequency coherency. In this paper, time-coherency has never been imposed, although this has already been considered in some modern algorithms (e.g. Tasse 2014).

• Several regimes are identifiable based on the observing frequency and the maximum baseline length (see Table 1). The adopted calibration strategy must vary accordingly, by identifying, given the desired dynamic range, which of the effects are relevant and which have direction-dependent properties.

• The multi-beam capability of modern phased-array interferometers needs to be tested extensively to see how large a reliable TEC screen can be. Furthermore, this capability can be used to develop novel techniques to transfer time-variable instrumental errors from the calibrator to the target fields.

• So far, GNSS data have been marginally used to calibrate radio interferometric observations. However, it would be important to check those data in real time to assess the ionospheric state and therefore decide whether to schedule ultra-low-frequency observations (or any interferometric observation at all).

• In sight of the large data volume produced by SKA, an automatic pipeline to reduce calibrator data could assess the quality of the observation in a short time period and provide a quantitative figure of merit to decide whether to archive the data or to re-schedule the observation (e.g. through plots similar to Fig. 3).

Fig. 8 First-and second-order ionospheric effects quantified by variation of differential TEC and differential Faraday rotation between the LOFAR core and all remote stations. Each panel is a separate observation towards 3C380. The gaps between observations are: 21, 21, 21, 21, 72, 21, 21, and 44 hr. Stations are colour-coded in alphabetical order. The similarity between the top and bottom panels is due to the TEC dependency of Faraday rotation.

Fig. 9 As in Fig. 8 but for 3C295. Gaps between observations are: 19, 19, 19, 19, 91, 19, and 42 hr.

Fig. 10 Absolute TEC derived from World Magnetic Model, dTEC and dRM of the first observation in Fig. 9 comparedto three different interpolated IGS model datasets.

Movies

Movie of Fig. 4 Access here

Movie of Fig. 6 (core) Access here

Movie of Fig. 6 (all) Access here

Acknowledgements

The authors want to thank Claudio Cesaroni for the many discussions and suggestions. FdG is supported by the VENI research programme with project number 639.041.542, which is financed by the Netherlands Organisation for Scientific Research (NWO). M.M. acknowledges support from the ERC (grant 339743, LOFARCORE). This paper is based on data obtained with the International LOFAR Telescope (ILT) under project code LC9_016. LOFAR 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. This work had made use of the Lofar Solution Tool (LoSoTo), developed by F. de Gasperin.

References

All Tables

All Figures

Fig. 1 Ionospheric-induced phase variations between the beginning and the end of a band of 1/2, 1 and 10 LOFAR sub bands (1SB ≃ 0.195 MHz). The dTEC is assumed to be 0.1, 0.3, and 0.5 TECU. These are typical values for distances of a few tens of kilometres. A phase variation larger than ~ 100° creates strong decorrelation. This plot can be used to estimate the maximum amount of averaging (or the maximum band usable to find a single solution) before decorrelating the signal. Coloured bands are the frequency ranges observed by LOFAR.

In the text

Fig. 2 Gain phase solutions (from + π: blue to − π: red) for the XX polarisation referenced to station CS 002 (located in the superterp centre). Each panel is one station. White pixels represents bad data that were removed. Stations CS 013, CS 026, CS 031, CS 101, and RS 310 were fully removed due to hardware issues or strong radio frequency interference (RFI) contamination. In the bottom image we show the residuals after subtracting all the effects shown in Fig. 3.

In the text

	Fig. 3 From top-left to bottom-right: instrumental clock delay (in s). Total electron content variation along the observation (in TECU). Faraday rotation (in rad m−2 ). Ionospheric third-order effect (in rad m−3). All values are differential between CS 001 (assumed constant at 0) and all remote stations (from blue to red in alphabeticalorder).
In the text

Fig. 4 Two examples of TEC-screen fits from the observation shown in Fig. 3. Values are in differential TECU withrespect to station CS 002. Each circle represents a LOFAR station; colour coded within the circle is the dTEC value for that station (we note that the markers are filled with the colour corresponding to the measured dTEC). The background colour is the TEC screen fitted across the array as described in the text. Stations crossed with an X are the same excluded from Fig. 2. The full movie is available in the online material.

In the text

	Fig. 5 As in Fig. 2, but here we show the differential phase solutions in circular polarisation (RR – LL). The plot isolates the effect of Faraday rotation. Since it depends on dTEC and dB, the effect is stronger on stations further away from the reference (CS 002). By fitting a 1∕ν2 dispersive delay we can recover the rotation measure displayed in the third panel of Fig. 3.
In the text

Fig. 6 Left to right: time evolution of amplitude solutions for LOFAR core stations for the observation shown in Fig. 3. Each plot is separated from the next by 10 s. Each point represents a core station and each is positioned in the plot according to its geographical location. Colour (blue to red) and size (small to large) are related to the value of the amplitude correction (from − 10% to + 10%). A wave travelling south–north is visible. Similar structures were also reported with single-station observations. The full movie is available in the online material.

In the text

	Fig. 7 Phasestructure function divided into five time chunks for the observation of Fig. 3. The initial part of the observations is visibly more affected by ionospheric disturbances also in Fig. 3. The phase variance is converted to the expected value at 150 MHz to compare it with other experiments such as LOFAR HBA and MWA.
In the text

Fig. 8 First-and second-order ionospheric effects quantified by variation of differential TEC and differential Faraday rotation between the LOFAR core and all remote stations. Each panel is a separate observation towards 3C380. The gaps between observations are: 21, 21, 21, 21, 72, 21, 21, and 44 hr. Stations are colour-coded in alphabetical order. The similarity between the top and bottom panels is due to the TEC dependency of Faraday rotation.

In the text

	Fig. 9 As in Fig. 8 but for 3C295. Gaps between observations are: 19, 19, 19, 19, 91, 19, and 42 hr.
In the text

	Fig. 10 Absolute TEC derived from World Magnetic Model, dTEC and dRM of the first observation in Fig. 9 comparedto three different interpolated IGS model datasets.
In the text