Interstellar Scintillation of an Extreme Scintillator: PKS B1144$-$379

In this paper, we use data from the COSMIC (Continuous Single-dish Monitoring of Intraday Variability at Ceduna) observing project to investigate the variable radio source PKS B1144$-$379. The COSMIC project was designed to provide high-cadence monitoring of the flux density variations of several strong sources of radio IDV using the 30-m Ceduna radio telescope. The observations started in early 2003 and continued until 2013. The Ceduna 30-m telescope was equipped with an uncooled receiver with two orthogonal circular polarizations operating over the frequency range of 6.4–6.9 GHz (centred at 6.7 GHz, covering a bandwidth of approximately 500 MHz). This frequency was selected as it corresponds to the approximate transition frequency between weak and strong scattering for extragalactic sources at high Galactic latitude and the transition frequency is where the variability amplitude is largest (Walker, 1998, 2001). More details about the COSMIC project, including the procedures used for the observations and an outline of the data reduction, can be found in McCulloch et al. (2005). The methods for reduction of COSMIC data were refined by Blanchard (2013) and we have utilised those methods here.

Here we give a brief summary of the observational method. A single measurement of the flux density of a source is obtained from two pairs of orthogonal scans. The first pair of scans has one taken with increasing Right Ascension with constant Declination, followed by another taken with decreasing Right Ascension. The second pair holds the Right Ascension constant and scans first in increasing, then decreasing Declination. All the scans are initially centred on the nominal source position, but the online observing system estimates pointing offsets from the orthogonal scans and updates the position for future observations. Updated pointing offsets are only utilised if less than four hours old, otherwise, the nominal pointing position is used. Scans are saved in FITS format files, along with appropriate metadata. Obtaining a flux density measurement from the two pairs of scans is a multi-step process, with each group of four scans toward a single source considered separately. The amplitude data is first normalised using the receiver noise diode and the amplitude as a function of position on the sky is then fitted with a Gaussian profile and a polynomial baseline. The next step is to average together the amplitude data for the pair of Right Ascension (constant Declination) scans and similarly for the pair of Declination scans. The two averaged scans are then compared for consistency and a group of four scans is only accepted for further processing if the amplitude for the two average scans is the same to within 10 $\%$. The next step is to estimate the pointing offset between the observed and actual source position and correct the measured amplitude for that offset. The observed amplitude for the Right Ascension average scan will be reduced by a Declination offset between the observed and measured source position. We measure this offset from the position of the peak in the Declination average scan, calculate the appropriate correction factor using the known shape of the beam profile and apply it to the Right Ascension scan. The same process is then repeated for the Declination average scan using the offset measured in the Right Ascension average scan. The two pointing-corrected averaged scans are then themselves averaged. Finally, a gain-elevation correction is made based on a gain curve calculated from one year of observations of PKS B1934$-$638, a southern calibrator known to be flux density stable, with an assumed flux density of 3.9383 Jy at 6.7 GHz (with an uncertainty of 2$\%$; Reynolds, 1994). The gain-elevation curve is applied to correct for small changes in the gain of the antenna with elevation.

To reduce the amount of observing time lost to antenna slewing the COSMIC project split its target sources into two groups, those that passed north and south of the zenith. Each group of sources was typically observed for periods of 10–14 days at a time, alternating between the two groups. PKS B1144$-$379 was an exception and from early 2004 was included in both groups (meaning near-continuous monitoring). Most sources (including calibrators) were only observed in one group and for northern sources, 3C 227 (assumed flux density of 1.9885 Jy) was the flux density calibrator (with an uncertainty of 2$\%$; Baars et al., 1977). The amplitude scale for the data was converted into flux density in Jansky (Jy) by calculating the noise diode amplitude through comparison with either PKS B1934$-$638 (southern calibrator) or 3C 227 (northern calibrator). The COSMIC project made observations of the target and calibrator sources only when they were at elevations greater than 10 degrees and the measured gain variations of the antenna are less than 5$\%$ (McCulloch et al., 2005). It was found that after scaling using the gain-elevation correction, the measured root mean square (RMS) of the flux density for three calibrator sources (PKS B1934$-$638, PKS B0945$+$076 (3C227) and PKS B1921$-$293) has a term that is some fraction of the source flux density. Fitting the observed RMS of these calibrator sources with the quadrature function gives a constant error term of 29 mJy and a fractional term of 0.84$\%$ of the source flux density (in Jy). All the scintillating sources observed in the COSMIC programme have mean flux densities at 6.7 GHz that are typically in the range 1–4 Jy.

