© The Authors 2025

1. Introduction

The synchrotron emission of active galactic nuclei (AGNs) observed in the radio regime is generated in the relativistic jets of AGNs. When the jet is pointing close to the line of sight, the AGN is called a blazar. Blazars have highly beamed emission, and they exhibit variability in frequencies ranging from radio to γ-rays (e.g. Hovatta & Lindfors 2019 and references therein).

Blazars and the causes for their variability have been studied intensely since their discovery in the 1960s. In blazars, the relativistic jet accelerates the synchrotron emitting electrons to such high frequencies that they are also visible in the optical range and X-rays (e.g. Urry & Padovani 1995).

Blazars can be divided into two categories: flat-spectrum radio quasars (FSRQs) and BL Lac objects (BLOs). The FSRQs have a spectral energy distribution (SED) where the peak is typically in the infrared regime (Hovatta & Lindfors 2019 and references therein). In BLOs, the defining features are a flat radio spectrum with a featureless optical spectrum (Stickel et al. 1991; Stocke et al. 1991). The SED peak is between the infrared and hard X-rays. A third radio-loud source type, a radio galaxy, differs from a blazar in that its relativistic jet does not point towards the line of sight.

A popular theoretical model for explaining the variability in the jet, the shock model, was proposed by Marscher & Gear (1985) and generalised to cover high-frequency radio variability in Valtaoja et al. (1992). In their model, the injected plasma flows outwards in a conical jet, and changes in the flow parameters cause the formation of shocks. The superluminal knots, observed in blazar jets and believed to be caused by these shocks (e.g. Blandford & Königl 1979), have been temporally related to flares in AGN radio light curves (e.g. Turler et al. 2000; Savolainen et al. 2002; Lindfors et al. 2006).

The emission observed from AGNs is noise produced by a stochastic process. The power spectral density (PSD) of a time series describes how its power is divided over temporal frequencies. In the simplest noise-process case, the PSD shape is a simple power law, and its slope determines what type of noise it is. Here, the most relevant noise types are white noise with slope β = 0, flicker noise with β = 1, and red noise with β = 2 (e.g. Press 1978; Uttley et al. 2002; Max-Moerbeck et al. 2014; Ramakrishnan et al. 2015). White noise is uncorrelated noise, whereas its integral, red noise, has uncorrelated increments. Flicker noise is found everywhere, from electronic components to biological signals and natural phenomena, but despite its prevalence, its generation mechanisms have not been successfully explained.

The AGN variability is often studied by analysing timescales, which can in some cases be extracted from the PSD. The term ‘timescale’ is used ambiguously in the literature and may mean any of the following things: the average time between flares, average time from the start to finish of a flare, a periodicity, or a quasiperiodicity. A simple power law may be fitted to the PSD of a source, and the best-fit slope then describes how rapid or slow the variability of the source is. Assuming that the simple power law is sufficient in describing the noise process, one can use it as a model to study potential deviations from it. In this case any deviations would reveal a timescale that is assumed to be caused by a separate process, that is, a periodic or quasiperiodic signal with some periodic properties not consistent with the noise process.

In general, though, a simple power law may not be sufficient for describing the long-term variability of AGNs. Using the method of looking for deviations from a simple power law has been challenged due to difficulties in estimating the correct underlying noise process as well as due to challenges in estimating the significance of these deviations (Vaughan et al. 2016). A deviation from the underlying noise process is not the only timescale of interest though: Especially in studies using data from the X-ray regime, a bending or broken power law has been used to obtain a specific timescale, which can be seen as a bend frequency in the PSD of some sources (e.g. Edelson & Nandra 1999; Uttley et al. 2002; McHardy et al. 2006). In fact, it is widely accepted that this type of a bend should be present in all frequency domains but at varying timescales and that this bend represents a type of characteristic timescale (Uttley et al. 2002). Variability timescales in the X-rays are typically on the orders of days, and thus if denser sampling can be achieved, relatively short monitoring campaigns are sufficient in characterising the X-ray variability. On the other hand, AGN radio variability is difficult to characterise reliably because long-term radio monitoring is required. Hovatta et al. (2007) showed that for some sources, accurate analysis of their variability in the radio regime can require over 25 years of observations. To obtain such extremely long light curves with dense-enough sampling, dedicated long monitoring programmes are required.

We generated the periodograms representing the PSDs of 123 sources from the Aalto University Metsähovi Radio Observatory (MRO) 37 GHz sample in order to probe for bend frequencies, which are physically motivated characteristic timescales. The used methods have been scrutinised thoroughly in order to make as accurate interpretations of the results as possible. An accompanying paper (Kankkunen et al. 2025, hereafter Paper II) has been written to clarify some of the common methodological pitfalls and to justify the decisions made for this analysis.

In Sect. 2, we describe the used data and the sample. In Sect. 3, we explain the methods and assumptions we used in our analysis. We dedicate Sect. 4 to some of the challenges that arise in this type of analysis, followed by the results in Sect. 5, and discussion in Sect. 6. We conclude the study in Sect. 7.

2. Data and the sample

The MRO 14 m telescope observes AGNs in both 22 GHz and 37 GHz, where the latter includes observations of up to 42 years. The entire MRO sample consists of approximately 1200 sources.

The long-term variability (1980–2005) of a large sample of MRO sources has been studied by Hovatta et al. (2007), and the five-year variability of a smaller sample in connection with gamma-ray variability has been studied by Ramakrishnan et al. (2015). The data in our study consist of observations extending until January 2023, thus bringing an additional 17 years of data to the long-term variability analysis by Hovatta et al. (2007).

We opted to only use the 37 GHz observations from the MRO sample, as they are more consistently sampled than the 22 GHz observations, and they include the longest monitoring periods. The observational and reduction methods are thoroughly described in Teräsranta et al. (1998). The detection limit of MRO is 200 mJy in ideal conditions.

The sample for this study was chosen from the 37 GHz MRO sample with the following criteria: a minimum of ten years of observations, maximum flux density of at least 1 Jy, and a minimum of N = 100 data points. With these criteria, the number of sources to be analysed was 123. The sources, their total number of data points, as well as their categories (Flat-Spectrum Radio Quasar, BL Lac Object, Radio Galaxy) are listed in Table A.1.

3. Methods

For this study, we tested different methods described thoroughly in Paper II. Below, we give the full description of the chosen analysis process with some references to Paper II, where further justification for these choices is given. The used methods are not new, but we have attempted to carefully analyse their caveats.

3.1. The power spectral density

The PSD of a signal describes its power content over the frequency domain. It can be obtained from the modulus squared of the discrete Fourier transform (DFT) defined as (Deeming 1975)

$$\begin{array}{c}\hfill |F_N{(\nu )|}^2={({\displaystyle \sum _{i=1}^N}f(t_i)cos(2\pi \nu t_i))}^2+{({\displaystyle \sum _{i=1}^N}f(t_i)sin(2\pi \nu t_i))}^2,\end{array}$$(1)

where f(t_i) is a time series, t is the time of the observation, N is the number of data points, and ν is the sampled frequency.

In the case of a simple power-law PSD, a smooth light curve corresponds to a steep power law because there is little power in short intervals (temporal distances of flares). Conversely, an erratic light curve of short temporal distances between flares likely corresponds to a flatter slope. However, if there are large-amplitude flares over long intervals, the power law steepens even if there are also variations over short intervals. The steepness of the power-law slope depends on the overall ratio in power content.

In general though, one should be careful in visually analysing light curves for their PSDs, as there are different ways of generating noise. Press (1978) shows that it is possible to generate very different-looking time series that still have the same PSDs.

3.2. The probability density function

The probability density function (PDF) gives the distribution and probability of a random variable having a certain value. In our case, the PDF describes the distribution and probability of flux densities in a light curve. It is still unclear whether blazars have one typical PDF and how it depends on the observing frequency. The PDF of any AGN cannot be purely Gaussian, as the flux density is always positive. We discuss the PDF in the context of simulations further in Paper II.

3.3. Initial data processing

