© ESO, 2011

1. Introduction

Active galactic nuclei (AGN) have been studied extensively in the wavelength range from cm-radio to γ-rays in the last decades (see, e.g., Kembhavi & Narlikar 1999; or Krolik 1999, and references therein for a review). There is overwhelming observational evidence that their emission originates from accretion onto supermassive black holes (SMBH) with masses M_• ≃ 10^6...9   M_⊙ (e.g., Ferrarese & Ford 2005, and references therein).

The viewing-angle unification standard model of AGN (Lawrence 1987; see also Urry & Padovani 1995) has long been established, although recent studies emphasize the impact of galaxy mergers on AGN fueling and luminosity (e.g. Hopkins et al. 2006). However, the details of AGN activity are still not well understood. One example is the nuclear emission; proposals for its origin range 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.

One possibility to constrain the AGN physics is offered by the analysis of characteristic timescales and the noise properties of the emission. Several studies aimed at identifying characteristic timescales, including possible quasi-periodic oscillations, from tens of minutes (using X-ray observations; e.g. Benlloch et al. 2001, and using mm-radio observations; e.g. Schödel et al. 2007) to tens of years (based on cm- and mm-radio monitoring data; Hovatta et al. 2007; 2008). Other works focused on the power spectral statistics of the emission which was found to globally follow red-noise type powerlaws (e.g. Lawrence & Papadakis 1993; Uttley et al. 2002).

The studies by Hovatta et al. (2007; 2008) used radio data in the frequency range 22−87 GHz spanning timelines of  ≈25 years. A key result of their study was that AGN show large-amplitude activity on timescales  ≳6 years. Therefore, analysing the long-term behaviour of AGN requires monitoring on timescales of  ≳10 years. Building up on this conclusion, we have analyzed archival flux monitoring data from the IRAM Plateau de Bure Interferometer (PdBI) of six AGN – 0923+392, 3C 111, 3C 273, 3C 345, 3C 454.3, and 3C 84 – in the (observatory frame) frequency range 80−267 GHz. Given the range of redshifts z ≈ 0.02−0.86 of our sample, this corresponds to a range of rest-frame frequencies of  ≈82−497 GHz. Our dataset comprises several thousand measurements for each source, covering timelines of  ≈14 years. We note that our study complements earlier AGN mm/radio surveys carried out with the IRAM Pico Veleta observatory (Steppe et al. 1988; 1992; 1993; Reuter et al. 1997).

2. Observations and data processing

The PdBI is able to observe in three atmospheric windows located around wavelengths of 1.3 mm, 2 mm, and 3 mm. Each of these bands covers a continuous range of frequencies that are available for observations; these ranges are 201−267 GHz for the 1.3 mm band, 129−174 GHz for the 2 mm band, and 80−116 GHz for the 3 mm band (Winters & Neri 2010).

Until the end of 2006, observations were made using dual-band, single-polarization receivers. One frequency in the 1.3-mm band and one frequency in the 3-mm band could be selected for simultaneous observations; the 2-mm band was not reachable. At the beginning of 2007, new single-band, dual-polarization receivers were installed. Since then it is possible to observe a single frequency located in a given band simultaneously in both linear polarizations. At the beginning of 2008, the 2-mm band became available for observations.

Flux monitoring data for AGN are available from October 1996 on as a by-product of regular observatory operations. The PdBI uses the most luminous (S_ν ≳ 3 Jy in the 3-mm band) active galactic nuclei as calibrators for the spectral response of the instrument. The sources we present here are the most luminous targets regularly available (especially in terms of declination) as spectral calibrators. Therefore at least one of these sources has been observed in almost any individual observing project, with each project lasting typically few hours. This leads to a dense sampling over the entire observation timeline of  ≳13 years; gaps in the lightcurves are due to technical works at the interferometer and seasonal sunavoidance constraints. For our final sample of sources to be analyzed we picked all objects for which we have at least 1000 measurements over the entire observing timeline without major (≳1 year) interruptions of the monitoring. These were 0923+392, 3C 111, 3C 273, 3C 345, 3C 454.3, and 3C 84.

Antenna temperatures are converted into physical flux densities by using empirical antenna efficiencies as conversion factors. Those factors are functions of frequencies and are located in the range from  ≈22 Jy/K (for the 3-mm band) to  ≈37 Jy/K (for the 1.3-mm band). In order to estimate the systematic uncertainties of our dataset, we analyzed observations of the radio continuum star MWC 349A (e.g. Tafoya et al. 2004) that were obtained and calibrated like the AGN measurements. Due to its well-known flux properties MWC 349 serves as a “standard candle” for the PdBI (Krips et al., in prep.). From the systematic scatter of the MWC 349 lightcurves we conclude that the systematic relative uncertainties of our AGN observations are  ≈10% for the 2-mm and 3-mm bands and  ≈16% for the 1.3-mm band. We note that the systematic accuracy of our data is limited by the automatic monitoring and calibration; when calibrating observations individually and interactively while using MWC 349 as reference, systematic errors can be reduced to  ≲5% (Krips et al., in prep.).

Before analysis, we had to remove erroneous values from the flux data. Errors are caused by either of two effects: (1) poor signal-to-noise ratio (S/N); or (2) failure of the flux calibration. We caught effect (1) by rejecting data with statistical flux errors larger than 0.1 Jy or S/N  <  10. Statistical errors of individual measurements used for analysis are thus  <0.1 Jy. Effect (2) is more difficult to treat as it leads to a systematic deviation of the observed from the true flux values. Especially, it leads to a clustering of data points in the range 0−1.5 Jy (details depending on the source), well separate from, and systematically below, the long-term lightcurves. In order to remove this effect, we rejected flux measurements with fluxes located below the “divides” that separate the long-term lightcurves from the clusters of erroneous points. Additionally, we checked the statistical behaviour of the data. For this, we analyzed the differences between the lightcurves and smoothed versions of the lightcurves, using smoothing by running medians. We found that the distributions of the residuals usually follow Gaussian distributions with extended tails. We regarded the data points located in the tails as “suspicious” outliers and rejected all data located more than 2.5σ away from the centers of the distributions. This conservative approach rejected between 2% and 16% of the data of a given lightcurve.

Fig. 1 The long-term lightcurves of our six target quasars spanning timelines from 1996 to 2010. Please note the different flux axis scales. For each source we give the fluxes vs. time separatly for the 1.3-mm (here labeled 1-mm for brevity), 2-mm, and 3-mm bands. For each source (except 3C 111) we approximate the most prominent, continuously sampled variations as exponential rises or decays according to Sν ∝ exp(± t/tx) (orange lines). Here t is the time, tx the characteristic rise/decay timescale of the curve labeled x   (x = a,b,c,...). The sign of t is  + t for exponential rises and −t for decays. All parameters are given in the corresponding diagrams.

The data processing as well as the subsequent analysis made use of the IRAM Grenoble data processing software GILDAS¹ and the MPE Garching data processing software DPUSER².

After processing the data, we were left with large numbers of measurements ranging from 974 to 3027 per target, corresponding to one measurements every 1.5−4.3 days. Our flux data base covers a timeline from September 1996 (for the 2-mm band: February 2008) to August 2010. In Table 1 we provide an overview over our target AGN and individual observation statistics. We collected the source types and redshifts from the NASA/IPAC Extragalactic Database (NED) and the kilo-parsec scale morphology from the MOJAVE database (Lister et al. 2009).