After applying the calibration procedures developed by Blanchard (2013), the examination of the flux density calibrators shows RMS residual variations of approximately 3$\%$. These variations have been identified as systematic errors due to a small temperature dependence of the noise diode amplitude (which is not in a temperature-controlled environment), and this, in turn, creates systematic errors in the gain-elevation curve. In the previous analysis using COSMIC data, the amplitude of the noise diode was assumed to be constant and was determined from the assumed flux density of the calibrator source (a constant value which does not vary with time). However, the temperature dependence of the noise diode amplitude means that daily temperature variations result in apparent diurnal variations in the flux density of the calibrator. Figure 1 shows the variations in the external temperature and the flux density of PKS B1934$-$638 measured assuming a constant amplitude for the noise diode, over a 6 day period. Diurnal variations in the apparent flux density of the calibrator are clearly visible in the lower panel of Figure 1. We investigated a number of different approaches to correct for the temperature dependence of the noise diode, but were unable to find any method which reliably reduced the systematic variations below 3$\%$. While we have good data on the external air temperature during the observations, the noise diode is located in the receiver room, where the ambient temperature differs from the external temperature (Senkbeil et al., 2008) and changes in the ambient temperature of the receiver room lag external temperature variations, but not by a constant nor predictable amount. The temperature dependence of the noise diode affects the flux density measurements of both the calibrator and the target sources. Fortunately, the amplitude of the typical ISS-induced variations in the target sources is much larger than 3$\%$, however, it is important during the data analysis that any methods applied to the target source must be applied to the calibrator as well. In this way, the effect of any systematic variations in the flux density of the target source can be estimated. It is vital to properly characterize the systematic variations in the flux density of the calibrator and target source, as it may impact the measurement of the characteristic timescale of variations, and through that subsequent time series analysis of the target source.

We filtered the flux density data for the calibrators PKS B1934$-$638 (southern calibrator) and 3C227 (northern calibrator) in order to remove extreme outliers. The mean and standard deviation ($\sigma$) of the flux density for one year of observations were calculated and any data that differed from the mean by more than 5$\sigma$ were eliminated. Removing the extreme outliers reduced the standard deviation in the flux density of the calibrators by $\sim$7$\%$. We then applied a smoothing procedure using a Simple Moving Average (SMA) technique to the flux density of both calibrators. This technique is a type of Finite Impulse Response (FIR) filter which is applied to a set of data points by creating an average of different subsets of the full dataset (Azami et al., 2012). The observations of the flux density calibrators for each year were smoothed with 20$\%$ of the sample size of the calibrator data as the weighted average. The next step is to measure the scale factor for the amplitude of the noise diode. The scale factor was determined from the ratio of the smoothed flux densities and the assumed flux density of the calibrators (3.9383 and 1.9885 Jy for the southern and northern calibrators respectively at a frequency of 6.7 GHz). Figure 2 shows the flux density data for the calibrator PKS B1934$-$638 (southern calibrator) after it has undergone all the procedures to reduce the systematic errors. After applying all the procedures, the standard deviation of flux density of the calibrator was reduced by $\sim$22$\%$. As mentioned before, it is important during the data analysis that any methods applied to the target source must be applied to the calibrator as well. We, therefore, measured the noise-diode scaling factor for the target source from the scale factor of the calibrator. This procedure was done by linearly interpolating scale factor of the target source from the scale factor of the calibrator for the particular period of the observing epoch (MJD). Figure 3 illustrates the flux density of the target source PKS B1144$-$379 prior and after we applied the scale factor. The scale factors for both calibrators and target sources are less than 5$\%$. The final estimate of the flux density for the calibrators and target source was obtained by multiplying the data with the scale factors. Prior to undertaking time-series analysis, the data for each source were split into “blocks”, which typically have a duration of $\sim$14 days. The blocks generally correspond to the periods of time for which we observed either the southern or northern source groups. For most blocks, we have flux density data for the calibrator and target sources at least once per hour for the times that they were above 10 degrees elevation for the Ceduna antenna and there are no breaks in the data for any block of more than 5 days. Any outliers in each block of data for target sources were removed before performing the time-series analysis. A previous analysis of COSMIC data used a different approach to removing the systemic variations that involved smoothing and resampling the data (e.g. Carter et al., 2009). Carter et al. (2009) combined the use of discrete autocorrelation functions and the related power spectral density to estimate the characteristic timescale of the variations. This approach was effective, however, it does intrinsically low-pass filter the timescale of variations. It was therefore not possible to reliably determine if there are multiple variability timescales present in the data with the method of Carter et al. (2009). In the present paper, we analyse the data using a different method and without smoothing or resampling, as we have reduced the level of systematic variations (3$\%$) well below the level of radio IDV observed for PKS B1144$-$379 which has high modulation index of 5–18$\%$. Our method potentially allows the detection of multiple variability timescales present in the data, although for the present paper we define a single characteristic timescale.

There are several methods that have been used to measure the variability timescale of radio IDV. For example Simonetti et al. (1985) and Dennett-Thorpe & de Bruyn (2003) use the structure function (SF) of the measured time series to define the characteristic timescale, while Rickett et al. (1995, 2002) estimate the variability timescale from the half-width at half maximum of the autocorrelation function (ACF), Kedziora-Chudczer et al. (1997) used the mean peak-to-peak time, Jauncey & Macquart (2001) adopted the peak-to-trough time and Gabányi et al. (2007b) estimated the timescale using 3 different methods; average peak-to-trough of the light curve, SF and the ACF. Due to the inherently stochastic nature of the flux density variations in AGN produced by ISS, it is difficult to define a characteristic timescale by simply examining the light curve (see for example Figure 4). For our analysis we chose to use the SF of the flux density variations (the light curve) to estimate the characteristic timescale. Prior to calculating the SF, the amplitude of the light curves of the calibrator and target source were normalised by dividing all the flux densities for each source light curve by the mean flux density of that source for the specific observing block. This ensures that the amplitudes of the SF are related to the fractional variation and can be directly compared for the calibrator and target sources. The SF method was first introduced by Simonetti et al. (1985) and further described by Spangler et al. (1989). The SF is defined to be the average squared difference between points of a data series, as a function of lag. The SF $D(\tau)$ is computed in a similar manner to the Discrete Correlation Function (DCF) by taking the difference between every point in the data series with every other point. The individual differences are then sorted according to the time between the differenced data points and binned accordingly

