The Micro-Arcsecond Scintillation-Induced Variability (MASIV) Survey II: The First Four Epochs

We report on the variability of 443 flat spectrum, compact radio sources monitored using the VLA for 3 days in 4 epochs at ~ 4 month intervals at 5 GHz as part of the Micro-Arcsecond Scintillation-Induced Variability (MASIV) survey. Over half of these sources exhibited 2-10% rms variations on timescales over 2 days. We analyzed the variations by two independent methods, and find that the rms variability amplitudes of the sources correlate with the emission measure in the ionized Interstellar Medium along their respective lines of sight. We thus link the variations with interstellar scintillation of components of these sources, with some (unknown) fraction of the total flux density contained within a compact region of angular diameter in the range 10-50 micro-arcseconds. We also find that the variations decrease for high mean flux density sources and, most importantly, for high redshift sources. The decrease in variability is probably due either to an increase in the apparent diameter of the source, or a decrease in the flux density of the compact fraction beyond z ~ 2. Here we present a statistical analysis of these results, and a future paper will the discuss the cosmological implications in detail.


Introduction
The discovery of centimeter-wavelength Intra-Day Variability (IDV) or "flickering" in Active Galactic Nuclei (AGN) by Heeschen (1984) initially raised concerns that some AGN possess brightness temperatures over six orders of magnitude above the 10 12 K inverse Compton limit for incoherent synchrotron emission (e.g. Quirrenbach et al. 1989). However, considerable evidence has now accumulated to demonstrate that interstellar scintillation (ISS) in the turbulent, ionized interstellar medium (ISM) of our Galaxy is the principal mechanism responsible for the IDV observed in AGN, as was proposed by Heeschen & Rickett (1987).
Two more recent observational techniques provide compelling evidence for the prevalence of ISS. Time delays of 1-8 min are observed in the arrival times of the flux density variations between telescopes on different continents for the three intra-hour variable sources B0405-385, B1257-326 and J1819+3845 (Jauncey et al. 2000;Bignall et al. 2006;Dennett-Thorpe & de Bruyn 2002). The delay arises due to the finite time required for the stochastic fluctuations associated with the ISM to drift across the Earth. A second observational signature of ISS relates to the modulation of IDV variability timescales with a period of exactly one year. This arises because the Earth's orbital motion about the Sun contributes to the effective velocity with which the interstellar scattering material moves relative to an Earth-bound observer; the variations are slow as the Earth moves parallel to the material and fast as it moves anti-parallel to it. Annual cycles in IDV variability timescales are reported in at least seven sources (e.g. Dennett- Thorpe & de Bruyn 2003;Rickett et al. 2001;Jauncey & Macquart 2001;Bignall et al. 2003;Jauncey et al. 2003), including several whose long variability timescales preclude detection of time delays in the scintillation pattern over intercontinental distances. For many lines of sight through the ISM the slowest variations are expected in September if the motion of the turbulent material is comparable to the local standard of rest (LSR).
The recognition of ISS as the dominant cause of IDV has not entirely alleviated the brightness temperatures problems posed by these sources. A source must be small to scintillate; in the weak scintillation case most frequently observed at frequencies near 5 GHz (Walker 1998), the source angular size must be comparable to or smaller than the angular size of the first Fresnel zone, θ F = c/2νπL. Here L is the distance to the scattering region, which we will refer to as the screen even though in some cases it may be better described as a slab extending from the Earth out to distance ∼ L. θ F is typically tens of microarcseconds for screen distances of tens to hundreds of parsecs, which implies source components with angular sizes two to three orders of magnitude finer than the scales probed by VLBI. The long time-scale over which IDV has been observed in some sources suggests that such scintillating components can be relatively long-lived despite their small physical sizes. This paper reports on the results of a Micro-Arcsecond Scintillation-Induced Variability (MASIV) survey for IDV in AGN. This year-long survey conducted observations of between 500-700 AGN over each of four epochs of three or four days duration in 2002 and 2003 at 4.9 GHz with the VLA. The aim of the survey was to provide a catalogue of at least 100 AGN which vary on timescales of hours to days to provide the basis of detailed studies of the IDV AGN population drawn from a well-defined sample. A description of the observations in epochs 2, 3 and 4 is presented in §2 as a supplement to descriptions of the first epoch observations and MASIV source selection in  (hereafter Paper 1). In §3 we describe how the time series of flux density for each source in each epoch was classified as variable or non-variable. In §4 we describe how we have quantified the amplitude and timescale for the variations from an analysis of the structure function combined from all 4 epochs and apply corrections for noise and other sources of flux density error. The basic hypothesis of the paper that the variations are predominantly due to interstellar scintillation is presented and examined in §5, including the influence of parameters of the interstellar medium ( §5.1 emission measure and galactic latitude) and of the sources themselves ( §5.2 mean flux density, spectral index). We now have redshifts for more than half of the sources and we present the dependence of the variability on redshift in §5.5 -a result which shows that the MASIV survey provides a new cosmological probe. In §5.6 we consider whether the variability is intermittent over the four epochs. Our data are listed in Table 4 and our conclusions are presented in §6.