We cleaned the observations by only choosing the data points with a signal-to-noise ratio (S/N) of over four, as data points with S/N ≤ 4 are considered non-detections in the MRO data. While it is uncommon, sometimes there are multiple observations of a given source during one day, often related to a high-activity state. Because we were not interested in intra-day variability in this study and the typical cadence of MRO observations would make such an analysis challenging, we only kept the first detection of the day.

In addition, we shifted each observation to the beginning of the day to allow our simulation grid to be daily, and we then binned the resulting data points to be weekly. This change in the simulation grid did not have any effect on the results, but it reduced the required computing time significantly, as we were only required to simulate one data point per day instead of hourly data points.

3.3.1. Interpolation

The periodogram introduced below in Sect. 3.4.1 expects evenly sampled data, which means that interpolation is required for our light curves. The commonly used generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009), an improved version of the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982), is a modification of the conventional periodogram and does not require evenly sampled data. However, we opted to not use it due to some of the issues underlined by VanderPlas (2018), which are also discussed in Paper II.

We used linear interpolation for our light curves, as more complex ways to interpolate did not seem beneficial. While interpolation appeared to work satisfactory in our analysis, it is not completely without drawbacks. We discuss some of the potential issues further in Paper II.

Uttley et al. (2002) observed their sources with different cadences, and they separated their PSD analysis into frequency ranges based on these observing schemes. They used the entire light curve length to determine the low-frequency portion of the PSD, the long-look portion of the light curve for the high frequencies, and the intensive monitoring segment to determine the intermediate frequencies. Because our sources have sometimes been observed with a less frequent observing cadence, specifically before the year 2000 (discussed further in Sect. 4.1), we investigated its potential effect on the obtained periodograms. In principle, our simulations should minimise the differences in the observing cadence, as we are applying the same sampling for our simulated light curves. However, some sources, notably 0716+714, have a long segment of significantly sparser sampling than the rest of the light curve. If the variability in this segment is not sufficiently well sampled, it is possible that the simulations cannot reproduce the sampling effects well enough.

We simulated data using the observing cadences of a selection of observed sources that were potentially affected and tested whether the simulations were able to replicate the behaviour. We wanted to understand whether the obtained fit for the entire light curve differed from using the shorter better-sampled segment for the higher frequencies below the bend frequency. Usually, there was no significant effect on the results. We also tested this on the observed light curves themselves to confirm they were not one of the rare extreme cases, where the variability was falsely interpreted. Due to the consistent results with both methods, we opted to use the entire light curve to determine the PSDs of each source for our analysis.

3.3.2. Direct current offset correction

The direct current (DC) offset occurs, when a time series has a non-zero mean. This zero-frequency power then corresponds to the offset from the expected zero-mean amplitude and results in a spike in the periodogram near zero frequency. This is undesirable, because it can smear the possible real features of the power spectrum. We removed the DC offset (zero-frequency power) by subtracting the mean of the light curve.

The subtraction of the mean may potentially cause some deviations in the position of the periodogram in the power axis. We discuss the implications and a possible solution in Paper II.

3.3.3. Red-noise leak and windowing

Red-noise leak is a well-known issue in PSD analysis, when the slope is expected to be steep β > 2. Red-noise leak causes the frequencies that are not the exact Fourier frequencies probed by the DFT to leak into the periodogram. Because we have a finite observing window of an on-going continuous process, we do not have the integer multiples of every possible frequency that the observed process contains. These incomplete frequencies are then spread out into our observed time frame in the frequency domain causing the periodogram to visually flatten in the higher frequencies.

This issue can be alleviated for example by using a windowing function, which reduces the side lobes of the convolution of the rectangular-shaped observing window. We chose the Hann window similarly to Max-Moerbeck et al. (2014). The potential downside is that the window function smears both ends of the light curve, and thus some low-frequency information is lost. This can sometimes be seen in the periodogram as a dip in low-frequency power, which affects the goodness of fit, albeit it should only increase the rejection confidence, not necessarily alter what the best-fit parameters are. We discuss both red-noise leak as well as methods for correcting it in Paper II, and show examples to demonstrate what causes the leak and the following flattening of the periodogram.

3.4. Monte Carlo simulations and the power spectral response method

3.4.1. The periodogram

In order to find the characteristic timescales (bend frequencies) for the sources, we analysed their PSDs with the periodogram. To obtain the periodogram, a normalisation to the modulus squared of the light curve Fourier transform is applied. One commonly used normalisation (Uttley et al. 2002) is

$$\begin{array}{c}\hfill P(\nu )=\frac{2T}{{\mu }^2{N}^2}{|F_N(\nu )|}^2,\end{array}$$(2)

where T is the length of the time series, μ is the mean of the light curve, and N is the number of data points. This normalisation is linked to the fractional root mean squared (RMS) variability (Uttley et al. 2002). The normalisation of the periodogram can be chosen quite freely when the power axis position is not of interest, which is the case in our analysis, where we are fitting only the shape of the periodogram. Vaughan et al. (2003) describe some of the commonly used periodogram normalisations in astronomy. We discuss some of the issues with normalisation in Sect. 4.

The periodogram is not a perfect descriptor of the PSD: The mean of the periodogram approaches the mean of the true power spectrum as the length of the time series increases. However, the variance of the true power spectrum is exaggerated by the periodogram and increasing the number of data points does not reduce it. Due to the consequential power fluctuations in the periodogram, binning of the periodogram frequencies is standard practice (Papadakis & Lawrence 1993).

3.4.2. The power spectral response method

The power spectral response method (PSRESP, Uttley et al. 2002) is a commonly used method for estimating the shape of the PSD of unevenly sampled AGN light curves. The sampling and other distortions of the original light curve are applied to the simulations to ensure that the simulated PSD contains all of the same non-physical effects as the PSD of the original observed light curve.

We used the Timmer & Koenig (1995) formulation for simulating the light curves. The algorithm randomises both the phase and the amplitude of the Fourier transform of the surrogate time series producing approximately Gaussian-distributed light curves, albeit the Gaussianity has recently been debated (Morris et al. 2019) and appears to be related to the slope steepness β. Emmanoulopoulos et al. (2013) formulated an improved version of this method including the matching of the PDF, which closely follows the Iterated Amplitude-Adjusted Fast Fourier Transform (IAFFT) method (Schreiber & Schmitz 2000). However, we did not use the Emmanoulopoulos et al. (2013) algorithm as the fitting of the PSD does not require a matching PDF. Instead, we exponentiated our simulated light curves to visually examine whether they appeared similar to the observed light curves when using the best-fit parameters. However, as mentioned above and discussed further in Paper II, this may not provide truly log-normal PDFs due to the possible non-Gaussianity of the original Timmer & Koenig (1995) simulated light curves. Exponentiating the simulated light curves did not affect the obtained results.

3.4.3. Models

In order to generate simulated light curves, we used both the simple power law and the bending power law. The simple power law is defined as

$$\begin{array}{c}\hfill P(\nu )=\frac{1}{{\nu }^{\beta }},\end{array}$$(3)

where β is the PSD slope.

The simple power law should provide a good fit for such sources whose characteristic timescale is long and not identifiable in the PSD. However, if the bend is already within the observing window or close to it, the simple power law may give a too flat best-fit slope. We discuss this issue further in Paper II and it is already noted in Uttley et al. (2002).

We used the following formulation for the bending power law similarly to Uttley et al. (2002):

$$\begin{array}{c}\hfill {\displaystyle P(\nu )=\frac{1}{{(1+{\left(\frac{\nu }{x_{\mathrm{b}}^{-1}}\right)}^2)}^{\frac{\beta }{2}}},}\end{array}$$(4)

where x_b is the inverse of the bend frequency (characteristic timescale). The slope in this formulation flattens to zero after the bend at low frequencies. The simple power-law and bending power-law models are shown side-by-side in Fig. 1.

Fig. 1. Comparison between the simple power law and the bending power law. The simple power law in black is a straight line of slope βspl = 2. The bending power law in orange follows the simple power law with the same slope β = 2 until the bend, after which the slope turns over to β = 0.