The SF of a source exhibiting variation due to ISS is usually described by a power law rising from a value of zero at zero lag to the peak plateau level, at a given lag. The amplitude of the SF is linked to the RMS difference at a given lag by twice the variance $(2\sigma^{2})$. In our case, the peak amplitude of the SF of all blocks exceeds the twice variance level (“overshoot”) and this can potentially lead to overestimating the variability timescale. The upper right-hand panel in Figure 5 shows the structure function calculated for the PKS B1144$-$379 light curve shown in Figure 4. The SF in the upper right-hand panel of Figure 5 shows a plateau commencing at around a lag of 0.5–1 day before an “overshoot” which peaks around 2.3 days. The overshoot is linked to quasi-periodic behaviour in the variability light curve and is often observed for radio IDV sources (e.g. Rickett et al., 2002). To better estimate the characteristic timescale, instead of trying to fit to determine the amplitude of the SF plateau, we use the theoretical estimate of the saturation level (plateau) of twice the variance (red line in upper right-hand panel of Fig. 5). We have used the time lag at half of the saturation value for our estimate of the characteristic timescale (blue line in upper right-hand panel of Fig. 5). The estimated characteristic timescale corresponds to the length of the average peak-to-peak modulation divided by 4.

As previously mentioned the Ceduna data suffers from $\sim 3\%$ systematic variations which are predominantly diurnal and manifest as a peak at 0.5 day lag (see grey points in upper left-hand panel of Fig. 5). Therefore, prior to estimating the characteristic timescale for target sources, we subtracted the SF of the calibrator (for which we assume the variations are residual systematic gain and diurnal variations) from the SF of target source. The resulting SF we assume to be a “clean” SF for the target source. We take the “clean” SF of the target source and performed a linear fit of the SF data with an amplitude less than twice the variance (saturation level). The intercept between this straight line fit and the variance (half the saturation level) we use as our initial estimate of the characteristic timescale. However, because there is intrinsic scatter in the SF the “final” (best) estimate of the characteristic timescale and its uncertainty is obtained by calculating the Gaussian probability distribution given by equation 2, which estimates the probability of a particular lag $\tau$ corresponding to the point at which the SF is equal to the light-curve variance. The bottom-left panel of Figure 5 shows the Gaussian probability distribution for the SF in the upper right-hand panel. We take the peak of the Gaussian probability distribution to be the characteristic timescale and its 1$\sigma$ uncertainty is computed from the Full Width Half Maximum (FWHM) of a Gaussian fitted to the Gaussian probability distribution. Quantile-Quantile (Q-Q) plots (lower right-hand panel) were constructed in order to quantify the reliability of Gaussian fitting using a Gaussian probability distribution (Fig. 5). A Q-Q plot is another graphic method for testing whether a dataset follows a given distribution, in which our case is the normal probability distribution. It differs from the probability plot in that it shows observed and expected values instead of percentages on the X and Y axes. If all the scatter points are close to the reference line, we can say that the dataset follows the given distribution.

We adopted the method of Rickett et al. (2002) to quantify the uncertainty in the estimated SF for each lag $\tau$ using their formula:

Limited sampling is also the dominant source of uncertainty in determining the amplitude of variability. One measure of flux density variability in a light curve is the modulation index and this has been used by various authors to characterise AGN variability (e.g. Rickett et al., 2001; Carter et al., 2009; Marchili et al., 2012). The modulation index, m = $\sigma_{\textnormal{S}}$/$\bar{S}$, is defined as the RMS variation $\sigma_{\textnormal{S}}$ normalised by the average flux density of the source $\bar{S}$. Although the modulation index has the advantage of being non-negative and more robust against outliers, it still consists of a convolution of intrinsic source variation and observational uncertainties. The large modulation index could be an indication of either a strongly variable source or a faint source with high uncertainties in the measurement of individual flux density. Due to this effect, the accurate measurement of modulation index requires one to take into account the measurement errors and the uncertainty in the modulation index caused by the finite number of flux density measurements. To examine the accuracy of our measurement of modulation index, we decided to measure the intrinsic modulation index, $\bar{m}$ = $\sigma_{\textnormal{o}}$/$\bar{S}_{\textnormal{o}}$, which is the intrinsic RMS variation $\sigma_{\textnormal{o}}$ normalised by the intrinsic average flux density of the source $\bar{S}_{\textnormal{o}}$. This method was introduced by Richards et al. (2011) where they used a likelihood analysis to obtain the intrinsic modulation index as well as its uncertainty. The term “intrinsic” is used to denote flux densities and variations as would be observed with a perfectly uniform sampling of adequate cadence and zero observational error. Thus, $\bar{m}$ is a measure of the true amplitude of variations in the source, rather than a convolution of true variability, observational uncertainties, and effects of finite sampling. Observational uncertainties and finite sampling will affect the accuracy of measurement of $\bar{m}$. Here, we decided to use the measured modulation index since the percentage difference between the modulation index (m) and the intrinsic modulation index ($\bar{m}$) only ranged between 1$\%$–17$\%$. The modulation index for the target source varied between 5$\%$ and 18$\%$ and its uncertainty was taken from the uncertainty of the intrinsic modulation index.

