Free Access
Issue
A&A
Volume 542, June 2012
Article Number A121
Number of page(s) 9
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201218977
Published online 19 June 2012

© ESO, 2012

1. Introduction

Radio variability on timescales from days to weeks was observed in AGNs as early as 1971 (see, e.g., Kinman & Conklin 1971; Wills 1971). In 1986, Witzel et al. (1986) reported the discovery of flux density variations occurring on timescales shorter than one day, which have been referred to as intraday variability (IDV). It is known that IDV affects 20% to 50% of the flat-spectrum radio source population (see Quirrenbach et al. 1992; Lovell et al. 2008). The origin of the variability, however, is not completely understood. Assuming that the variability is intrinsic to the source, causality arguments constrain the maximum size of the emission region using the observed variability timescales. For IDV sources, this frequently implies very high brightness temperatures, exceeding by far the inverse-Compton limit (1012 K; see Kellermann & Pauliny-Toth 1969).

Propagation effects, such as interstellar scintillation (ISS), were suggested as a likely alternative scenario. In this model, the variability is caused by the scintillation of very compact components (of the size of tens of μas) in the IDV source by a scattering screen placed between the source and the observer. The radiation can be either focused or defocused by the screen, causing the measured flux density to be a function of the line of sight to the observer. Owing to the Earth’s revolution around the Sun, the observer’s line of sight toward the screen changes continuously, causing the detected rapid variations in the flux density measurements. Consequently, the faster the Earth’s movement with respect to the scattering screen, the shorter the variability timescale. Since the Earth’s velocity follows a one-year periodic cycle, the relative velocity between the scattering screen and the observer should change accordingly, resulting in a seasonal cycle of the variability timescale.

The detection of annual modulation patterns has been essential in demonstrating the source-extrinsic origin of IDV for the most extreme sources (timescales of hours, amplitudes of  ~100%), the so-called fast scintillators, namely PKS 0405-385, J1819+384, and PKS 1257-326 (see Kedziora-Chudczer et al. 1997; Dennett-Thorpe & de Bruyn 2000; Bignall et al. 2003). For type-II IDV sources (characterized by less extreme variability, with amplitudes of  ~1–20% and timescales  ≲ 2 days, see, e.g., Quirrenbach et al. 1992), the situation is less clear. It may partially depend on their relatively long variability timescales, which make seasonal cycles harder to identify, given the unrealistically long observing time that would be required to measure accurate timescale estimates. Seasonal cycles have been detected in J1128+592 (see Gabányi et al. 2007) and claimed for other sources, such as S4 0917+62 (see Rickett et al. 2001). The MASIV survey (Lovell et al. 2003, 2008) showed that a significant part of the flux density variations observed in compact radio sources are due to ISS. On the other hand, the discovery of correlated optical-radio IDV in the flux density of the type-II source S5 0716+71 suggested that its origin may be intrinsic to the source (see Quirrenbach et al. 1991).

The detailed analysis of the variability characteristics of IDV sources plays an essential role in understanding the ISS contribution to their variability pattern. From 2005 to the end of 2009, the Max-Planck-Institut für Radioastronomie and the Urumqi Observatory carried out a program of intensive monitoring of type-II IDV sources with the Nanshan 25 m radio telescope, located about 70 km south of Urumqi city (China). The main aims of the project are: (i) monitoring the long-term IDV characteristics of a larger sample of IDV sources, looking for seasonal cycles in their variability patterns, and (ii) detecting new IDV sources. The main target of the monitoring program were the type-II sources AO 0235+164, S5 0716+71, S4 0917+62, S4 0954+65, J1128+592, and HB89 1156+295. During the observations, about sixty additional sources were sporadically observed while attempting to find new IDV candidates.

In the present paper, we focus on our results obtained for S4 0954+65, which is a well-studied BLLac object at redshift z = 0.368 whose IDV features have been investigated since the early ’90s (Wagner et al. 1993). The origin of the variability is still controversial – S4 0954+65 is the second source for which a radio-optical correlation has been claimed (Wagner et al. 1990). However, the detection of extreme scattering events (Fiedler et al. 1987, 1994; Cimò et al. 2002) and the possible presence of a seasonal cycle in the variability timescales of the source, which was revealed by a previous IDV monitoring project (see Fuhrmann 2004), seem to indicate that at least part of the variability is caused by extrinsic mechanisms, and that more than one scattering screen may be interposed between the source and the Earth.

In the next section, we present an overview of the observations and the data calibration procedure. The analysis methods used to extract the variability characteristics of the light curves are presented in Sect. 3. The main results are shown in Sect. 4. In Sects. 5 and 6, we discuss our results and present a summary of the paper, respectively.

2. Observation and data calibration

The Nanshan 25 m telescope is equipped with a 4.8 GHz single-horn dual-polarization receiver provided, along with the new telescope driving program, by the MPIfR (for more details, see Sun et al. 2006). The receiver’s bandwidth is 600 MHz. The flux density measurements were performed in cross-scan mode. Each scan consists of eight sub-scans in perpendicular directions across the source position in order to estimate the pointing offsets in the two scanning directions and correct for them. The raw data were processed at the MPIfR via TOOLBOX, a Python-based software package that subtracts baseline drifts and fits a Gaussian profile to each cross-scan. The amplitude of the profile provides an estimate of the source’s flux density expressed in units of antenna temperature (K).