3. Analysis

3.1. Lightcurves

In Fig. 1 we present the lightcurves of our six target AGN. For each source we give flux densities S_ν vs. time t separately for the 1.3-mm, 2-mm, and 3-mm bands as far as available.

In order to quantify the variability timescales, we approximated the most prominent – referring to smooth, continuously sampled variations in flux by factors  ≳2 – variations of each lightcurve as exponential rises or decays like (1)with t_e being the evolution timescale and the sign of the exponent being positive for exponential rises and negative for decays. We attempted those approximations only for parts of lightcurves that were sufficiently smooth for being described by a single analytical function and well sampled. These limitations prevented us from applying this approach to the lightcurve of 3C 111 which shows a combination of irregular sampling and fast variability. For the other sources, we give the results in the corresponding diagrams of Fig. 1.

3.2. Flux distributions

In Figs. 2 and 3 we present for each source the flux distributions for the 1.3-mm and 3-mm bands. The 2-mm lightcurves do not contain enough data for calculating corresponding distributions. In order to prevent distortions introduced by irregular sampling, we binned our data in time. The size of the time bins for a specific lightcurve is given by (2)Here T is the total observation period, N is the number of data points. For the special case of regular sampling, Δt corresponds to the inverse of the Nyquist frequency.

Fig. 2 Distributions of the 3-mm (left column) and 1.3-mm band (right column) fluxes of 0923+392, 3C 111, and 3C 273. Please note the different axis scales and bin sizes. In order to minimize the impact of irregular sampling, we binned all lightcurves (compare Fig. 1) in time using N bins with sizes Δt. The parameters N, Δt, and the sizes of the flux bins ΔSν are given in the respective diagrams. Errorbars indicate binomial errors. Continuous grey lines indicate approximations to the data composed of individual Gaussians (dashed grey lines). In three cases (0923+392 at 1.3 mm, 3C 111 at 1.3 mm and 3 mm) we also indicate the best-fitting log-normal distributions (dotted black curves).

Fig. 3 Same as Fig. 2 for 3C 345, 3C 454.3, and 3C 84.

As the resulting histograms show somewhat complicated morphologies, we investigated the minimum number of smooth components necessary for a reasonable description of their profiles. For this we approximated each distribution as a superposition of two (0923+392, 3C 111, 3C 345) or three (3C 273, 3C 454.3, 3C 84) Gaussian distributions (3)with S denoting the flux, S₀ being the center position, a₀ being the amplitude at S₀, and σ² denoting the variance. The choice of Gaussian distributions was governed by the fact that a priori such distributions are most likely to occur in case of uniform stochastic emission processes.

In three cases (0923+392 at 1.3 mm, 3C 111 at 1.3 mm and 3 mm) we also attempted to describe the flux histograms by log-normal distributions (4)where ln denotes the natural logarithm and σ and μ are the scale parameter and the location parameter, respectively. As log-normal distributions are characteristic for multiplicative stochastic processes³, they may be expected in AGN flux distributions (compare, e.g., Dodds-Eden et al. 2011).

3.3. Spectral indices

As we observed our target AGN over a range of frequencies located in three spectral bands, we were able to analyze their mm-spectroscopic properties. This analysis was however complicated by fast intrinsic variability as well as by the fact that the spectral information partially needed to be derived from non-simultaneous observations. We therefore inserted all measurements within selected time windows into flux-frequency diagrams and derived time-averaged spectral indices α via means of least-squares fits. Spectral indices were defined via the relation (5)For each source we present the spectral index as function of time in Fig. 4. We selected the time windows such that a given window covers a phase of similar flux levels, meaning a “quiet” phase, an outburst, or a low-flux regime between two outbursts.

In order to test possible correlations between spectral index and source flux, we present α as function of the 3-mm flux for each source in Fig. 5. For this diagram we use the same binning in time as for Fig. 4. In order to preserve some time information, we connect and label the data points according to their order in time.

Fig. 4 Evolution of the spectral indices from 1996 to 2010. Please note the different axis scales. The spectral index α is defined via Sν ∝ ν−α. Errorbars along the time axes denote the bin sizes, error bars along the α axes denote the statistical 1σ errors. Horizontal dashed lines (where visible) mark the α = 0 lines. For each source we chose the bins such that each bin covers a phase of similar flux levels (quiescent phases, flares).

Fig. 5 Evolution of the spectral index α as function of the 3-mm flux. The data points are connected by dotted lines according to their order in time. Labels 1−3, ... mark the first, second, third, ... measurement in time. The time bins used here are those defined for Fig. 4. Error bars along the flux axis denote the flux variability (standard deviation), error bars along the α axis denote statistical 1σ errors.

3.4. Power spectra

As the lightcurve of our sources show variability on many timescales, we performed a power spectrum analysis in order to quantify this variability. For each lightcurve we computed a normalized Scargle periodogram (6)(Scargle 1982) from our data. Here A_f is the periodogram amplitude evaluated at frequency f, S_j and t_j are the flux value and the time of the jth data point, respectively, and σ² denotes the variance of the data. The base frequencies are f = f_min, 2f_min, 3f_min, ..., f_max with f_min = 1/T, f_max = N/(2T), T being the total observation period, and N being the number of data points. For the special case of regular sampling, f_max corresponds to the Nyquist frequency. The timelines and sampling of our observations (see Table 1) lead to typical values of f_min ~ 2  ×  10^-4 days^-1 and f_max ~ 0.2 days^-1. Before calculating a periodogram we subtracted the average source flux from the data in order to remove the zero-frequency power. In order to exclude systematic errors due to spectral effects, we performed the periodogram analysis separately for the 1.3-mm and 3-mm bands of each source (for the 2-mm band the amount of data is insufficient for a comparison to the other bands).

The power spectra we find all show a global decrease of spectral power with increasing frequency, i.e. they correspond to red-noise spectra. In order to quantify this behaviour, we derived powerlaw indices by fitting straight lines to the data in log-space. The spectral index β of the power spectra is defined via A_f ∝ f^−β.

In several cases we find that the low-frequency ends of the power spectra follow the average powerlaw slopes but are located systematically above the best-fitting powerlaw models. In those cases we approximated the affected spectral ranges as upscaled (by factors  ≈2−3) versions of the overall power spectra. We present our results in Figs. 6 and 7.

Fig. 6 Scargle periodograms (spectral power Af vs. frequency f) of the 3-mm (left column) and 1.3-mm band (right column) lightcurves of 0923+392, 3C 111, and 3C 273. Please note the different axis scales. The black curves denote the measured power spectra, the continuous grey lines marks the best-fitting powerlaws. The spectral index β of the power spectra is defined via Af ∝ f−β. Dashed grey lines (where present) correspond to the best-fiting powerlaws (continuous grey lines) upscaled by factors s as given in the corresponding plots.

Fig. 7 Same as Fig. 6 for 3C 345, 3C 454.3, and 3C 84.

3.5. Time offsets among spectral bands