An investigation of the reliability of Gaussian fitting to the Gaussian probability distribution of the target source SF was performed through Monte Carlo simulations. In this simulation, the initial model was a simple sinusoid with mean flux density of 10 Jy, a modulation amplitude of 3 Jy, a period of 0.16 days (equivalent to $T_{\textnormal{period}}$/4), a light-curve duration of 3.5 days and constant Gaussian white noise of amplitude 1 Jy. The light-curve duration of 3.5 days was chosen as it represents the period of first quarter of each observation block of the target source ($\sim$14 days) where the first “scints” can be observed. We adopted the first “scints” in estimating the variability timescale to reduce the uncertainties due to less sampling. The average power of the signal ($P_{\textnormal{s}}=\textnormal{amplitude}^{2}/{2}$) and average power of uniform white noise ($P_{\textnormal{n}}=\textnormal{amplitude}^{2}/\sqrt{3}$) were 4.5 and 0.6 respectively, and the measured SNR (signal to noise ratio) of the signal is 8 ($\textnormal{SNR}=10.\log_{10}\left\{(P_{\textnormal{s}}~{}-~{}P_{\textnormal{n}})/P_{\textnormal{n}}\right\}$) (Sijbers et al., 1996). In order to examine how changes in data sampling frequency affect the light curves on the SF pattern, the peak of Gaussian fitting and the value of its FWHM (assumed as the uncertainty of the estimated characteristic timescale), 5 simulated datasets of light curves were created. These simulated light curves consisted of 991, 772, 496, 234 and 178 data points and had the same amplitude of Gaussian white noise and modulation period. The selected values represent the range of our observed data in each observation block of the target source. We then measured the SF for each of the simulated light curves and estimated the characteristic timescale as well as its uncertainty from the Gaussian fitting of the Gaussian probability distribution. Gaussian fitting of the Gaussian probability distribution was undertaken 1000 times through Monte Carlo simulations. Figure 6 shows 5 of these simulated light curves, the measured SF and the Gaussian fitting of the Gaussian probability distribution. We measured the peak of Gaussian fitting for each simulated light curves data to be 0.16, 0.17, 0.17, 0.22 and 0.24 days, for the light curves with 991, 772, 496, 234 and 178 data points, respectively. The uncertainties were 0.03 days for simulated light curves with 496 or more points and 0.02 for simulated light curves with 234 or 178 points. Our results suggest that a lack of data points located at the peaks and troughs of a light curve significantly changes the SF pattern, increasing the estimated characteristic timescale and decreasing its uncertainty by up to 30$\%$. Most of the light curves from the Ceduna dataset were similar in appearance to the 496 point simulated data sets; relatively well sampled but lacking some data at light curve peaks and troughs. This implies that the estimated characteristic timescale of Ceduna datasets are likely to be slightly overestimated compared to a fully sampled light curve.

Temporal variations in the flux density are caused by the relative transverse velocity of the scintillation pattern and the observer, $\textbf{{V}}_{\textnormal{rel}}$. For an isotropic scattering medium, the characteristic timescale that is observed (as a function of day of year (DOY)) is determined by $\textbf{{V}}_{\textnormal{rel}}$ and the characteristic spatial scale of the ISS pattern, $s_{\textnormal{o}}$:

The relative transverse velocity, $\textbf{{V}}_{\textnormal{rel}}$, has components that are both parallel and perpendicular to the equatorial plane, $\textbf{{V}}_{\alpha}$ and $\textbf{{V}}_{\delta}$. These components are in the directions $\hat{\textbf{{e}}}_{\alpha}$ and $\hat{\textbf{{e}}}_{\delta}$ respectively, where:

The relative velocity of the Earth and the scattering screen in the ISM, $\textbf{{V}}_{\textnormal{Earth-ISM}}$ depends on the relative velocities of the Earth and the Sun, $\textbf{{V}}_{\textnormal{Earth-Sun}}$ (a function of the DOY due to the orbital motion of the Earth), as well as the Sun and the LSR $\textbf{{V}}_{\textnormal{Sun-LSR}}$:

When projected onto the plane perpendicular to the line of sight, equation 7 gives the transverse relative velocity of the scintillation pattern:

The predicted annual cycle depends on the coordinates of the source ($\hat{\textbf{{e}}}_{\textnormal{source}}$). Figure 7 shows the predicted timescale as a function of DOY for PKS B1144$-$379, assuming an isotropic scattering screen with a variable velocity offset compared to the LSR. For a scattering screen moving with the LSR, the relative transverse velocity drops to almost zero near DOY 210, which result in a rapid increase in variability timescale (a near “stand-still” in the variations) around that time of the year. There are 3 fitted parameters obtained from the isotropic annual cycle model; the characteristic spatial scale of the phase variations on the scattering screen ($s_{\textnormal{o}}$) and both components of relative transverse velocity of the scintillation pattern and the observer ($\textbf{{V}}_{\alpha}$ and $\textbf{{V}}_{\delta}$). A non-linear least squares fit that assumed an isotropic (but non-stationary compared to the LSR) scattering screen was performed. We then implemented the Levenberg-Marquadt technique for non-linear least squares regression. Initially, each variability timescale data point was weighted with the inverse of the corresponding estimated uncertainty. However, as explained in §4, we also fit the data with uniform weighting for each data point. The $\chi^{2}$ values were calculated to quantify the accuracy of the model using the equations:

We take the square of the difference between the timescale estimate, $y_{\textnormal{i}}$, and the expected timescale according to the model, $f_{\textnormal{i}}$, scaled by the uncertainty in the measured timescale, $\sigma_{\textnormal{i}}$, for all $N$ timescale estimates. The reduced chi-squared value, $\chi_{\textnormal{Red}}^{2}$, takes into account the number of degrees of freedom (here m is the number of free parameters in the model). For the isotropic scattering screen case, $\textit{m}=3$.

Next, an anisotropic model of the scattering screen in the ISM was considered. An anisotropic interference pattern is assumed to be in the form of ellipse with axial ratio R, and its major axis inclined at angle $\beta$ to the direction of the relative transverse velocity (Carter, 2008). Figure 8 shows the predicted timescale as a function of DOY for PKS B1144$-$379, assuming an anisotropic scattering screen with a variable velocity offset compared to the LSR. Similar to the isotropic model, the anisotropic model also predicts a sudden increase in the variability timescale around DOY of 200. The predicted timescale as a function of the DOY is:

Several tests to examine the reliability of the annual cycle fitting for both isotropic and anisotropic models were also performed using Monte Carlo simulations. The first test was to examine the effect of the distribution of data points on the annual cycle fitting (variations of the characteristic timescales). In this simulation, we sampled an annual cycle curve at 55, evenly spaced, points over the year. The characteristic timescales ranged from 0.01–1.60 days. The longest timescales were around DOY of 200 and 248 for isotropic and anisotropic fitting models. We then added noise to the sampled characteristic timescales. For each simulation the amplitude of the random noise was fixed at a value which was varied between $\pm$0.1 and 0.5 day in steps of 0.1 days. Average annual cycle fittings were determined through 1000 iterations of Monte Carlo simulations for each noise amplitude. Figure 9 shows the average annual cycle fittings for both isotropic and anisotropic models. We found that for both isotropic and anisotropic annual cycle models, the estimate of the slowest period in the annual cycle was in error by 6$\%$. However, the anisotropic annual cycle model, with more free parameters, was much more susceptible to the noise in the characteristic timescale estimates for random noise more than 0.5 days. These results indicate that the uncertainty in the characteristic timescale estimates plays a major role in the accuracy of the annual cycle fitting. We also investigated how the accuracy of the annual cycle fitting depended on the number of characteristic timescale measurements. Five sets of data of the characteristic timescale measurements with constant values of uncertainties were created (Fig. 10). We then determined the average annual cycle fitting for each of the dataset through 1000 iterations of Monte Carlo simulations. Figure 10 illustrates the average annual cycle fitting for both isotropic and anisotropic models for all five simulated datasets. For the isotropic model, the annual cycle fits were very similar and the estimate of the slowest period of the annual cycle was only 4$\%$ off as the number of timescale estimates was decreased. On the other hand, the anisotropic model showed much greater sensitivity, with the annual cycle fittings changing remarkably. The estimated slowest period of the annual cycle had a 19$\%$ error as we reduced the timescale values. Our results suggest that in order to accurately fit the annual cycle, at least 20 independent measurements of characteristic timescale covering the slow-down period are required. Estimation of fitted parameters from the annual cycle fitting can be improved by more continuous monitoring.

In Table 1, we summarize results of the observations of PKS B1144$-$379 taken as part of the COSMIC observing project between 2003 to 2011. The mean flux density of PKS B1144$-$379 at 6.7 GHz over the period 2003–2011 was found to vary by a factor of 4. The long-term monitoring of the flux density of this source shows that there were two flares, which commenced in 2005 November and 2008 August. Further investigation of the long term variability and high resolution imaging is part of a future publication. The characteristic timescale of the source (equivalent to $T_{\mbox{period}}/4$) was measured to vary from 0.22–3.37 day, with the longest timescale occurring during the first flare (2005 November). We found the modulation index of PKS B1144$-$379 ranged from 5–18 $\%$.

Figure 11 shows the variability timescale versus day of the year of PKS B1144$-$379 from 2003–2011 fitted with both isotropic and anisotropic annual cycle models. We also performed a comparison between weighted and unweighted fitting. The annual cycle in the variability timescale was found to be more prominent when the fitting does not take into account the estimated uncertainties. The estimated uncertainty in the timescale depends on the duration of the observing session compared to the characteristic timescale. The majority of our observing sessions were 10–15 days in duration, so the result is that the sessions corresponding to longer characteristic timescales (slow-down periods) have larger estimated uncertainties and hence have lesser weight in the weighting fitting process (i.e. the fit is down-weighting the slow-down period). Our Monte Carlo simulation shows that the slow-down periods are most important to constrain the parameters of the annual cycle as the absence of these crucial data points increases the fitted value for the spatial scale of the screen by a factor of 6, inflating the $\textbf{{V}}_{\alpha}$, $\textbf{{V}}_{\delta}$ and R to unrealistic values and decreases the $\beta$ by a factor of 2. Therefore, in this paper, we chose to determine the fitted parameters from the unweighted annual cycle fitting. Although for some years there are fewer data points than for others, the plots show consistent peaks of longer characteristic timescales (slow scintillation) around DOY 50–100 and 300–350. Shorter timescales (fast scintillation) are observed between DOY 150–250. These patterns are seen also by combining the data in 3 year groupings and combining the entire 9 years of observations (Fig. 12). The persistent detection of an annual modulation pattern is a strong evidence that interstellar scintillation is the cause of the intraday variability observed in PKS B1144$-$379.