The data calibration was performed in semi-automatic mode (for more details, see Marchili 2009; Marchili et al. 2010). The quality of the individual sub-scans was checked by considering their pointing offset, the full-width at half-maximum (FWHM) of their Gaussian profile, and the error in their flux density measurement. Sub-scans that showed pointing offsets of more than 100′′, deviations from the average FWHM (i.e.  ~600′′) larger than 150′′, or uncertainties in the Gaussian fit estimate larger than 10% of the amplitude were automatically discarded. Subsequently, a weighted average of the remaining sub-scans was obtained with a weight proportional to the inverse of the squared errors, separately for the two scanning directions. After the pointing-error correction, we averaged the mean sub-scan amplitudes in the two scanning directions, providing a single flux-density measurement per scan.

In each observing session, about 50% of the time was dedicated to measurements of primary and secondary calibrators, i.e., compact radio-sources with steep spectra and no evidence of IDV. The 5 GHz flux densities of the so-called primary calibrators are fairly constant in time and well-known. The secondary calibrators, instead, show variability on timescales of months. We used the primary and secondary calibrators to estimate the variations in the gain of the antenna due to changes in both its parabolic shape and the position of the focus. We used the secondary calibrators – which were selected to be at small angular separations from the target sources – to correct for possible changes in the collected flux densities caused by instabilities in the antenna-receiver system or changing weather conditions. The data calibration was completed by the conversion of measured antenna temperatures (in K) to flux densities (in Jy). The time-interpolated conversion factors were derived by means of the known flux densities of the primary calibrators, namely NGC 7027, 3C 286, and 3C 48. The complete data-calibration procedure guarantees a high level of accuracy, which, based on the residual variability in the calibrators, can be estimated in the range between 0.2 and 0.8%, under normal weather conditions.

thumbnail Fig. 1

Variability curves of S4 0954+65 (black dots) in April 2006 (upper panel) and April 2008 (lower panel). The turquoise squares show the residual variability in the fluxes of the secondary calibrators S5 0836+71 and S4 0951+69.

During the monitoring campaign, between August 2005 and December 2009, 41 IDV datasets for S4 0954+65 were collected. Two examples of light curves are shown in Fig. 1. In three cases, the observing sessions were too short (≤1.7 d) to provide a reasonable estimation of the variability characteristics of the sources, were consequently discarded from the analysis we describe in the following. An overview of the main characteristics of the remaining 38 sessions is reported in Table 1. In the few cases in which the duration of the observations is longer than 5 days or shorter than 2.5 days, a bias of the estimated timescales toward, respectively, the slowest and the fastest variability components cannot be excluded.

Table 1

Basic information about the 38 observing sessions at the Urumqi Observatory during which S4 0954+65 was observed.

3. Variability analysis

Following the scheme proposed by Quirrenbach et al. (2000) and Kraus et al. (2003), two different quantities were used to evaluate the significance and amplitude of the variability, namely the reduced chi-squared test statistic χr2 and the variability amplitude Y.

The reduced chi-square (see, e.g., Bevington & Robinson 1992) was used to test for given datasets the likelihood of the hypothesis of a constant flux density χr2=1N1i=1N(SiSδSi)2,\begin{equation} {\chi_{\mathrm{r}}}^2=\frac{1}{N-1} \sum_{i=1}^N{\left(\frac{S_{{i}}-\langle S\rangle}{\delta S_{{i}}}\right)}^2, \end{equation}(1)where N is the number of datapoints,  ⟨ S ⟩  the average flux density, Si the ith flux density measurement, and δSi its error. A dataset was regarded as variable if the χr2 test returned a probability of constant flux below 0.1%, which is equivalent to a 99.9% probability for significant variability.

An estimation of the variability in a dataset was provided by the modulation index, mi [%]  = 100   σS/ ⟨ S ⟩ , where σS is the root mean square (rms) in flux density. However, this parameter does not take into account the residual noise in the calibrated data, which varies from epoch to epoch, mainly depending on weather conditions. We defined the variability amplitude Y as an equivalent to a 3-σ confidence interval for Gaussian noise, modified by the quadratic subtraction of the rms noise floor, m0, to take into account the noise bias Y=3mi2m02.\begin{equation} Y=3\sqrt{m_{{i}}^2-m_{\mathrm{0}}^2}. \label{eq:y} \end{equation}(2)

3.1. Variability timescale

The characteristic timescales of the variable datasets were estimated using two different methods of time series analysis – a first-order structure function (see Simonetti et al. 1985) and a wavelet-based analysis.

The structure function of a given time series  {S(ti)} i is calculated as SF(τ)=1nij[S(ti)S(tj)]2,\begin{equation} {\it SF}(\tau)=\frac{1}{n}\sum_{ij}[S(t_i)-S(t_j)]^2, \end{equation}(3)where the sum is extended to the n pairs (ti,tj) for which  ||titj ||  ≃ τ. Ideally, when a structure function is applied to a red-noise-like signal, it shows a monotonic increase that can be described as a power law. At the time lag τsf corresponding to the maximum coherency time of the signal, the structure function shows a plateau, for which τsf provides an estimate of the characteristic timescale of the signal. Owing to both the uneven sampling and the noise in the analyzed data, the position of the plateau is subject to some degree of uncertainty. For each plateau, it is possible to determine a minimum and a maximum timescale – respectively τmin and τmax. The characteristic timescale of the structure function τsf was calculated to be (τmax + τmin)/2. The uncertainty in the timescale estimation was evaluated as the maximum between two terms, namely (τmaxτmin)/2 (which reflects the significance of the plateau) and 0.35(τsf3/2)/(obsd)1/2 (see Marchili et al. 2011), where obsd indicates the duration of the observation. The latter term accounts for the significance of the timescale, which is limited by the small number of cycles that it is possible to observe, given the limited duration of the observations. To estimate the proportionality factor, we calculated all the time-intervals between a local minimum (or maximum) and the following maximum (or minimum) for a set of light curves characterized by rapid variability; for each light curve, we calculated the average of the observed time-intervals (which provides an estimate of the timescale that is consistent with τsf) and their standard deviation (which provides an estimate of its uncertainty). A linear regression of the uncertainties, plotted against the values of (τsf3/2)/(obsd)1/2 for the corresponding light curves, returned the proportionality factor 0.35.