From previous studies, especially in the X-ray regime (e.g. Uttley et al. 2002, and references therein), we know that the light curve PSDs contain a bend. This should be the case in every frequency regime simply because we do not expect a light curve to exhibit infinite variance. Mathematically, this would mean that the behaviour is divergent, whereas convergence is required.

We chose to use the simple power-law and bending power-law models instead of more complex models, such as one with multiple bends, as we are trying to understand whether AGNs observed in the radio domain already contain a visible bend. Cyg X-1 is an X-ray binary that has been shown to contain two bends (Belloni & Hasinger 1990) in the hard state and thus is possibly a motivation for such behaviour in AGNs as well (Uttley et al. 2002). However, due to the slow variability of AGNs in the radio domain we refrain from doing such analysis now and leave it for a future study.

An important thing to note is that the bending power law is equal to the simple power law at higher frequencies. This is visible in Fig. 1, where the higher frequencies converge to the same slope for both models. Because the bend is not a sharp turnover and due to sampling-induced effects, the timescale needs to be sufficiently longer than the monitoring period length for it to not have an effect on the periodogram.

3.4.4. Parameter space

From earlier research (e.g. Ramakrishnan et al. 2015) and by visually comparing light curves between different frequency regimes, it has been established that radio variability is generally slower than in the optical, gamma-ray and X-ray frequencies.

The typical parameters for radio variability are a slope β of around 2 (e.g. Max-Moerbeck et al. 2014; Ramakrishnan et al. 2015) and variability on timescales on the order of years (Hovatta et al. 2007). Because of prior information on radio variability and considering the length of each monitoring period, the searched parameter space was chosen to be

$$\begin{array}{c}\hfill \beta =1:0.1:3.5\end{array}$$

with a bend timescale

$$\begin{array}{c}\hfill x_{\mathrm{b}}(\mathrm{days})=500:500:7000.\end{array}$$

We additionally probed for higher timescales of x_b (days) = 100, 200, 400, and 800 to cover more of the high-frequency range. The chosen timescales are arbitrary. We decided on this sparse slope and timescale grid due to the unevenly sampled data and to avoid the risk of overanalysing our results. The final timescale, 7000 days, is approximately half of 40 years. Of course, not all sources have 40 years of observations, but this is not an issue as there are no restrictions for fitting a bending power law with the inverse timescale (bend frequency) longer than the monitoring period. This can be understood from Fig. 1, as the bending power law will simply resemble the simple power law, when moving to longer timescales.

3.4.5. Mitigating for aliasing and red-noise leak

We simulated the light curves with a daily resolution and then binned the data points weekly exactly in the same way as the observed light curves to mimic sampling-induced deviations from the true PSD. We ensured with additional simulations that more dense sampling was not required.

Aliasing is caused by discrete sampling of a continuous process below the Nyquist frequency. Due to this insufficient sampling, power from higher frequencies leaks into the periodogram causing flattening of the high-frequency portion. Simulations with a daily resolution as well as identical binning and sampling with the observed data sufficiently mimicked this effect.

Red-noise leak causes a similar flattening effect as aliasing, but it only seriously affects larger slope values (β > 2). We corrected for red-noise leak by simulating both 10 and 10 000 times longer light curves compared to the original observed light curves. We confirmed by tests on simulated data that 10 times longer light curves provided similar red-noise leak effects as the 10 000 times longer light curves, and we opted for the saving of computational time.

3.4.6. Observational errors

We drew a random Gaussian number with μ = 0 and σ = 1 for each observation. We then multiplied the numbers with the standard error associated with each observation to scale the Gaussian random numbers to the correct standard deviations. These errors were then added to the simulated data to mimic observational errors.

Observational errors can be seen as a flattening of the PSD slope β in the highest frequencies. This limited our ability to constrain the steepest PSD slopes and it is described in more detail in Sect. 5.

3.5. The goodness of fit

The goodness of fit is defined as (Uttley et al. 2002)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{dist}}^2={\displaystyle \sum _{\nu ={\nu }_{\min }}^{{\nu }_{\max }}}\frac{{(\overline{P_{\mathrm{sim}}}(\nu )-P_{\mathrm{obs}}(\nu ))}^2}{\mathrm{\Delta }\overline{P_{\mathrm{sim}}}{(\nu )}^2},\end{array}$$(5)

where P_obs is the periodogram of the observed data, and $\overline{P_{\mathrm{sim}}}$ the mean of the simulated periodograms. The χ_dist² distribution is otherwise similar to the chi-square statistic, but it does not assume Gaussianity.

We constructed the periodograms for both the real observations and simulations in the exact same way. We chose to use logarithmic binning with a factor of 1.3 similarly to Uttley et al. (2002). We created $\overline{P_{\mathrm{sim}}}$ by taking the mean of all simulated light curves.

In the PSRESP formulation, both the periodogram of the observed data as well as the periodograms of each simulated light curve are compared to the mean periodogram of the simulations. The best-fit parameters for the periodogram are found by minimising χ_dist² and calculating the percentile of the simulated χ_dist² distribution that is below the χ_dist² calculated for the observed data. This is the p-value for each parameter combination. The errors of the observed power spectrum can be estimated from the spread of the simulated periodograms around the constructed mean. The full procedure is explained in Uttley et al. (2002).

4. Challenges

4.1. Flare bias and uneven sampling

The MRO sample consists of unevenly sampled data. Observations are made based on source priority determined by different factors such as prior variability in the radio regime, and multi-frequency campaigns. Weather affects the observations introducing gaps and altering the telescope detection limit. Interest in high-activity sources causes the observations to be biased towards flaring events. This is challenging to apply in the Monte Carlo simulations as the PSRESP method requires replicating the sampling in order to reproduce any sampling artefacts. Flare bias may occur randomly even in the simulations, but depending on the variability, the biased sampling will also randomly occur during low-activity states. An additional issue is the sparseness of the observations before the year 2000, when a larger fraction of telescope time was used in the 22 GHz band.

Bias towards flaring causes changes in the variance and mean, which may not be similar to the variance or mean in an evenly sampled signal. This poses some problems when simulating light curves and normalising them to match the observed light curve mean and variance. Interpolation alleviates the issue but does not completely remove it.

A potential issue is also the artificial steepening of the periodogram slope due to interpolation. Interpolation of the light curve was shown to not substantially affect the periodogram when the light curve variability was expected to be slow and the bend timescale at an intermediate or low frequency. If the bend timescale was at a very high frequency or the variability was very fast, then the steepening of the spectrum due to the effect of interpolation became more prominent. This steepening is caused by the smoothness that interpolation causes, magnified by long gaps between observations. Flare bias complicates this issue as the effects of interpolation cannot be completely reproduced with simulations. We discuss this issue with interpolation further in Paper II.

4.2. Periodogram normalisation

Normalisation is also complicated by flare bias, as the variance and mean tend to be overestimated. It is also an issue due to the nature of red noise, where the parameters can vary due to statistical fluctuations (e.g. Vaughan et al. 2003).

We attempted to normalise the periodograms to have zero mean in logarithmic space. This way, any inconsistencies in the estimate of mean or variance would not alter the position of the PSD on the power axis. However, we concluded that the steepness was often overestimated with this shifting of the periodogram and benefits over the usual normalisation were unclear. We discuss the normalisation and our results regarding it further in Paper II.

4.3. P-value and constraining the characteristic timescale

We used a significance level of 0.1 for the bending power law and the simple power law, and we calculated the p-values using the PSRESP formulation. We report a ‘rejection confidence level’, that is 1 − p, using the p-value of the best-fit parameter combination for sources for which we were able to constrain a characteristic timescale. Thus, if we have a p-value = 0.1, then the rejection confidence is 90%. We obtained the 90% confidence regions similarly to Uttley et al. (2002) by including all of the parameter combinations with a p-value p ≥ 0.1. While this method of obtaining the confidence region is not incorrect, it does not consider the uncertainty in both the used model and the PSD shape as the periodogram is prone to be erratic. Mueller & Madejski (2009) discuss the issues in generating confidence levels, and they suggest a different method, which is not tied to the obtained p-value. In our analysis, if the best-fit p-value is close to the chosen significance level, the confidence region may not exist, and thus it may appear as if the result is more constrained than it is.