The observed pattern in the timescale for PKS B1144$-$379 is quite different to those shown in Figures 7 and 8, for isotropic and anisotropic scintillation models with screens moving at velocities close to the local standard of rest (§3). Those models predict “slow scintillation” to be observed around DOY 210. In contrast, the observed variability timescales were found to be shortest during this particular period. We have examined the reduced $\chi^{2}$ values for both the isotropic and anisotropic models for each of the years. The reduced $\chi^{2}$ values for the isotropic fit vary by a factor of 4 (excluding 2009 and 2011, where we have insufficient sampling of the change in characteristic timescale over the year to obtain a meaningful fit). For the anisotropic fit the reduced $\chi^{2}$ varies by a factor of 6. For most years both the isotropic and anisotropic fits have reduced $\chi^{2}$ greater than 2, demonstrating that the model does not fit the data particularly well. The most likely reason for the relatively poor fit of the annual cycle models is that the source structure evolves over the year as indicated by changes in both source flux density and structure index (Shabala et al., 2014). Further evidence for this is that the highest reduced $\chi^{2}$ values are seen for 2005 and 2008, the two years corresponding to the outbursts. It should be noted that the longest timescale observed over the 8 years covered in the current work coincides with the start of the 2005 November flare. At this point, the modulation index is large in addition to the timescale being long. This particular one period might be just an anomalously large "scintle" (large and slow modulation) due to the stochastic nature of the scattering. Another possibility is that the scintillating component has both significantly brightened (compact fraction is large, hence large modulation) and expanded (longer timescale).

In Table 2, we summarize results of the fitted parameters for both isotropic and anisotropic annual cycle models for PKS B1144$-$379. We list the fitted parameters for each year from 2003 through 2011, inclusive. As previously mentioned, due to the larger formal uncertainty when the source exhibits longer characteristic timescales, we chose to determine the fitted parameters from the unweighted annual cycle fitting. For the year 2011, the anisotropic model could not be fitted due to insufficient data. For some years, the fitted parameters are not at all well constrained (larger by a factor of $>>$ 100), whereas for other years (2006 and 2007) they appear to be reasonably well constrained (by a factor of 3). The large range of fitted parameters is due to some unrepresentative outliers and small amount of data points. We also found that for some observing years (2005 and 2010), the screen velocity has unrealistic values and the R has gone to zero (or infinity), demonstrating that the fitted parameters were not well constrained, likely due to intrinsic source evolution which is entangled with the interstellar scintillation. Both these periods coincide with new flares which are likely to correspond to ejection of new plasma components. The emergence of new components during the source evolution leads to lengthening of the characteristic timescales. This will cause ambiguity in the fitted parameters as the shape of the slow down period is critical in determining the ISM velocity and anisotropy in the scintillation pattern (Bignall et al., 2003). The lower reduced $\chi^{2}$ for the anisotropic models shows that it is able to provide a better fit to the data.

We investigated the effect of fixing some of the model parameters, as for some of the years the fitted values for the screen velocity and the axial ratio R of the scattering screen produce unrealistic values. We chose to fix $\textbf{{V}}_{\alpha}$ and $\textbf{{V}}_{\delta}$ to be 30 and 7.7 km s^-1, respectively. These values come from the fitted parameters of observing year 2006 and 2007 where they appear to be reasonably well constrained. A persistent slow-down period near DOY 50 and towards the end of the year (Fig. 11 and Fig. 12) is a strong indication that there is a repeating annual cycle in the variability timescale of PKS B1144$-$379. Therefore, it seems reasonable to assume that the velocity and other properties of the scattering screen remain unchanged over the entire observing period. Figure 13 shows the variability timescales versus day of the year for PKS B1144$-$379 from 2003 through 2011 fitted with fix values of $\textbf{{V}}_{\alpha}$ (30 km$s^{-1}$) and $\textbf{{V}}_{\delta}$ (7.7 km$s^{-1}$). Figure 14 displays plots of the variability timescale combining data in three year groups and for the entire 9 years of observations. Both Figure 13 and Figure 14 show the unweighted annual cycle fittings. Overall, we observed an improvement in the annual cycle fitting as the optimum fittings now pass through almost all the data points, especially at the peak of “slow-scintillation”. The fitted parameters values were also found to be better constrained. In Table 3, we summarize results of the fitted parameters for both unweighted isotropic and anisotropic annual cycle models of PKS B1144$-$379, with the screen velocity fixed. Overall, the fitted parameters appear to be reasonably well constrained. The 2011 data can now be fitted with the anisotropic annual cycle model as we only have 3 parameters that are allowed to freely vary ($s_{\textnormal{o}}$, R and $\beta$).