The structure function in some cases was found to have more than one plateau, which is indicative of multiple variability timescales. In these cases, we identified the characteristic timescale with the shortest one.

The wavelet-based algorithm that we developed is based on the Ricker (“Mexican hat”) mother wavelet Ψ(t)=C(1t2)exp(t22)\begin{equation} \Psi(t)=C\,\left(1-t^2\right)\exp\left(-\frac{t^2}{2}\right) \end{equation}(4)where C is a normalization factor and t is the time parameter. An extensive description of this analysis method can be found in Appendix A. The timescales obtained by applying the wavelet analysis are indicated with τwt.

Examples of structure functions and wavelet results are provided in Fig. 2 for the light curves shown in Fig. 1.

thumbnail Fig. 2

Upper panels: the structure function plots of the April 2006 (left panel) and April 2008 (right panel) datasets; the blue vertical lines show the estimated timescales. Lower panels: the wavelet-based analysis for the same epochs; the green vertical lines show the estimated timescales. The existence of two bumps in the right panels is the signature of two distinct variability components in the dataset.

4. Results

The variability characteristics of S4 0954+65 are summarized in Table 2. Our study had two main outcomes, which are discussed below: (i) we found strong evidence of an abrupt change in the variability characteristics of S4 0954+65 between February and March 2008; (ii) the IDV timescales of the source seem to undergo an annual modulation.

Table 2

Variability characteristics of S4 0954+65 in 38 epochs of observation with the Urumqi Telescope.

4.1. Change in the overall variability characteristics

In the long-term light curve of the source (Fig. 3, upper panel), two considerably different phases can be recognized. From August 2005 up to February 2008, the flux density variations were generally moderate in intensity (the average peak-to-peak variation is  ~0.2 Jy) and occurred on timescales of three months; March 2008 marked the starting point of a strong outburst phase, which up to December 2009 was still ongoing, characterized by large variations (~0.4 Jy) on a timescale of about seven months. For simplicity, we refer to the first and the second timespan, respectively, as the pre-outburst and the outburst phase of S4 0954+65.

thumbnail Fig. 3

Long-term variability of S4 0954+65 (upper panel). In the middle panel, the wavelet timescales (red dots) and a three-point running average (blue line), showing the sudden increase starting from March 2008. The structure function timescales are shown in the lower panel. The red vertical line divides the observing sessions into two phases characterized by remarkably different variability characteristics.

The change in the long-term activity of the source was also reflected in the remarkable variation in its IDV characteristics. There are three lines of evidence for such a variation:

  • During the pre-outburst phase of the source, the IDV timescalesprovided by the wavelet-based method and the structurefunction (see Fig. 3, middle and lower panels)were highly correlated with each other. A linear regressionapplied to the timescales (Fig. 4, left panel)results in a correlation coefficient of 0.82 and a slope of 1.05,which proves the consistency in the definition of timescalesused for the two methods. During the outburst phase, a simi-lar analysis provides a correlation coefficient of –0.11 witha slope of –0.19 (Fig. 4, right panel), imply-ing that there is a clear lack of correlation. It is useful to apply alinear regression to a comparison sample comprising all thetimescales estimated for the other main sources in the Urumqimonitoring program (i.e., AO 0235+164,S5 0716+71, S4 0917+62,J1128+592, and HB89 1156+295). This re-turns a correlation coefficient of 0.86 with a slope of 1.01,which proves the high degree of compatibility between thestructure function and the wavelet-based method results.This demonstrates that the variability characteristics of thesource are atypical during the outburst phase, as illustrated bya Kolmogorov-Smirnov test for the comparison of datasets.There is a 74% probability that the distribution of the values ofτsf/τwt during the pre-outburst phase is equal to the one estimated for the comparison sample. A similar test for the τsf/τwt values during the outburst phase returns a probability below 0.1%.

  • The wavelet timescales during the pre-outburst phase of S4 0954+65 follow a slowly decreasing trend, pinpointed by the three-point running average in the middle panel of Fig. 3 (blue line). A sudden increase (from  ~0.5 to  ~1.5 d) occurs in March 2008, after which the estimated timescales remain considerably longer than before. It is again helpful to apply a Kolmogorov-Smirnov test when comparing datasets; the probability that the two distributions of timescales are equal is below 1%.

  • We calculated an average structure function from the light curves obtained during the pre-outburst phase, and compared it with the average one for the outburst phase1. The difference between the two can be seen in Fig. 5. In the structure function of the pre-outburst phase, there is considerably more power on short timescales than in the structure function of the outburst phase. On longer timescales, we find the opposite situation, suggesting that there is some source of slower variability that is more dominant during the outburst phase.

Possible explanations of the large changes that occur in the IDV characteristics of S4 0954+65 at the beginning of its outburst phase, are discussed in Sect. 5.

thumbnail Fig. 4

Comparison between the structure function and wavelet results for the comparison sample (black dots) and the S4 0954+65 light curves. In the left panel, the source’s timescales estimated up to February 2008 (cyan squares); in the right panel, the timescales from March 2008 to December 2009 (green squares). The black, cyan, and green lines show the results of applying linear regressions to the respective datasets.

thumbnail Fig. 5

Average structure function calculated from the normalized light curves obtained during the pre-outburst phase (turquoise diamonds) and the outburst phase (red squares).

4.2. Annual modulation of the timescales

