Issue 
A&A
Volume 533, September 2011



Article Number  A97  
Number of page(s)  18  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201015558  
Published online  07 September 2011 
The longterm millimeter activity of active galactic nuclei^{⋆}
^{1}
Seoul National UniversityDepartment of Physics and Astronomy,
599 Gwanakro, Gwanakgu,
Seoul
151742,
South Korea
email: trippe@astro.snu.ac.kr
^{2}
Institut de Radioastronomie Millimétrique (IRAM),
300 rue de la Piscine, 38406 Saint Martin d’Hères,
France
^{3}
LERMA, Observatoire de Paris, UMR 8112 du CNRS,
75014
Paris,
France
^{4}
Observatoire de Paris, LESIA, 92195
Meudon,
France
^{5}
Istituto di Radioastronomia – INAF, via Gobetti 101,
Bologna,
Italy
^{6}
ESO, Karl Schwarzschild Str. 2, 85748
Garching bei München,
Germany
^{7}
INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi
5, 50125,
Firenze,
Italy
Received: 10 August 2010
Accepted: 22 July 2011
We analyze the longterm evolution of the fluxes of six active galactic nuclei (AGN) – 0923+392, 3C 111, 3C 273, 3C 345, 3C 454.3, and 3C 84 – in the frequency range 80−267 GHz using archival calibration data of the IRAM Plateau de Bure Interferometer. Our dataset spans a long timeline of ≈14 years with 974−3027 flux measurements per source. We find strong (factors ≈2−8) flux variability on timescales of years for all sources. The flux density distributions of five out of six sources show clear signatures of bi or even multimodality. Our sources show mostly steep (α ≈ 0.5−1), variable spectral indices that indicate outflow dominated emission; the variability is most probably due to optical depth variations. The power spectra globally correspond to rednoise spectra with five sources being located between the cases of white and flicker noise and one source (3C 111) being closer to the case of random walk noise. For three sources the lowfrequency ends of their power spectra appear to be upscaled in spectral power by factors ≈2−3 with respect to the overall powerlaws. In two sources, 3C 454.3 and 3C 84, the 1.3mm emission preceeds the 3mm emission by ≈55 and ≈300 days, respectively, probably due to (combinations of) optical depth and emission region geometry effects. We conclude that the source emission cannot be described by uniform stochastic emission processes; instead, a distinction of “quiescent” and (maybe multiple) “flare” states of the source emission appears to be necessary.
Key words: galaxies: active / quasars: general / radiation mechanisms: nonthermal
© ESO, 2011
1. Introduction
Active galactic nuclei (AGN) have been studied extensively in the wavelength range from cmradio 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).
Physical properties and observation journal for our six target AGN.
The viewingangle 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 quasiperiodic oscillations, from tens of minutes (using Xray observations; e.g. Benlloch et al. 2001, and using mmradio observations; e.g. Schödel et al. 2007) to tens of years (based on cm and mmradio monitoring data; Hovatta et al. 2007; 2008). Other works focused on the power spectral statistics of the emission which was found to globally follow rednoise 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 largeamplitude activity on timescales ≳6 years. Therefore, analysing the longterm 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 restframe 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 dualband, singlepolarization receivers. One frequency in the 1.3mm band and one frequency in the 3mm band could be selected for simultaneous observations; the 2mm band was not reachable. At the beginning of 2007, new singleband, dualpolarization 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 2mm band became available for observations.
Flux monitoring data for AGN are available from October 1996 on as a byproduct of regular observatory operations. The PdBI uses the most luminous (S_{ν} ≳ 3 Jy in the 3mm 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 3mm band) to ≈37 Jy/K (for the 1.3mm 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 wellknown 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 2mm and 3mm bands and ≈16% for the 1.3mm 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 signaltonoise 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 longterm lightcurves. In order to remove this effect, we rejected flux measurements with fluxes located below the “divides” that separate the longterm 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 longterm 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.3mm (here labeled 1mm for brevity), 2mm, and 3mm 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/t_{x}) (orange lines). Here t is the time, t_{x} 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^{1} and the MPE Garching data processing software DPUSER^{2}.
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 2mm 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 kiloparsec 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.3mm, 2mm, and 3mm 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.3mm and 3mm bands. The 2mm 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 3mm (left column) and 1.3mm 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 bestfitting lognormal distributions (dotted black curves). 
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_{0} being the center position, a_{0} being the amplitude at S_{0}, and σ^{2} 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 lognormal distributions (4)where ln denotes the natural logarithm and σ and μ are the scale parameter and the location parameter, respectively. As lognormal distributions are characteristic for multiplicative stochastic processes^{3}, they may be expected in AGN flux distributions (compare, e.g., DoddsEden 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 mmspectroscopic 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 nonsimultaneous observations. We therefore inserted all measurements within selected time windows into fluxfrequency diagrams and derived timeaveraged spectral indices α via means of leastsquares 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 lowflux regime between two outbursts.
In order to test possible correlations between spectral index and source flux, we present α as function of the 3mm 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 3mm 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 σ^{2} 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 zerofrequency power. In order to exclude systematic errors due to spectral effects, we performed the periodogram analysis separately for the 1.3mm and 3mm bands of each source (for the 2mm 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 rednoise spectra. In order to quantify this behaviour, we derived powerlaw indices by fitting straight lines to the data in logspace. The spectral index β of the power spectra is defined via A_{f} ∝ f^{−β}.
In several cases we find that the lowfrequency ends of the power spectra follow the average powerlaw slopes but are located systematically above the bestfitting 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 A_{f} vs. frequency f) of the 3mm (left column) and 1.3mm 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 bestfitting powerlaws. The spectral index β of the power spectra is defined via A_{f} ∝ f^{−β}. Dashed grey lines (where present) correspond to the bestfiting powerlaws (continuous grey lines) upscaled by factors s as given in the corresponding plots. 
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.3mm 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 3mm lightcurve preceeds (follows) the 1.3mm lightcurve. In each diagram, we indicate the ranges (as greyshaded 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.3mm and 3mm 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.3mm and 3mm 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 3mm lightcurve preceeds (lags behind) the 1.3mm 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 3mm (1.3mm) 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 presentday 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 3mm flux distribution shows a double peak profile that can be described by two Gaussian components separated by ≈1.6 Jy. The 1.3mm flux distribution shows an asymmetry indicative of a second component but no significant secondary peak. We therefore attempted to describe the histogram by a lognormal distribution according to Eq. (4). However, we did not find a satisfactory solution, with the best fitting curve showing a χ^{2}/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 3mm 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.3mm and 3mm lightcurves are simultaneous within the 3σ confidence interval. The power spectra of the 3mm and 1.3mm flux measurements can be well described as simple rednoise 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 lowactivity 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 3mm (1.3mm) 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) lognormal distributions. In both cases we find bestfit solutions with χ^{2}/d.o.f. ≈ 1.6. Although this is a relatively high value, the lognormal distribution appears to be an acceptable description of the data. Within errors, the 1.3mm and 3mm 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 3mm (1.3mm) flux data globally follows a rednoise powerlaw spectrum with β ≈ 1.0 (β ≈ 0.6). The lowfrequency end (f ≈ 0.002−0.01 days^{1}) of the 3mm power spectrum is located systematically above the bestfitting overall powerlaw; this region can be described as an upscaled (by a factor ≈2) version of the global rednoise spectrum.
4.3. 3C 273
At the beginning of the monitoring, the 3mm (1.3mm) 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, doublepeaked 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.3mm and 3mm 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 3mm and 1.3mm lightcurves follow rednoise powerlaws with identical (within errors) indices β ≈ 0.45. For frequencies f ≲ 0.04 days^{1} the spectral power is located systematically above the bestfitting 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 3mm (1.3mm) 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 mid1998. 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 3mm (1.3mm) flux distributions follow doublepeaked 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.3mm and 3mm 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.3mm lightcurve is ahead in time of the 3mm lightcurve. 
Fig. 10 A visualization of the time delay between the 1.3mm and 3mm lightcurves of 3C 84. Left panel: the 1.3mm and 3mm 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.3mm 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.3mm and 3mm flux data only, with the 1.3mm 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 mmwavelengths is consistent with the classification of 3C 345 as FSRQ: the classification as FSRQ is based on cmobservations, whereas at mmwavelengths 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 rednoise spectra, with powerlaw indices β ≈ 0.5 for the 3mm and β ≈ 0.4 for the 1.3mm 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 3mm (1.3mm) 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 3mm (1.3mm) 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.3mm emission preceeds the 3mm 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.3mm emission preceeds the 3mm 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.3mm and 3mm flux density data only. We compensated the spectral time delay by artificially retarding the 1.3mm 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.3mm 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 3mm (1.3mm) 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 doublepeak 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 rednoise spectra, with powerlaw indices β ≈ 0.7 for the 3mm and β ≈ 0.4 for the 1.3mm lightcurve, respectively.
The discrete correlation function computed from the 1.3mm and 3mm 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.3mm emission preceeds the 3mm 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.3mm and 3mm 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.3mm and 3mm 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.3mm 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 ≈14year 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 lognormal distribution, i.e. an asymmetric distribution with an extended wing towards high values (see Uttley et al. 2005, for a discussion of AGN Xray lightcurves).
The flux distributions of 3C 111 are (at least marginally) consistent with lognormal 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 mmemission 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 kiloparsec 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 sourceintegrated 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 powerspectral 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 Xray (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 2mm and 3mm fluxes and ≈16% for the 1.3mm flux data (see Sect. 2). The presence of white measurement noise leads to the observation of powerspectral 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 rednoise 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.3mm 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.
True spectral indices β_{true} vs. observed indices β_{obs}.
From our test we can draw various conclusions:

The presence of Gaussian measurement noise flattens thepowerspectral indices toβ_{obs} < 1.5 even if β_{true} ≈ 2.

As β_{obs} ≈ 0.4−0.7 for all sources except 3C 111, the intrinsic powerspectral indices should be flatter than 1 in general.

Increasing the amplitude of measurement noise from 0.1 S_{ν} (like for our 3mm flux data) to 0.16 S_{ν} (like for our 1.3mm flux data) leads to a further flattening of the powerspectral index by about 0.1. This agrees with our observations: we find a systematic trend towards flatter (by ≈0.1) indices for the 1.3mm data compared to the 3mm 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^{4}.
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.0050.01 days^{1}. In combination with the information from the flux distributions (Sect. 5.1; Figs. 2, 3) this strongly suggests that “quiescent” highfrequency and “flaring” lowfrequency 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 rednoise distributions. However, a physical dichotomy of quiescent vs. flaring source states has been suggested for the nearinfrared emission observed from the Galactic nuclear black hole, Sagittarius A* (e.g. Gillessen et al. 2006; Trippe et al. 2007; DoddsEden 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 wellknown 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 Xray spectral properties (e.g., Tananbaum et al. 1972; Mirabel & Rodriguez 1998; 1999; Takeuchi & Mineshige 1998). A possible explanation for this transition (KleinWolt et al. 2002) is the difference between the ejection of fairly localized knots from the black hole vicinity (leading to Xray and radio flares) and the ejection of quasicontinuous, extended jets that are initially optically thick (leading to “plateaus” in the Xray and radio lightcurves). For microquasars, 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 blackhole mass scaling relation – in order to be consistent with our observations.
Notably, the “flare” and “highfrequency” 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.3mm and 3mm 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 cmradio range. Hovatta et al. (2008) conclude that those limits on time delays, both in the cm and mmregimes, can be understood in the frame of the semiempiric “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.3mm emission preceeds the 3mm 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 highfrequency 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 shockinjet 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^{5}. 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. Highresolution radiointerferometric 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 shockinjet 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 lightyear.
6. Conclusions
We have analyzed the longterm 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, nonperiodic 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 rednoise 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) powerspectral 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 lowfrequency 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.3mm and 3mm lightcurves are simultaneous within 3σ confidence levels. This is in agreement with shockmodels 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.3mm 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, kiloparsec 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. Longterm, multiwavelength monitoring of AGN appears to be the most important tool towards a deeper understanding of this phenomenon.
Developed and maintained by the IRAM GILDAS team, see http://www.iram.fr/IRAMFR/GILDAS/
Developed and maintained by Thomas Ott (MPE Garching), see http://www.mpe.mpg.de/~ott/dpuser/index.html
This is a simple consequence of the central limit theorem which states that sums of random numbers are normally distributed. As the logarithm of a product of random numbers equals the sum of the logarithms of the individual random numbers, it follows (under certain weak conditions) that the logarithms of products of random numbers are normally distributed (e.g., Montroll & Shlesinger 1982, and references therein).
Although there seems to be some confusion in the literature; e.g., Do et al. (2009) conclude that the infrared emission from Sagittarius A* follows a red noise power spectrum with β ≈ 2.5, without providing a physical interpretation of this unexpected result.
This observation changes only marginally when artificially removing the spectral time delay; see Fig. 11. The most notable difference between the two calculations is the time when inverted indices are observed.
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 DoddsEden (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
 Abramowicz, M. A., Bao, G., Lanza, A., & Zhang, X.H. 1991, A&A, 245, 454 [NASA ADS] [Google Scholar]
 Asada, K., Kameno, S., Shen, Z.Q., et al. 2006, PASJ, 58, 261 [Google Scholar]
 Benlloch, S., Wilms, J., Edelson, R., Yaqoob, T., & Staubert, R. 2001, ApJ, 562, L121 [NASA ADS] [CrossRef] [Google Scholar]
 Cuadra, J., Nayakshin, S., & Martins, F. 2008, MNRAS, 383, 458 [NASA ADS] [Google Scholar]
 Do, T., Ghez, A. M., Morris, M. R., et al. 2009, ApJ, 691, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 DoddsEden, K., Gillessen, S., Fritz, T. K., et al. 2011, ApJ, 728, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrarese, L., & Ford, H. 2005, SSR, 116, 523 [Google Scholar]
 Fromm, C. M., Perucho, M., Ros, E., et al. 2011, A&A, 531, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillessen, S., Eisenhauer, F., Quataert, E., et al. 2006, ApJ, 640, L163 [NASA ADS] [CrossRef] [Google Scholar]
 Ginzburg, V. L., & Syrovatskii, S. I. 1965, ARA&A, 3, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1 [Google Scholar]
 Hovatta, T., Tornikoski, M., Lainela, M., et al. 2007, A&A, 469, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hovatta, T., Nieppola, E., Tornikoski, M., et al. 2008, A&A, 485, 51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hufnagel, B. R., & Bregman, J. N. 1992, ApJ, 386, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Jorstad, S. G., Marscher, A. P., Larionov, V. M., et al. 2010, ApJ, 715, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, S. 2000, PASJ, 53, 1 [NASA ADS] [Google Scholar]
 Kaulakys, B., & Alaburda, M. 2009, JSMTE, 2, 51 [NASA ADS] [Google Scholar]
 Kelly, B. C., Sobolewska, M., & Siemiginowska, A. 2011, ApJ, 730, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Kembhavi, A. K., & Narlikar, J. V. 1999, Quasars and Active Galactic Nuclei (Cambridge University Press), ISBN 0521479894 [Google Scholar]
 KleinWolt, M., Fender, R. P., Pooley, G. G., et al. 2002, MNRAS, 331, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K., Inui, T., Matsumoto, H., & Tsuru, T. G. 2008, PASJ, 60, S201 [NASA ADS] [CrossRef] [Google Scholar]
 Krolik, J. H. 1999, Active Galactic Nuclei (Princeton University Press), ISBN 0691011524 [Google Scholar]
 Lawrence, A. 1987, PASP, 99, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Lawrence, A., & Papadakis, I. 1993, ApJ, 414, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Lister, M. L., Aller, H. D., Aller, M. F., et al. 2009, AJ, 137, 3718 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Y., Jiang, D. R., & Gu, M. F. 2006, ApJ, 637, 669 [NASA ADS] [CrossRef] [Google Scholar]
 Makse, H. A., Havlin, S., Schwartz, M., & Stanley, H. E. 1996, Phys. Rev. E, 53, 5445 [NASA ADS] [CrossRef] [Google Scholar]
 Marchesini, D., Celotti, A., & Ferrarese, L. 2004, MNRAS, 351, 733 [CrossRef] [Google Scholar]
 Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, L., Schödel, R., Eckart, A., et al. 2007, A&A, 473, 707 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meyer, L., Do, T., Ghez, A., et al. 2008, ApJ, 688, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Mirabel, I. F., & Rodriguez, L. F. 1998, Nature, 392, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Mirabel, I. F., & Rodriguez, L. F. 1999, ARA&A, 37, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Montroll, E. W., & Shlesinger, M. F. 1982, Proc. Natl. Acad. Sci. USA, 79, 3380 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, H., Suzuki, K., Asada, K., et al. 2010, PASJ, 62, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Pang, N.N., Yu, Y.K., & HalpinHealy, T. 1995, Phys. Rev. E, 52, 3224 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H. 1978, Comm. Astrophys., 7, 103 [Google Scholar]
 Reuter, H.P., Kramer, C., Sievers, A., et al. 1997, A&ASS, 122, 271 [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Schawinski, K., Evans, D. A., Virani, S., et al. 2010, ApJ, 724, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559 [NASA ADS] [CrossRef] [Google Scholar]
 Schödel, R., Krips, M., Markoff, S., Neri, R., & Eckart, A. 2007, A&A, 463, 551 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steppe, H., Salter, C. J., Chini, R., et al. 1988, A&ASS, 75, 317 [NASA ADS] [Google Scholar]
 Steppe, H., Liechti, S., Mauersberger, R., et al. 1992, A&ASS, 96, 441 [NASA ADS] [Google Scholar]
 Steppe, H., Paubert, G., Sievers, A., et al. 1993, A&ASS, 102, 611 [NASA ADS] [Google Scholar]
 Tafoya, D., Gómez, Y., & Rodríguez, L. F. 2004, ApJ, 610, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuchi, M., & Mineshige, S. 1998, ApJ, 505, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Tananbaum, H., Gursky, H., Kellogg, E., Giacconi, R., & Jones, C. 1972, ApJ, 177, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Timmer, J., & König, M. 1995, A&A, 300, 707 [NASA ADS] [Google Scholar]
 Tornikoski, M., Valtaoja, E., Terasranta, H., et al. 1994, A&A, 289, 673 [NASA ADS] [Google Scholar]
 Trippe, S., Paumard, T., Ott, T., et al. 2007, MNRAS, 375, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Urry, C. M., & Padovani, P. 1995, PASP, 99, 309 [Google Scholar]
 Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Valtaoja, E., Terasranta, H., Urpo, S., et al. 1992, A&A, 254, 71 [NASA ADS] [Google Scholar]
 van der Laan, H. 1966, Nature, 211, 1131 [NASA ADS] [CrossRef] [Google Scholar]
 Winters, J. M., & Neri, R. 2010, An Introduction to the IRAM Plateau de Bure Interferometer, public IRAM document, version 4.101 [Google Scholar]
All Tables
All Figures
Fig. 1 The longterm 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.3mm (here labeled 1mm for brevity), 2mm, and 3mm 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/t_{x}) (orange lines). Here t is the time, t_{x} 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 3mm (left column) and 1.3mm 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 bestfitting lognormal 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 3mm 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 A_{f} vs. frequency f) of the 3mm (left column) and 1.3mm 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 bestfitting powerlaws. The spectral index β of the power spectra is defined via A_{f} ∝ f^{−β}. Dashed grey lines (where present) correspond to the bestfiting 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.3mm and 3mm 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 3mm lightcurve preceeds (lags behind) the 1.3mm 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.3mm lightcurve is ahead in time of the 3mm lightcurve. 

In the text 
Fig. 10 A visualization of the time delay between the 1.3mm and 3mm lightcurves of 3C 84. Left panel: the 1.3mm and 3mm 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.3mm 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.3mm and 3mm flux data only, with the 1.3mm 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.