The long timeline of observations and the good time sampling of our flux data make it possible to test if the lightcurves of the various spectral bands are simultaneous or if time lags are present. In order to quantify the presence or absence of time offsets, we computed for each source a discrete correlation function (DCF) according to the scheme proposed by Edelson & Krolik (1988). For two discrete datasets  {a_i} ,  {b_j}  with averages , one computes the unbinned discrete correlations (7)for all pairs (a_i,b_j) with time differences Δt_ij. The are the variances, the δ_a,b are the mean measurement errors of  {a_i} ,  {b_j} . Accordingly, the terms  are the variances corrected for artificial broadening due to finite measurement errors; this is necessary in order to preserve the correct normalization in noisy data. The actual DCF for a time offset τ is obtained via averaging over all N UDCF_ij(Δt_ij) for which Δt_ij is located within a selected bin Δτ (i.e., τ − Δτ/2 ≤ Δt_ij < τ + Δτ/2): (8)By definition, DCF(τ) is located in the range  [−1, + 1] , with (−1) +1 indicating perfect (anti-)correlation and 0 corresponding to the absence of any correlation. The position of the maximum of the DCF indicates the time offset between the lightcurves. The statistical uncertainty of DCF(τ) is given by the standard error of mean (9)For each of our six target AGN, we examined ranges of τ sufficient for covering the largest time offsets to be expected, either τ  ±  365 days or τ  ±  730 days. In the special case of 3C 454.3, we additionally probed two subsets of data with τ  ±  180 days (see Sect. 4.5). In order to preserve high time resolution while excluding the risk of creating oversampling artefacts, we used sampling times Δτ = T_1   mm/N_1   mm, with T_1   mm (N_1   mm) being the total observing time (number of measurements) obtained in the 1.3-mm band which is for all sources the wavelength band with the sparser sampling (see Table 1). We note that the effective resolution depends only weakly on the choice of bin sizes: smaller Δτ increase the formal resolution but also increase the errors as they cover smaller numbers of UDCF_ij per bin (see Eq. (9)), and vice versa.

We present the results of our analysis in Fig. 8. In our convention, a positive (negative) time lag means that the 3-mm lightcurve preceeds (follows) the 1.3-mm lightcurve. In each diagram, we indicate the ranges (as grey-shaded areas) covered by the range [max(DCF),max(DCF) − 3σ_max(DCF)], thus covering all values of DCF that are in agreement with the maximum of the DCF on a 3σ confidence level. We probe the null hypothesis “the 1.3-mm and 3-mm lightcurves are simultaneous”. Accordingly, we regard the lightcurves as simultaneous as long as DCF(τ = 0) ≥ max(DCF) − 3σ_max(DCF).

Fig. 8 Discrete correlation functions (DCF) as functions of time lags τ between the 1.3-mm and 3-mm lightcurves (see Fig. 1) of each AGN. The black graphs indicate the measured DCF, bin sizes Δτ are indicated in the corresponding diagrams. The maxima of the curves correspond to the most probable time lag between the spectral bands. Grey shaded areas mark the DCF ranges  [max(DCF),max(DCF) − 3σ]  with σ being the statistical error of max(DCF). Thus all values located within the shaded area are consistent with the maximum of the DCF curve within 3σ confidence. A positive (negative) time lag means that the 3-mm lightcurve preceeds (lags behind) the 1.3-mm lightcurve. Vertical dashed lines mark the null positions of the time lag axes, corresponding to the absence of time lags.

4. Results

4.1. 0923+392 (4C +39.25)

The 3-mm (1.3-mm) flux densities of this flat spectrum radio quasar (FSRQ) dropped from  ≈6 Jy (≈3 Jy) in 1996 to  ≈3 Jy (≈1.5 Jy) in 2003 and have been rising again since then to present-day levels of  ≈4.5 Jy (≈3 Jy). For the flux density decrease from 1996 to 2003 we find characteristic (exponential) timescales of  ≈5−7 years. The 3-mm flux distribution shows a double peak profile that can be described by two Gaussian components separated by  ≈1.6 Jy. The 1.3-mm flux distribution shows an asymmetry indicative of a second component but no significant secondary peak. We therefore attempted to describe the histogram by a log-normal distribution according to Eq. (4). However, we did not find a satisfactory solution, with the best fitting curve showing a χ²/d.o.f.  =  13.8.

The spectral index α is  ≈0.8 in the first years of monitoring before 2000 and steepens to  ≈1 around 2002. The index then flattens again until reaching  ≈0.6 in the time window 2007−2010. There appears to be no correlation between spectral index and 3-mm flux. Even though we find the steepest spectral indices for the lowest source fluxes, we find similar variations in α (by  ≈0.2) within the same flux ranges.

The 1.3-mm and 3-mm lightcurves are simultaneous within the 3σ confidence interval. The power spectra of the 3-mm and 1.3-mm flux measurements can be well described as simple red-noise spectra with fairly flat powerlaw indices β ≈ 0.4. Within the errors, the values for β are in agreement for both frequency bands.

4.2. 3C 111

This Seyfert 1 galaxy shows episodic strong variability on timescales of years. The source seems to spend most of the time in a low-activity or “quiescent” state with fluxes in the range  ≈1−3 Jy. In addition, we can identify at least five pronounced spikes in the lightcurves in 1997, 2001, 2006, 2008, and 2009. These outbursts reach flux densities of  ≈6−15 Jy, thus exceeding the quiescent level by factors 3−6. The outbursts appear to rise and decay on fast (≲1 yr) timescales; unfortunately, the irregular sampling of the lightcurves prevents a reliable quantitative analysis. The 3-mm (1.3-mm) flux distributions can be described by smooth, asymmetric profiles composed of two Gaussian components for fluxes  ≲8 Jy (≲5 Jy) plus sparsely populated tails ranging out to  ≈14 Jy (≈10 Jy). Except for the tails, it is also possible to describe the histograms with (unimodal) log-normal distributions. In both cases we find best-fit solutions with χ²/d.o.f. ≈ 1.6. Although this is a relatively high value, the log-normal distribution appears to be an acceptable description of the data. Within errors, the 1.3-mm and 3-mm lightcurves are simultaneous.

For most of the monitoring time the spectral index remains at α ≈ 0.5 (within errors). The only exception appears around 1997 when the spectral slope is α ≈ 0.8; this phase coincides with the first outburst. In the flux-α diagram the various points do not indicate a trend.

The power spectrum of the 3-mm (1.3-mm) flux data globally follows a red-noise powerlaw spectrum with β ≈ 1.0 (β ≈ 0.6). The low-frequency end (f ≈ 0.002−0.01 days^-1) of the 3-mm power spectrum is located systematically above the best-fitting overall powerlaw; this region can be described as an upscaled (by a factor  ≈2) version of the global red-noise spectrum.

4.3. 3C 273

At the beginning of the monitoring, the 3-mm (1.3-mm) lightcurves of this FSRQ rose fast towards a strong outburst peaking in 1997 at levels of  ≈38 Jy (≈30 Jy). From then on, the flux densities steadily decreased down to minimum values of  ≈6 Jy (≈4 Jy) in 2004 with an exponential decay timescale of  ≈3 yr (measured for the decrease in 1997−2000). Since about 2006, 3C 273 has been in a state of high variability, with fluxes fluctuating in ranges  ≈10−27 Jy (≈6−20 Jy). The flux distributions follow complex, double-peaked profiles; they can be approximated by superpositions of three Gaussian distributions peaking at  ≈8 Jy,  ≈18 Jy, and  ≈28 Jy (≈4 Jy,  ≈10 Jy, and  ≈16 Jy). Within errors, the 1.3-mm and 3-mm lightcurves are simultaneous.

