PROBING FINE-SCALE IONOSPHERIC STRUCTURE WITH THE VERY LARGE ARRAY RADIO TELESCOPE

The ionosphere causes a phase delay to incoming radio waves which is that for radiation passing through a plasma and can be approximated by

$$\begin{equation} \phi \approx \frac{e^2}{c\,m_e\nu }\int {N_e(l)\,dl}, \end{equation} \tag{ 1 }$$

where ϕ is the phase delay in radians, e is the electron charge, m_e is the electron mass, c is the speed of light in vacuum, ν is the light frequency, and the integral is the column density of free electrons (N_e(l)) along the line of sight (Rybicki & Lightman 1979). For a source at zenith and referencing to the VLSS frequency of ν = 74 MHz, this can be simplified to

$$\begin{equation} \phi = 1.14\times 10^2\,\left(\frac{\nu }{74\;\hbox{MHz}}\right)^{-1}\, \left(\frac{\hbox{TEC}}{1\;\hbox{TECU}}\right), \end{equation} \tag{ 2 }$$

where TEC stands for total electron content and a TECU (or TEC unit) is a unit of vertical electron column density equal to 10¹² electrons cm⁻².

In the limit of small refraction angles (total refraction ⩽10⁻³ rad), the refracted path passes through very nearly the same ionosphere as the unrefracted path, and so the total phase delay of a single ray can be approximated by integrating the total ionospheric delay along the original straight-line path. Therefore, a proper ray tracing is unnecessary. We follow this "small refraction angle" approximation throughout this paper.

As a radio interferometer only measures differences between signals received at pairs of antennas, a constant phase delay across an interferometric array will have no effect. If the phase delay is different from one antenna to another, the phase difference between the antennas produces a geometric phase delay due to the wave front coming from a different direction as seen in Figure 1. This phase delay between the paths arriving at two antennas causes an apparent position shift in the source location by an angle

$$\begin{equation} \theta = \frac{\lambda \,\Delta \phi }{2\pi B}, \end{equation} \tag{ 3 }$$

where λ is the wavelength of the light ray, B is the baseline length, and Δϕ is the difference in ionospheric phase delays between the two antennas. If the ionospheric phase delay varies linearly across an array of antennas, all baselines "see" the same source shift, and the only effect on the resulting image of a celestial source is that it appears shifted from its true sky position. Combining Equations (2) and (3), we can relate the gradient of the TEC along a direction x on the ground to the observed angular shift in that direction, θ_x, by

$$\begin{equation} \hspace{-5pt}\frac{d}{dx}\hbox{TEC} = 6.59\times 10^{-5}\;{\hbox{TECU km}^{-1}}\, \left(\frac{\nu }{74\;\hbox{MHz}}\right)^2\,\left(\frac{\theta _x}{1^{\prime \prime }}\right). \end{equation} \tag{ 4 }$$

Note the strong dependence on frequency, ν, which explains why the ionosphere is mainly a problem for low-frequency radio observing. If the ionospheric phase delay does not vary linearly across the array, but also has higher order terms (curvature), the image is not only shifted, but also distorted. However, in the regime of 74 MHz VLA observations with the B-configuration (11 km array), this distortion is rarely significant except under unusual ionospheric conditions.

Zoom In Zoom Out Reset image size

In practice, absolute positions of sources are difficult to measure because that would require differentiating between the ionospheric phases and instrumental phases during the initial calibration of the array. The instrumental phases are determined by observing a bright calibrator source, but these solutions are degenerate to a phase gradient across the array corresponding to the unknown ionospheric position offset of that calibrator source. That phase gradient does not affect astronomical observations as it is corrected in later calibration steps. But it does corrupt the measurement of ionospheric position offsets of other observed sources. Unless this effect is determined and removed, which has yet to be done, it corrupts the measurement of absolute ionospheric gradients.

What can be measured with great accuracy is the change of ionospheric phase gradients across a single field of view. While the absolute ionospheric offsets of sources in this field of view are corrupted by an arbitrary phase gradient applied in the initial calibration, the differences in the position shifts of two different sources in the image can be determined very accurately because the initial corruption applies equally to both sources and so cancels out in the difference. From now on we will refer to these differences in position shifts as the "differential refraction" between pairs of sources measured simultaneously in the same field of view. Source pair differential refraction measures the difference in the ionospheric phase gradient between the paths leading to the two source locations.