We investigated the possible annual cycle in the variability timescales of S4 0954+65 by analyzing the evolution of the timescale τ versus day of the year (DoY), separately for the wavelet-based and the structure function results. Following the scheme proposed by Qian & Zhang (2001), updated to the case of anisotropic scattering (see, e.g., Bignall et al. 2006; Gabányi et al. 2007), we expressed τ in terms of the vector s – which defines the orientation of the elliptical scintillation pattern –, the relative velocity between the scattering screen and the observer, v, the distance to the screen, D, and the anisotropy factor, rτ(DoY)Drv2(DoY)+(r21)(v(DoY)×s)2·\begin{equation} \tau(\mathrm{DoY}) \propto \frac{{D} \sqrt{{r}}}{\sqrt{{{v}^2(\mathrm{DoY})+({r}^2-1)\,({{\vec v}(\mathrm{DoY}) \times {\vec s}})^2}}}\cdot \end{equation}(5)We note that the anisotropy may either be induced by an anisotropic scattering medium or be intrinsic to the source (caused by an anisotropic emitting component).

We developed a code, based on a least squares fitting algorithm, to calculate the scintillation pattern parameters that most closely describe the detected variability timescales (for details, see Marchili 2009). We note that the distance cannot be estimated without making “ad hoc” assumptions about the effective angular size of the source; in the following, we adopt the model proposed by Qian & Zhang (2001), which assumes an intrinsic angular size of the scintillating component of 70 μas. In Sect. 5.2, we modify this size estimate, to obtain an improved fit to the detected variation in the IDV pattern.

Table 3

Scintillation pattern parameters deduced by the distribution of the timescales as a function of the day of the year.

thumbnail Fig. 6

Annual modulation plots based upon the wavelet results (left panels), and the structure function results (right panels). The upper panels show all the estimated timescales, while lower panels show monthly averages. The green and turquoise lines show the annual modulation patterns which best fit the data.

thumbnail Fig. 7

Annual modulation plots, as in Fig. 6; in the upper panels, the timescales obtained for the epochs between August 2005 and February 2008 are highlighted in magenta squares. In the middle and lower panels are displayed the annual modulation plots (green and turquoise lines) and the 40-day timescales averages for the August 2005–February 2008 and the March 2008–December 2009 epochs. The red vertical lines indicate the peaks in the annual modulation patterns obtained from the complete set of timescales.

The plots of the timescales versus DoY – the so-called annual modulation plots –, comprising all the estimated timescales between August 2005 and December 2009, are shown in Fig. 6, separately for the results provided by the two analysis methods. The wavelet results (left panels in the figure) suggest that two slow-down phases occur during the year, centered on May (DoY 140) and November (DoY 312), corresponding to a cycle of approximately six months. The structure function timescales (right panels in the figure) show a less clear cycle, with only two data-points being considerably higher than the others, corresponding to the observing sessions of December 2005 and November 2006. Nevertheless, they appear to confirm that there is an annual cycle in the timescales that has a prominent slow-down phase peaking in November (DoY 320), and a much weaker one peaking in May (DoY 136) – in fair agreement with the wavelet results. The best-fit scintillation pattern parameters are reported in Table 3. The observed annual modulation cycle is compatible with the presence of a scattering screen at a distance between 180 pc and 270 pc. The scintillation pattern is slightly anisotropic, with an anisotropy degree around 2 and an anisotropy angle around 90 deg. Were we to assume that the scattering screen is associated with the Ursa Major cloud complex, the estimated distance would be compatible with the value of 240 pc adopted by Heithausen et al. (1993), but less consistent with the distance of only 100–120 pc provided by Penprase (1993) for the molecular clouds MBM 29–MBM 31 in the complex. A distance to the screen of between 120 and 290 pc can also be deduced from the Hipparcos catalog (Perryman & ESA 1997; van Leeuwen 2007), by looking at the parallaxes of the stellar objects within a two-degree circle around S4 0954+658.

5. Discussion

The existence of a break in the variability characteristics of the source, described in Sect. 4, raises the question of whether this break affects the search for an annual variation in the characteristic timescales. The consequences of the March 2008 break can be evaluated by separately analyzing the variability timescales before and after it occurs. The corresponding annual modulation plots and scintillation pattern parameters are presented in Fig. 7 and Table 3 for the wavelet and structure functions separately. In the upper panels of Fig. 7, the August 2005 – February 2008 timescales are shown as magenta squares, while the March 2008 – December 2009 timescales are shown as black circles. In the middle and lower panels, we show 40-day averages of the timescales during the pre-outburst and the outburst phases, respectively. In all plots, the green and turquoise lines display the best annual modulation fits to the data. Apparently, the yearly cycles found before and after March 2008 in the wavelet-based results are quite similar, with slow-down phases at DoY 144 and 320 for the former and at DoY 137 and 309 for the latter. In the case of the structure function results, instead, the difference in the cycles detected in the two timespans is significant. To understand the nature of the change in the variability characteristics of S4 0954+65 between the pre-outburst and the outburst phases, it is essential to find the reason why the wavelet and the structure function results are affected differently.

5.1. Multiple timescales in the light curves after February 2008

The general increase in the wavelet timescales after February 2008 could be explained in terms of a general slow-down in the variability of the source. This, however, could not account for the systematic discrepancies between the results of the two analysis methods, highlighted in Fig. 4. These discrepancies are most likely the consequence of the appearance, starting from February 2008, of multiple variability contributions in the light curves, of different timescales and comparable strengths. The wavelet and structure function plots in the right panels of Fig. 2, which show the results for one of the light curves collected during the outburst phase, confirm this interpretation. The wavelet plot shows two bumps, corresponding to timescales of 2.2 and 0.5 days, which indicate that there are two distinct variability contributions; the two plateaus in the structure function plot, which have similar timescales, can be explained in the same way.