The obtained p-value was tied to for example the factor chosen for the logarithmic binning of the periodogram, that is it varied depending on different decisions made in the data processing steps. Thus, the p-values should not be considered to be too exact. In some cases, especially with a low number of observations, many of the fits were exceptionally good. This was likely caused by interpolation as well as errors, where a low number of data points caused the light curves to be similarly smooth and errors limited the steepness of the slope. The magnitude of this effect was also dependent on the choice of both t_bin of the light curve and the periodogram bin size. We discuss the effect of sampling on our ability to differentiate the different PSDs further in Paper II.

Due to Fourier frequencies, the lowest frequency bins have the lowest number of data points. This in part affects the smoothness of the periodogram around the expected characteristic timescale causing the fits to be susceptible to lower significance levels (Uttley et al. 2002).

We constructed the PSDs of all sources using the periodogram, and we conducted Monte Carlo simulations to find the best-fit parameters for the PSDs of each source. We ran the simulations with 1000 iterations per parameter combination. We took a smaller sub sample and tested whether running the simulations with 10 000 iterations would alter the results, but they mostly just refined the fits. Thus, we decided to keep the n = 1000 simulations for most sources and only used the n = 10 000 iterations in the cases where a timescale appeared to be constrainable.

We can directly compare the obtained best-fit p-values from the simple power law and the bending power law. As discussed in Sect. 3.4.3, when we move to longer bend timescales, we see that the bending power law starts to resemble the simple power law as expected. Although we fit the simple power law to the PSDs, it is equivalent to fitting a bending power law with a much longer timescale compared to the observing window. Thus, the p-values are directly comparable without using methods such as the Akaike Information Criterion to compare the models (due to a different number of parameters). If we get a reasonable fit for the simple power law, it in principle means that we cannot set an upper limit for our timescale.

5. Results

Our results suggest that in the radio regime, many of the bend frequencies are either at such low frequencies that they are not well identifiable, or that the number of data points is not enough for an accurate PSD estimate. This does not imply that the characteristic timescale is necessarily longer than the probed length of the observing period but rather that observations must be continued for longer for the bend frequency to start dominating in the periodogram, or that more data points are required to overcome issues caused by sparse sampling. This is also apparent upon visual examination of such periodograms that are mostly consistent with a simple power law.

Observational errors also limited our ability to constrain the slope β. This occurs because steeper slopes correspond to smoother light curves: The white noise of the observational errors start to dominate the light curve more and more changing the shape of the periodogram back into a flatter apparent high-frequency slope. We visually examined the fits and if we saw a parameter space where the fit first worsened when going to steeper slopes and then started to improve again, we discarded those steepest β values.

All of the obtained fits are listed in Table A.1. The confidence limits for the bending power-law fits are derived directly from the p-values: All fits with a p-value ≥0.1 are considered acceptable. However, as we present them in a tabular form, we are unable to convey the fact that not all parameter combinations are possible. This can be seen below in Sect. 5.2 and Fig. 2, where the heat map shows how the parameters are dependent on each other. In addition, the heat maps do not convey the uncertainties in the periodogram fits themselves as discussed in Sect. 4.3. The simple power-law confidence limits are similarly given based on the fits, where p ≥ 0.1.

From Table A.1, we can see that for all sources, except for 1036+054, we were unable to find an upper limit for the 90% confidence region. This is a direct consequence of the simple power law being a good fit for most sources where the bending power law with the bend approaching infinity equals the simple power law as discussed in Sect. 3.4.3 and demonstrated below. The upper limit of 1036+054 is the same as its best-fit timescale because the best-fit p-value is p = 0.1. This is caused by the confidence regions being determined from the p-values. We were able to constrain a lower limit for the bend timescale for most sources.

5.1. Interpretation of results

The characteristic timescale we are interested in is the bend frequency in the PSD. This timescale should be intrinsic to the source noise process, that is, it is not a separate periodic/quasiperiodic component, and the timescale should be related to the duration of the flares in the time series. The assumption of a relationship between flare duration and the timescale follows from our understanding of the noise process: The PSD bend frequency shows a bend after which a flare should not affect any future flares. Thus, as long as a flare is ongoing, it will affect all future flares increasing the overall variance. This matter is discussed further in Paper II and the reader is referred to for example Press (1978) and Halford (1968) for more accurate descriptions of power-law noise and its generation.

The slope of a PSD shows the distribution of power over temporal distances in the time domain, that is, in principle, the distance between flares and what their amplitude relations are. If we have a large amount of high-amplitude short-timescale variability, the PSD slope should be flatter. If the variability appears smooth, that is the variations between long temporal distances are clearly dominant, then the PSD slope should be steeper. However, the characteristic timescale can make visual analysis more challenging as this power-law relationship continues only until a characteristic timescale, and given long-enough monitoring periods, the low-frequency portion in the PSD will become dominant. This means that given long enough monitoring periods in comparison to the characteristic timescale, all data will resemble white noise.

For our results, we generated heat maps which show the goodness of fit of all fitted parameter combinations. The x-axis shows the slope values, and the y-axis the timescales in days. We used the heat map for our results because we found it to be the best way to convey the large number of good fits. If we only showed confidence intervals attached to singular best-fit values for the β and x_b, we would be unable to convey how the limits in both parameters are related to each other.

5.2. Example source 3C120

We chose source 3C120 as an example to demonstrate some of the concepts described earlier regarding the behaviours of the two models. This was the only source for which we analysed a shorter segment of five years in addition to the long-term analysis. We chose the shorter segment quasi-randomly, that is the segment was randomly chosen from the better-sampled portion of the light curve.

Fig. 2 shows the results for fitting the periodogram of source 3C120 with both five years of data (2010–2015) and with using the entire light curve. The heat maps show all of the fitted parameter combinations with colour indicating the goodness of fit between 0 and 1, that is the obtained p-value. The white graph overlaid with the heat map shows the simple power-law fit with the x-axis being equivalent to the heat map x-axis. The y-axis for the simple power-law fit approximately equals the heat-map legend labels with p-value from 0 to 1, that is the higher the white graph reaches, the better the fit.

Fig. 2. Comparison between the results for source 3C120, when using a shorter segment of data and the entire monitoring period. The first row of plots uses data between 2010 and 2015, and the middle row uses the entire light curve shown in bottom row, where the red portion corresponds to the five-year segment. The upper left-side heat map shows the PSRESP fits for each bending power-law parameter combination using five years of data. The white graph overlaid with the heat map shows the simple power-law fit with the x-axis being equivalent to the heat map x-axis. The y-axis for the white graph equals the heat map legend range scaled to agree with the label positions. The right-side plot shows the source periodogram (black) and each simulated mean periodogram with cyan indicating the best fit and magenta the worst fit. The heat map and PSD graph below show the equivalent results using the entire light curve of 3C120 shown on the bottom row. The dotted vertical lines on both periodogram graphs show the best-fit bend frequency obtained by using the entire light curve, and the dashed horizontal lines show the noise floor.

We should get a good fit with the bending power law even if the bend is not yet visible in the PSD of the source. The issue is constraining the characteristic timescale: As can be seen in the upper left-side heat map in Fig. 2, for five years of data, we obtain a similar best fit as for the full 40 years of data below it. However, the parameter space of near-equal fits is much larger. The periodogram of a source with a long characteristic timescale compared to the observing window length will resemble a simple power law and naturally so will the periodograms of the simulated light curves, thus causing the observed large parameter space of good fits. The p-value is 0.75 for the simple power-law fit, that is the rejection confidence is only 25%.