At this point it is useful to relate differential refraction to the geometry of the ionosphere. The diffraction-limited field of view at 74 MHz for the 25 m VLA antenna dishes is about 15° across. While this field of view probes a cone-shaped region through the ionosphere, most of the electrons are usually located within a relatively thin layer near the altitude of maximum electron density in the F-layer (Kelley 1989; Schunk & Nagy 2000). It is therefore instructive to consider a thin-shell model of the ionosphere, in which the ionosphere is approximated as a very thin layer at a constant height above the Earth. For a shell height of 400 km (see Section 3.3), the field-of-view cone intersects the shell in a circle of diameter ∼100 km at zenith. In the 11 km B-configuration, the lines of sight from individual VLA antennas toward any single source define a set of pierce points in an 11 km region corresponding to the VLA projected to the F-layer altitude (Figure 2). The observed position offset of that source is determined by the average ionospheric phase gradient in that set of pierce points. The assumption that for individual sources ionospheric distortions are primarily refractive is equivalent to assuming that in any 11 km region of pierce points, the ionospheric variations can be approximated as linear. However, this refraction, or phase gradient, may be different at the location of another 11 km region toward a different source in that field of view, and thus the refraction across the field of view changes and the image becomes "warped." Differential refraction measures the magnitude of the vector difference in ionospheric phase gradients between the ionospheric pierce points toward two celestial sources.

Zoom In Zoom Out Reset image size

Measuring differential refraction rather than absolute refraction may seem like a major limitation. However, the differential refraction of source pairs is actually more relevant to the unique problems of high-resolution low-frequency radio imaging than the absolute refraction of an image. This is because any constant refraction across the field of view causes nothing more than a constant position shift of the entire image, which is easy to correct and remove using existing calibration techniques (such as self-calibration) commonly used at higher frequencies (⩾1 GHz). The differences in the ionospheric phase gradient cause the refraction to vary throughout the field of view causing a stretching or "warping" of the image which standard self-calibration cannot correct. Additionally, this warping is time-variable on a timescale of minutes or less at subarcminute resolution. It is this warping of the field that will be the major challenge for future low-frequency instruments to handle, and differential refraction measurements probe this effect directly.

As for studying the ionosphere, the absolute refraction corresponds to an overall TEC gradient across the entire ∼100 km field of view. Differential refraction measures higher order structure in addition to that overall gradient. To be limited to differential refraction is to be sensitive only to finer-scale structural variations than the roughly 100 km field of view. Larger scale TEC variations can be explored by current techniques such as line-of-sight TEC measurements from ground- or space-based GPS receivers, UV limb scans of recombination emission, radar, and ionosonde. Radio interferometers therefore can extend the measurement of ionospheric structure to much smaller spatial scales than are currently possible.

Ionospheric calibration for the VLSS was done by measuring the relative position shifts of all sources in the field of view that are detected within a 1–2 minute time interval. That time interval has been determined empirically as the maximum length of time during which time variations in the source shifts are small compared to the synthesized beam. The synthesized beam for the VLSS has an FWHM of 80'' (3.9 × 10⁻⁴ rad) which corresponds to a physical size of about 160 m at a zenith height of 400 km. A two-dimensional "phase screen" of ionospheric variations is fit to the position shifts of the detected sources for each time interval. The corrections are applied in the image plane. This method is called "field-based" calibration (Cotton et al. 2004; Cotton 2005) and its specific implementation for the VLSS is described in detail in Cohen et al. (2007). This calibration method works remarkably well for the 74 MHz VLA in most circumstances, producing high-quality images that are astrometrically accurate within about 1/4 of a synthesized beamwidth and with a typical ionospheric smearing that averages only about 1/10 of a synthesized beamwidth. For some fraction of the time, the ionosphere was so turbulent that either it was impossible to image the bright calibrator sources or it was impossible to fit a second-order Zernike polynomial to the source offsets. This was typically 10%–20% of the total observing time.

The data set used for this ionospheric study consists of the position offset measurements for bright (⩾3 Jy⁷) and compact sources in the field of view that were taken and recorded during the field-based ionospheric calibration of all VLSS data acquired in the 11 km B-configuration. This resulted in a list of source position offsets for roughly 3–8 sources for each 2 minute interval in over 500 hr of observing time. In addition to position offsets, the fitted peak brightnesses and integrated flux densities for each calibrator source were also recorded. The time of day was also recorded for each time interval, and from this and the source positions the elevation above the horizon was calculated.