The spectral index varies substantially over the 14 years of observation, starting at α ≈ 0.3 around 1997, steepening up to α ≈ 0.8 in 2000−2003 and flattening again to α ≈ 0.4 after 2007. There is a global trend towards flatter indices at higher flux levels, though no clear correlation: on the one hand α remains at similar levels over wide flux ranges, on the other hand strong variations in α occur within the same flux ranges.

The power spectra of the 3-mm and 1.3-mm lightcurves follow red-noise powerlaws with identical (within errors) indices β ≈ 0.45. For frequencies f ≲ 0.04 days^-1 the spectral power is located systematically above the best-fitting overall powerlaw; these ranges can be described as varieties of the overall power spectra upscaled by factors  ≈3 (3 mm) and  ≈2 (1.3 mm), respectively.

4.4. 3C 345

The 3-mm (1.3-mm) lightcurves of this FSRQ start at flux densities of  ≈4 Jy (≈2 Jy) in 1997 and reach their highest values of  ≈7 Jy (≈4 Jy) in mid-1998. From then the fluxes decrease overall down to  ≈3 Jy (≈2 Jy) in the time range 2004−2008. Around 2009 an outburst occurs, reaching flux levels up to  ≈6.5 Jy (≈4.5 Jy) and slowly decreasing since then. The 2009 outburst raises in 2008−2009 with an exponential timescale of  ≈1 yr. The 3-mm (1.3-mm) flux distributions follow double-peaked profiles we approximated as superpositions of two Gaussians peaking at  ≈3 Jy and  ≈5.5 Jy (≈2 Jy and  ≈3.5 Jy). Within errors, the 1.3-mm and 3-mm lightcurves are simultaneous.

Fig. 9 Like Fig. 8, for 3C 454.3 before 2005.2 (left panel) and after 2005.2 (right panel). During the “quiescent” phase before 2005.2, the peak of the DCF is in good agreement with a time delay of zero. We note the low absolute value of the correlation, DCF(0) ≈ 0.6. For the “flaring” phase after 2005.2, we find a significant time delay in the range (as defined by the 3σ confidence interval) of  ~ −15...−80 days, peaking at τ ≈ −55 days. The negative delay indicates that the 1.3-mm lightcurve is ahead in time of the 3-mm lightcurve.

Fig. 10 A visualization of the time delay between the 1.3-mm and 3-mm lightcurves of 3C 84. Left panel: the 1.3-mm and 3-mm lightcurves. Each lightcurve is normalized to zero mean and unity standard deviation, in analogy to the calculation of the unbinned discrete correlations according to Eq. (7). Here Sν indicates the flux density,  ⟨ Sν ⟩  is the average flux density, and σ denotes the standard deviation of each lightcurve. Right panel: same as left panel, but with the 1.3-mm lightcurve artificially retarded by 300 days. Comparison of the two panel shows improved agreement of the two lightcurves when taking into account the time delay τ ≈ −300 days found from the discrete correlation function; see also Fig. 7, bottom right panel.

Fig. 11 Spectral index α vs. time based on 1.3-mm and 3-mm flux data only, with the 1.3-mm lightcurves artificially retarded. Left panel: 3C 454.3, using a spectral time delay τ = −55 days. Right panel: 3C 84, using a spectral delay τ = −300 days. The time delays are given by the correlation analysis illustrated in Fig. 7. This figure should be compared to Fig. 4.

Within errors, the spectral index remains on very similar levels in the range α ≈ 0.6−0.8 for most of the monitoring time. After 2007 – the time frame of the most recent outburst – the spectrum flattens somewhat towards α ≈ 0.4. We note that this somewhat steep index at mm-wavelengths is consistent with the classification of 3C 345 as FSRQ: the classification as FSRQ is based on cm-observations, whereas at mm-wavelengths the spectrum bends towards α > 0. In a S_ν − α diagram the data points move with time on an almost circular trajectory, leaving no room for a correlation between flux and spectral index.

The power spectra can be described as simple red-noise spectra, with powerlaw indices β ≈ 0.5 for the 3-mm and β ≈ 0.4 for the 1.3-mm lightcurve, respectively.

4.5. 3C 454.3

This flat spectrum radio quasar shows the most spectacular lightcurves among the sources of our sample. For the first nine years (1996−2005) of monitoring the activity of this source was restricted to variability on levels of  ≈2−8 Jy. This changed dramatically in 2005 when the first of at least three major outbursts occured, peaking at flux densities of  ≈35−40 Jy. The series of adjacent outbursts has continued until the present day, making 3C 454.3 one of the brightest extragalactic mm/radio sources currently observed (see also, e.g., Jorstad et al. 2010). All outbursts show short exponential rise and decay times in the range  ≈0.1−0.5 years. The 3-mm (1.3-mm) flux distribution is dominated by a Gaussian peak centered at  ≈6 Jy (≈4 Jy) with an extended tail ranging out to  ≈40 Jy (≈35 Jy).

The spectral index has been fluctuating around α ≈ 0, with the most extreme values being α ≈ −0.4 and α ≈ 0.5, respectively. A somewhat steep value of α ≈ 0.5 can be observed for the time 1996−2000. Notably, the most extreme change (from α ≈ −0.4 to α ≈ 0.4) occurs at the time of the peak of the first bright outburst in 2005. The data almost fill the entire S_ν − α area observed, leaving the impression that flux and spectral index are uncorrelated.

The 3-mm (1.3-mm) power spectrum globally follows a red noise powerlaw with index β ≈ 0.7 (β ≈ 0.5). However, for frequencies f ≲ 5  ×  10^-3 days^-1 the power spectrum seem to be an upscaled – by a factor  ≈3 (≈2) – version of the average powerlaw.

The discrete correlation function points to time delays τ significantly different from zero in terms of the 3σ confidence threshold, located in the range τ ≈ −20...−110 days. The negative delay indicates that the 1.3-mm emission preceeds the 3-mm emission. However, as 3C 454.3 shows a dramatic transition from “quiescent” to “flaring” around epoch 2005.2, the curve presented in Fig. 8 (bottom left) might be misleading as the timing of the emission might be different in the two states. We therefore computed in addition one DCF for each of the two activity phases, one covering all data from 1996.8 to 2005.1, one including all measurements obtained between 2005.2 and 2010.6. The results are summarized in Fig. 9. We note that for the “quiescent” phase (before 2005.2) the most probable time delay – on a low absolute level, DCF(0) ≈ 0.6 – is consistent with zero, whereas we find for the “flaring” phase (after 2005.2) τ ≈ −15...−80 days, with the maximum of the DCF located at τ ≈ −55 days. Again, the negative delay indicates that the 1.3-mm emission preceeds the 3-mm emission.

