Radio Variability and Random Walk Noise Properties of Four Blazars^{1}^{1}affiliation: Based on observations obtained by the University of Michigan Radio Astronomy Observatory (UMRAO)
Abstract
We present the results of a time series analysis of the longterm radio lightcurves of four blazars: 3C 279, 3C 345, 3C 446, and BL Lacertae. We exploit the data base of the University of Michigan Radio Astronomy Observatory (UMRAO) monitoring program which provides densely sampled lightcurves spanning 32 years in time in three frequency bands located at 4.8, 8, and 14.5 GHz. Our sources show mostly flat or inverted (spectral indices ) spectra, in agreement with optically thick emission. All lightcurves show strong variability on all time scales. Analyzing the time lags between the lightcurves from different frequency bands, we find that we can distinguish highpeaking flares and lowpeaking flares in accord with the classification of Valtaoja et al. (1992). The periodograms (temporal power spectra) of the observed lightcurves are consistent with randomwalk powerlaw noise without any indication of (quasi)periodic variability. The fact that all four sources studied are in agreement with being randomwalk noise emitters at radio wavelengths suggests that such behavior is a general property of blazars.
Subject headings:
Galaxies: active — radiation mechanisms: nonthermal — methods: statistical1. Introduction
The strong and complex temporal flux variability of Active Galactic Nuclei (AGN; see, e.g., Beckmann & Shrader 2012 and references therein for a recent review) provides valuable information on the internal conditions of accretion zones and plasma outflows. Various characteristic variability patterns have been associated with a wide range of physical phenomena, from shocks in continuous (e.g., Marscher & Gear 1985) or discontinuous (e.g., Spada et al. 2001) jets to orbiting plasma “hotspots” (e.g., Abramowicz et al. 1991) or plasma density waves (e.g., Kato 2000) in accretion disks. Accordingly, multiple studies have aimed at quantifying the properties of AGN variability on all time scales and throughout the electromagnetic spectrum. At radio wavelengths, variability time scales probed by observations range from tens of minutes (Schödel et al. 2007, studying the mm/radio lightcurve of M 81*; see also Kim & Trippe 2013 for a discussion of the detectability of intraday variability) to tens of years (Hovatta et al. 2007, in a statistical analysis of the longterm flux variability of 80 AGN). Of particular interest is the possible presence of quasiperiodic oscillations (QPO) which has been reported by several studies of blazar lightcurves (e.g. Rani et al. 2009, 2010; Gupta et al. 2012).
Fourier transform, period folding, power spectrum, and periodogram methods (cf. Priestley 1981 for an exhaustive review of time series analysis) have been used extensively for quantifying the statistical properties of AGN variability and for the search for possible QPOs (e.g. Benlloch et al. 2001 for Xray, Webb et al. 1988 for optical, Fan 1999 for near infrared, and Aller et al. 2003 for radio observations). As already noted by Press (1978), power spectra of AGN lightcurves globally follow power laws with , corresponding to red noise;^{1}^{1}1In the context of time series analysis, the term “noise” refers to stochastic emission from a source of radiation, not to measurement errors or instrumental noise. here denotes the power spectral amplitude as function of sampling frequency . Lightcurves composed of pure Gaussian white noise have flat power spectra (). Other important special cases are random walk noise () – which is the integral of white noise – and flicker noise (, “ noise”) as intermediate case between white noise and random walk noise (see also Park & Trippe 2012 for a detailed technical discussion). Lawrence et al. (1987) found that the power spectrum of the Seyfert galaxy NGC 4051 can be described by flicker noise. Red noise power spectra were also observed by Lawrence & Papadakis (1993) who used 12 highquality “long look” Xray lightcurves of AGN.
Object  RA  DEC  Type  Redshift  [yr]  