Some radio sources are bright enough that their side-lobes can contaminate the detection of calibrator sources. For this reason, we did not use field pointings within 15° of Cassiopeia A or Cygnus A, or within 8° of Hydra A, the Crab, Virgo A, or Hercules A. In order to remove false source detections and improve the accuracy of the position offset measurements used, we restricted our study to sources within 9° of the field center, with fitted peak brightnesses above 4 Jy beam^−1 8 that were less than 300'' from their known locations. After applying these restrictions, there were 16,032 usable 2 minute time intervals in the survey. During these time intervals, there were 96,575 source detections, and 299,359 pairs of simultaneous source detections.

We can estimate the position error along each axis by using the results of Condon (1997) applied to the case where we have (1) a point source, (2) a circular synthesized beam, and (3) the FWHM of the Gaussian noise correlation function equal to that of the synthesized beam. This results in a position error estimate of

$$\begin{equation} \sigma _{\rm pos} = \frac{\theta _b}{r\,\sqrt{8\,\ln (2)}}, \end{equation} \tag{ 5 }$$

where r is the signal-to-noise ratio of the detection and θ_b is the FWHM of the synthesized beam. Because the average rms noise level in the 2 minute snapshot images is about ∼0.6 Jy beam⁻¹, the 4 Jy beam⁻¹ brightness cutoff ensures that we used sources detected with signal-to-noise ratios of 6.67 or better. Applying Equation (5) then indicates that the predicted position errors are 51 or less.

For each pair of simultaneous source detections, we calculated the magnitude of their differential refraction as

$$\begin{equation} r_{ij} = |\vec{r}_{ij}| = |\vec{r}_i - \vec{r}_j|, \end{equation} \tag{ 6 }$$

where $\vec{r}_i$ and $\vec{r}_j$ are the measured vector position offsets of sources i and j with respect to their known positions taken from the higher frequency (1.4 GHz) Northern VLA Sky Survey (NVSS; Condon et al. 1998) and $\vec{r}_{ij}$ is the vector difference between $\vec{r}_i$ and $\vec{r}_j$. If the two components of the vector $\vec{r}_{ij}$ were independent and normally (Gaussian) distributed with an rms of σ, the magnitude, r_ij, would have a Rayleigh distribution described by the following probability density function:

$$\begin{equation} P(r_{ij}|\sigma) = \frac{r_{ij}}{\sigma ^2}\exp \left(\frac{-r_{ij}^2}{2\sigma ^2}\right). \end{equation} \tag{ 7 }$$

To measure the actual probability distribution of r_ij, we needed to avoid variations in r_ij caused by source pair angular separation and elevation. Therefore, we considered source pairs with angular separations in the narrow range between 5° and 55 and only elevations higher than 60° above the horizon, where elevation effects should be small. There were 6387 such pairs of sources, and their distribution of r_ij is shown as the histogram in Figure 3. Also plotted in Figure 3 as the smooth solid line is the best-fitting Rayleigh distribution. Clearly, the Rayleigh distribution is a poor fit, mostly because of the relatively greater number of very high values of r_ij which extend into a very long "tail" which persists at a non-negligible level for well past 10 times the peak of the distribution. (This, of course also results in a corresponding deficit of values within the main peak of the distribution.) This indicates that the components of the two-dimensional vector $\vec{r}_{ij}$ are not normally (Gaussian) distributed, but rather have a statistical overabundance of very high values, presumably corresponding to times of unusually disturbed ionospheric conditions. These very high values of r_ij constitute a small fraction of the data and so have relatively little effect on the median value of r_ij. However, because those values of r_ij are so high, they have a more significant effect on the mean value of r_ij. Indeed the mean value of r_ij for the distribution in Figure 3 is 331, which is a full 72% higher than the median value of 192. For a Rayleigh distribution the mean and median are $\sqrt{\pi /2}\,\sigma$ and $\sqrt{\ln {4}}\,\sigma$, respectively, which differ by only 6%. Our goal is to measure typical ionospheric conditions. Therefore the median value of r_ij, which is less affected by a small number of extreme events, is the more robust measurement, and that is what we use for all further analysis. We define the term D(θ_ij) to be the median magnitude (r_ij) of differential refraction for sources with angular separation of θ_ij.

Zoom In Zoom Out Reset image size