Finally, we probed the impact of the spectral time delay on the calculation of the spectral index α. For this, we repeated the calculation outlined in Sect. 3.3 using 1.3-mm and 3-mm flux density data only. We compensated the spectral time delay by artificially retarding the 1.3-mm lightcurve by 55 days. We present the resulting modified spectral index evolution in Fig. 11 (left panel). Time bins are the same as in Fig. 4 (bottom left panel) which should be compared to Fig. 11 (left panel). We note that retarding the 1.3-mm data removes the inverted values for α that are visible in Fig. 4 (bottom left panel) during the outburst 2005/2006.

4.6. 3C 84

The 3-mm (1.3-mm) lightcurve of this Seyfert 2 galaxy varies smoothly during the entire monitoring time between  ≈3.5 Jy (≈2 Jy) in 2002 and  ≈11 Jy (≈7 Jy) in 2010. Since about 2004 the fluxes raise continuously with exponential rise time scales of  ≈6 years (≈4 years). The flux distributions show somewhat complicated double-peak profiles that can be approximated as superpositions of three Gaussians peaking at  ≈4 Jy,  ≈5.5 Jy, and  ≈9 Jy (≈2 Jy,  ≈4 Jy, and  ≈6 Jy).

The spectral index is steep, with α ≈0.9 (within errors) until about 2005 and then drops smoothly to α ≈ 0.5 after 2007. We find the flattest spectral index for the highest flux value; however, as this is based on a single data point that occupies a large range in flux, this observation does not establish a significant trend.

The power spectra can be described as simple red-noise spectra, with powerlaw indices β ≈ 0.7 for the 3-mm and β ≈ 0.4 for the 1.3-mm lightcurve, respectively.

The discrete correlation function computed from the 1.3-mm and 3-mm lightcurves finds a spectral time delay significantly different from zero. The most probable time delay is located at τ ≈ −270...−330 days, meaning that the 1.3-mm emission preceeds the 3-mm emission by almost one year. We note the high absolute value of the maximum correlation, DCF(τ = −277   d) = 1.04  ±  0.05, which is in agreement with the theoretical maximum value of 1. However, as already DCF(0) ≈ 0.7, one might wonder how large the difference between τ ≈ 0 and τ ≈ −300 actually is in absolute terms. We therefore tested the impact of the time offset graphically. For this we normalize the 1.3-mm and 3-mm lightcurves to zero mean and unity standard deviation, in analogy to the computation of the unbinned discrete correlations according to Eq. (7). In Fig. 10 we accordingly present the quantity (S_ν− ⟨ S_ν ⟩ )/σ for each of the two frequency bands; here S_ν is the flux density,  ⟨ S_ν ⟩  denotes the time average of the flux density, and σ indicates the standard deviation of the flux density. In Fig. 10 we present the normalized data separately for the two cases τ = 0 days (left panel) and τ = −300 days (right panel). Indeed the case τ = −300 days shows a better agreement between the 1.3-mm and 3-mm lightcurves than the case τ = 0 days. This agrees with the result obtained from the correlation analysis. However, we note that the lightcurves show clear trends that are only partially covered by our observations. This can mislead a correlation analysis, meaning that the time delay we find has to be taken with care.

Like for 3C 454.3, we probed the impact of the spectral time delay on the calculation of the spectral index. We compensated the spectral time delay by artificially retarding the 1.3-mm lightcurve by 300 days. We present the resulting modified spectral index evolution in Fig. 11 (right panel). Time bins are the same as in Fig. 4 (bottom right panel) which should be compared to Fig. 11 (right panel). Comparison of the two diagrams indicates that removing the spectral time delay systematically shifts α towards steeper values by  ≈0.1.

5. Discussion

5.1. Lightcurves and flux distributions

Our analysis reveals substantial – by factors  ≈2−8 – flux variability on timescales of several years in all sources of our sample. Even though all lightcurves (see Fig. 1) show strong variations, the characteristic (exponential rise/decay) timescales differ by more than one order of magnitude, ranging from  ≲0.5 years (3C 454.3) to  ≈7 years (0923+392). Accordingly, we find lightcurves with smooth morphologies during the entire  ≈14-year monitoring time (0923+392, 3C 84) as well as lightcurves that are “spiky” (3C 111, 3C 454.3).

The flux distributions (Figs. 2, 3) of five out of six (the exception being 3C 111) AGN show bi- or even multimodality. This indicates that the flux variability of the concerned sources cannot be described by uniform stochastic processes:

• In case of Gaussian (white) noise processes the distribution ofamplitudes would be (obviously) Gaussian.

• In case of red noise processes (see also Sect. 3.4), the distribution of amplitudes is Gaussian for a given sampling frequency interval  [f,f + df]  (e.g. Timmer & König 1995). The distribution of amplitudes resulting when summing over all f would thus a priori be a superposition of Gaussians of varying widths (i.e. a unimodal, symmetric distribution). However, for physical lightcurves one needs to take into account that fluxes are constraint to positive values. This introduces a bias towards higher flux values, meaning that the resulting flux histogram shows an extended wing towards high values.

• In case of multiplicative stochastic processes the distribution of amplitudes would follow a unimodal log-normal distribution, i.e. an asymmetric distribution with an extended wing towards high values (see Uttley et al. 2005, for a discussion of AGN X-ray lightcurves).

The flux distributions of 3C 111 are (at least marginally) consistent with log-normal distributions, pointing towards an underlying multiplicative process. Such processes can be understood as chains of events where each event is triggered by the previous one with a certain probability  <1 (e.g., Montroll & Shlesinger 1982). Candidates are outbursts that trigger secondary emission, or shocks in jets that trigger secondary shocks or fragmentation of emission regions.

For the remaining five sources, it appears to be necessary to distinguish separate states of source activity like “quiescent” and (maybe multiple) “flare states” where each state is characterized by an individual flux distribution. The final, observed flux distribution would then be a superposition of the individual distributions.

5.2. Spectral indices

All sources of our sample show (mostly) fairly steep spectral indices within ranges α ≈ 0.5−1 (Fig. 4). All sources show a significant, smooth variability of their indices on timescales of years. Only 3C 454.3 provides an exception from both of the aforegiven statements: it shows a faster variability on timescales of  ≈1−2 years and reaches flat and even inverted indices α ≈ −0.5−0. As our sample is composed of very luminous AGN, we can safely assume that their mm-emission is dominated by synchrotron radiation. We can thus classify the sources of our sample (with the possible exception of 3C 454.3) as outflow (i.e. halo or jet) dominated sources for which typically α ≈ 0.5−1 (e.g. Krolik 1999). This is consistent with the fact that for all sources luminous halos or jets on kilo-parsec scales have been observed (see Table 1).

The variations of the spectral indices point towards variations of intrinsic source properties on timescales of years. A priori, there are three potential mechanisms:

• 1.

  Variations of the optical depth. Synchrotron emission fromoptically thin emission regions shows steep spectral indices withα ≳ 0.5; strong absorption in emission regions (i.e. optically thick emission) leads to approximately flat (α ≈ 0) or even inverted, complex spectra (see, e.g., Kembhavi & Narlikar 1999, for a review). Accordingly, variations of the optical depth can lead to variations of the spectral index inbetween these extreme cases.