The simple power-law fit for the five years of data and 40 years of data are significantly different in both their slope-values β and their best-fit p-values. With five years of data, the fit is steeper (β_spl = 1.9) and closer to the value obtained by the bending power law. The best fit of the simple power-law slope is not exactly where the best-fit slopes of the bending power law are: This may occur even if the characteristic timescale is outside the observing window due to for example sampling effects or shifting and scaling of the light curve. What this means is that the bend may need to be at a significantly lower frequency than what the inverse of the observing window length is for the bend to not have any effect on the simple power-law fit. We discuss the effects of sampling and normalisation further in Paper II. The simple power law nevertheless shows a good fit also closer to the hot spot of the heat map.

For 40 years of data, the best-fit simple-power law slope is flatter, β_spl = 1.6, consistent with our expectations of seeing a flatter fit with the periodogram, when the bend is clearly visible in the PSD. The simple power-law fit significance is also much reduced compared to the bending power-law fit: The rejection confidence for the best fit simple power-law fit β_spl = 1.6 is now 70%.

These effects are also evident in both periodograms in the right-side column. There appears to be no flattening of the slope with 5 years of data, whereas with the entire light curve of over 40 years, the slope clearly flattens to a power-law slope β = 0.

5.3. Sources with well-constrained timescales and slopes

We considered the timescale and slope constrained, if the heat map of best fits showed an even distribution of values, that is there were no multiple hot spots in the heat map suggesting uncertainty in the fit. We also required for the hot spot to be sufficiently within the parameter limits so that the hot spot did not continue beyond the longest probed timescale. We also required for the simple power-law fit to be clearly worse than the bending power-law fit (minimum of 20% difference if both fits were poor). These were conservative decisions, and six borderline cases were left out from in-depth analysis (PKS0422+0036, OJ248, 1324+224, PKS1725+044, S52007+77, and PKS2022−077). For PKS0422+0036 the flattening of the periodogram in the low frequencies appeared to be an artefact; thus we omitted it from our analysis in this paper.

We analysed the sources by dividing them according to their best-fit timescales. Due to the often large parameter spaces of good fits and wide bins, this is only used as an aid for comparisons. The chosen grid of timescales is not a very accurate one but sufficient for a broad estimate. For example, a timescale of 1200 days would likely have a best fit centred around x_b (days) = 1000. We chose these timescales for simplicity as they are enough to give a general idea of the source fits and should help in choosing candidates for future work with finer parameters grids. As discussed in Sect. 3.4.4, we also did not want to overanalyse our results because it is unclear how accurately the timescales can be determined.

A common feature of the well-constrained sources is that their baseline flux-density levels remain fairly constant. If we assume the bend frequency to represent some function of flare duration, it would be expected for the flux density levels to return to the baseline level within the process duration which we see with these sources. A very high occurrence of flares may temporarily increase the baseline, but it should not increase infinitely as that would indicate either infinitely growing flare amplitudes or infinitely increasing occurrence of flares. These aspects are discussed in more detail in Paper II.

Below we have divided the well-constrained sources by their bend timescales. For these 11 sources a well-constrained best fit for the bending power law was found. We provide some visual analysis of the light curves in relation to the obtained results, but wish to remind the reader that there are some caveats in this and no quantification of results should be attempted visually.

5.3.1. 500 days

0716+714. Examining the light curve in Fig. A.1 shows short temporal distances between flares, as well as an apparent short duration of flares consistent with a short characteristic timescale. The baseline flux density does not remain completely constant, but it increases during the onset of a period with high-amplitude flares, which could be caused by their overlap.

The periodogram visually matches one with a flat slope but the bend is less evident. The bending power-law fit with a best-fit x_b = 500 days and β = 1.8 gives a rejection confidence of 63% and the simple power-law fit with β_spl = 1.2 gives a rejection confidence of 86%.

5.3.2. 800 days

2144+092. The light curve in Fig. A.2 has a baseline flux density slightly above the MRO nominal detection limit with fast apparent variability but low overall flare amplitudes. This source has one of the largest parameter spaces of good fits for the constrained sources which can be, at least partially, explained by the fewer number of data points as discussed in Sect. 4.3. The periodogram of the source shows a visually prominent bend and flattening into white noise at low frequencies. The bending power-law fit with a best-fit x_b = 800 days and β = 2.6 gives a rejection confidence of 8 % and the simple power-law fit β_spl = 1.5 gives a rejection confidence of 66%.

In the heat map, the effect discussed in the beginning of this section is visible, that is the increase of fits at steeper slopes caused by observational errors. This occurs especially when the PSD slope is steep, as the observational errors are more prevalent in smoother light curves. Thus, we disregard the steepest values.

5.3.3. 1000 days

PG0007+106. The light curve in Fig. A.3 shows a baseline flux density close to the nominal MRO detection limit with several flares that appear to have slight overlap due to multiple peaks. The best-fit power-law slope is β = 2.4, and it may correspond to a light curve, where higher amplitude flares occur relatively rarely. The bending power-law fit with a best-fit x_b = 1000 days and β = 2.4 gives a rejection confidence of 19%. The simple power-law fit is rejected with over 90% confidence.

Also in this source, the observational errors affect the bending power-law fit in the steeper slopes.

0235+164. In the heat map of Fig. A.4, we see that the best-fit power-law slope is β = 2.4, same as the best-fit slope of the source PG0007+106. However, the light curve of this source exhibits higher amplitudes for its strongest flares and the baseline flux density is closer to 1 Jy. There is a slight dip in the low-frequency part of the periodogram, which is likely a sampling artefact or potentially caused by the Hann window. Nevertheless, visual examination of the periodogram shows an apparent bend. The bending power-law fit with a best-fit x_b = 1000 days and β = 2.4 gives a rejection confidence of 49%. The simple power-law fit is rejected with over 90% confidence.

0736+017: Looking at Fig. A.5, this source has similar slope values (β = 2.5) to the previous two sources but to the eye its light curve looks strikingly different with more erratic variations. The amplitudes of the flares are similar to PG0007+106 but its baseline flux density is higher at approximately 1 Jy. It is possible that due to the different mechanisms in which the same spectra can be generated, the flares have a smaller amplitude but a higher occurrence. Of course, it is also possible that the true timescales for the sources differ by up to 500 days or potentially more within error margins. The bending power-law fit with a best-fit x_b = 1000 days and β = 2.5 gives a rejection confidence of 1% and the simple power-law fit β_spl = 1.5 gives a rejection confidence of 75%.

Visually, we see that the periodogram is very smooth explaining the low rejection confidence of the bending power-law fit. That is, with a smooth periodogram, partially obtained by chance, the fits will be high as the periodograms on average have more variations. This does not necessarily mean that the periodogram fit of this source gives a better description of the true PSD than the fits of the other sources.

Similarly to 2144+092 and PG0007+106, the observational errors affect the fit in the steeper slope values, and we have disregarded them in our analysis.

PKS1749+096. The baseline flux density of the light curve varies around 2 and 3 Jy (Fig. A.6) and its best-fit slope is β = 2.1. The light curve exhibits high-amplitude flares with short temporal separations between them. Visually, the higher amplitude flares appear to occur within shorter intervals than in 0235+164 supporting the slightly flatter PSD slope. As can be seen in the heat map, this source has one of the best-constrained parameter spaces for the bending power law. The bending power-law fit with a best-fit x_b = 1000 days and β = 2.1 gives a rejection confidence of 10% and the simple power-law fit β_spl = 1.5 gives a rejection confidence of 81%.

2230+114. The light curve has a baseline flux density of approximately 2 Jy and its best-fit slope is β = 2.3 (Fig. A.7). We see multiple large flares in fairly large temporal distances from each other, supporting the steeper periodogram slope. The bending power-law fit with a best-fit x_b = 1000 days and β = 2.3, gives a rejection confidence of 67%. The simple power-law fit is rejected with over 90% confidence.

5.3.4. 1500 days

0415+379. The light curve of this source is strikingly smooth with rare high-amplitude flares, and this is also in agreement with the best-fit slope β = 2.8 (Fig. A.8). The baseline flux density level is approximately at 2–3 Jy but appears to be decreasing slightly after the year 2010. The heat map of this source includes multiple good fits and it is not as constrained as for some of the other sources. The periodogram shows a flattening in the low frequencies preceded by a uniquely steep part in the middle of the spectrum. The bending power-law fit with a best-fit x_b = 1500 days and β = 2.8 gives a rejection confidence of 32% and the simple power law fit β_spl = 1.9 gives a rejection confidence of 69%.