There is no analytical formula for the standard deviation in the median value of a nonstandard distribution. However, if we assume that the standard deviation is proportional to N^−1/2 (where N is the number of source pairs) and to the median differential shift, we have

$$\begin{equation} \sigma (D(\theta _{ij})) = A\,\frac{D(\theta _{ij})}{\sqrt{N}}, \end{equation} \tag{ 8 }$$

where A is a constant of proportionality. We can determine A by further examining the sample in the previous paragraph. We can estimate σ(D(θ_ij)) for this case by taking random samples of N = 1000 source pairs from the overall sample of 6387 source pairs and measuring the variance in the resulting medians. Taking 10,000 such subsamples, we get an rms variation in the median differential refraction of σ(D(θ_ij)) = 076. We already know that D(θ_ij) = 192, and can solve for the remaining variable, A = 1.25.

We estimate D(θ_ij) itself by binning the source pairs according to angular separation, taking the median value of the magnitude of the differential shifts. The errors are then calculated according to Equation (8). The resulting function is shown in Figure 4. Because of the field-of-view limitations, there are relatively fewer pairs of points with large angular separations and so the error bars are larger for the largest separations. There were pairs with separations of more than 16°, but not enough for meaningful statistics, so we cut off the plot at this value. Note that the median differential source shift generally increases with source separation, though the increase is slower at larger separations. Also note that the differential refraction in the limit of zero angular separation, D(θ_ij → 0), does not appear to be zero. This can be seen more clearly in Figure 5 which plots the median differential refraction for small angular separations with smaller bins so that the trend toward zero angular spacings can be seen more clearly. Figure 5 only shows source pairs that are at high elevation (⩾60°) to remove most elevation effects. A linear least-squares fit of the relation between angular separation and median differential refraction estimates by extrapolation a value of D(θ_ij → 0) = 6467.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are two possible reasons why D(θ_ij → 0) is nonzero. First is that there are significant variations in ionospheric spatial gradients on an angular scale well below about 03. Such an angular scale would correspond to a spatial separation of about 2 km at a characteristic ionospheric height of 400 km (see Section 3.3). But that is much smaller than the 11 km size of the VLA in B-configuration, and therefore would violate our assumption of ionospheric variations being roughly linear within the array size. If that were the case, sources would typically be significantly distorted rather than simply being refracted. This occurs in unusually bad conditions, but is not seen typically.

The second, and more likely, explanation for a nonzero value of D(θ_ij → 0) is map-noise-induced errors in the position measurements of the sources. Thus, rather than tending toward zero, D(θ_ij → 0) tends to the measurement error. The measurement error can be estimated again with Equation (5). The median peak brightness for a source in our sample is 6.88 Jy beam⁻¹. The typical noise level for the 2 minute snapshot observations used to measure source positions is about 0.6 Jy beam⁻¹, for a signal-to-noise ratio of ρ = 11.5. For a resolution of 80'', this implies a measurement error of 296 for each axis. A differential refraction measurement consists of four such measurements (two axes each for two sources) and so the combined error (added in quadrature) is 592. This estimate is less than 10% lower than what the data show, and so this explanation seems reasonable. Ultimately we are interested in the median differential refraction caused by the ionosphere, not measurement errors. Therefore, for all future analysis, we will subtract the value of 6467 in quadrature from all median differential refraction estimates.

We examined the dependence of differential refraction on elevation by splitting the data into three elevation bins as shown in Figure 6. For each elevation bin, we fit the data to a power-law function of the form

$$\begin{equation} D(\theta) = \beta \,\left(\frac{\theta }{1^\circ }\right)^\alpha, \end{equation} \tag{ 9 }$$

where θ is the source-pair angular separation and D(θ) is the median differential refraction within a bin centered on θ. The exponent, α, and scale factor, β, were determined by fitting a linear model to the logarithms of the data points using the least-squares method. The resulting power-law fits are shown as the solid lines in Figure 6 and their parameters are listed in Table 1. As Figure 6 demonstrates, the data are fairly well described by a power-law fit. It is also clear that the differential refraction increases significantly as elevation decreases.

Zoom In Zoom Out Reset image size