To investigate whether there are correlations between the radio IDV and long-term variations in the total flux density, we plotted the temporal changes in the mean flux density, modulation index and characteristic timescales ($t_{\textnormal{char}}$) against time (Fig. 15). No correlation was found between the mean flux density and the modulation index as well as the mean flux density and the characteristic timescale. Figure 16 shows the RMS of the flux density and the mean flux density on the same temporal scale and a clear correlation with correlation coefficient of 0.7 between the two is evident (refer to Figure 17). This demonstrates that the amplitude of the radio IDV phenomenon is larger when the mean flux density is high and hence that the radio IDV is entangled with the temporal evolution of the source. Higher values of the modulation index imply that the source becomes more compact and scintillates more, or that the fraction of the flux density in the compact component increases. Our results are in good agreement with Carter et al. (2009). They reported a loose correlation between the mean flux density and modulation index for PKS B1519$-$263 which led them to suggest a significant fraction of the change in the source mean flux density was due to evolution of the compact scintillating core. Bignall et al. (2003) reported a large increase in the RMS variations is expected if the increase in the total flux density is due to an increase in the flux density of the scintillating component. In addition, if the scintillating component expands, this will manifest as an increase in the characteristic timescale, which can also be detected (Liu et al., 2012).

The long-term variations in the mean flux density of PKS B1144$-$379 are most likely due to intrinsic source evolution. The amplitude of the short timescale variability due to ISS depends on the fraction of the source flux density in the scintillating components, the angular size of the source, the observed frequency, and the ratio of the observed frequency to the transition frequency between weak and strong scattering (Walker, 1998). PKS B1144$-$379 is located at Galactic coordinates of $289.24^{o}$ (longitude), $22.95^{o}$ (latitude) and redshift $z$ =1.048 (Kedziora-Chudczer et al., 1997, 2001b). For this line of sight, the transition frequency ($\nu_{\textnormal{o}}$) estimated by the Cordes & Lazio (2001) model of the Galactic electron distribution is 14.4 GHz. Thus, the Ceduna 6.7 GHz observing frequency is expected to be in the refractive scattering regime ($\nu$ < $\nu_{\textnormal{o}}$) (Turner et al., 2012).

In this paper, we adopted two techniques to determine the physical characteristics of the source; the “Earth Orbital Synthesis” approach which utilizes the fitted parameters from annual cycle fitting (Macquart & Jauncey, 2002) and a simple model of an extragalactic point source introduced by Walker (1998). Modelling using “Earth Orbital Synthesis” allows us to estimate the source angular size from the fitted minor axis length scale (i.e. $a_{\textnormal{min}}$ = L$\theta_{\textnormal{s}}$, where L is the distance of the scattering screen in kpc and $\theta_{\textnormal{s}}$ is the source angular size in microarcsec). We consider this relationship as we assumed that the source size is “resolved” by the ISM and it determines the minor axis length scale as changes of both the minor axis length scale and the modulation index appear to be related to the total flux density changes. Therefore, we would expect that the source is most of the time larger than the Fresnel scale regardless of the scattering strength, and larger than refractive scale if we are in strong scattering. Walker (1998) give the scaling of point source modulation index in the regime of refractive scintillation as:

Several possibilities of the screen distance (L) and strength of the scattering screen ($\xi$) were considered (refer to Table 4). In our annual cycle fitting, the minor axis scale of the scintillation pattern ($a_{\textnormal{min}}$ = $s_{\textnormal{o}}$ $\sqrt{R}$, where $s_{\textnormal{o}}$ is the characteristic spatial scale of the ISS pattern in km and R is the axial ratio) ranged from 0.48–1.95 ($\times 10^{6}$ km). If instead we use $a_{\textnormal{min}}$ taken from the modelling with the velocities of the screen fixed to be $\textbf{{V}}_{\alpha}$ = 30 km $s^{-1}$ and $\textbf{{V}}_{\delta}$ = 7.7 km $s^{-1}$, the $a_{\textnormal{min}}$ was found to range from 0.71–2.00 ($\times 10^{6}$ km). For better constraint of fitted parameters, here we only considered the $a_{\textnormal{min}}$ measured from the fixed velocity fits. If the scattering screen is at 0.84 kpc, the Fresnel scale is 3.37 $\mu$as, which is only slightly smaller than the smallest angular size from the annual cycle fitting ($\theta_{\textnormal{s}}$ = 5.65–15.90 $\mu$as). Assuming the source angular size corresponds to the refractive scale and the equation is valid close to the transition frequency, the scattering strength would be 1.68. This would suggest that the scattering is not very strong if the screen is at 0.84 kpc or beyond it. If we move the screen closer, at 500 pc, we have the smallest angular size of 9.49 $\mu$as ($\theta_{\textnormal{s}}$ = 9.49–26.73 $\mu$as), the Fresnel scale would be 4.37 $\mu$as while the scattering strength is 2.17. At screen distance of 100 pc, the smallest angular size is expected to be 47.45 $\mu$as ($\theta_{\textnormal{s}}$ = 47.45–133.66 $\mu$as), the Fresnel scale and the scattering strength would be 9.78 $\mu$as and 4.85, respectively. The screen distance is not well constrained, but a closer screen implies a lower source brightness temperature. We assumed that the screen distance is probably closer than 0.84 kpc or otherwise the scattering strength will be less than 1.68 which seems to be a lot weaker. On the other hand, if the screen distance is closer than 100 pc, this would imply a very strong scattering ($>$ 4.85) in a very nearby screen which is also not expected and the compact fraction of the source would possibly be higher than what we infer from the VLBI data. Variations of timescales longer than a day for PKS B1144$-$379 also indicate that more distant screen ($>$ 100 pc) is expected as measured by Turner et al. (2012). In addition, a handful of rapid IDV sources scintillating through nearby screens have been argued to be scintillating in weak scattering at frequencies as high as 6.7 GHz (e.g. Kedziora-Chudczer et al., 1997). For PKS B1144$-$379, the long timescale of variation and the annual cycle fit constrains the characteristic length scale of the scintillation pattern to $>$ $10^{6}$ km. For the more rapid scintillators that have had an annual cycle measured (PKS 1257$-$326 and J1819$+$3845), the scintillation pattern scale is more than the order of magnitude smaller, and the scattering screen is inferred to be within tens of parsec. If PKS B1144$-$379 is scintillating through such a nearby screen, the scintillation pattern scale must be much larger than the Fresnel scale, implying either the source angular size is much larger or the scattering is very strong (leading to a larger refractive length scale). Such strong scattering on the nearby screen is not expected (or at least would have been unprecedented). Figure 18 shows the source angular size estimated using the fitted parameters of the annual cycle model when we constrained the values of both velocities component of scattering screen, $\textbf{{V}}_{\alpha}$ (30 km $s^{-1}$) and $\textbf{{V}}_{\delta}$ (7.7 km $s^{-1}$) for the years 2003–2011. Here we used 0.84 kpc, the upper limit screen distance to measure the source angular size. The source angular size of PKS B1144$-$379 determined using the “Earth Orbital Synthesis” ranged from 5.65–15.90 $\mu$as.