4C29.45. The baseline flux density of the light curve is approximately 1 Jy (Fig. A.9). The flares have clear features and appear smooth. The slope steepness β = 2.5 indicates slower variability, and the light curve matches this visually. The bending power-law fit with a best-fit x_b = 1500 days and β = 2.5 gives a rejection confidence of 24% and the simple power-law fit β_spl = 1.7 gives a rejection confidence of 65%.

5.3.5. 2500 days

3C120. The baseline flux density of the light curve is approximately 2 Jy, although there is a slight increase after year 2015 (Fig. 2). The already non-zero baseline of the flux density increases during the strongest flare and appears to remain approximately the same or slightly decreased for multiple years. The source has a best-fit slope β = 2.0 consistent with fairly fast variability compared to some of the other sources. The periodogram of the source shows a clear flattening in the low frequencies. The source shows almost equally good fits for timescales of 1500 days and 2000 days. The bending power-law fit with a best-fit x_b = 2500 days and β = 2.0 gives a rejection confidence of 29% and the simple power-law fit β_spl = 1.6 gives a rejection confidence of 70%.

5.3.6. 3000 days

3C454.3. The light curve is visually very striking as it has a high baseline flux density of approximately 8 Jy with moderate levels of variability, after which three prominent high-amplitude flares occur back-to-back, strongest of which exceeds 49 Jy (Fig. A.10). After a decay period, there is another flare that is likely followed by multiple flares with high temporal frequency. The cooling down of the initial flare has visually lasted for nearly 10 years. The source has several good fits even beyond the best-fit timescale of 3000 days and the periodogram of the source shows a flattening in the lower frequencies. The bending power-law fit with with a best-fit x_b = 3000 days and β = 2.3 gives a rejection confidence of 21% and the simple power-law fit β_spl = 1.8 gives a rejection confidence of 68%.

5.4. Examples of sources with unconstrained timescales

5.4.1. 3C84

A strikingly different source is 3C84, which already upon visual inspection of the light curve includes extremely long timescales (Fig. A.11). It is also a secondary flux calibrator in MRO due to its slow variability, especially prominent before the year 2012.

The best simple power-law fit gives a slope β_spl = 2.6 and the best bending power-law fit a slope β = 3.5. This might in some cases suggest the presence of a bend due to the flattening effect caused by fitting a simple power law to a bending power law. However, as can be seen in Fig. A.11, the cause of this difference is more likely that the bend is in fact at such a long timescale that a better bending power-law fit would be found if longer timescales were included.

The periodogram is also consistent with a simple power law. The flat region in the higher frequency part occurs when the slope itself is steep, that is the light curve is smooth. Then the observational errors dominate the region causing it to flatten. This effect is more prominent in the simple power law, where the smoothness can continue infinitely, whereas in the bending power law the timescale forces a discontinuity.

5.4.2. BL Lac

For BL Lac, the best-fit bending power-law parameters are β = 1.9 and x_b = 7000 days (Fig. A.12). The simple power-law fit gives a best-fit slope β_spl = 1.7, where β_spl = 1.9 is still within error margins. For this source the slope appears fairly well constrained in the heat map but more data are required to confirm it and especially to confirm the timescale.

6. Discussion

As discussed briefly in Sect. 5, it is challenging to quantify an upper limit for the characteristic timescale. A simple power law will in many cases provide an adequate fit even when the bend frequency is clearly visible in the PSD of a source. This is because the high-frequency portion of the PSD equals a simple power law. Only when the power-law portion is under-represented compared to the zero-slope low-frequency portion (after the bend) should the simple power-law fit become possibly such that it can be disregarded with the usual p-value limits. In the orange graph in Fig. 1, the bend is at a high frequency; thus a simple power law would clearly provide a poor fit to the periodogram. If the bend was moved to a lower frequency, there would be a larger portion of the periodogram consistent with a simple power law. This effect is magnified if no logarithmic binning is used.

For some of the sources with a clearly better fit with the bending power law, an upper limit for the timescale may be reasonable even if it is difficult to define using the p-values. For example, the source PKS1749+096 is well sampled, and it has a clearly constrained hot spot in the heat map (Fig. A.6). The bending-power law fit gives a rejection confidence of only 10% and the simple power-law fit gives a rejection confidence of 81%. With the maximum fitted timescale of 7000 days, the rejection confidence is 73% and it is consistently increasing starting from the hot spot. Given the explanation in Sects. 4.3 and 3.4.3, the rejection confidence should continue increasing until such a timescale where the bending power law and simple power law converge.

Even for a source with a clear timescale, the result will always be a parameter space. The relationship between the slope β and the timescale x_b causes there to be a correlation between them: Given a simulated light curve of for example β = 2, and x_b = 500 days, the shape of the bending power law will be sufficiently similar for a good fit if we move to slightly longer timescales and a flatter slope. We discuss this further in Paper II and show examples of how sampling affects our ability to differentiate between fits.

6.1. Flattening to a non-zero slope

In addition to the simple and bending power laws, Uttley et al. (2002) also used a high-frequency break model, where they replaced the direct flattening to β_low = 0 with a slope β_low = 1 using a sharp break. Markowitz (2010) additionally used a smooth turnover model, which is more similar to the (smoothly) bending power law used in this analysis. The motivation behind using a non-zero low-frequency slope comes from studies of black hole X-ray binaries in the high state, where for example Cyg X-1 has been analysed to contain a low-frequency slope β_low = 1 (e.g. Cui et al. 1997). Here we discuss the effect of non-zero low-frequency slope on our results.

As mentioned in Uttley et al. (2002), there is no physical reason for a sharp break in the PSD, and thus we decided to simulate evenly sampled light curves without observational errors using the smoothly bending power law similarly to Markowitz (2010):

$$\begin{array}{c}\hfill {\displaystyle P(\nu )=\frac{{\nu }^{-{\beta }_{\mathrm{low}}}}{{(1+\left(\frac{\nu }{x_{\mathrm{b}}^{-1}}\right))}^{{\beta }_{\mathrm{high}}-{\beta }_{\mathrm{low}}}},}\end{array}$$(6)

where β_low is the low-frequency slope and β_high the high-frequency slope corresponding to β in Eq. (4).

In our simulations, we were unable to differentiate the significance of the goodness of fit between the two models irrespective of which of the two models were used as the underlying model. We used varying parameter combinations with our lowest bend timescale being 500 days. Uttley et al. (2002) obtained similar results using the two different break models for their observed data: Both models provided a good fit with good constraints to the p-values for their sources.

Our simulations showed a clear connection between the bend frequencies: The model with a low-frequency slope β_low = 1 would always fit a higher bend frequency to the mock light-curve PSD than if the underlying model had β_low = 0. Thus, using Eq. (4) gives the most conservative (longest) bend-frequency estimate and we can be more certain that the constrained sources in fact do contain a detectable bend. If the model we use does not describe the reality well, we may leave out sources with a detectable bend in their PSDs; however, we do not obtain false positives that we consider a more serious concern. Our results also suggest that using a different low-frequency slope, β_low < 1, would affect the bend position similarly, moving it to higher frequencies.

The high-frequency slopes obtained in our simulations were often consistent between both models, only the bend frequency position varied significantly. Thus, the slope values obtained with Eq. (4) are likely close to the ones we would obtain using Eq. (6) with a low-frequency slope β_low = 1. Due to this, we used the high-frequency slope we obtained in the simulations using Eq. (4) to find the best-fit bend frequency when the low-frequency slope was set to β_low = 1. This should provide a reasonable approximation for the bend frequency. The results for the bend timescale with β_low = 1, using the best-fit high-frequency slope to find the bend frequency, are also reported in Table A.1.

Because we used a fixed high-frequency slope, it is expected that the obtained 90% confidence regions would be more constrained than when the high-frequency slope is allowed to vary freely. However, this only seems to be the case for the 11 constrained sources as well as for some of the borderline cases. This result shows how poorly we are able to constrain the timescale for most sources.