Because elevation is mainly a geometrical effect, it should be possible to analytically "remove" its effects on the data. We do this with the commonly used "thin-shell" model of the ionosphere in which the ionosphere is approximated as very thin layer encircling the Earth at a constant height above the Earth's surface (Lanyi & Roth 1988; Ma & Maruyama 2003). In this model, the entire ionospheric phase delay of incoming radiation occurs at the "pierce-point" where the observing path to the celestial source intersects the ionospheric shell. Generally, the shell-height above the Earth's surface, h, is taken to be the height of maximum electron density, with values between about 300 and 500 km typically used. In particular, Nava et al. (2007) found that a shell height of h = 400 km minimized errors when using the thin-shell model to convert from slant to vertical GPS TEC measurements, and we use that value for all upcoming calculations.

According to this formulation, elevations affects the measurement of differential refraction in two ways. First, it changes the relationship between the angular separation of a source pair and the spatial separation of the ionospheric pierce-points of the observing paths to those sources. We correct the observed source-pair angular separation for this by calculating the spatial pierce-point separations for each source pair based on the altitude and azimuth of each source, and converting this to the equivalent angular separation if the sources were centered at zenith. Second, elevation determines the effective path length through the ionosphere. For an angle of incidence, (where = 90° at zenith), the effective path length through the ionosphere is increased by a factor of sec(). The ionospheric phase delays also increase by this factor (see Equation (1)), as does the effective gradient of the phase delays. Because sources within each pair are a few degrees apart, they have slightly different values of sec(), and we used the quadratic average to convert the measured differential refraction to its predicted zenith value.

These two adjustments attempt to convert all differential refraction and angular separation measurements to what they would be if the measurement had been taken for source pairs centered at the zenith. Figure 7 shows the same data as Figure 6, but with the differential refractions converted to their zenith equivalent as described above. As can be seen by the data points and the fitted power-law curves, the applied elevation adjustments have removed nearly all differences between the three elevation bins. We note that our data cover a wide range of elevations, but generally avoid elevations below 30° (see Figure 8).⁹ So it is unclear if these corrections continue to work at lower elevations however, as Figure 7 indicates, the applied corrections appear sufficient for the data used in this study. Therefore, we conclude that elevation effects can be successfully removed analytically, and we adjust all data to their zenith equivalent for all further analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Correcting for elevation effects allows the examination of other effects with far greater accuracy because data from a variety of elevations can now be combined to reduce greatly the statistical measurement errors. In this section, we examine variations due to the time of day. It is known that the ionosphere is thicker during the day, but it is not obvious how that should affect differential refraction. To determine this, we have split the data into day (8:00 to 16:00 local solar time) and night (20:00 to 4:00) observations. The resulting plots of median differential refraction versus source pair angular separation are shown in Figure 9. Clearly, the median differential refraction is significantly greater during the day. For reference, on the right axis of this plot we have converted the differential refraction into the equivalent TEC gradient difference using Equation (4). Also, on the top axis, we give the source-pair angular separation as a spatial separation of the source-pair pierce points projected to the ionospheric F-layer altitude, taken to be 400 km.

Zoom In Zoom Out Reset image size

The error bars generally reflect the quantity of data available. Less observing time was scheduled during the day than at night¹⁰ as can be seen by the histogram in Figure 10. Also, the error bars are greater for the largest angular separations because these sources pairs tend to be at the edges of the fields of view where the primary antenna beam is increasingly attenuated and so fewer source pairs were detected. We also note that the time of year was not evenly sampled, with most of our observations occurring during either early spring or late fall (Figure 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We fit the data to a power law with the form given by Equation (9) in the same manner described in Section 3.3. For night observations, the scale factor is β = 726 and the exponent is α = 0.497, while for the day observations the scale factor is β = 101 and the exponent is α = 0.728. These power-law fits to the data are also plotted in Figure 9. Once again, the power-law relation describes the data extremely well, with the possible exception of the largest separations during the day. In this case, the median differential refraction does not continue increasing quite as fast as the power-law relation would predict. This could be caused by the fact that at such large separations, the ionosphere becomes completely uncorrelated and so further increases in the separation no longer result in increased differential refraction.

To get a better idea of how the ionospheric phase distortions vary throughout the day, we have binned the data into 2 hr intervals of local solar time and produced power-law fits to the data in each bin, the parameters of which are listed in Table 2 along with fits to the day and night observations. From these power-law fits, we predict the median differential refraction at four representative source pair angular separations of 2°, 4°, 8°, and 16° (corresponding to F-layer spatial separations of 14, 28, 56, and 112 km, respectively) as a function of local solar time (Figure 12).

Zoom In Zoom Out Reset image size