Observations and Data Reduction
The VLA observations took place over four periods during January, May and September 2002 and January 2003 (see Table 1). All epochs were 72 hours in duration except September 2002, which included an additional 24 hours. The additional time in this epoch was added as an attempt to detect the slower variation expected in September due to interstellar scintillation caused by material moving at a velocity comparable to the LSR.
All observations took place during array reconfiguration and in each case the array was being moved to a more compact configuration (except for epoch #1). In each epoch the VLA was divided into five independent subarrays. In the first epoch each subarray observed a subset of the 710 source sample. Subarrays 1 through 4 observed the "core" 578 sources of our compact, flat spectrum sample , namely the weak (105 < S 8.6 GHz < 130 mJy) and strong (S 8.6 GHz > 600 mJy) sub-samples, while the fifth subarray observed a sample of intermediate flux density sources in two regions of the sky. For the next three epochs all sources previously observed in subarrays 1 through 4 were reobserved, thus providing full coverage for the core sample. However, subarray 5 was rededicated to observing a smaller number of sources comprising all objects found to be variable in subarray 5 during the first epoch, as well as the most rapid variables found during the first epoch which required faster sampling in order to determine time scales. In this paper we compare results for the core sample from subarrays 1 through 4 only.
Shadowing was a serious concern at low elevations in the more compact VLA configurations. The impact of shadowing was minimized by assigning 5 or 6 antennas to each subarray in such a way that at any given time a source could be observed with at least 3 unshadowed antennas. Any data from baselines containing shadowed antennas were flagged. Antennas that had either not moved during reconfiguration or were moved but had pointing solutions already applied were assigned to subarrays 1 to 4 where possible.
The core sample was divided into four roughly equal parts by Declination and a subarray assigned to each. The Declination ranges for the four subarrays were 0 • ≤ δ < 14 • , 14 • ≤ δ < 34.08 • , 34.08 • ≤ δ < 49 • , δ ≥ 49 • . The Declination boundary between the second and third subarray was set to the latitude of the VLA to avoid long slews in azimuth when changing between sources transiting north and south of the zenith.
Each subarray was scheduled so that every source was observed for one minute every ∼2 h while it was above an elevation of 15 • . We used the standard VLA frequency configuration for continuum 4.9 GHz observations (dual polarisation and two 50 MHz bandwidth IFs per polarisation) and a 3.3 s integration time. For flux density calibration, each subarray observed B1328+307 (3C286) and J2355+4950 every ∼ 2 hours. B1328+307 is the primary flux density calibrator for the VLA and J2355+4950 is a GPS source, not likely to vary over short timescales, and is monitored regularly at the VLA as part of a calibrator monitoring program.
Following the observations we calibrated the data in AIPS using the standard technique for continuum data. The task FILLM was used to load the data where corrections were made for known antenna gain variations as a function of elevation and for atmospheric opacity. Upon inspection of the data it was clear that there were residual time-dependent amplitude calibration errors. We ascribe this to the fact that in each epoch, some VLA antennas had recently been moved and their pointing calibration observations were not complete. The residual pointing errors may depend on azimuth and elevation so we chose several bright, non-variable sources in each subarray at a range of RA as gain calibrators for surrounding sources. Precautions were taken to ensure that the calibrators themselves were not variable: if a given calibrator caused the majority of sources against which it was applied to vary, then another calibrator was chosen. These sources, typically recommended VLA calibrators, were drawn from our source sample. The numbers of sources used as secondary calibrators in the four epochs were 42, 33, 20 and 36. On the epochs that a particular source was used as a calibrator, it is by definition non-variable and was excluded from the structure function analysis discussed below.
Following calibration, the data for each source were inspected and occasional outlying samples were flagged. The data were then incoherently averaged on a one-minute timescale over all baselines. For all sources, the formal errors obtained were less than those estimated due to the residual constant and fractional errors discussed later in this paper. Incoherent averaging was chosen because, on the assumption that all sources are unresolved, the phase should be zero and any residual phase errors in the data would artificially reduce the average flux density. In low signal-to-noise regimes an incoherent average can induce an upward bias, as visibility amplitudes follow a Rice distribution, which has a mean that is systematically higher than the true mean (Thompson et al. 2001). In the case of our observations though, the signal-to-noise ratio of our observations is sufficiently high that this effect is negligible As the VLA array configurations became more compact, our data became more sensitive to extended structure in the source and in nearby objects. The changing response of the VLA with time due to this structure can appear as variability in a light curve. Fortunately, as our observations were scheduled in sidereal time and were repeated at the same times every day, false variability due to structure appears as a repeating pattern with a period of one day. In general this was easy to recognise but in some cases imaging was carried out to verify contaminating structure. For the purposes of our analysis, such sources were removed from our sample. A total of 102 sources were removed due to structure or confusion, and one was removed due to an error in our initial sample selection process. Following these removals, we are left with a sample of 475 point sources common to all four epochs.

Classification of Variables
We used two separate approaches to analyze the flux density variations. This section describes how we classified the sources based on the apparent modulation index of their intensity. The following section describes a structure function analysis. The virtues of these two techniques are complementary: the first is based on a conservative but robust criterion of source variability, while the second uses a simple statistical estimator which allowed us to quantify the amplitude and time scale of the variations.
In order to ascertain whether an individual source is variable it is necessary to understand the sources of error inherent to the measurements. There are uncertainties in the individual measurements due to calibration errors causing a fractional error, p, and additive errors, s, due to thermal noise and confusion. Calibration errors include contributions from antenna pointing errors, system gain variation between the observations of flux calibrators and variable atmospheric absorption. Since they are a small percentage of the mean source flux density,S, they can be approximated as additive and added in quadrature to the noise as given by σ err,s,p = (s/S) 2 + p 2 . (1) Here σ err,s,p is the rms error in each flux density estimate normalized by the mean for each epoch. In our initial analysis (Paper I) we estimated s = 1.5 mJy and p = 0.01.
An initial inspection classified sources based on a their modulation index, defined as the rms of the three days observations divided by the mean flux density, computed for each epoch. A source was identified as variable if its modulation index exceeded twice the expected contribution from the measurement errors, 2σ err,s,p , as in a χ 2 test. However, direct inspection of the data revealed that many of the slower variables, that is sources with variability time-scales longer than three days, were not detected as variable using this criterion.
We therefore introduced an alternative variability criterion, which classified a source as variable if the modulation index of its daily average flux densities exceeded 2σ err . This process yielded more detections, but visual inspection again revealed that many sources that clearly exhibited variability were undetected by either test. In particular, examination revealed that our two criteria do not detect variability in those sources with low-level monotonic flux density changes during the 3-4 d duration of our observations. This is because the χ 2 statistic used in our two selection criteria is not an ordered statistic and is therefore sub-optimal in detecting such a low level trend.
Our selection criteria were therefore augmented with a visual inspection of the remaining light curves. Inspection was performed by two of us independently and the results compared. Each source was considered non-variable unless otherwise demonstrated, and any source where there was disagreement on its classification was reviewed; we adopted the conservative approach of classifying as non-variable any source is which no final agreement was reached on the classification. In the left panel of Figure 1 we show a scatter plot of the raw modulation index m raw of all sources against their mean flux density from epoch 1, using differing symbols for variables and non-variables. The predicted error in a single flux density measurement from equation (1) is shown by the dashed line, which roughly separates the variable from non-variable classifications. This emphasizes that our classification is determined largely by m raw . Table 4 lists all the raw modulation indices and variability classifications of all the sources and Figures 2 and 3 show some sample light-curves from the survey.
In Table 2 we list the number of sources classified as variable in zero, one, two, three or all four observing epochs and associated percentages. In this table and in all of what follows we define a set of 443 sources from the 475 sources that excludes those used as secondary calibrators in two or more epochs. Note that 12% of the sources were seen to vary during all four epochs. With any analysis of a large number of observations, however, false positives, i.e. sources that are incorrectly classified as variable, are a significant concern. Since the visual classification was very close to the two sigma criterion, our variability classification is reliable with ∼ 95% confidence. Thus we give in columns 3 and 4 the predicted fraction of sources misclassified based on the null hypothesis that all 443 sources are non-variable. Evidently these are negligibly small relative to the observed fractions of variables in two or more epochs.
The large numbers of sources classified as variable on multiple epochs firmly establishes that our classification process is reliable. If all 443 sources were non-variable, then the number expected to be classified variable on 2, 3 or 4 epochs is a mere 6.7 or 1.4%, compared with the observed 192 (43%). Even if our classification was 90% reliable, then the number of false positives with 2 or more detections remains less than 5%.
The fraction of sources that can be reliably classified as variable can be deduced using Table 2. Denoting N as the actual number of non-variable sources, T the actual number of one-time variables and F the number of non-variables misclassified as one-time variables, then the apparent number of non-variables 161 = N − F and from Table 2, F/N = 0.172. This yields N = 194, F = 33 and hence T = 90 − F = 56, and the corrected total number of variables is 192+56=248. Therefore the fraction of sources that exhibited variability in one or more epochs in our survey is 56%. This value is significantly higher than the 15-20% found in previous IDV surveys (Quirrenbach et al. 1992;Kedziora-Chudczer et al. 2001a). In comparing with other surveys one must be careful to consider the selection criteria applied. As described in §2 we started with 710 sources selected from spectral index and mean flux density criteria which was reduced to 443, when we excluded those used as calibrators in more than one epoch and those exhibiting resolution effects, raising the percentage of variable sources. In addition the large number of one-, two-and three-time variable sources raises the percentage. We believe that this is due to intermittency in the IDV phenomenon, which we describe by a simple model in §5.6.
For the purposes of subsequent analysis we conservatively define as "non-variable" those 161 sources that showed no variability in any of the four epochs, and we define as "variable" those 192 sources that showed variability on two or more of the four epochs. With these definitions we have two large and reliable samples each of approximately 200 sources, where the non-variables act as a control sample for the variables. Each was drawn from the same selection criteria and cover the same overall area of sky.

Structure Function of the Variations
Here we discuss how we quantify the flux density variation for each source using the structure function (SF) of each time series, where S j is a flux density measurement normalized by the mean flux density of the source over all four epochs and N τ is the number of pairs of flux densities with a time lag τ binned in 2-hour increments. The SF is a statistically reliable estimator which can be modeled even for short data spans. It is defined independent of any variability classification. Examples are shown in the lower panels of Figure 2 and 3. The SF is preferable to the autocorrelation function for short data spans, which can be badly biased by a poor estimation of the mean.
For an idealized observation of stationary stochastic variations spanning a time much longer than their characteristic time τ char , D(τ ) rises with time lag and tends to saturation at twice the true variance. However for our observations τ char is typically more than two days and saturation is rarely seen. We thus chose D(τ = 2 d) as a standard characterization of the intra-day variations because τ = 2 d is the maximum lag out to which our structure functions contain reliable information for the 72 hour observing sequences. Note that we estimate D(τ = 2 d) by combining data from all four epochs, which makes it a more stable statistic under the hypothesis of stationary stochastic variations. Hence we do not use this statistic to determine whether a source's variability differs over the four epochs.

Measurement errors
Additive flux density errors that are independent of the ISS (or intrinsic) variations contribute an additive term to the SF. If the errors are independent from one sample to the next (i.e. "white"), the error contribution to SF would have a mean value D noise = 2σ 2 err independent of time lag. While random pointing errors will be white, there may be non-white errors caused by systematic pointing errors due to recent antenna relocations and residuals of gain and atmospheric variation, which may contribute a term that is a function of time lag. Nevertheless we now proceed by assuming that all the errors are white and thus the errors add a constant (D noise ) (plus random variations about the constant due to estimation error) and postpone discussion of possible non-white measurement errors.
Our goal is to subtract D noise from the raw SF D(τ ) in order to estimate the SF of any true variations in a source's flux density. But first we must consider how best to estimate D noise . One estimate is 2σ 2 err from equation (1), which depends on the constants s and p. Another estimate of D noise comes from the SF itself -evaluated at its shortest time lag (2h). This would be an unbiased estimate for D noise if all real flux density variation had time scales much longer than 2 h. Hence we examined single-epoch estimates of D(t = 2h), which we plot as an equivalent modulation index m 2h = 0.5D(t = 2h) in the right panel of Figure 1. The sources classified as variable include a substantial number with large values of m 2h , which are due to rapid flux density variations stronger than expected for noise. Note in particular the highest point which is quasar J1819 + 3845 which shows large amplitude rapid ISS. Its timescale is typically shorter than our 2 h sampling and so its SF has already saturated at 2 h and its variations are "white" in our sampling. Similarly the other sources with elevated m 2h are probably due to ISS with short time scales. Now consider the non-variables, which are plotted as pluses and provide a set of sources with low or zero variation in epoch 1, which are useful in studying the noise processes. In the right panel of Figure 1 the mean of the plus symbols lies significantly below the dashed line (equation (1) with s = 0.0015 Jy and p = 0.01) particularly for the higher flux density sources. The solid line with s = 0.0013 Jy and p = 0.005 provides a better model for the noise as discussed in §4.2. In the absence of any real variations the estimates m 2h should be scattered equally above and below their mean value. Comparing the right and left hand panels we see that m raw is typically higher than m 2h even for the non-variables, which suggests that m raw is increased due to low level variations with a timescale longer than 2 h.
The discussion above reduces the choice for estimating D noise to either D(2 h) or the value corresponding to the solid line D err,s,p = 2σ 2 err with p = 0.005. Whereas it is appealing to use the observable D(2h), any rapid real flux density variations can contribute to D(2 h) which thus overestimates D noise . Thus we adopted D err,s,p as an estimator for D noise , which states that the measurement errors are well described by equation(1) but with calibration errors contributing at a 0.5% rather than a 1% level. However we discuss a slight revision of this in the next section.

Structure Function Correction and Fitting
For each source at each epoch we computed the raw SF as defined in equation (2). As the time lag increases the number of available lagged products (N τ ) drops and so increases the error in the SF. A threshold for plotting an estimate for D(τ ) was set at N τ > 20% of the total number of data samples in that epoch. As the example in Figure 2 shows, this gives estimates clustered at lags 2-8 h and near 24 h and near 48 h. In order to characterize the variability amplitudes we initially estimate the SF at a lag of 2 d from the mean SF and its estimated error, calculated using the values at lags in the range τ = 48 ± 2 hours. The overplotted model is described in §4.3.
First consider the equivalent 2-day modulation index m 2 d = 0.5D(2 d) (without noise correction) plotted against the mean source flux densityS for epoch 1 in Figure 4. It should be compared with the left and right panels of Figure 1 for the non-variable sources. The comparison makes it clear that m 2 d > m 2 h for almost all of the non-variable sources. Indeed there is a close correspondence between m 2 d and m raw . Thus it is clear that there are low level flux density variations on a timescale longer than 2 h in the sources classified as non-variable. In order to investigate what causes these we examined the difference ∆D = D(2 d) − D(2 h) for the non-variable sources from each epoch.
We examined how ∆D depends on mean flux density, Galactic latitude and Hα emission -quantities which we find in §5 influence D(2 d) for the variable classifications. The results showed that the mean of ∆D in each epoch was significantly higher for the sources weaker than 0.4 Jy than for those stronger than 0.4 Jy. However, it showed no significant dependence on Galactic latitude or Hα emission. Hence the process responsible for the lowlevel variability in the sources classified as non-variable is not ISS and is unlikely to have an astrophysical origin.
We consider the most likely cause to be confusion, which can be due either to extended source structure which is partially resolved on the baselines of each sub-array or to low level confusing sources in the primary beam. We had already eliminated obvious cases of confusion by removing sources whose light curves showed clear daily patterns of variability which repeated in each day of a 3 day sequence. Since our SF estimation is from the light curves normalized byS an increase in ∆D at lowS is consistent with the effect of low level confusion characterized by a certain rms in Jy. This is supported by the VLA documentation that gives 2.3 mJy as the brightest source expected in a single antenna beam at 5 GHz.
We examined the shape of the mean D(τ ) for high and lowS sources in each epoch. Evidence for the effects of confusion was found in the consistent minima in D(τ ) near time lags τ = 1d and 2d. However, there was also a component rising from the noise floor at 2h and typically saturating between 12 and 24h. All three components were substantially larger for the weak group of sources than for the strong sources. When averaged over four epochs the noise floor was, respectively, at ∼ 3.0 × 10 −4 and ∼ 0.61 × 10 −4 and the averages for ∆D were, respectively, ∼ 2.8 × 10 −4 and ∼ 1.1 × 10 −4 for these groups (mean flux density 0.13 Jy and 1.4 Jy). By equating these average SF amplitudes to 2σ 2 err using equation (1) we estimated s = 0.0013 Jy and p = 0.0073. The s values for each epoch ranged from 0.001-0.002 Jy and for p 0.003-0.01, but no consistent patterns were seen versus epoch.
We also considered the effect of long-term variations as characterized by the variation between epochs of the mean flux density from each epoch. These slower variations are mostly intrinsic and might contribute to a trend within the 3 or 4 days of each epoch. So we calculated the rate of change of flux density from the average of the magnitude of the differences in flux density between neighbouring epochs, divided by the number of days between epochs and the resulting magnitude of D(2d) due to such a trend. Surprisingly, for all sources this was smaller than the noise-corrected D(2d) except for those which were negative as a result of noise subtraction. Further the highest of them was 0.0001, which is one quarter of the threshold value. So we conclude that long-term variations did not make a significant contribution to the variations observed within the epochs. We note, however, that we can usefully estimate the structure function on time lags of 3 and 6 months from the MASIV survey data, which may be useful in studies of intrinsic variability. Thus we include the epoch averaged flux density for each epoch in the data table 4.
The foregoing studies reveal that the apparently non-variable light curves include not only white noise but also a low level contaminating process whose rms values can be approximately characterized by equation (1). While we can quantify the white noise process by this equation, the non-white contamination is not well enough understood to be reliably characterized by an SF which could be subtracted from the estimates for all sources at each epoch. So in order to reduce the effect of these non-white variations, we increased the estimated D noise to to be subtracted by using the values s = 0.0013 Jy and p = 0.007, and we also set a threshold on D(2 d) = 4 × 10 −4 below which the noise-corrected SF values may be contaminated and therefore should only be interpreted as upper limits to the true SF of flux density variations.
It is of interest that, since confusion effects will be precisely repeated at 24 h intervals, the samples of the structure function at 24 and 48 hr will be unaffected by confusion. The results we present in subsequent sections are from fitting the SF over all time lags, since the fit makes better use of the data. However, we also analyzed the single sample estimates of D(24h) and D(48h) and found very similar results, though with somewhat worse statistical errors. As an extra precaution we re-reviewed all the light curves and SF plots for 24 h periodic patterns and found 34 sources that might be contaminated at a level near the threshold. For these sources we used the lower of D(2 d) from the fit and that estimated from 48 ± 2 h.

Variability Timescales
Although we estimated the SF as described above for each source at each epoch, the single-epoch D(2 d) is only based on about one independent sample of a 2 d variation and so it has a large statistical error -such that its rms error is about equal to its mean (Rickett et al. 2000). Thus we do not attempt to evaluate the variability amplitude on a source by source basis for each of the epochs. Rather we average the SF for each source over all 4 epochs (after subtracting the D noise as defined in the previous section). Hence the SF results are insensitive to any intermittency or annual changes in ISS due to the effects of the Earth's velocity.
We then fitted the following simple model to the SFs: where 2m 2 is the amplitude at which the function would saturate and τ char is the characteristic timescale where the SF reaches half of its saturation. The motivation for this form is described in Appendix A. It approximates the form expected from ISS caused by a turbulent interstellar medium uniformly distributed through a thick scattering region, as opposed to turbulence confined to a thin layer (see equation (A4)). We estimated two parameters only: the timescale τ char and the value of D(τ = 2 d) for each source. It should also be noted that a light curve that is dominated by a linear trend in flux density gives rise to a parabolic SF, which is not well fitted by equation (3). The value of D(τ = 2 d) in such a case will be somewhat underestimated.
An example of fitted SF is shown in Figure 2. The points show D(τ ) increasing (noisily) with time lag τ but not reaching saturation. Since the timescale is defined at half the saturation value it is poorly constrained in this example: τ char = 1.0 ± 0.5 days. However, from the same fit D(t = 2 d) = (44 ± 3) × 10 −4 , which is quite well constrained. With observations limited to 3 days (4 days for epoch 3) it is not possible to estimate timescales longer than about 3 days. However, in those cases it was possible to recognize that the characteristic timescale is longer than 3 days from the shape of the structure function. Two other examples are shown in Figure 3, in which there is evidence for faster variations. For source J0949+5819 the variations are very strong in epoch 1 and much weaker in epoch 3. From the epoch-average we find τ char = 0.02 ± 0.05 days and D(t = 2 d) = (19 ± 2) × 10 −3 . For J1328+6221 the variations are more consistent over the epochs with τ = 0.2 ± 0.2 days and D(t = 2 d) = (10 ± 1) × 10 −3 . Given our 2 hr sampling and the typically large fractional errors in the timescale we have simply classified the timescale into fast τ char < 0.5 d, medium 0.5 < τ char < 3 d and slow 3 d < τ char . We also looked for any correlation between the timescale and D(t = 2 d) but found no consistent pattern.
During the visual examination of the light curves for each source at each epoch, the timescales were estimated by counting the number of inflection points (i.e. change in sign of the derivative) for those epochs classified as variable. Since in the visual examination there was an effective smoothing, inflection points due to noise-like deviations were not counted. The majority of sources were found to show none or at most one inflection point indicating variability timescales that are predominantly longer than 3 days. The observed distribution of inflection points is shown in Figure 5. Only a small number of sources (20%) showed 2 or more inflection points. A comparison of the distribution of inflection points for the weak and strong sources revealed no significant difference between the two classes. Overall, the distribution of timescales was statistically the same for each epoch, remembering that epoch 3 was four days rather than the three days of the other epochs. An important conclusion from the timescale study is our 3 or 4 d lightcurves commonly underestimate both the timescale and true modulation index for many of the sources.
The annual cycle reported in a number of IDV sources is due to the changing relative velocities of the Earth and the ISM responsible for the scattering . If the ISM velocities follow the Local Standard of Rest, many sources would be expected to exhibit slower variations in the third quarter of the year (i.e. during the third epoch), and hence may more easily be missed because of the lengthened time-scales. Figure 5 shows the numbers of variables found at each of the four epochs. A Chi-squared contingency test shows no evidence that the numbers differ from the mean in any epoch, even though epoch 3 lasted for four rather than three days. The uniformity of variable numbers in each epoch suggests a lack of evidence for a third-quarter slow-down, and it follows that the majority of the scattering material is not moving at the LSR. This is perhaps not unexpected; both PKS 1257-326 and J1819+3845, the two sources for which reliable screen velocities have been measured, have measured screen velocities that differ significantly from the LSR (Bignall et al. 2006;Dennett-Thorpe & de Bruyn 2003;Linsky et al. 2007).
In summary we used the visual analysis to classify each source at each epoch as variable or not variable. We computed SFs for all sources and then examined the SFs of those classified as non-variable in order to quantify the measurement errors. We were able to correct the SFs by subtracting a constant versus time lag due to errors that are independent over the 2 h sampling and are characterized by equation 1. In addition we found a low level contaminating process with a timescale of 1-2 d, which we suggest is due to low-level confusion with an rms of 1-2 mJy. The SF of the slower contamination could not be reliably estimated and so sets a limit on the minimum detectable variation in flux density. By fitting a simple curve to the epoch-averaged and noise-corrected SFs, we estimated D(t = 2 d) for each source. The contamination is minimized by requiring this quantity to be above a threshold of 4×10 −4 for a useable estimate of the timescale τ char . Most sources were classified as slow variables.

Comparison of SF with the visual variability classification
We now compare the SF analysis with the variability classification from §3. Figure 6 shows the number of epochs in which a source was classified as variable plotted against its value of D(2 d) obtained as described from a fit to the structure function of the cumulative data from all 4 epochs. The large circle shows the mean values for each group of sources. While it is clear that the sources with higher D(2 d) were classified as variable more frequently, there is a very wide distribution in the rms level of the variation over 2 d among the sources. The vertical line marks the threshold D(2 d) = 4 × 10 −4 above which we have made a timescale estimate. Values below this should be regarded as upper bounds in view of the possibility of low level confusion. Table 3 lists the source counts sorted by the number of "variable" epochs for the 443 sources (as always excluding those used as calibrators in more than one epoch.) It also shows the mean values of D(2 d) and the numbers of sources above and below the SF threshold. In total 37% of them are above the threshold versus 45% from the variability classification on 2 or more epochs. In the latter classification process we do not characterize the rms amplitude but attempt to quantify any intermittency in the phenomenon. In contrast the SF analysis quantifies the rms amplitude over 2 d averaged over all 4 epochs.

Interpretation as Interstellar Scintillation
We now examine our basic hypothesis, stated in the introduction, that the variations in flux density detected in the MASIV survey are caused by interstellar scintillation (ISS). The extremely small diameters of pulsars revealed the ISS phenomenon almost as soon as pulsars were discovered and they still provide the best information on the distribution of small scale structure in the ionized interstellar medium -indeed they provide a calibration of the ISS phenomenon. The pulsar observations have been combined into a model for the distribution of electron density in the interstellar medium by Taylor & Cordes (1993) and revised by Cordes and Lazio (2005).
Because the refractive index of an ionized medium varies with radio frequency, there is a transition frequency (f w ) above which the scintillation of a point source, like a pulsar, is weak in the sense that its scintillation (modulation) index (m pt ) is less than one. This frequency is on the order of 5 GHz but depends on the strength of the turbulent fluctuations in electron density on the given line of sight (Walker 1998;Cordes and Lazio 2005). Above f w the ISS of a point source has a single timescale approximately given by t F = r F /V , where V is the effective velocity of the Earth through the ISS diffraction pattern and r F = Lλ/2π and L is the typical distance to the scattering region.
In our observations the angular diameters of the extra-galactic sources are considerably larger than those of pulsars and so their ISS is heavily quenched. See Rickett et al. (2006) for a discussion of how the "low-wavenumber approximation" can be applied for quenched scintillation even below the transition frequency. If we approximate the scattering as taking place in a thin region of the Galactic plane, we obtain simple expressions for the reduction in scintillation index and lengthening of timescale (e.g. Rickett (1986)). In Appendix A we apply the same simple "screen" model to the structure function analysis and obtain expressions for how the observable D(2 d) might vary with angular size of each source and on distance to the scattering region and the level of turbulence on each line of sight. Of course the level of ISS also depends on properties of the source -in particular the fraction of its flux density in its most compact component and on the effective diameter of that component.

Galactic Dependence of ISS
We start by comparing our D(2d) results with the emission measure (column density of the square of the electron density) as estimated from observations of Hα emission. We find the intensity of Hα emission (in Rayleighs) from the WHAM Northern sky survey on a 1 degree grid (Haffner et al. 2003) nearest to each source. We use the intensity summed over all velocities, which Haffner et al. (2003) interpret as proportional to the ISM emission measure on that line of sight, assuming the temperature of the emitting gas does not vary by a large percentage. We expect the level of ISS to be related to the emission measure on that line of sight, as described by Spangler & Cordes (1998)   We stress the complete independence of the two data sets in this figure and that the bin averages are independent of any threshold set on D(2d). We exclude the extreme IHV source J1819+3845 from this and subsequent bin average plots because its D(2d) value is 0.25 which is so much higher than the next highest at 0.015 that it distorts the mean and the variance within its bin.
The bottom panel shows that the fraction of slowly scintillating sources clearly increases with emission measure and vice-versa for the fast scintillators. This finding that the longer timescales occur when seen through greater column density of electrons is consistent with enhanced ISS from strongly ionized regions of the ISM, which are typically at low Galactic latitudes and at greater distances L. An increase in L increases the scale of the scintillation pattern which slows the scintillation time -see Appendix A.
Since the strength and effective distance of the scattering layer depends on Galactic latitude, we also expect a dependence of ISS on Galactic latitude. For our visual classification of sources (as variable or non-variable) we asked the simple question "are their latitude distributions the same?". A Chi-squared contingency test dividing the sources into two samples, a low latitude sample, |b| < 40 degrees, and high latitude sample, |b| > 40 degrees, shows that the two distributions differ at the 98% confidence level. There are fractionally more variables at low latitudes than there are at high latitudes, supporting ISS as the origin of the intra-day variability.
For the structure function analysis we simply plot D(2d) against the Galactic latitude of each source, in a fashion similar to that of Heeschen & Rickett (1987). The upper panel of figure 8 is a scatter plot, differentiated by the timescale group (fast, medium or slow). The middle panel averages D(2d) into 30 degree wide bins for all sources, which as already noted is independent of the threshold on D(2d). There is a low level of scintillation above 60 degrees, increasing in the mid range (30-60 deg), in both northern and southern hemispheres. However in the low latitudes (0-30 deg) the ISS increases in the south of the plane but decreases north of the plane. While the figure is in reasonable agreement with the latitude dependence found by Rickett et al. (2006) in their analysis of the modulation index of 146 flat-spectrum sources observed at 2 GHz with the Green Bank Interferometer, there are competing effects in our MASIV survey since data were sampled for no more than 4 days.
As the latitude b decreases both the distance L to and the path length through the scattering medium increase (∝ cosec b). The increased path length makes the scintillation stronger at lower latitudes so that the modulation index should increase. However, the increased distance L increases the scale of the scintillation pattern which slows the scintillation time so that the structure function will saturate at times longer than 4 days, causing a decrease in D(2d). The combination of these two effects requires careful modeling.
The model for the structure function described in Appendix A is a starting point for analyzing the effects of Galactic latitude. Equation (A5) shows that as the distance L increases D(2d) should decrease. However, the equation assumes that the scintillation index of a point source m pt = 1, omitting any increase due to the longer scattered path length at lower latitudes. In a more realistic model m pt may be less than one looking out of the Galactic plane (|b| ∼ 90 deg) and should increase with decreasing latitude, reaching unity at ∼ ±45 deg (Walker 1998). At lower latitudes still it will increase only slowly due to effects of refractive ISS. The increase in m pt between 90 and 45 deg partially compensates for the reduction in D(2d) due to increasing distance L. At still lower latitudes D(2d) might be expected to decrease. In the observations one sees a difference of low latitude behavior between the Northern and Southern hemispheres. This asymmetry can be understood by looking at the center panel of Figure 8 which shows that the emission measure is commonly higher for Southern latitudes.
In a complete model the (unknown) compact fraction f c of the source flux density in the scintillating component must also be considered. Thus the expected variation of the ISS level with Galactic latitude must be combined with the probability distributions for the flux fraction and for the diameters of the compact components, which will further dilute the variation of D(2d) with latitude.
In the Green Bank Interferometer observations cited above, the data were sampled daily (or on alternate days) over many years and so provided estimates of the actual modulation index, which were not reduced by the lengthening ISS time scale at low latitudes. However, even for these data the latitude dependence is not a strong effect. Note that the asymmetry about the Galactic plane is very similar to the asymmetry in the typical Hα emission as shown by the circles in the center plot. The lower panel of figure 8 plots the fraction of sources in the three timescale groups versus latitude and shows clear evidence that the fast scintillators dominate at high latitudes and that slow scintillators dominate at low latitudes. This agrees with the expected increase in ISS timescale at low latitudes outlined above.

Dependence of ISS on Source Spectral Index and Flux Density
As stated earlier the ISS of extragalactic sources is expected to be strongly suppressed, relative to that of pulsars, by the smoothing effect of their larger angular diameters. Consequently, we expect that the more compact sources should show higher levels of ISS. Here we examine the influence of mean flux density (S) and spectral index α (S ∝ ν α ), since we expect synchrotron emitting sources to be more compact for larger α and lowerS, due to the effects of synchrotron self-absorption and inverse Compton losses.
In Heeschen (1984)'s initial survey of short-term variability of a large sample of both steep-spectrum and flat-spectrum sources he found that the flat-spectrum sources varied, "flicker", but the steep-spectrum sources did not. This can be understood as the steepspectrum sources are dominated by optically thin synchrotron emission with low brightness temperatures, while the flat spectrum sources are dominated by synchrotron self-absorbed components with very high brightness temperatures (Scheuer & Williams 1968), making them compact enough to show ISS.
The MASIV sources were selected to have flat spectral indices α > −0.3  and so are predominantly quasars with compact cores. However, it is possible that there are also some very compact galaxies in the sample. Figure 9 shows the spectral index distributions separately for the sources with the visual classification as variable or non-variable. The spectral indices shown are those used to form the sample: 1.4 GHz NVSS (Condon et al. 1998) to 8.5 GHz JVAS (Patnaik et al. 1992Browne et al. 1998;Wilkinson et al. 1998) or CLASS (Myers et al. 1995) flux densities. This shows a slight increase in the fraction of sources that are variable with increasing α, in agreement with the expectation that the flatter (and inverted) spectrum sources are more compact. Though in Figure 10 the mean D(2d) shows no significant trend with spectral index, the bottom panel shows a slight increase in the fraction of variable sources for α > 0. This is mostly due to an increase in the fraction of slow variables which constitute the largest timescale group. It is worth pointing out that the surveys from which the flux densities were drawn to obtain spectral index were not coeval. It is likely then that any change in sample properties as a function of spectral index will be blurred as many of the sources vary intrinsically.
Turning to the influence of mean flux density, we first discuss the visual variability classification and then plot D(2d) against flux density. The selection of sources for the MASIV survey divided them into a high mean flux density group (strong) and a low mean flux density group (weak S < 0.3 Jy). As reported in Paper 1 there was a greater fraction of variable sources in epoch 1 from the low flux density group than from the high flux density group. Combining all four epochs the numbers of weak sources that varied in 0-4 epochs is (94,46,36,39,33) and for strong sources the numbers are (62,45,49,22,23). Thus there are significantly more 3 & 4 time variables among the weak sources than among the strong ones, though this trend is not supported in the 2 or 1 time variables. Figure 11 shows D(2d) versus mean flux density. The center panel shows a clear downward trend with increasing flux density in the lowest three bins, while in the fourth bin it increases but with fewer sources the mean has a large error. The interpretation is an increasing angular diameter of the compact source components with increasing total mean flux density. Furthermore, the bottom panel shows a decrease in the fraction of fast and medium scintillators with increasing mean flux density. These are exactly the trends expected if their effective angular diameters are constrained by a maximum brightness temperature due to self-absorption or inverse Compton losses (θ ∝ (S/T B ) 0.5 ).

Source Models
The foregoing analysis establishes that about half of the 443 compact flat spectrum radio sources in the MASIV survey show ISS at an rms level above 1% over times of 2 d. We now consider a simple model for the compact source structure based on Appendix A. Equation (A5) gives an approximate relation between D(2d) and parameters of the source (f c and θ src ) and of the interstellar medium (L and V ). We proceed by assuming a basic model for the latter parameters L = 500 pc and V = 50 km s −1 and finding constraints on the source. Figure 12 shows (solid) contours of D(2d) versus the source parameters: compact component diameter and compact flux fraction -defined in the observer's frame. Also shown are (dashed) contours of T b /S Jy . The majority of sources are in the range 0.0004 < D(2d) < 0.01 which, for sources of 0.1 Jy mean flux density, maps to maximum brightness temperatures 10 12 K to 10 14 K. These figures require substantial Doppler factors in the AGN jets comparable with those estimated from VLBI. However, we note that the plot of Figure 12 only provides a guide since it is based on 500 pc as the distance to an interstellar scattering screen.
The distribution of scattering electrons along the line of sight to each source is likely to be much more complex and can extend from tens of pc to a few kpc. Scattering at tens of pc has been shown to be important for the rare rapid scintillators (IHV) and scattering from more than 1 kpc is responsible for the slower timescale ISS associated from sight lines with large emission measures. Since the implied angular diameters scale roughly with the distance, the uncertainty in L corresponds to an extremely large range in implied brightness temperatures. We note, however, the fast ISS sources in our sample are mostly scintillating at levels of only 1-5%, unlike the large amplitude variations in the well studied IHV sources B0405-385, B1257-326 and J1819+3845. This suggests that the nearby scattering regions responsible are relatively rare covering only a small fraction of the sky. Lazio et al. (2008) discuss the relative importance of sparsely distributed "clumps of scattering material" and a more uniformly distributed interstellar scattering plasma, suggesting that the former could be more important for ISS and the latter for angular broadening in the ISM. These authors find minimum diameters of ∼ 1−2 mas at 1 GHz, which they suggest is caused by interstellar scattering, which predicts 40 − 80µas when scaled to 5 GHz. It will be important to use the full MASIV data set to re-examine these questions, but since our emphasis here is on presenting the data, we postpone them to a later paper.

Combined effects of Intrinsic Variations and ISS
Here we examine the competing contributions that scintillation and intrinsic variability would potentially make to the measured lightcurves. A source at an angular diameter distance D A undergoing intrinsic variations on a timescale τ has an implied maximum intrinsic angular size: where D is the Doppler factor. Further, a variation in flux density ∆S implies an observed brightness temperature of When mapped into the emission rest frame the brightness temperature is then reduced by a factor D 3 (1 + z) −3 , under the hypothesis of intrinsic variation.
As an example consider a source that undergoes 100 mJy intrinsic fluctuations in 100 days, as observed between epochs for several of our sources. At a typical distance D A ∼ 1Gpc the implied maximum source size would be ∼ 17 Dµas. Suppose further that it does not show ISS within an epoch, which implies it must be larger than 80 µas in the observed frame and so D > 4.6. Hence mapping equation (5) into the emission frame gives T B,em < ∼ 5 × 10 10 K. We conclude that sources showing intra-epoch (intrinsic) variation and no ISS have relatively low emission frame brightness and, conversely, higher brightness sources that show intra-epoch variation have to show ISS.

Dependence of ISS on Source Redshift
We found redshifts for about half of the 443 sources in the survey from the published literature and we have subsequently measured another 69 (Pursimo et al. 2008) for a total of 275 redshifts. This constitutes the largest sample of ISS measurements versus redshift. Figure 13 plots D(2d) versus redshift and reveals a highly significant decrease in the prevalence of ISS as redshift increases. In particular the middle panel shows that when binned in redshift the mean level of ISS drops steeply above redshift 2. As in other plots in this format the binned averages are independent of any SF threshold. However, note that the exact value of the lowest binned averages of D(2d) depend on the details of the noise subtraction. We note that, since the effect of confusion is to slightly increase D(2d), our estimates become upper limits when below the threshold of about 0.0004.
The bottom panel suggests that the fraction of fast variables drops more quickly with redshift than the fraction of slow and medium variables. However, the error bars show that this is only marginally significant, and in view of the importance of this question we list the source counts in Table 5. If true it would imply that the drop in mean ISS level is due to an increase in angular diameter with redshift, which also lengthens the ISS timescale. While the drop in ISS level seen in the middle panel could be due either to an increase in diameter of the compact core of emission from these sources or a decrease in the compact fraction of flux density in that core, the latter interpretation would not explain a steeper decrease in fast ISS with redshift than in slow ISS. The most likely conclusion from this analysis is that the extremely compact emitting regions responsible for the ISS in over half the quasars studied appear broader in angular diameter with redshift above 2. The interpretation of this result involves either an evolution in the emitting regions with redshift or an angular broadening phenomenon due to propagation. We caution that a full consideration of selection effects must be made when interpreting this result. For example, redshifts are currently more complete in the strong (84%) than the weak (43%) sub-samples. Therefore, although this new cosmologically important result is the major finding from the MASIV survey, we postpone a full discussion of the interpretation to a paper in preparation (Macquart et al. 2008) pending a thorough investigation of selection effects (Pursimo et al. 2008). Interested readers can consult preliminary discussions of the interpretation by Rickett et al. (2007). We also note that Lazio et al. (2008) plotted angular diameter at 1 GHz against redshift from a much smaller sample of scintillating and non-scintillating sources but could draw no firm conclusions.

Intermittent Variability
We now ask whether the sources that were only classified as variable in one to three of the four epochs are varying episodically or are the result of statistical uncertainty and a fixed threshold for variability on the raw modulation index. As a result of the measurement uncertainties there can be both false positives and false negatives, whose probabilities we can estimate. Concentrating on the 91 sources that were classified as variable in only one epoch and correcting for false positives leaves us with 61 sources ∼ 13% of the total. These may be a category of short-lived, episodic scintillators revealed by our regular sampling over a full year.
We note the strong case for intermittent ISS in the rapid variable PKS 0405-385 (Kedziora-Chudczer et al. 2001b), and so consider a simple model for the intermittency in terms of the longevity of each IDV episode and the duration between episodes. Intermittency may arise either due to fluctuations in the level of turbulence responsible for ISS, or it may arise if the lifetimes of the bright, microarcsecond source components that undergo ISS are short. In the latter case one might expect sources to undergo IDV in conjunction with each outburst of the central engine.
Consider a simple model in which the IDV episodes have a finite duration ∆T and an interval, T f , between outbursts. Obviously in reality both quantities will have a distribution of possible values, but given the infrequency of our time sampling we restrict ourselves to this simple assumption here. For any given source the IDV episode commences at some random time t i ∈ [0, T f ), with the probability distribution of episodes distributed uniformly: If we make a single observation the probability that the source will be exhibiting IDV is so for a survey that examines N intermittent sources the expected number showing IDVs at any one time is N∆T /T f . One can similarly calculate the probability of detecting IDV in a source during one or more of multiple observing epochs of a multi-epoch survey (see Appendix B). In particular, the probability of detecting IDV in a source in one or more epochs of a four epoch survey, with epochs separated evenly in time by t obs , is The number of IDV sources detected is a maximum when the interval between observing epochs exceeds the duration of IDV episodes because for shorter epoch intervals, after one merely discovers few new IDVs after the first epoch, one only keeps reobserving all the IDVs that were present in the first epoch. (Obviously, in the limit when t obs is small, multiple observations discover the same number of sources as a single-epoch survey.) Now, the mean detection rate of IDV sources in each epoch is 30%, whereas the fraction of sources that exhibited IDV in one or more of our four epochs is 58%. These two numbers imply a typical burst duration ∆T = 1.2 yr and a burst period of T f = 3.8 yr.
We can also calculate the corresponding probability that IDV is observed in a source in all four epochs: Based on the foregoing estimates of ∆T and T f one estimates that only 4% of all sources should be common to all four epochs. However, the actual detection fraction is 12%.
It should be remembered that this model does not take into account several effects which are likely to be important, including: (i) annual cycle effects influence the number of sources one detects at any one epoch (ii) there is likely a distribution of episode durations and repetition rates, (iii) the repetition is likely irregular and (iv) not all IDVs are likely to be episodic. We favor (ii), the importance of which is illustrated by the fact that many one-time IDVs were seen in epochs 2 & 3 that were (obviously) not detected in epoch 4. However, our model implies that sources that commenced IDV in these epochs should have been detected subsequently because the predicted burst duration exceeds the interval between observing runs (i.e. ∆T > 3t obs ).

Conclusions
We have reported results from the four epochs of the MASIV survey. There were 710 sources with flat spectra (α < −0.3) near 5 GHz selected in weak and strong flux density groups surveyed for variability in four epochs over a year. These flat spectrum sources are predominantly quasars with compact emission probably from a core and jet, many with effective diameters small enough to show interstellar scintillation (ISS). In each epoch the flux density was measured using sub-arrays of the VLA every 2 hrs for about 12 hrs each day for 3-4 days. Sources were removed from the study if they showed evidence for changing correlated flux density due to confusion or resolution of their more extended structure, leaving 443 sources which were analyzed for variability in two ways.
The first was a binary classification based on the raw modulation index (visual method) in which 43% of the sources were classified as variable in 2, 3 or 4 of the epochs. The second was a fit to the epoch-averaged structure function parameterized by D(2d) and a timescale τ char . In view of the uncertainties in the latter we classified sources as slow (> 3 days), medium (0.5-3 days) or fast (< 0.5 days) if D(2d) exceeded 4 × 10 −4 . By this criterion 37% of the sources varied with more than 1.4% modulation index over 2 days, which is similar to the 43% variables by the visual classification.
We found that D(2d) and timescale varied both with coordinates in the Galaxy and also with source-based quantities. This confirms that the variations are dominated by ISS, which depends on both the strength of scattering and the distance to the scattering region and also on the fraction of flux density in its most compact component and its effective angular diameter. The following is a summary of our findings: • The amplitude of 2-d variability increases with increasing emission measure estimated from Hα intensity for each line of sight. Emission measure is the column density for the square of the electron density which is expected to be strongly correlated with inhomogeneity in the ionized medium that causes ISS. This result provides observational evidence that ISS is the dominant cause of the variations. We find fast variations dominate for low emission measure, as expected since such regions will be seen out of the plane and closer to the Earth, and slow variations dominate for high emission measure which are typically seen at greater distances toward the Galactic plane especially for southern latitudes where the Hα intensity is high.
• The amplitude of 2-d ISS variability varies significantly with Galactic latitude but differs substantially between positive and negative latitudes. The expected behaviour is complicated; greater path lengths at low latitudes, where the scattering should be stronger, cause the scattering to be slower which should reduce the rms over 3 d. However, the observed timescales show that there are more sources with fast variations at high latitudes and more sources with slow variations at low latitudes in both hemispheres, in clear support of ISS as the dominant cause.
• The ISS modulation index tends to decrease with increasing mean flux density, as expected if the compact emission is limited by synchrotron self absorption or inverse Compton losses to have a maximum brightness temperature. In that case the expected angular diameter ∝ √S which will quench the ISS of the stronger sources.
• There is little change in the ISS amplitude with spectral index for our sample with α > −0.3.
• There is evidence that the ISS can be intermittent on times of 3-6 months for some sources, but this is hard to quantify from the 3 day observing sequences, when the time scale of the variations is of the same order.
• We model D(2d) as a function of compact source component fractional flux density and angular diameter, from which we find compact diameter to lie in the range 0.005 -0.15 milli arcseconds and brightness temperatures in the range 10 12 − 10 14 K.
• The most far-reaching result reported here is the discovery of a decrease in both the fraction of sources that scintillate and in their scintillation amplitude beyond redshifts around 2. We conclude that there is an increase in the typical angular diameter of the most compact radio-emitting regions of the quasars beyond reshift 2. The possible interpretations of this exciting result will be presented in a companion paper (Macquart et al. 2008).
• A further surprise (at least to us) was the apparent absence of the very rapid variables (IHV). J1819+3845 fell in our sample, but it was the only source to show remarkable large amplitude variability. J0929+5013 showed rapid variability in the January 2002 epoch ) but, although monitored closely, revealed only slower, manyhour variability in the three later epochs. We had expected to find more of these rapid variables especially given that two of the three known, J1819+3845 and PKS1257-326 were discovered serendipitously. This strongly suggests that the IHV sources lie behind discrete local interstellar clouds which cover a small fraction of the sky.
The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We are extremely grateful for the technical support provided by NRAO staff at Socorro, in particular we would like to thank Ken Sowinski, Miller Goss, Mark Claussen and Jim Ulvestad for helping to implement the five subarrays at the VLA. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We also made use of data from the Wisconsin H-Alpha Mapper, which is funded by the National Science Foundation. BJR thanks the US-NSF for funding under grant AST 0507713. He also thanks both the Cavendish Astrophysics group at Cambridge and the ATNF, Epping for hospitality.