6.2. Comparison to 43 GHz very long baseline interferometry data

Light curves from one frequency band alone are limited in how much information they convey. As discussed further in Paper II, unless the process is Gaussian-Markov, the PSD does not reveal all statistical information of the random process (Press 1978). In general, light curves should be compared between different frequency bands or for instance against very long baseline interferometry (VLBI) images to be able to better analyse the underlying physical conditions.

To lay some ground work for future analysis, we compared a selection of light curves and their obtained characteristic timescales to the 43 GHz VLBI sample by Boston University (Jorstad et al. 2017; Weaver et al. 2022)

6.2.1. Large-scale emission

We tested our hypothesis of how flare superposition affects the baseline flux density levels by comparing the 37 GHz emission to the 43 GHz VLBI emission of sources 3C120 and 0716+714. AGNs emit on different scales and with a single-dish telescope all emission from parsec-scale knots and kiloparsec-scale diffuse emission is visible. With VLBI, one samples only the smaller compact structures; it is not sensitive to the diffuse emission and the flux density in VLBI images is only from the parsec-scale knots. Comparing between the two light curves should allow for identification of how much the large-scale emission increases the baseline of the MRO observations. Here we assume that the large-scale emission does not exhibit any significant variability, and thus increases in the baseline flux density compared to the combined 43 GHz core+jet emission from VLBA images would then be caused by flare superposition.

We added both the core and jet flux density into a combined flux density per observation epoch and then compared it to the MRO curve in Fig. 3. The fast variability observed in 3C120 complicates the analysis as the flux density can increase over 1 Jy in one day and observations rarely coincide perfectly but rather have a few days of difference. The feature in January 2015, where the 43 GHz combined flux density exceeds 37 GHz, is likely a manifestation of this as the MRO data point is two days earlier than the VLBI 43 GHz data point.

Fig. 3. Light curves of 3C120 and 0716+714 using both the MRO 37 GHz data and VLBI 43 GHz (core + jet) data between 2013 and 2019.

There are a few data points, where the flux densities were measured within 24 hours of each other. We took the average difference of these as the approximate emission excess considering that intra-day variability may cause some inaccuracies. Due to there being a limited number of such points, the error margins may be higher.

Considering this excess, the baseline of 3C120 in 37 GHz should be at around 1 Jy, as the average difference of 10 data points was 0.92 Jy. We see it briefly decrease to that position during 2013, but otherwise it appears to remain close to 2 Jy. This may potentially be an effect caused by emission from the core and overlapping of flares.

In the case of 0716+714, the baseline drops to near the MRO nominal detection limit at multiple points despite the approximate mean-baseline level being above it (Fig. A.1). One would then assume that the 43 GHz emission would be very similar to the 37 GHz light curve as it appears that there is little diffuse emission visible to the single dish observations. Indeed, both flux densities follow each other closely in the 2013 and 2019 segment with an average difference of 0.04 Jy for 8 data points (Fig. 3).

6.2.2. Knot timescales

In the radio regime, multiple authors have been able to temporally relate flaring and VLBI knots (e.g. Marscher & Gear 1985; Turler et al. 2000; Savolainen et al. 2002; Lindfors et al. 2006). In order to examine the connection between the obtained characteristic timescales and jet knots, we decided to find the maximum timescale an individual knot is visible from the 43 GHz sample containing approximately 10 years of data between 2008 and 2018. The characteristic timescale in the PSD should in principle tell us for how long a flare affects the entire light curve. Thus, it seems logical that a knot with the longest appearance would have some correlation with the bend frequency.

We chose our sample for the comparison based on the constrained sources and our additional example sources 3C84 and BL Lac. Three sources (PG0007+106, 0736+174, and 2144+092) were not in the VLBI 43 GHz sample. We determined the approximate knot timescales by using Table 8 from Weaver et al. (2022). We compared the first observation for each knot to their last observation. The maximum knot timescale was then the maximum value found by this method. These should be considered as minimums for the longest knot timescales, that is the longest actual knot timescale may exceed these values.

Table 1 lists the compared sources with their characteristic timescales and VLBI knot timescales. The reported 90% confidence regions in Table 1 differ from the values in Table A.1 as we are now considering only the uncertainty in the timescale associated with the best-fit slope β obtained from the best bending power-law fit. As discussed in Sect. 5, not all combinations of slope and timescale are possible and correspondingly the confidence regions vary. In Fig. 4 the timescale constraints reported in Table A.1 are plotted in grey and they depict the minimum and maximum values for the timescales considering all possible slope-timescale combinations within the 90% confidence region. The best-fit characteristic timescales with error bars corresponding to the 90% confidence regions using only the obtained best-fit slope β are plotted in colour.

Fig. 4. Scatter plot showing the results from Table 1. The knot timescale is on the x-axis and the characteristic timescale on the y-axis. The coloured error bars are based on the 90% characteristic timescale confidence regions associated with the best-fit slope β obtained from the best bending power-law fit. The grey error bars are the confidence regions of the timescale when the slope is allowed to vary. For each source, except 2230+114, the knot timescale and characteristic timescale associated with the fixed best-fit slope β are consistent within uncertainties given the sparse parameter grid.

Sources 0235+164, 0415+379, 0716+714, and 4C29.45 have a best-fit timescale (near) equal to the maximum knot timescale, considering the sparse parameter grid. For all sources, except 2230+114, the timescales are close to the knot timescale within error margins. Examining the heat map of 2230+114 (Fig. A.7), we see that a timescale of 2500 days – close to the knot timescale – still gives a relatively good fit using a different slope β. However, this is the only source where the obtained timescale is not close to the knot timescale within error margins when using the best-fit slope.

These preliminary results are encouraging and warrant further investigation into the potential connection between the characteristic timescale and the maximum time a jet knot is visible. We discuss the reasoning behind the connection further in Paper II.

6.3. Similar studies in the radio domain

The 37 GHz data consists of emission arriving from closer to the central supermassive black hole compared to the two other long-term radio monitoring programmes by the University of Michigan Radio Astronomy Observatory (UMRAO, 14.5 GHz) and Owens Valley Radio Observatory (OVRO, 15 GHz). Park & Trippe (2017) note that for higher observing frequencies the slope of the PSD should be flatter but the characteristic timescale at a higher frequency (shorter timescale). The reason for this are the shorter cooling times of higher energy particles (e.g. Marscher & Travis 1996).

6.3.1. Comparison to 15 GHz

Max-Moerbeck et al. (2014) used the OVRO 15 GHz data to analyse correlations between radio domain and gamma-ray data observed with the Fermi LAT Gamma-Ray Space Telescope. They used four years of radio domain data from 2008 to 2012. Park & Trippe (2017) used UMRAO data (4.8, 8, and 14.5 GHz) to analyse the long-term variability of AGNs. Due to the significantly shorter light curves used by Max-Moerbeck et al. (2014), we only compare our results to the analysis with UMRAO data.

For some of the sources analysed by Park & Trippe (2017) using the lower 14.5 GHz data, the results are similar to what we obtained in our analysis. For our results, the mean of the PSD slope, when using the simple power law, was $\overline{{\beta }_{\mathrm{spl}}}=1.6$ which is the same result that Park & Trippe obtained. They noted that their mean result of 1.6 was flatter than expected and that it might have to do with them using the Lomb-Scargle periodogram. However, Max-Moerbeck et al. (2014) and Ramakrishnan et al. (2015), to which they compare their results to, used much shorter observing periods because they were interested in comparing the 37 GHz variability to gamma-ray variability. Thus, it is likely that the flatter mean slope is caused by the differences in the monitoring periods and that it is in fact result from the bend frequency influencing the simple power-law fit for some of the sources (see Sect. 6.3.2).