• 2.

  Variations of the relative luminosities of different source components. Emission from (optically thick) nuclei shows usually fairly flat spectra (α ≈ 0) whereas the emission from (optically thin) outflows (jets, lobes) typically shows steep spectra (α ≈ 0.5−1; e.g. Krolik 1999). Variations in the relative luminosities of the various components therefore modify the spectral index of the source-integrated emission. The spectrum becomes steeper if the relative contribution by outflows increases, it becomes flatter if the relative contribution by the nucleus increases.

• 3.

  For optically thin synchrotron emission from an ensemble of electrons with a powerlaw energy spectrum with index γ the spectral index of the radiation is given by α = (γ − 1)/2 (e.g. Ginzburg & Syrovatskii 1965); for α = 0.5−1 this corresponds to γ = 2−3. Accordingly, variations in the electron energy spectrum express themselves as variations of the spectral index.

For three of our six target sources (namely 3C 111, 3C 273, and 3C 454.3) the spectral indices are – at least occasionally – significantly flatter than the range permitted for the assumption of optically thin emission (i.e. we observe α < 0.5) for parts of the monitoring time. This suggests that scenario 3 cannot be applied regularly as this scenario is valid only for optically thin emission.

Additional information is provided by the relation between source fluxes and their spectral indices (Fig. 5). Especially, the S_ν − α relation reveals that scenario 2 can play only a minor role for the variability of the spectral index: if the α variations were due to (or at least dominated by) luminosity increases or decreases of parts of the sources, there should be a tight correlation between source fluxes and spectral indices – which we do not observe. From the combined information we therefore conclude that variations of the optical depth with time (scenario 1) are the most likely mechanism for variations of the spectral index.

5.3. Power spectra

The power spectra of our sources (Figs. 6, 7) globally follow red noise distributions, meaning powerlaws with spectral indices β > 0 (cf. Sect. 4.3). This, as well as the apparent absence of periodic signals, shows that the emission processes are intrinsically stochastic. The power-spectral indices β ≈ 0.4−1 are located between those of Gaussian (white) noise (β = 0) and flicker noise (β = 1) processes (e.g. Timmer & König 1995). They are significantly flatter than the indices found from optical and cm/radio (e.g. Hufnagel & Bregman 1992) as well as X-ray (e.g. Lawrence & Papadakis 1993) studies of AGN emission which are located in the range 1−2.

When discussing the power spectral indices we have to take into account that we observe lightcurves that are affected by white (Gaussian) measurement noise. This noise is dominated by a systematic scatter that amounts to  ≈10% for the 2-mm and 3-mm fluxes and  ≈16% for the 1.3-mm flux data (see Sect. 2). The presence of white measurement noise leads to the observation of power-spectral indices β that are actually flatter (i.e. closer to the case of white noise, β = 0) than the indices intrinsic to the source emission. In order to estimate the intrinsic spectral indices, we performed Monte Carlo tests using simulated red-noise lightcurves with known ideal spectral indices. For this test we applied the algorithm by Timmer & König (1995). We calculated artificial lightcurves corresponding to red noise lightcurves with β = 2, β = 1.5, β = 1, and β = 0.5, respectively. For each case, we examined three scenarios: the absence of measurement noise, Gaussian measurement noise with σ = 0.1   S_ν (i.e. a Gaussian width corresponding to 10% of the observed fluxes), and Gaussian measurement noise with σ = 0.16   S_ν (i.e. a Gaussian width corresponding to 16% of the observed fluxes). We simulated 100 lightcurves with 500 data points each for each of the 12 cases examined. In order to assess the impact of gaps in the data, we applied a sampling approximately corresponding to the worst case actually encountered which is the sampling of the 1.3-mm lightcurve of 3C 111 (compare Fig. 1, top right panel). For this, we inserted – by simply leaving out data – a total of eleven artificial gaps, each covering  ≈2% to  ≈9% of the total observing time. For each artificial lightcurve we calculated a Scargle periodogram (Eq. (6)) and derived the powerlaw index. We present the results in Table 2. Due to the stochastic nature of the noise signals, the “observed” values for β show substantial scatter around the averages; the 1σ widths of the corresponding distributions are quoted as uncertainties in Table 2.

From our test we can draw various conclusions:

• The presence of Gaussian measurement noise flattens thepower-spectral indices toβ_obs < 1.5 even if β_true ≈ 2.

• As β_obs ≈ 0.4−0.7 for all sources except 3C 111, the intrinsic power-spectral indices should be flatter than 1 in general.

• Increasing the amplitude of measurement noise from 0.1   S_ν (like for our 3-mm flux data) to 0.16   S_ν (like for our 1.3-mm flux data) leads to a further flattening of the power-spectral index by about 0.1. This agrees with our observations: we find a systematic trend towards flatter (by  ≈0.1) indices for the 1.3-mm data compared to the 3-mm power spectra of the corresponding sources.

This also shows that the difference between our results and those reported by Hufnagel & Bregman (1992) or Uttley et al. (2002) cannot be attributed to measurement noise effects. Instead, we have to attribute this difference to the different frequency regime (and thus emission region) we observe. This indicates that at millimeter wavelengths, adjacent AGN emission events tend to be less correlated than in other regimes.

The fact that intrinsically β < 1 (except for 3C 111) points towards stochastic emission processes located between the cases of white noise and flicker noise, with 3C 111 being closer to the case of random walk noise (β = 2). Formally, the various cases are connected by different levels of correlation among adjacent data. In case of white noise (β = 0) adjacent data are completely uncorrelated. Random walk noise (β = 2) corresponds to the integral of white noise over an infinite interval in time, meaning that adjacent data are fully correlated. Intermediate cases (0 < β < 2) can be expressed as “half” integrals of white noise (again, over infinite intervals in time). Such a “half” integration is performed by convolving a white noise time series with a powerlaw kernel G ∝ t^−ξ with ξ > 0 and t being the time. For the specific case of flicker noise (β = 1) this was discussed in detail by Press (1978), leading to ξ = 1/2. From this, one may a priori expect to find any 0 ≤ β ≤ 2 in AGN time series driven by stochastic emission processes⁴.

Even though there is a straightforward mathematical recipe for the creation of red noise (in the range 0 ≤ β ≤ 2) time series, there is not yet a consistent understanding which physical mechanisms create those time series. This is a problem of importance not only in astrophysics but in many fields of physics, and various approaches towards a general explanation have been proposed (e.g., Press 1978; Pang et al. 1995; Makse et al. 1996). More recently, Kaulakys & Alaburda (2009) and Kelly et al. (2011) proposed time series models based on superpositions of exponentially decaying outbursts occuring at random times and with random amplitudes. Both studies are able to reproduce the occurence of powerlaw power spectra. However, as already noted by Press (1978) who discussed (and eventually rejected) similar proposals, those approaches might only shift the problem from the question “What process creates those power spectra?” to the question “What process generates those outburst sequences?” without actually solving the problem.