The discrepancies between the estimated timescales reveal the differences between the two analysis methods. The wavelet timescales are generally longer than the structure function ones, suggesting that the slowest variability contributions are the most powerful. For the structure function, it was our choice to identify the characteristic timescale with the one corresponding to the first significant plateau of the function. This rigid approach, which was forced by the necessity to minimize the arbitrariness of the definition of the characteristic timescales, should not interfere with the detection of possible annual modulation cycles. Assuming that a fast and a slow timescales are both detected in the majority of the light curves – which is the case for the observing sessions during the outburst phase – and that they are caused by the same scattering screen, they should follow the same annual modulation cycle. Therefore, theoretically, this can be identified by systematically looking at the fastest variability component – which is the one that can be estimated with the highest accuracy, given the increase in the uncertainty with lengthening timescale. In practice, however, the situation is more complicated; the large difference in the annual modulation plots of the structure function timescales before and after March 2008 demonstrates that a proper interpretation of the results cannot be achieved without simultaneously taking into account the information provided by both of the time analysis methods.

thumbnail Fig. 8

Wavelet (black dots) and structure function (red squares) timescales calculated during the outburst phase of S4 0954+658 are indicative of multiple variability contributions. In the left panel, the variability is interpreted as the superposition of a source-extrinsic effect (blue line) and a source-intrinsic process (green line). In the right panel, both the variability contributions are attributed to ISS and the two annual modulation patterns have the same scattering screen parameters, the only difference being the size of the scintillating components (90 μas for the blue line, 30 μas for the green line).

5.2. The nature of the different variability components

We present below two possible models for explaining the annual modulation plots of the outburst phase of S4 0954+658 (see Fig. 7, lower panels). Both assume that a source of variability is the propagation effects caused by a scattering screen between the source and the observer, and that a reasonable estimate of its main parameters is provided by the pre-outburst wavelet results (see Table 3 and Fig. 7, middle-left panel). These assumptions are based on the existence of an annual cycle in the pre-outburst timescales and the detection of a similar cycle in the wavelet results during the outburst phase of the source.

In the first model, we hypothesize that the second contribution to the variability is provided by a source-intrinsic mechanism; the corresponding variability timescale should show no annual modulation effect. The frequent detection of timescales of about one day in the wavelet and structure function results suggests that the source-intrinsic process could be described in terms of a variability contribution of constant timescale. The present scenario is represented in Fig. 8, left panel. The wavelet and structure function results are plotted along with a one-day constant value (green line) and an annual modulation pattern (blue line), which assumes an intrinsic angular size of the scintillating component of 50 μas. To improve the fit to the data, the anisotropy degree of the scattering pattern was increased from three to seven. The most reasonable way of explaining this change is to hypothesize that the origin of the anisotropy is the emitting component of the source; structural changes in this component would also explain the observed break in the variability characteristics of S4 0954+658.

In the second model, both the variability contributions are caused by ISS through the same scattering screen, parameterized according to the fit to the pre-outburst wavelet results. The corresponding annual modulation curves, plotted in the right panel of Fig. 8, assume intrinsic angular sizes of the scintillating components of 90 (blue line) and 30 μas (green line). Following Marscher et al. (1979), and considering an amplitude of the IDV of  ~0.1 Jy, from the size of the smaller component we can calculate its brightness temperature. This turns out to be higher than 1012 K; in order to avoid a violation of the inverse-Compton limit, we have to assume that the emitting component is moving relativistically with a Doppler boosting factor  ≥ 3.

The characteristic timescales estimated during the outburst phase of S4 0954+658 seem compatible with both models. The source-intrinsic scenario provides a simple explanation of the 11 occurrences of timescales of about one day detected in the 17 observing sessions, but its fit to the wavelet timescales is worse than the one obtained through the fully source-extrinsic scenario. None of the models explains the long timescale detected by the wavelet method in July 2008.

The limited amount of information provided by the variability analysis does not allow us to reach any final conclusion. Whatever the origin of the additional variability component, its appearance is in all cases remarkable, especially when associated with the outburst phase of the source. It seems unlikely that the simultaneity of the two events is accidental. Since the outburst phases of blazars are often explained by changes in the structure of the sources, such as the emission of new emitting components (see, e.g., Mutel et al. 1990), and since the size and luminosity of the emitting regions are supposed to play a fundamental role in IDV, the existence of a correlation between outburst phases and IDV characteristics of a source appears pretty reasonable. A possible connection between changes in the variability patterns of IDV sources and variations in their structures was reported by Fuhrmann et al. (2002) and Macquart & de Bruyn (2007), and our findings seem to confirm this picture. We note that many AGNs show intrinsic structural variability on timescales of months to years. These varying source structures complicate the search for annual modulation patterns, which may be detectable in AGNs only during sufficiently long periods of quiescence.

6. Summary

We have reported on the evolution of the intraday variability characteristics of S4 0954+65 during the period from August 2005 to December 2009. For this source, 41 IDV observing sessions have been collected at the Urumqi Astronomical Observatory at a frequency of 4.8 GHz; 38 of these sessions were long enough to provide a reliable estimation of the variability characteristics.

The data analysis has been performed using standard tools of statistics and time series analysis, such as the variability amplitude, the chi-square test for evaluating the probability of a constant flux density and the first-order structure function. A new method for the estimation of the characteristic timescale – a modified wavelet analysis based on the Ricker mother wavelet – has been developed and successfully applied to the data. The results provided by this new algorithm show a high degree of correlation with the ones obtained using a structure function analysis. We collated a comparison data set of more than 120 light curves, comprising the light curves of all the main targets of the monitoring program except S4 0954+65; for this data set, a linear regression between wavelet and structure function timescales recovers a correlation coefficient of 0.86 and a slope of 1.01, which demonstrates the consistency between the definitions of timescales in the two methods. The wavelet estimation of the timescales offers the important advantages of providing an unambiguous definition of the variability timescale and the possibility of full automation.