Our results are consistent with Park & Trippe (2017) in the number of sources that differ from the simple power law. From their smaller sample of 43 sources, only 4 were inconsistent with a simple power law. These sources are 0235+164, 3C120, 4C29.45, and 3C454.3. The sources for which we have differing results are 0415+379, 0716+714, PKS1749+096, and 2230+114. For most of these sources their best-fit simple power-law slope β estimate is consistent with ours varying between β = 1.4 and β = 1.6 in 14.5 GHz (or 8 GHz if the higher frequency is not available). Only for 2230+114, their best-fit slope β = 2.0 is clearly steeper than our simpler power-law estimate β = 1.6.

Park & Trippe (2017) were also concerned about how well a timescale in the radio frequencies can be discovered. They simulated pulses of varying length, amplitude and arrival times and concluded that a higher degree of superposition influences the timescale, shifting it to lower frequencies. They simulated pulses with rise and decay timescales which they multiplied by a factor τ to represent the timescale that would control the degree of overlap. Their conclusion was that the PSDs of AGN light curves would follow broken power laws intrinsically but that a high degree of superposition would move the timescale to lower frequencies.

Their analysis results are in line with what we assume the characteristic timescale to be. That is, the timescale describes how long a flare can potentially influence the light curve due to superpositioning. At some point a flare must decay sufficiently to not contribute to the light curve mean and variance anymore and that in our understanding is the characteristic timescale, which in the PSD appears as a bend frequency.

6.3.2. Comparison to five years of MRO data

Ramakrishnan et al. (2015) studied the Fermi LAT Gamma-Ray Space Telescope γ-ray and 37 GHz radio correlations. They defined the simple power-law fits for 55 sources using five years of data observed by MRO. The short segment was chosen according to the coincidence with γ-ray observations. We would expect the five-year results to be close to the slope of the bending power law in our long-term data with constrained timescales. Indeed, the matching sources 0235+164, 0716+714, PKS1749+096, 2230+114, and 3C454.3 have the same estimated slopes (within error margins) for both five years of data and for their full monitoring period when using the bending power law, though the parameter space for the slopes is large. Interestingly though, one would expect the bend to already affect the fit especially with 0716+714, whose best-fit characteristic timescale is 500 days, well within five years of observations. The reason for this discrepancy may be in the individual decisions discussed earlier, such as bin size or scaling of light curve.

7. Conclusions

We analysed the long-term radio variability of 123 sources observed by the Metsähovi Radio Observatory in 37 GHz. We utilised the longest monitoring periods in the 37 GHz band to date, maximum of which extended to 42 years of observations with an average monitoring period of 34.5 years. Our sample of 123 sources was also exceptionally large for such analysis, allowing unique insights into the long-term radio variability of AGNs.

We fitted the periodograms of each source with both the bending and simple power-law models. In addition, we compared the obtained timescales to the VLBI 43 GHz sample and the maximum knot timescales. Our main results were:

1. We found a well-constrained timescale and PSD slope for 11 sources. The obtained timescales varied between 500 and 3000 days and the slopes between 1.8 and 2.8.

2. Constraining the parameters was challenging and the result was always a large parameter space. We used a sparse grid for the parameters due to the uncertainty of the benefits of using a finer grid without prior knowledge of how well a timescale can be constrained.

3. In addition to the bending power law, where the low-frequency slope flattens to zero, we used a high-frequency bend model with β_low = 1. We were unable to differentiate which model fits the data better. This appears to be a feature of the models and it may not be possible to differentiate between the two for observed radio data.

4. We found a preliminary correlation with the characteristic timescale and the maximum knot timescale of the jet. This suggests that the characteristic timescale is connected to the duration of the flares.

5. The results from the PSD fits suggest that radio-frequency observations require long monitoring periods, sometimes exceeding 40 years of monitoring especially when the data are sparse and unevenly sampled.

In future work, we will analyse the constrained sources further with finer parameter grids and proper error estimates for the PSD bin powers. We will also explore the VLBI knot maximum-timescale connection with the characteristic timescale using more sources in order to confirm the found correlation statistically.

Acknowledgments

The authors would like to thank Dr. Kari Nilsson for his valuable comments on the draft. S.K. was supported by Jenny and Antti Wihuri Foundation, Väisälä Fund and Academy of Finland project 320085. T.H. was supported by Academy of Finland projects 317383, 320085, 322535, and 345899. This publication makes use of publicly available data from the Metsähovi Radio Observatory (https://www.metsahovi.fi/opendata/), operated by Aalto University in Finland.

References

Appendix A: Heat maps and periodogram fits for the constrained sources and example comparison sources, with a result table

Fig. A.1. Results from the periodogram analysis of source 0716+714. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.2. Results from the periodogram analysis of source 2144+092. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.3. Results from the periodogram analysis of source PG0007+106. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.4. Results from the periodogram analysis of source 0235+164. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.5. Results from the periodogram analysis of source 0736+017. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.6. Results from the periodogram analysis of source PKS1749+096. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.7. Results from the periodogram analysis of source 2230+114. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating best fit and magenta worst. The bottom row shows the source light curve.

Fig. A.8. Results from the periodogram analysis of source 0415+379. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.9. Results from the periodogram analysis of source 4C29.45. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.10. Results from the periodogram analysis of source 3C454.3. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.11. Results from the periodogram analysis of source 3C84. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

Fig. A.12. Results from the periodogram analysis of source BL Lac. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.

All Tables

All Figures

	Fig. 1. Comparison between the simple power law and the bending power law. The simple power law in black is a straight line of slope βspl = 2. The bending power law in orange follows the simple power law with the same slope β = 2 until the bend, after which the slope turns over to β = 0.
In the text

Fig. 2. Comparison between the results for source 3C120, when using a shorter segment of data and the entire monitoring period. The first row of plots uses data between 2010 and 2015, and the middle row uses the entire light curve shown in bottom row, where the red portion corresponds to the five-year segment. The upper left-side heat map shows the PSRESP fits for each bending power-law parameter combination using five years of data. The white graph overlaid with the heat map shows the simple power-law fit with the x-axis being equivalent to the heat map x-axis. The y-axis for the white graph equals the heat map legend range scaled to agree with the label positions. The right-side plot shows the source periodogram (black) and each simulated mean periodogram with cyan indicating the best fit and magenta the worst fit. The heat map and PSD graph below show the equivalent results using the entire light curve of 3C120 shown on the bottom row. The dotted vertical lines on both periodogram graphs show the best-fit bend frequency obtained by using the entire light curve, and the dashed horizontal lines show the noise floor.

In the text

	Fig. 3. Light curves of 3C120 and 0716+714 using both the MRO 37 GHz data and VLBI 43 GHz (core + jet) data between 2013 and 2019.
In the text

Fig. 4. Scatter plot showing the results from Table 1. The knot timescale is on the x-axis and the characteristic timescale on the y-axis. The coloured error bars are based on the 90% characteristic timescale confidence regions associated with the best-fit slope β obtained from the best bending power-law fit. The grey error bars are the confidence regions of the timescale when the slope is allowed to vary. For each source, except 2230+114, the knot timescale and characteristic timescale associated with the fixed best-fit slope β are consistent within uncertainties given the sparse parameter grid.

In the text

	Fig. A.1. Results from the periodogram analysis of source 0716+714. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.2. Results from the periodogram analysis of source 2144+092. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.3. Results from the periodogram analysis of source PG0007+106. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.4. Results from the periodogram analysis of source 0235+164. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.5. Results from the periodogram analysis of source 0736+017. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.6. Results from the periodogram analysis of source PKS1749+096. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.7. Results from the periodogram analysis of source 2230+114. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating best fit and magenta worst. The bottom row shows the source light curve.
In the text

	Fig. A.8. Results from the periodogram analysis of source 0415+379. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.9. Results from the periodogram analysis of source 4C29.45. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.10. Results from the periodogram analysis of source 3C454.3. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.11. Results from the periodogram analysis of source 3C84. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text

	Fig. A.12. Results from the periodogram analysis of source BL Lac. The upper left-side plot shows the heat map overlaid with the simple power-law fit in white. The plot on its right is the periodogram of the source including all of the mean periodograms of the simulated bending power laws with cyan indicating the best fit and magenta the worst. The bottom row shows the source light curve.
In the text