Although the power spectra globally follow red noise laws, three out of six sources (3C 111, 3C 273, and 3C 454.3) are hardly consistent with a single stochastic emission process. Instead, we observe excess (by factors  ≈2−3) spectral power with respect to the average powerlaw for sampling frequencies f ≲ 0.005-0.01 days^-1. In combination with the information from the flux distributions (Sect. 5.1; Figs. 2, 3) this strongly suggests that “quiescent” high-frequency and “flaring” low-frequency emission have to be treated as separate source states. This conclusion deviates from the more simple picture drawn by other studies (e.g., Hufnagel & Bregman 1992; Uttley et al. 2002) that treat AGN time series and power spectra as uniform red-noise distributions. However, a physical dichotomy of quiescent vs. flaring source states has been suggested for the near-infrared emission observed from the Galactic nuclear black hole, Sagittarius A* (e.g. Gillessen et al. 2006; Trippe et al. 2007; Dodds-Eden et al. 2011), and has been hotly debated since then (compare, e.g., Meyer et al. 2007 vs. Meyer et al. 2008).

We also note that a segregation of distinct source states is well-known for the case of Galactic microquasars. These sources mirror many properties of AGN (taking into account scaling with the black hole mass) and their emission shows distinct states that are ususally classified as “hard” and “soft” according to their X-ray spectral properties (e.g., Tananbaum et al. 1972; Mirabel & Rodriguez 1998; 1999; Takeuchi & Mineshige 1998). A possible explanation for this transition (Klein-Wolt et al. 2002) is the difference between the ejection of fairly localized knots from the black hole vicinity (leading to X-ray and radio flares) and the ejection of quasi-continuous, extended jets that are initially optically thick (leading to “plateaus” in the X-ray and radio lightcurves). For micro-quasars, these events happen on timescales from minutes (minimum time between the ejection of knots) to tens of days (maximum length of “plateau” phases observed); for AGN, corresponding phenomena would have to occur on timescales of years – thus deviating strongly from a simple linear black-hole mass scaling relation – in order to be consistent with our observations.

Notably, the “flare” and “high-frequency” regimes of the power spectra are consistent with having the same powerlaw index though on different (by factors  ≈2−3 in spectral power) energy levels. This indicates that both emission regimes, although they appear to be physically disconnected, are caused by intrinsically stochastic processes with similar correlation properties. Remarkably, the three sources in question – 3C 111, 3C 273, and 3C 454.3 – are also those for which we find the (temporarily) flattest spectral indices, indicating optically thick emission.

5.4. Spectral time delays

We tested the presence or absence of time lags between the 1.3-mm and 3-mm fluxes by means of a correlation analysis (cf. Eq. (8)). We can exclude the presence of significant – referring here to a 3σ confidence level – offsets in time for four of our targets, namely 0923+392, 3C 111, 3C 273, and 3C 345 (see Fig. 8). This statement holds within accuracies (referring to the ranges of timelags consistent with a zero offset within the 3σ confidence band) of few tens of days. Our results are consistent with those obtained by the time delay analyses of Tornikoski et al. (1994) and Hovatta et al. (2008) in the cm-radio range. Hovatta et al. (2008) conclude that those limits on time delays, both in the cm- and mm-regimes, can be understood in the frame of the semi-empiric “generalized shock model” by Valtaoja et al. (1992).

For the remaining two sources, 3C 454.3 (after 2005.2; cf. Fig. 9) and 3C 84, we do find indication that the 1.3-mm emission preceeds the 3-mm emission by  ≈55 days (3C 454.3) and  ≈300 days (3C 84). We note that the result for 3C 454.3 is consistent with the results by Tornikoski et al. (1994) and Hovatta et al. (2008) who do not find a time delay in their radio lightcurves; the data of Hovatta et al. (2008) end in 2005.3, just at the beginning of the “flare phase” of 3C 454.3.

For 3C 454.3, we suspect variations of the optical depth with frequency as the principal cause for the time delay. Adiabatically expanding plasmas become optically thin at higher energies first, introducing a systematic time delay between two lightcurves obtained at two frequencies with the high-frequency emission leading. Details depend on various parameters, notably the electron energy distribution, the optical depth at a reference frequency, the magnetic field strength, the electron density, and the emission region size (van der Laan 1966). This mechanism has been incorporated into the shock-in-jet model by Marscher & Gear (1985); proper quantitative treatment takes into account the combination of Compton, synchrotron, and adiabatic cooling (e.g., Fromm et al. 2011, and references therein). Our interpretation is supported by the fact that the spectral index shows fast variability from inverted to almost optically thin spectral indices in 2005−2007 and remains in the optically thick regime after 2007⁵. This is in agreement with the results by Jorstad et al. (2010) who find time delays up to  ≈215 days between 37 GHz and 230 GHz lightcurves, with the 230 GHz emission leading.

For 3C 84, the answer is less clear. High-resolution radio-interferometric observations let Asada et al. (2006) conclude that the lightcurve, as observed until 2001, is in agreement with adiabatic cooling of an expanding radio lobe. The rise of the emission around 2005 nicely agrees with observations of the ejection of a new radio jet component in 2005 (Nagai et al. 2010). Those observations are a priori in agreement with optical depth variations in the context of the shock-in-jet model. However, we note the unusually large time delay of  ≈300 days – assuming the result of the DCF analysis is correct. Such a large delay suggests additional travel time effects given by the geometry of the emission region. If the emission is composed of various components emitted from spatially separate zones, according time lags may occur. Given the observed time offsets of  ≈300 days, this would suggest a spatial extension of the emission region of the order of one light-year.

6. Conclusions

We have analyzed the long-term activity of six luminous active galactic nuclei using mm/radio lightcurves spanning  ≈14 years in time. From our analysis we draw the following principal conclusions:

• 1.

  All sources show substantial, non-periodic variability byfactors  ≈2−8 on timescales of years. This is in good agreement with the results reported by Hovatta et al. (2007; 2008) and shows that source monitoring times on scales of decades are necessary in order to obtain a realistic overview on the source activity. We also note that at least for the special case of Sagittarius A* both observational (e.g. Koyama et al. 2008) and theoretical (e.g. Cuadra et al. 2008) evidence suggests activity timescales of several hundred years. An even more extreme example is provided by recent observations of “Hanny’s Voorwerp” which suggest flux variations of the galactic nucleus of IC 2497 by factors  >100 on timescales  <70 000 years (e.g., Schawinski et al. 2010). Therefore one should be aware that AGN monitoring studies, like our work, might observe some of their targets at special times even when covering timelines of decades.

• 2.

  The flux density distributions of five out of six sources show clear signatures of bi- or even multimodality. This indicates that the emission processes cannot be treated as uniform stochastic processes. Instead, we have to assume a segregation of “quiescent” and (maybe multiple) “flare” source states that correspond to different physical regimes. Each state would be characterized by an individual (probably red-noise type) flux distribution.

• 3.

  Our sources show mostly steep (α ≈ 0.5−1) spectral indices that indicate outflow (halos, jets) dominated synchrotron emission. The spectral indices show substantial variability that is not correlated with the source fluxes. From the combined information we conclude that optical depth variations are the dominant mechanism for the spectral variability.

• 4.

  All sources show power spectra that are globally consistent with red noise (i.e. β > 0). We find that the intrinsic (i.e. measurement noise corrected) power-spectral indices are flatter than 1 for five of our sources, placing their emission between the cases of flicker noise (β = 1) and white noise (β = 0). The only exception is 3C 111 whose power spectrum is consistent with being located between flicker noise and random walk noise (β = 2). For three sources we find that the low-frequency ends of their power spectra appear to be upscaled in spectral power by factors  ≈2−3 with respect to the overall powerlaws, presenting another indication towards a segregation between “quiescent” and “flaring” source states.