The combined use of a wavelet and structure function analysis proved to be a powerful tool for tracing the evolution of the IDV timescales of S4 0954+65. Both methods reveal that there is an annual modulation cycle in the variability timescales, suggesting an essential contribution of propagation effects to the detected variability. The seasonal cycle is not as strong as the ones detected in fast scintillators, such as J1819+3845 (Dennett-Thorpe & de Bruyn 2000), or in the type-II source J1128+5925 (Gabányi et al. 2007), but the agreement between the cycles observed in the wavelet timescales before and after February 2008, with slow-down phases around DoY 140 and 310, and the confirmation of a slow-down phase around DoY 320 in the structure function results are strong arguments in favor of a source-extrinsic origin of the variability.

The comparison of the results from the two methods contributed to the discovery of a substantial change in the IDV characteristics of the source, occurring between February and March 2008 at the beginning of a strong outburst phase. The signature of this change is a sudden lengthening in the average timescales detected by the wavelet-based analysis, an increase in the variability amplitude measured through the structure function, and a dramatic weakening in the correlation among the results provided by the two methods. The variability detected in the light curves is likely the result of the superposition of two independent contributions. The data seem to favor a source-extrinsic origin for both the contributions, but the results cannot be unambiguously interpreted. Nevertheless, the clear correlation between the long-term flux density variations in S4 0954+65 – likely associated with changes in the source structure – and its IDV features shows that the intrinsic properties of the source exert a significant influence on the characteristics of the variability on the shortest timescales.


1

Note that, in the case of ISS-induced variability, the rms of the flux density is expected to increase linearly with the average flux density (see Narayan 1992). This implies that the structure functions calculated during the outburst phase are expected to be systematically higher than the ones obtained during the pre-outburst phase. To make the results comparable, we applied the structure function to the light curves normalized by their average flux density.

Acknowledgments

We would like to thank Peter Müller, who provided the Python-based software package TOOLBOX for the processing of the Urumqi raw data. N.M. is funded by an ASI fellowship under contract number I/005/11/0. K.É.G. was supported by the Hungarian Scientific Research Fund (OTKA, K72515). This paper made use of data obtained with the 25 m Urumqi telescope of the Xinjiang Astronomical Observatory of the Chinese Academy of Sciences (CAS). Liu X. is supported by the National Natural Science Foundation of China under grant No.11073036 and the 973 Program of China (2009CB824800).

References

  1. Bevington, P. R., & Robinson, D. K. 1992, 2nd edn. (New York: McGraw-Hill) [Google Scholar]
  2. Bignall, H. E., Jauncey, D. L., Lovell, J. E. J., et al. 2003, ApJ, 585, 653 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bignall, H. E., Macquart, J.-P., Jauncey, D. L., et al. 2006, ApJ, 652, 1050 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cimò, G., Beckert, T., Krichbaum, T. P., et al. 2002, PASA, 19, 10 [NASA ADS] [CrossRef] [Google Scholar]
  5. Dennett-Thorpe, J., & de Bruyn, A. G. 2000, ApJ, 529, L65 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  6. Fiedler, R. L., Dennison, B., Johnston, K. J., & Hewish, A. 1987, Nature, 326, 675 [NASA ADS] [CrossRef] [Google Scholar]
  7. Fiedler, R., Dennison, B., Johnston, K. J., et al. 1994, ApJ, 430, 581 [NASA ADS] [CrossRef] [Google Scholar]
  8. Fuhrmann, L. 2004, Ph.D. Thesis, Bonn University, Germany [Google Scholar]
  9. Fuhrmann, L., Krichbaum, T. P., Cimò, G., et al. 2002, PASA, 19, 64 [Google Scholar]
  10. Gabányi, K. É., Marchili, N., Krichbaum, T. P., et al. 2007, A&A, 470, 83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Heeschen, D. S. 1984, AJ, 89, 1111 [NASA ADS] [CrossRef] [Google Scholar]
  12. Heeschen, D. S., Krichbaum, T., Schalinski, C. J., & Witzel, A. 1987, AJ, 94, 1493 [NASA ADS] [CrossRef] [Google Scholar]
  13. Heithausen, A., Stacy, J. G., de Vries, H. W., Mebold, U., & Thaddeus, P. 1993, A&A, 268, 265 [NASA ADS] [Google Scholar]
  14. Kedziora-Chudczer, L., Jauncey, D. L., Wieringa, M. H., et al. 1997, ApJ, 490, L9 [NASA ADS] [CrossRef] [Google Scholar]
  15. Kellermann, K. I., & Pauliny-Toth, I. I. K. 1969, ApJ, 155, L71 [NASA ADS] [CrossRef] [Google Scholar]
  16. Kinman, T. D., & Conklin, E. K. 1971, Astrophys. Lett., 9, 147 [NASA ADS] [CrossRef] [Google Scholar]
  17. Kraus, A., Krichbaum, T. P., Wegner, R., et al. 2003, A&A, 401, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Lovell, J. E. J., Jauncey, D. L., Bignall, H. E., et al. 2003, AJ, 126, 1699 [NASA ADS] [CrossRef] [Google Scholar]
  19. Lovell, J. E. J., Rickett, B. J., Macquart, J.-P., et al. 2008, ApJ, 689, 108 [NASA ADS] [CrossRef] [Google Scholar]
  20. Macquart, J.-P., & de Bruyn, A. G. 2007, MNRAS, 380, L20 [NASA ADS] [Google Scholar]
  21. Marchili, N. 2009, Ph.D. Thesis, Bonn University, Germany [Google Scholar]
  22. Marchili, N., Martí-Vidal, I., Brunthaler, A., et al. 2010, A&A, 509, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Marchili, N., Krichbaum, T. P., Liu, X., et al. 2011, A&A, 530, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Marscher, A. P., Marshall, F. E., Mushotzky, R. F., et al. 1979, ApJ, 233, 498 [NASA ADS] [CrossRef] [Google Scholar]
  25. Mutel, R. L., Phillips, R. B., Su, B., & Bucciferro, R. R. 1990, ApJ, 352, 81 [NASA ADS] [CrossRef] [Google Scholar]
  26. Narayan, R. 1992, Roy. Soc. London Philos. Trans. Ser. A, 341, 151 [CrossRef] [Google Scholar]
  27. Qian, S.-J., & Zhang, X.-Z. 2001, Chin. J. Astron. Astrophys., 1, 133 [NASA ADS] [CrossRef] [Google Scholar]
  28. Quirrenbach, A., Witzel, A., Wagner, S., et al. 1991, ApJ, 372, L71 [NASA ADS] [CrossRef] [Google Scholar]
  29. Quirrenbach, A., Witzel, A., Kirchbaum, T. P., et al. 1992, A&A, 258, 279 [NASA ADS] [Google Scholar]
  30. Quirrenbach, A., Kraus, A., Witzel, A., et al. 2000, A&AS, 141, 221 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Penprase, B. E. 1993, ApJS, 88, 433 [NASA ADS] [CrossRef] [Google Scholar]
  32. Perryman, M. A. C., & ESA 1997, ESA Spec. Publ., 1200 [Google Scholar]
  33. Rickett, B. J., Coles, W. A., & Bourgois, G. 1984, A&A, 134, 390 [NASA ADS] [Google Scholar]
  34. Rickett, B. J., Witzel, A., Kraus, A., Krichbaum, T. P., & Qian, S. J. 2001, ApJ, 550, L11 [NASA ADS] [CrossRef] [Google Scholar]
  35. Shapirovskaya, N. Y. 1978, Sov. Ast., 22, 544 [NASA ADS] [Google Scholar]
  36. Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46 [NASA ADS] [CrossRef] [Google Scholar]
  37. Sun, X. H., Reich, W., Han, J. L., Reich, P., & Wielebinski, R. 2006, A&A, 447, 937 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Wagner, S., Sanchez-Pons, F., Quirrenbach, A., & Witzel, A. 1990, A&A, 235, L1 [NASA ADS] [Google Scholar]
  40. Wagner, S. J., Witzel, A., Krichbaum, T. P., et al. 1993, A&A, 271, 344 [NASA ADS] [Google Scholar]
  41. Wills, B. J. 1971, ApJ, 169, 221 [NASA ADS] [CrossRef] [Google Scholar]
  42. Witzel, A., Heeschen, D. S., Schalinski, C., & Krichbaum, T. P. 1986, Mitteilungen der Astronomischen Gesellschaft Hamburg, 65, 239 [NASA ADS] [Google Scholar]