The simple model of an extragalactic point source introduced by Walker (1998) assumes that all the flux density is in the scintillating component. However, observations typically show that only a fraction of the total flux density is time variable with some of the emission arising in a more extended milliarcsecond-scale jet. Furthermore, the source angular size determined from this approach is derived using only the modulation index and the assumption that the transition frequency is 14.4 GHz. This approach does not utilise information obtained from the modelling of the annual cycle in the modulation timescale (known as “Earth Orbital Synthesis” ), which also yields information on the source angular size. The modulation index has been used by various authors to characterise the linear scale of the region responsible for the intraday variability (Quirrenbach et al., 1992; Peng et al., 2000; Kedziora-Chudczer et al., 2001b). Here, we altered the equation of refractive scintillation given by Walker (1998) and to include the compact fraction of the source ($f_{\textnormal{c}}$ = $S_{\textnormal{core}}$/$S_{\textnormal{total}}$) as:

In addition to the COSMIC dataset presented in this paper, information on the physical characteristics of the source can also be determined from VLBI radio images. Figure 20 compares the source angular size estimated using the “Earth Orbital Synthesis” of the COSMIC data at 6.7 GHz and major axis of core component measured from the 8.6 GHz VLBI radio images. The VLBI radio images were obtained from The Radio Reference Frame Image Database (RRFID) (https://www.usno.navy.mil/USNO/astrometry/vlbi-products/rrfid) and the Astrogeo VLBI FITS image database (http://astrogeo.org/vlbi-images). The upper limit of the VLBI core size was measured from an elliptical Gaussian model fit to the image data made using the DIFMAP software package (Shepherd, 1997). In the period 2001 to 2013, the major axis of the core component was estimated to vary between 0.16 and 0.81 mas. We found quite similar pattern between these two datasets with the source angular size estimated from the COSMIC data modulation index lagging the value obtained from the major axis of the core component of VLBI radio images. The core appears to expand first at the scales of ISS, and then later at VLBI scales implying a new component moving out in the jet. Figure 21 compares the flux density of both the compact and extended emission regions for PKS B1144$-$379 estimated through two independent techniques. From the COSMIC data we can measure the flux density of both the scintillating and non-scintillating (extended) components, while from the VLBI images we have the flux density of the core component compared to the total flux density of the source.

The modulation index of the scintillating component ($m_{\textnormal{scint}}$) is expected to be:

As previously mentioned, a closer screen implies a lower brightness temperature. We deduced the brightness temperature and inferred Doppler factor are extremely high if the screen is at 0.84 kpc. Table 4 shows several possibilities of the screen distance (L) that we considered to estimate the brightness temperature and the Doppler factor. Equation 17 used to measure the flux density of scintillating component and we assumed the source angular size corresponds to the refractive scale. If the screen is at 0.84 kpc, the brightness temperature would range from $9.48\times 10^{13}$ to $2.74\times 10^{14}$ K and this would suggest Doppler factor of $\sim$ 600. If we move the screen closer, at 500 pc, the brightness temperature could be lowered down to 3.36–9.72 ($10^{13}$ K) while the Doppler factor is expected to be 253. At a screen distance of 100 pc, the brightness temperature and the Doppler factor would be 1.34–3.89 ($10^{12}$ K) and 17, respectively. A reasonable value of brightness temperature and inferred Doppler factor imply that the screen distance of PKS B1144$-$379 is probably closer than 0.84 kpc.