• 5.

  For four of our targets – 0923+392, 3C 111, 3C 273, and 3C 345 – we find that the 1.3-mm and 3-mm lightcurves are simultaneous within 3σ confidence levels. This is in agreement with shock-models like the “generalized” model by Valtaoja et al. (1992); given the measurement accuracies of tens of days however these constraints are fairly loose. For 3C 454.3 (in its “flare state(s)” after 2005.2) and 3C 84 we find spectral time delays of  ≈55 days and  ≈300 days, respectively, with the 1.3-mm emission leading. The behaviour of 3C 454.3 can be traced back to variations of the optical depth with time and frequency, in agreement with previous work. In case of 3C 84, we suspect a combination of optical depth and emission region geometry effects as cause for the time lag.

• 6.

  We do not see a correlation between the observed source properties – variability, spectral index, power spectrum – with physical source parameters – redshift, black hole mass, kilo-parsec morphology – as given in Table 1.

The discovery of distinct source emission states challenges the common assumption that AGN emission can be described by uniform stochastic processes. A promising road for further studies is indicated by the cases of microquasars, Sagittarius A*, and IC 2497 where similar segregations of emission states have been seen. Long-term, multi-wavelength monitoring of AGN appears to be the most important tool towards a deeper understanding of this phenomenon.

Acknowledgments

We are grateful to the entire IRAM and PdBI staff whose hard work over many years made this study possible. We would like to thank Katie Dodds-Eden (MPE Garching) for pointing out the power of flux distribution analyses. Our work has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California

Institute of Technology, under contract with the National Aeronautics and Space Administration. We have also made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2009). Last but not least, we are grateful to an anonymous referee whose comments helped to improve the quality of this paper.

References

All Tables

All Figures

Fig. 1 The long-term lightcurves of our six target quasars spanning timelines from 1996 to 2010. Please note the different flux axis scales. For each source we give the fluxes vs. time separatly for the 1.3-mm (here labeled 1-mm for brevity), 2-mm, and 3-mm bands. For each source (except 3C 111) we approximate the most prominent, continuously sampled variations as exponential rises or decays according to Sν ∝ exp(± t/tx) (orange lines). Here t is the time, tx the characteristic rise/decay timescale of the curve labeled x   (x = a,b,c,...). The sign of t is  + t for exponential rises and −t for decays. All parameters are given in the corresponding diagrams.

In the text

Fig. 2 Distributions of the 3-mm (left column) and 1.3-mm band (right column) fluxes of 0923+392, 3C 111, and 3C 273. Please note the different axis scales and bin sizes. In order to minimize the impact of irregular sampling, we binned all lightcurves (compare Fig. 1) in time using N bins with sizes Δt. The parameters N, Δt, and the sizes of the flux bins ΔSν are given in the respective diagrams. Errorbars indicate binomial errors. Continuous grey lines indicate approximations to the data composed of individual Gaussians (dashed grey lines). In three cases (0923+392 at 1.3 mm, 3C 111 at 1.3 mm and 3 mm) we also indicate the best-fitting log-normal distributions (dotted black curves).

In the text

	Fig. 3 Same as Fig. 2 for 3C 345, 3C 454.3, and 3C 84.
In the text

Fig. 4 Evolution of the spectral indices from 1996 to 2010. Please note the different axis scales. The spectral index α is defined via Sν ∝ ν−α. Errorbars along the time axes denote the bin sizes, error bars along the α axes denote the statistical 1σ errors. Horizontal dashed lines (where visible) mark the α = 0 lines. For each source we chose the bins such that each bin covers a phase of similar flux levels (quiescent phases, flares).

In the text

Fig. 5 Evolution of the spectral index α as function of the 3-mm flux. The data points are connected by dotted lines according to their order in time. Labels 1−3, ... mark the first, second, third, ... measurement in time. The time bins used here are those defined for Fig. 4. Error bars along the flux axis denote the flux variability (standard deviation), error bars along the α axis denote statistical 1σ errors.

In the text

Fig. 6 Scargle periodograms (spectral power Af vs. frequency f) of the 3-mm (left column) and 1.3-mm band (right column) lightcurves of 0923+392, 3C 111, and 3C 273. Please note the different axis scales. The black curves denote the measured power spectra, the continuous grey lines marks the best-fitting powerlaws. The spectral index β of the power spectra is defined via Af ∝ f−β. Dashed grey lines (where present) correspond to the best-fiting powerlaws (continuous grey lines) upscaled by factors s as given in the corresponding plots.

In the text

	Fig. 7 Same as Fig. 6 for 3C 345, 3C 454.3, and 3C 84.
In the text

Fig. 8 Discrete correlation functions (DCF) as functions of time lags τ between the 1.3-mm and 3-mm lightcurves (see Fig. 1) of each AGN. The black graphs indicate the measured DCF, bin sizes Δτ are indicated in the corresponding diagrams. The maxima of the curves correspond to the most probable time lag between the spectral bands. Grey shaded areas mark the DCF ranges  [max(DCF),max(DCF) − 3σ]  with σ being the statistical error of max(DCF). Thus all values located within the shaded area are consistent with the maximum of the DCF curve within 3σ confidence. A positive (negative) time lag means that the 3-mm lightcurve preceeds (lags behind) the 1.3-mm lightcurve. Vertical dashed lines mark the null positions of the time lag axes, corresponding to the absence of time lags.

In the text

Fig. 9 Like Fig. 8, for 3C 454.3 before 2005.2 (left panel) and after 2005.2 (right panel). During the “quiescent” phase before 2005.2, the peak of the DCF is in good agreement with a time delay of zero. We note the low absolute value of the correlation, DCF(0) ≈ 0.6. For the “flaring” phase after 2005.2, we find a significant time delay in the range (as defined by the 3σ confidence interval) of  ~ −15...−80 days, peaking at τ ≈ −55 days. The negative delay indicates that the 1.3-mm lightcurve is ahead in time of the 3-mm lightcurve.

In the text

Fig. 10 A visualization of the time delay between the 1.3-mm and 3-mm lightcurves of 3C 84. Left panel: the 1.3-mm and 3-mm lightcurves. Each lightcurve is normalized to zero mean and unity standard deviation, in analogy to the calculation of the unbinned discrete correlations according to Eq. (7). Here Sν indicates the flux density,  ⟨ Sν ⟩  is the average flux density, and σ denotes the standard deviation of each lightcurve. Right panel: same as left panel, but with the 1.3-mm lightcurve artificially retarded by 300 days. Comparison of the two panel shows improved agreement of the two lightcurves when taking into account the time delay τ ≈ −300 days found from the discrete correlation function; see also Fig. 7, bottom right panel.

In the text

	Fig. 11 Spectral index α vs. time based on 1.3-mm and 3-mm flux data only, with the 1.3-mm lightcurves artificially retarded. Left panel: 3C 454.3, using a spectral time delay τ = −55 days. Right panel: 3C 84, using a spectral delay τ = −300 days. The time delays are given by the correlation analysis illustrated in Fig. 7. This figure should be compared to Fig. 4.
In the text