Appendix A: Wavelet analysis

The wavelet-based analysis we used to estimate the timescales in our datasets is based on the Ricker mother wavelet Ψ(t)=C(1t2)exp(t22),\appendix \setcounter{section}{1} \begin{equation} \Psi(t)=C\,\left(1-t^2\right)\exp\left(-\frac{t^2}{2}\right), \end{equation}(A.1)where C is a normalization factor and t is the time. Given the scale s and the time-shift l, we define the wavelet Ψs,l(t)=iN(1(til)2s2)exp((til)22s2)\appendix \setcounter{section}{1} \begin{equation} \Psi_{s,\,l}(t)=\sum_{i}^{N}\,\left(1-\frac{(t_\mathrm{i}-l)^2}{s^2}\right)\exp\left(-\frac{(t_\mathrm{i}-l)^2}{2\,s^2}\right) \end{equation}(A.2)and the wavelet transform ψ(s,l)=iN(SiS)(1(til)2s2)exp((til)22s2)·\appendix \setcounter{section}{1} \begin{equation} \psi(s, l)=\sum_{i}^{N}\,(S_\mathrm{i}-\langle S\rangle)\,\left(1-\frac{(t_\mathrm{i}-l)^2}{s^2}\right)\exp\left(-\frac{(t_\mathrm{i}-l)^2}{2\,s^2}\right)\cdot \label{eq:wav} \end{equation}(A.3)For each variability curve, we calculated ψ(s,l) on a range of scales s, from twice the average sampling to half the total length of the curve, with l varying between (t0 + s/2) and (tNs/2); we defined the variability timescale as 3sM\hbox{$\sqrt{3}\,s_M$}, where sM is the scale for which ψ(s,l) is maximum, noting that, given a scale s, 3s\hbox{$\sqrt{3}\,s$} is the time lag between the minimum and the maximum of Ψs,   l(t).

The results provided by the method described above are not satisfactory. To allow a proper comparison between the wavelet transforms, iNΨs,l(ti)\hbox{$\sum_i^N \Psi_{s,\,l}(t_i)$} should be close to zero for each set of parameters (s,l). This condition is often unfulfilled, because of the uneven sampling of the variability curves and their limited duration. To correct for this, we introduced a modified wavelet transform φ(s,l)=Cs,lψ(s,l),\appendix \setcounter{section}{1} \begin{equation} \phi(s, l)=C_{s,\,l}\,\psi(s, l), \label{eq:modwav} \end{equation}(A.4)where Cs,   l is defined as Cs,l=iN((1(til)2s2)exp((til)22s2))2·\appendix \setcounter{section}{1} \begin{equation} C_{s,\,l}=\sqrt{\sum_{i}^{N}\left(\left(1-\frac{(t_\mathrm{i}-l)^2}{s^2}\right)\exp\left(-\frac{(t_\mathrm{i}-l)^2}{2\,s^2}\right)\right)^2}\cdot \label{eq:cls} \end{equation}(A.5)As for the standard wavelet analysis, the variability timescale is defined as 3sM\hbox{$\sqrt{3}\,s_M$}, where sM is the scale that corresponds to the maximum value of φs,   l. The standard wavelet analysis is generally more sensitive to slow variability components than to fast ones, often detecting unreasonably long characteristic timescales. The results provided by the modified wavelet, instead, appear to be much more reliable.