A. Structure function for ISS
In all of the MASIV data the intrinsic source diameters (θ src ) are large enough to suppress the scintillations relative to those of a point source (such as a pulsar), which at 6 cm would be expected to have a true modulation index m pt near unity. In this section we describe a model for the structure function for such extended sources; this is easier to interpret than the apparent modulation index m i .
The source diameter smoothes the ISS of a point source and so reduces the modulation index to m θ and lengthens the timescale to τ θ . If the scattering is concentrated at a distance L from the Earth and we are near the transition from weak to strong scintillation, a useful approximate relation is: Here θ F = λ/(2πL) is the angular size subtended by the Fresnel scale (r F ) and θ src is the angular radius of the source. The exponents 7/6 and 7/12 apply for a Kolmogorov spectrum of interstellar turbulence (Coles et al. 1987). To the same accuracy the ISS timescale for a point source would be τ F = Lθ F /V and the source smoothing would increase it to: where V is the effective transverse velocity of the observer relative to the layer of scattering plasma. Note that when the source diameter is sufficiently large to suppress the scintillations we have where we have set m pt = 1 (Rickett 1986).
The theoretical form of the auto-correlation function for an extended source that substantially suppresses the scintillation index is given by the low wavenumber approximation of Coles et al. (1987). This in turn gives the theoretical structure function, whose detailed shape depends on both the spectrum of the plasma density and on its distribution along the line of sight through the Galaxy. Figure 14 of Rickett et al. (2006) shows the form for sources with Gaussian brightness with peak temperature 10 11 − 10 13 K, when the medium is modeled by the Cordes and Lazio (2005) electron distribution at a Galactic latitude of 45 • (away from the Galactic center). The form of the structure function at small time lags depends strongly on the density distribution in the local ISM. A useful approximation to these results is given by: where 1 ≤ a ≤ 2 is a constant that depends on the density distribution in the local ISM. Here a ∼ 2 for a local bubble with low turbulence such that the effective scattering distance is beyond the bubble (> 100 pc) and alternatively a ∼ 1 if the medium is uniformly turbulent out to a scale height (as in the "disk" of the CL05 electron density model). We have also introduced an extra variable f c that is the fraction of the source flux density in the bright (core) component.
Equations (A4) and (A3) thus provide a simple interpretation for our estimates of D(t = 2d). Inserting m θ and τ θ from equation (A3) we obtain: An important feature of this result is that D(2d) decreases steeply with increasing θ src and so provides a sensitive measure of source diameter. In estimating D(2d) for every source-epoch we set a = 1, which allows us to estimate the scintillation timescale τ θ without having to also estimate the exponent a. D(2d) can be converted to an effective 2-day modulation index by m 2d = 0.5D(2d). We note that m 2d can substantially exceed the apparent modulation index m i when the time scale is longer than about 2 days.