3C 279  12:56:11  05:47:22  FSRQ  0.536  32.52  1086  1337  1473 
3C 345  16:42:59  +39:48:37  FSRQ  0.593  32.54  1323  1315  1415 
3C 446  22:25:47  04:57:01  BL Lac  1.404  32.12  680  902  1088 
BL Lac  22:02:43  +42:16:40  BL Lac  0.069  32.54  1256  1315  1755 
Note. – J2000 coordinates, source types, and redshifts are taken from the NED. We also give the total monitoring time (in years) and the numbers of flux data points for 4.8, 8.0, and 14.5 GHz, respectively.
A multitude of studies illustrates the difficulties of determining the statistical significance of supposed QPO signals in the power spectra of AGN lightcurves. The analysis of Benlloch et al. (2001) concluded that a previously reported quasiperiodic signal in Xray lightcurves of the Seyfert galaxy Mrk 766 was actually statistically insignificant. Uttley et al. (2002) pointed out the importance of sampling effects leading to rednoise leaks and aliasing. Vaughan (2005) gives an analytical approach to derive significance levels for peaks in rednoise power spectra. Do et al. (2009) demonstrated the power of MonteCarlo techniques for deriving significance levels by comparing the rednoise power spectra of actual and simulated flux data.
The temporal flux variability of AGN can be exploited for elucidating the physical conditions within active galaxies especially at radio frequencies where monitoring observations of hundreds of targets have been conducted over several decades by various observatories. Remarkably, many studies aimed at analyzing longterm AGN radio variability do not take into account the intrinsic rednoise properties of the lightcurves. The incorrect assumption of constant (as function of ) significance levels in power spectra (following from the assumption of whitenoise dominated lightcurves) has lead to reports of “characteristic” time scales which are actually not special at all (cf., e.g., Ciaramella et al. 2004; Hovatta et al. 2007; Niepolla et al. 2009).
Blazars, characterized by violent flux variability across the entire electromagnetic spectrum, are a subset of AGN which include BL Lacertae (BL Lac) objects and Flat Spectrum Radio Quasars (FSRQs). In accordance with the standard viewing angle unification scheme of AGN (Urry & Padovani, 1995), it is commonly assumed that their observed emission is generated by synchrotron radiation – dominating from radio to optical frequencies – and inverse Compton emission – dominating at frequencies higher than optical – from relativistic plasma jets (almost) aligned with the line of sight. In order to perform a thorough study of the statistical properties of blazar emission, we analyze the lightcurves of four radiobright blazars with strong flux variability – 3C 279, 3C 345, 3C 446, and BL Lac – provided by the University of Michigan Radio Astronomy Observatory (UMRAO) monitoring program of AGN. The data set comprises data spanning 32 years in time and covering three frequency bands located at 4.8 GHz, 8.0 GHz, and 14.5 GHz.
2. Target Selection and Flux Data
For our study we exploited the AGN monitoring data base of the 26meter University of Michigan Radio Astronomy Observatory (UMRAO); the instrument, observations, and calibration procedures are described in detail by Aller et al. (1985). Our analysis required the use of densely sampled highquality lightcurves spanning several decades in time and obtained at several observing frequencies. Accordingly, we selected targets with (i) data available for all three UMRAO bands (4.8, 8, and 14.5 GHz); (ii) continuously spanning at least 30 years in time; (iii) dense – faster than monthly at each frequency – sampling over the entire monitoring time; (iv) a minimum flux (at all frequencies) of 2 Jy; and (v) strong flux variability by factors 2. Our very strict selection criteria left us with a sample of four blazars: 3C 279, 3C 345, 3C 446, and BL Lac. Table 1 provides an overview over the key properties of our targets (partially taken from the NASA/IPAC Extragalactic Database (NED)^{2}^{2}2http://ned.ipac.caltech.edu/) and data. The median statistical error of a flux measurement was 0.09 Jy. The lightcurves cover a time line from 1980 to 2012, slightly more than 32 years.
3. Analysis
3.1. Lightcurves
We selected our data for the purpose of time series analysis which can be misled by irregular sampling. In order to minimize such sampling effects, we binned our lightcurves in time such that the bin size is the time interval
(1) 
where is the total observing time and is the number of data points; for our sources, is on the order of three weeks typically. (For the special case of regular sampling, corresponds to the inverse of the Nyquist frequency.) The final lightcurves are shown in Fig. 1; evidently, all four sources show strong variability on various time scales.
3.2. Spectral indices
The fast – but not simultaneous – sampling of our targets at three frequencies made it possible to study their spectral evolution. A combination of (i) nonsimultaneous sampling and (ii) rapid intrinsic flux variability made it necessary to group our data into time windows; we eventually chose a time window of one year. Within each time window, we jointly described all data (covering all frequency bands) with a standard powerlaw model
(2) 
where is the observing frequency, is the flux density, and is the spectral index. We present the spectral index as function of time in Fig. 1.
3.3. Time offsets among spectral bands
The long overall time line and good sampling of our data made it possible to probe the data for time lags among the fluxes observed at different frequencies. We applied the discrete correlation function proposed by Edelson & Krolik (1988) to each pair of spectral bands for each source. In a first step, we computed the unbinned discrete correlation function (UDCF) for two discrete datasets {} and {} with ,
(3) 
Here and are the averages of {} and {}, respectively; denotes the difference of the observing times of the data pair (); are the variances; and denote the mean measurement errors of , .
The actual discrete correlation function (DCF) for a given time offset results from averaging over all for which falls into a selected bin (i.e., ):
(4) 
with () +1 corresponding to perfect (anti)correlation and 0 indicating the absence of any correlation. The position of the maximum of the DCF corresponds to the time offset between the lightcurves. The statistical uncertainty of DCF() is given by the standard error of mean
(5) 
An important effect to be considered is the interplay between (i) variations of the spectral index and (ii) a time delay between the lightcurves belonging to different frequency bands: an offset between spectral bands causes variations of the observed values for even if the intrinsic (i.e., corrected for the offset) spectral index is constant. Inspection of Fig. 1 suggests that this might indeed be the case at several occasions in 3C 279, 3C 345, and 3C 446. Turning this argument around, this implies that we can divide our data into activity phases defined by the times when (i.e. the reference spectral index for flatspectrum AGN). Accordingly, we divided our lightcurves into two or three phases (A, B, C – except for the case of BL Lac), and computed the DCF for each phase separately. We chose ranges sufficient for covering the largest time offsets expected between lightcurves, eventually adopting years (except for the lightcurve pair 4.8/14.5 GHz of 3C 279 where it was necessary to extend the range to years). In order to preserve a good time resolution, we usually used a bin size for two lightcurves “1” and “2”. In case of the 4.8/8 GHz lightcurve pair of 3C 345 it was necessary to increase by factors up to six in order to suppress sampling artifacts.
We present the resulting DCF in Fig. 2. In our convention, a positive (negative) time lag implies that the flux at the higher frequency precedes (follows) the flux at the lower frequency. As all values of a given DCF curve are computed from the same time series, adjacent DCF points are highly correlated. Depending on the bin size, this can cause the scatter of the DCF values to be substantially smaller than their individual errors given by Eq. 5. This implies that (i) the scatter of the DCF values is not a proxy for their uncertainties and (ii) averaging over or fitting of parametric models to multiple DCF values does not improve the accuracy of the estimate for the location of the peak of the DCF curve. With this in mind, we consider the null hypothesis “the lightcurves are simultaneous” as rejected if and only if DCF() is located below the line defined by (with denoting the statistical error of the maximum value of the DCF).
3.4. Periodograms
For a quantitative analysis of flux variability we employed the normalized Scargle periodogram
(6) 
where is a time offset satisfying
(7) 
(Scargle, 1982). Here is the angular frequency, is the amplitude of the periodogram evaluated at sampling frequency , is the th flux value, denotes the time when was obtained, and is the variance of the data. The base frequency is , the sampling frequencies are . Here is the total observing time and is the number of flux data points. The Scargle periodogram is preferable over standard Fourier transform methods because it can be applied to data with arbitrary sampling and has a wellunderstood statistical behavior (Priestley, 1981; Scargle, 1982). We present the periodograms of our lightcurves in Fig. 3.
The power spectra of AGN lightcurves are known (cf. § 1) to follow rednoise powerlaws. However, when dealing with lightcurves that are sampled irregularly and show gaps, a complication arises in form of aliasing: the power at frequencies above the Nyquist frequency is transferred to lower frequencies because variations on time scales shorter than the sampling period cannot be distinguished from variations on longer time scales. Aliasing introduces an approximately constant offset which adds to the power spectrum (Uttley et al., 2002). As (i) the amplitudes of power spectra tend to span several orders of magnitude and (ii) are affected by multiplicative noise (Scargle, 1982; Vaughan, 2005), periodograms have to be treated in logarithmic space. Accordingly, we assumed the functional form
(8) 
3.5. Simulated lightcurves and significance levels
The detection of deviations from a rednoise powerlaw periodogram, especially of (quasi)periodic signals at specific sampling frequencies, requires the establishment of reliable significance levels – a problem that has been notoriously difficult (cf., e.g., Vaughan 2005). Arguably the most straightforward ansatz is provided by MonteCarlo simulations that compare simulated periodograms and lightcurves to actual data (e.g., Benlloch et al. 2001; Do et al. 2009), and this is the approach we adopted.
We simulated rednoise lightcurves using the method suggested by Timmer & König (1995). For each sampling frequency , we drew two random numbers from Gaussian distributions with zero mean and unit variance for the real part and the imaginary part, respectively. We multiplied both numbers with to generate powerlaw noise with slope . The result was an array of complex numbers which corresponded to the complex Fourier transform of the artificial lightcurve. We constructed each complex array such that the values of the Fourier transform obey , with denoting complex conjugation, to obtain a real valued time series. We computed artificial lightcurves by taking the inverse Fourier transform of the complex arrays.
Each artificial lightcurve consisted of 4096 data points initially. We then identified the time between two adjacent data points with the time scale defined in Eq. 1 and omitted values from the artificial lightcurve at the locations of gaps in the observed lightcurves, thus mapping the sampling pattern of the UMRAO observations into the simulated lightcurve. Eventually, we computed a periodogram from each remapped artificial lightcurve. As the observed indices are for all blazar lightcurves (cf. § 4), we used values of 1.5, 1.75, and 2 for in the simulations. In order to decide which value for to adopt for a given observed periodogram, we (i) computed 10 000 simulated periodograms for each of the three choices of , and (ii) compared the average of the simulated periodograms to the observed periodogram via a weighted leastsquares test.
From the set of 10 000 artificial periodograms for a given blazar lightcurve we determined, separately for each sampling frequency , a 99.9% significance level (corresponding to in Gaussian terms) for the periodogram derived from the UMRAO data. Our simulation procedure is based on the null hypothesis “the observed periodogram originates from a rednoise lightcurve”. Accordingly, the spectral power of a deviation from a rednoise periodogram, especially a candidate periodic signal, needs to exceed the aforementioned significance levels in order to be potentially significant. The results of our analysis are illustrated in Fig. 3.
4. Results
4.1. 3c 279
Using our criteria outlined in § 3.3 we divided the lightcurves of 3C 279 into three activity phases (see also Fig. 1). Phase A, ranging from 1980 to 1990, is characterized by relatively weak variability with fluxes ranging approximately from 10 Jy to 14 Jy. Flux densities in all three frequency bands are very similar for most of the time, leading to spectral indices except of the very beginning of phase A. The DCF analysis (§ 3.3 and Fig. 2) finds that the 14.5GHz lightcurve precedes the 4.8GHz lightcurve by years ( uncertainty interval); time lags between the other lightcurve pairs are consistent with zero. Phase B, ranging from 1990 to 2003, is characterized by a strong increase in emission from 10 Jy to 30 Jy (at 14.5 GHz). This outburst occurs the earlier the higher the frequency; for the frequency pair 4.8/14.5 GHz, the time lag between the lightcurves is years. The observed spectral index becomes inverted, with as low as about . Phase C, starting in 2003, is characterized by multiple flux density fluctuations in the range 10–20 Jy for most of the time, with a strong increase – up to 35 Jy at 14.5 GHz – since 2010. The spectral index remains mildly inverted () for most of the time but reaches in 2012, coinciding with the observed flux maximum at the end of the monitoring. The time lags are consistent with zero for all frequency pairs.
The periodograms (§ 3.4, Fig. 3) of all three lightcurves decrease toward increasing sampling frequencies – as is characteristic for red noise – and show a flattening at the highest sampling frequencies – as expected in case of notable aliasing. By fitting the model given by Eq. 8 to the spectra we find approximate powerlaw indices of 2.5, 1.8, and 1.5 for the 4.8GHz, 8GHz, and 14.5GHz periodograms, respectively, with statistical errors between 0.2 and 0.4 ( confidence intervals). Comparison of observed periodograms to the ones found from MonteCarlo simulations (averages of 10 000 realizations of rednoise periodograms; cf. § 3.5) leads to the conclusion that all three observed power spectra are consistent with being generated by randomwalk noise () lightcurves. None of the observed periodograms shows statistically significant excess power with respect to a rednoise power spectrum.
4.2. 3c 345
We divided the lightcurves into two activity phases. Phase A, ranging from 1980 to 1994, is characterized by strong variability, with fluxes (at 14.5 GHz) moving between 5 and 17 Jy. This variability is also expressed in the spectral index which fluctuates between and 0.2. The DCF analysis finds a time lag of years ( confidence intervals) between the 4.8 GHz and 14.5 GHz lightcurves, with smaller, marginally significant, lags between the other frequency pairs. Phase B, beginning in 1994, is characterized by strong flux variability between 5 and 13 Jy and a mostly flat () spectrum. Time lags between the lightcurves are (marginally) consistent with zero.
The periodograms of all three lightcurves are consistent with pure rednoise spectra. The bestfitting parametric model solutions (Eq. 8) show quite extreme powerlaw slopes with statistical ( confidence) uncertainties of about 0.3. When comparing the data to the results of MonteCarlo simulations, we find that all three observed periodograms are consistent with randomwalk noise spectra.
4.3. 3c 446
During the entire observing time this source shows strong variability in both flux density – with values ranging from 3 Jy to 10 Jy (at 14.5 GHz) – and spectral index – with values fluctuating between and 0.3. According to our criteria (§ 3.3) we divided the lightcurves into three phases, ranging from 1980 to 1986 (phase A), 1986 to 1995 (phase B), and from 1995 onward (phase C). Whereas in phases A and B the lightcurves of all three frequencies are consistent with being simultaneous, we find a time lag of years ( confidence interval) for the pair 4.8/14.5 GHz.
The periodograms of all three lightcurves are consistent with being pure rednoise spectra. Our parametric model fit (Eq. 8) finds powerlaw slopes for all three periodograms, with statistical () errors between 0.2 and 0.4. From our MonteCarlo simulations we find that the observed periodograms are consistent with randomwalk noise lightcurves ().
4.4. BL Lac
The lightcurve of BL Lac is characterized by rapid variability throughout the entire monitoring time of 32 years. Flux densities vary between 2 Jy and 7 Jy with the notable exception of a flare that reaches about 15 Jy (at 14.5 GHz) in 1980. The spectral index varies only slowly (except at the time of the 1980 flare) – on time scales of years to decades – between and . Accordingly, we did not attempt to identify separate activity phases – qualitatively, the behavior of BL Lac is actually rather uniform. The lightcurves at the three observing frequencies follow each other closely not only in amplitude but also in time: all time lags identified by the crosscorrelation analysis are consistent with zero.
The periodograms of all three lightcurves are in agreement with being pure rednoise power spectra. The parametric model (Eq. 8) finds identical (within the errors of 0.1–0.2) for all periodograms. Comparison to the simulation results shows the best agreement with a theoretical, intrinsic powerlaw slope of .
5. Discussion
5.1. Spectral indices
For all sources the spectral index remains at values that are close to zero or even negative (). The low values of imply that the emission originates from optically thick synchrotron sources, leading to approximately flat () or even inverted () spectra (Ginzburg & Syrovatskii, 1965; Pacholczyk, 1970; Kembhavi & Narlikar, 1999; Krolik, 1999). This is in agreement with blazars being AGN with jets pointing (almost) toward the observer, resulting in a high column density of matter (potentially belonging to multiple individual plasma clouds) along the line of sight – we do not find any indication for a deviation from this (expected) behavior even over a time line of three decades.
Taking a closer look at individual flux maxima (outbursts or flares of radiation), we can distinguish two types of behavior: (i) events with (a) fluxes at higher frequencies preceding the ones at lower frequencies and (b) fluxes at higher frequencies being significantly higher than the fluxes at lower frequencies (implying ); and (ii) events with all lightcurves being simultaneous and of approximately equal amplitude (implying ). A noteworthy example is provided by 3C 345, where both types of events occur within about 20 years (cf. phase A vs. phase B in Fig. 1). An interpretation is readily provided by the “generalized shock model” of Valtaoja et al. (1992) which is based on the assumption that outbursts of radio emission in AGN are caused by shocks propagating through jets and which distinguishes two scenarios: (i) in highpeaking flares (“high” with respect to the observing frequency), the maximum luminosity is reached at frequencies well above the observing frequency. This implies that the flare is decaying at the time of observation, resulting in an observational signature equivalent to the shockinjet model by Marscher & Gear (1985). This model describes shocks in AGN jets as adiabatically expanding plasmas that become optically thin at higher frequencies first, thus causing a systematic time delay between the spectral bands with the higher frequency leading with higher amplitude of flux. In case of (ii) lowpeaking flares, the maximum luminosity is reached at frequencies well below the observing frequency. Lightcurves at different frequencies are almost simultaneous and have almost identical amplitudes. Applying this framework to our sources, highpeaking flares are present in phase B and C of 3C 279, phase A of 3C 345, and the entire lightcurves of 3C 446. The other flux outbursts can be described as lowpeaking flares.
5.2. Spectral time delays
We examined the presence or absence of time lags between lightcurves at different frequencies for each phase of each source via discrete correlation functions. The first feature we note is the large range of time delays – on the order of months – permitted by the confidence limits. The different temporal evolutions of the outbursts in different spectral bands, in combination with long variability time scales on the order of years, lead to broad, asymmetric DCF curves.
As discussed already partially in the context of the spectral index analysis, we find (i) significant positive time lags, as well as (ii) phases of no or weak positive time lags between different spectral bands. Phase B of 3C 279 ( years) and phase A of 3C 345 ( years) correspond to case (i); the amplitudes reach their maxima first at 14.5 GHz, with the 8 and 4.8 GHz lightcurves trailing – as expected for highpeaking flares (Valtaoja et al., 1992). The long evolution time scales of the outbursts, on the order of years, suggest physical sizes of the expanding emission regions on the order of lightyears. The remaining activity phases correspond to case (ii), with time delays close to or in agreement with zero. Here we find both low and highpeaking flares: phase C of 3C 279 and all phases of 3C 446 show very fast spectral index variability ranging from to , in agreement with the behavior of highpeaking flares; BL Lac however shows simultaneous lightcurves with approximately identical amplitudes throughout the entire monitoring time of three decades – in agreement with the behavior expected for lowpeaking flares.
5.3. Power spectra
The periodograms of all four blazars are in good agreement with lightcurves generated by powerlaw noise with index – i.e., random walk noise – and being affected by aliasing caused by irregular sampling. Furthermore, our statistical tests (§ 3.5) show that all periodograms are consistent with being pure rednoise power spectra without significant (quasi)periodic signals (cf. Fig. 3). The powerlawnoise nature of their lightcurves implies that none of our target blazars shows any “characteristic” activity time scale.
As already outlined in § 1, the rednoise nature of AGN lightcurves is observationally well established (albeit this discussion is complicated by the potential presence of multiple states of emission; cf. e.g. DoddsEden et al. 2011; Trippe et al. 2011; Park & Trippe 2012). Empirically, the slopes of the power spectral density of different AGN tend to scatter over a wide range of values, roughly from (e.g., Trippe et al. 2011) to (e.g., Do et al. 2009), with “typical” values (e.g., Press 1978). Indeed, the presumed flickernoise nature of AGN lightcurves triggered a search for an underlying physical mechanism which has lasted for more than three decades (e.g., Press 1978; Lyubarskii 1997; Kelly et al. 2011), without any clear picture emerging as yet.
Given the incoherent picture of the statistical properties of temporal AGN variability, the clarity of our results comes as a surprise: we find the lightcurves of all four blazars to be consistent with being randomwalk signals (). Within the obvious limits of lownumber statistics, this suggests that randomwalk noise radio lightcurves are characteristic for blazars. We note the importance of a careful treatment of data (§§ 3.1, 3.4) as well as a careful modeling of red noise lightcurves (§ 3.5): only the combination of good data quality, periodogram analysis, awareness of sampling effects, and Monte Carlo simulations of powerlaw noise lightcurves unveils the intrinsic randomwalk noise behavior. Evidently, this raises the question if randomwalk noise lightcurves could be a general feature of blazars that is frequently masked by limited data quality, irregular sampling, inappropriate modeling of power spectra, et cetera.
6. Conclusions
We studied highquality radio lightcurves of four luminous blazars – 3C 279, 3C 345, 3C 446, and BL Lac – spanning 32 years in time and covering the frequencies 4.8, 8, and 14.5 GHz. We analyzed the temporal evolution of fluxes and spectral indices. Our work leads us to the following principal conclusions:

Our sources show mostly flat or inverted () spectral indices, in agreement with optically thick synchrotron emission. The lightcurves of different frequencies are either simultaneous (within errors) or shifted relative to each other such that the highfrequency emission leads the lowfrequency emission by up to 1.5 years. We are able to distinguish highpeaking and lowpeaking flares according to the classification of Valtaoja et al. (1992).

All lightcurves show variability on all time scales. Their periodograms (power spectra) are in agreement with being pure rednoise powerlaw spectra without any indication for (quasi)periodic signals. When taking into account the sampling patterns via dedicated Monte Carlo simulations, we find that all lightcurves are consistent with being random walk noise signals with powerlaw slopes . Given that we find this behavior in all four sources under study, this suggests that random walk noise lightcurves are a general feature of blazars.
Our results imply that careful time series analysis of highquality blazar lightcurves provides information on the source structure even if a target is not resolved spatially. Obviously, it will be necessary to systematically study much larger blazar samples in order to decide if the trends we have uncovered are indeed general.
References
 Abramowicz et al. (1991) Abramowicz, M.A., et al. 1991, A&A, 245, 454
 Aller et al. (1985) Aller, H. D., Aller, M. F., Latimer, G. E., & Hodge, P. E. 1985, ApJS, 59, 513
 Aller et al. (2003) Aller, M. F., Aller, H. D., & Hughes, P. A. 2003, ApJ, 586, 33
 Beckmann & Shrader (2012) Beckmann, V., & Shrader, C. 2012, Active Galactic Nuclei, Weinheim: WileyVCH
 Benlloch et al. (2001) Benlloch, S., Wilms, J., Edelson, R., et al. 2001, ApJ, 562, L121
 Ciaramella et al. (2004) Ciaramella, A., Bongardo, C., Aller, H.D., et al. 2004, A&A, 419, 485
 Do et al. (2009) Do, T., Ghez, A.M., Morris, M.R., et al. 2009, ApJ, 691, 1021
 DoddsEden et al. (2011) DoddsEden, K., et al. 2011, ApJ, 728, 37
 Edelson & Krolik (1988) Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646
 Fan (1999) Fan, J. H. 1999, MNRAS, 308, 1032
 Fan et al. (2007) Fan, J. H., Liu, Y., Yuan, Y. H., et al. 2007, A&A, 462, 547
 Ginzburg & Syrovatskii (1965) Ginzburg, V.L. & Syrovatskii, S.I. 1965, ARA&A, 3, 297
 Gupta et al. (2012) Gupta, A. C., Krichbaum, T. P. Wiita, P. J. et al. 2012, MNRAS, 425, 1357
 Hovatta et al. (2007) Hovatta, T., Tornikoski, M., Lainela, M., et al. 2007, A&A, 469, 899
 Kato (2000) Kato, S. 2000, PASJ, 53, 1
 Kelly et al. (2011) Kelly, B.C., et al. 2011, ApJ, 730, 52
 Kembhavi & Narlikar (1999) Kembhavi, A. K., & Narlikar, J. V. 1999, Quasars and Active Galactic Nuclei, Cambridge: Cambridge University Press
 Kim & Trippe (2013) Kim, J. Y., & Trippe, S. 2013, JKAS, 46, 65
 Krolik (1999) Krolik, J. H. 1999, Active Galactic Nuclei: Princeton University Press
 Lawrence & Papadakis (1993) Lawrence, A., & Papadakis, I. 1993, ApJ, 414, 85
 Lawrence et al. (1987) Lawrence, A., Watson, M. G., Pounds, K. A., & Elvis, M. 1987, Nature, 325, 694
 Lyubarskii (1997) Lyubarskii, Yu. E. 1997, MNRAS, 292, 679
 Markowitz et al. (2003) Markowitz, A., Edelson, R., Vaughan, S., et al. 2003, ApJ, 593, 96
 Marscher & Gear (1985) Marscher, A. P. & Gear, W. K. 1985, ApJ, 298, 114
 Montroll & Shlesinger (1982) Montroll, E.W. & Shlesinger, M.F. 1982, Proc. Natl. Acad. Sci. USA, 79, 3380
 Niepolla et al. (2009) Nieppola, E., Hovatta, T., Tornikoski, M., et al. 2009, AJ, 137, 5022
 Pacholczyk (1970) Pacholczyk, A.G. 1970, Radio Astrophysics, San Francisco: W.H. Freeman & Co.
 Park & Trippe (2012) Park, J. H., & Trippe, S. 2012, JKAS, 45, 147
 Press (1978) Press, W. H. 1978, Comment. Astrophys., 7, 103
 Priestley (1981) Priestley, M.B. 1981, Spectral Analysis and Time Series, London: Elsevier
 Rani et al. (2009) Rani, B., Wiita, P. J., & Gupta, A. C. 2009, ApJ, 696, 2170
 Rani et al. (2010) Rani, B., Gupta, A. C., Joshi, U. C. et al. 2010, ApJ, 719, L153
 Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
 Schödel et al. (2007) Schödel, R., Krips, M., Markoff, S., et al. 2007, A&A, 463, 551
 Spada et al. (2001) Spada, M., et al. 2001, MNRAS, 325, 1559
 Timmer & König (1995) Timmer, J., & König, M. 1995, ApJ, 300, 707
 Trippe et al. (2011) Trippe, S., Krips, M., Piétu, V., et al. 2011, A&A, 533, A97
 Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
 Uttley et al. (2002) Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231
 Valtaoja et al. (1992) Valtaoja, E., Tersranta, H., Urpo, S., et al. 1992, A&A, 254, 71
 Vaughan (2005) Vaughan, S. 2005, A&A, 431, 391
 Webb et al. (1988) Webb, J. R., Smith, A. G., Leacock, R. J., et al. 1988, AJ, 95, 374