The modified wavelet transform has a strong advantage over the structure function analysis in dealing with light curves with variability on multiple timescales. While the identification of the characteristic timescale in the structure function requires a subjective interpretation of which timescale is the most predominant, the wavelet analysis always provides a unique value, because the set (s,l) that maximizes φs,   l is always univocally defined. This completely removes the ambiguity when interpreting the analysis results – regardless of who performs the analysis, the same timescale is found. It also implies that the wavelet analysis can be fully automated, which can be very useful in dealing with large datasets. On the other hand, one should always keep in mind that the timescale obtained by the wavelet method, strictly speaking, does not provide an estimation of the characteristic variability timescale. It refers, instead, to the strongest variability feature in the light curve. When the variability timescale in a light curve is much shorter than the duration of the observations, the structure function should be regarded as a more suitable tool for the data analysis.

We note that the IDV of blazars – such as S4 0954+65 – at radio frequencies can be generally described as red noise. The amplitude of the variability is stronger on longer timescales, which implies that the variability timescale is typically on the same order of magnitude as the duration of the observations. In these conditions, the identification of the characteristic timescales with the strongest variability feature in the light curves appears quite reasonable. This conclusion is convincingly supported by the high degree of correlation between the structure function and the wavelet results (see Sect. 4).

The error in the timescale estimates provided by the modified wavelet transform can be ascribed to two main contributions. The first consists in the combined effect of the uneven sampling and the uncertainty in the flux density measurements; the second has to do with the limited duration of the observations, which causes a gradual decrease in the reliability of the wavelet

transforms on larger scales s. Both effects are reflected in the width of the peak ψ(sM,l). Wider peaks are indicative of a low contrast among the strengths of the variability on different scales s, hence larger errors in the estimation of the timescale. We performed tests with synthetic light curves of different durations, sampled both evenly and unevenly. We found that a reasonable evaluation of the error in the timescales is given by (s2s1)/2, where s2 and s1 are respectively the maximum and the minimum scales for which ψ(s,l) equals 0.95ψ(sM,l).

All Tables

Table 1

Basic information about the 38 observing sessions at the Urumqi Observatory during which S4 0954+65 was observed.

Table 2

Variability characteristics of S4 0954+65 in 38 epochs of observation with the Urumqi Telescope.

Table 3

Scintillation pattern parameters deduced by the distribution of the timescales as a function of the day of the year.

All Figures

thumbnail Fig. 1

Variability curves of S4 0954+65 (black dots) in April 2006 (upper panel) and April 2008 (lower panel). The turquoise squares show the residual variability in the fluxes of the secondary calibrators S5 0836+71 and S4 0951+69.

In the text
thumbnail Fig. 2

Upper panels: the structure function plots of the April 2006 (left panel) and April 2008 (right panel) datasets; the blue vertical lines show the estimated timescales. Lower panels: the wavelet-based analysis for the same epochs; the green vertical lines show the estimated timescales. The existence of two bumps in the right panels is the signature of two distinct variability components in the dataset.

In the text
thumbnail Fig. 3

Long-term variability of S4 0954+65 (upper panel). In the middle panel, the wavelet timescales (red dots) and a three-point running average (blue line), showing the sudden increase starting from March 2008. The structure function timescales are shown in the lower panel. The red vertical line divides the observing sessions into two phases characterized by remarkably different variability characteristics.

In the text
thumbnail Fig. 4

Comparison between the structure function and wavelet results for the comparison sample (black dots) and the S4 0954+65 light curves. In the left panel, the source’s timescales estimated up to February 2008 (cyan squares); in the right panel, the timescales from March 2008 to December 2009 (green squares). The black, cyan, and green lines show the results of applying linear regressions to the respective datasets.

In the text
thumbnail Fig. 5

Average structure function calculated from the normalized light curves obtained during the pre-outburst phase (turquoise diamonds) and the outburst phase (red squares).

In the text
thumbnail Fig. 6

Annual modulation plots based upon the wavelet results (left panels), and the structure function results (right panels). The upper panels show all the estimated timescales, while lower panels show monthly averages. The green and turquoise lines show the annual modulation patterns which best fit the data.

In the text
thumbnail Fig. 7

Annual modulation plots, as in Fig. 6; in the upper panels, the timescales obtained for the epochs between August 2005 and February 2008 are highlighted in magenta squares. In the middle and lower panels are displayed the annual modulation plots (green and turquoise lines) and the 40-day timescales averages for the August 2005–February 2008 and the March 2008–December 2009 epochs. The red vertical lines indicate the peaks in the annual modulation patterns obtained from the complete set of timescales.

In the text
thumbnail Fig. 8

Wavelet (black dots) and structure function (red squares) timescales calculated during the outburst phase of S4 0954+658 are indicative of multiple variability contributions. In the left panel, the variability is interpreted as the superposition of a source-extrinsic effect (blue line) and a source-intrinsic process (green line). In the right panel, both the variability contributions are attributed to ISS and the two annual modulation patterns have the same scattering screen parameters, the only difference being the size of the scintillating components (90 μas for the blue line, 30 μas for the green line).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.