B. A simple model for IDV intermittency
Consider a model in which a source outbursts every duration T f in time and each IDV episode lasts ∆T . The initial outburst time is unknown, but its probability is evenly distributed in the interval , which assumes the value one whenever the source shows no IDV but the value zero when it is on.
Thus the fraction of the time in which the source is off for initial burst durations between t i and t i + dδt i isf (t i )p(t i )δt i . Thus the probability that IDV is absent is Thus the probability that the source is observed to exhibit IDV is 1 − p 1off = ∆T /T f . Now suppose we look for IDV at times t = 0, t obs , 2t obs , 3t obs . The fraction of the burst times between t i and t i + δt i for which IDV is absent in all four observations isf (t i )f (t i + t obs )f(t i + 2t obs )f(t i + 3t obs )p(t i )δt i . If we assume that the repetition period exceeds the total duration of our observations (i.e. T f > 3 t obs ), the probability of observing no IDV over all four epochs is Thus the probability that IDV is observed in any one or more of these four epochs is We can similarly consider the probability of observing IDV in multiple observations by defining the function f (t) = H(t) − H(t − ∆T ) which takes the value one whenever the IDV is on and zero otherwise. The probability that IDV is observed in all four epochs is thus This preprint was prepared with the AAS L A T E X macros v5.2.      Table 3).       Table 4. MASIV source variability characteristics. Column 4 gives Hα in Rayleighs for WHAM. Columns 5-8 give mean flux density for epochs 1-4; Columns 9-12 give raw modulation index in % for epochs 1-4; Column 13 σ err as %; Column 14 gives the visual classification per epoch for variable (y), non-variable (n) and non-variable source used as a secondary calibrator (c); Column 15 gives D(2d) from fit to epoch average SF where error -1 flags a non-convergent fit; Column 16 gives timescale classifications, 1:τ char < 0.5d, 2:0.5d< τ char < 3d, 3:τ char > 3d, 0: D(2d) below threshold of 4 × 10 −4 .