© The Authors 2026

1. Introduction

Active galactic nuclei (AGNs) are highly luminous centres of active galaxies and are known for their multi-wavelength variability. In the radio domain, the variability has been attributed to shocks in jets (Marscher & Gear 1985), where the observed luminosity and rate of variability are boosted by relativistic flow velocities and the jet orientation being close to the line of sight. The continuum observations of flares have been connected to very long baseline interferometry (VLBI) observations of ejected features in the jet (e.g. Türler et al. 1999; Savolainen et al. 2002; Lindfors et al. 2006).

The time variability of AGNs has been studied extensively in all wavelength domains (e.g. Hughes et al. 1992; Uttley et al. 2002; McHardy et al. 2006; Hovatta et al. 2007; Kelly et al. 2009; Ramakrishnan et al. 2015). Typically, variability studies attempt to look for so-called characteristic timescales of variability. A characteristic timescale describes the typical variability in the light curves, such as how often the source flares and for how long the flaring lasts. These timescales may be connected to the physical conditions in AGNs, and they are an important piece in understanding the drivers of AGN variability over different wavelengths.

One typical way to estimate characteristic variability from observed light curves is to analyse their power spectral densities (PSDs). The PSD gives the power content of a signal (or noise) in frequency space, and it can be estimated using the periodogram. The PSD of an AGN is typically modelled as a red noise power law with a potential bend after which the slope turns over to a flatter power-law slope or converges to white noise in low frequencies. This observed low-frequency flattening is expected in real physical processes to avoid infinite integrated power (e.g. Van der Klis 1989).

A bend in the bending power law is considered to be the inverse of a physically meaningful characteristic timescale; in black hole X-ray binaries (BHXRBs), the emergence of jets is associated with changing X-ray spectral states (e.g. Fender et al. 2004), and their light-curve PSDs have either one or two observable bends (e.g. Belloni & Hasinger 1990; McHardy et al. 2006; Wilkinson & Uttley 2009). The location of the bend has been connected to central engine parameters of both BHXRBs and AGNs (e.g. Uttley et al. 2002; Papadakis 2004; McHardy et al. 2006; Gonzalez-Martin & Vaughan 2012; Paolillo et al. 2023). Specifically, McHardy et al. (2006) found the X-ray PSDs of radio-quiet AGNs to be similar to the soft-state PSDs of BHXRBs with a single bend. They found the location of the PSD bend to scale with the black hole mass and mass accretion rate. This scaling was also found to apply to the PSD of the blazar 3C273 (McHardy 2009). Thus, of interest is to analyse whether the PSD bend of radio light curves can also be associated with the central engine.

In Kankkunen et al. (2025) (hereafter Paper I), we were able to identify 11 sources for which there was a constrainable characteristic timescale. Due to the small sample of sources with a PSD bend, we attempt to extract the characteristic timescale from the light curves via light-curve decomposition. We test the hypothesis that the timescale we obtain from the PSD, the inverse of the bend frequency, is connected to either the durations of flares or to their temporal separations. The idea that the longest timescales are related to the flare durations is supported by the preliminary analysis in Paper I, where we compared the PSD timescales to the longest times a knot component was visible in very long baseline interferometry (VLBI) 43 GHz data (Weaver et al. 2022). Similarly, Park & Trippe (2017) found that introducing longer-duration flares to their simulated light curves moved the bend location to lower frequencies. On the other hand, Mukherjee et al. (2019) used simulated data and light curves to analyse the connection between the PSD timescale and the light-curve characteristics and found the PSD bend to be the inverse of the intervals between consecutive shocks moving down the jet. Thus, it is not clear, which timescale observed in the light curve corresponds best to the PSD bend, and different timescales may also be blended. Nevertheless, if a connection is found, we can use the corresponding timescale as a proxy for the PSD bend and test for a connection with the central engine with a larger sample of sources.

In addition to the connections between the PSD timescale and the central engine, we analyse whether the jet and central engine parameters can be related to the emission region size. For example, a study by Potter & Cotter (2015) found compelling evidence of a connection between the bulk Lorentz factor of the jet and the radius of the jet transition region. The transition region is defined as the region where the jet shape changes from parabolic to conical and is in equipartition. A transition between the jet shapes has been observed in multiple AGNs (e.g. Asada & Nakamura 2012; Kovalev et al. 2020; Boccardi et al. 2021; Okino et al. 2022). Kovalev et al. (2020) suggest that the transition region likely coincides with the location of a recollimation shock. In millimetre wavebands, the radio core is often considered to be a standing recollimation shock (Daly & Marscher 1988), although in some sources it may still be the τ = 1 surface (e.g. Kutkin et al. 2014; Fromm et al. 2015). With simplifications and assuming a direct relationship between the length of the parabolic region and the width of the jet shape transition region (e.g. Vlahakis & Königl 2003; Beskin & Nokhrina 2006; Lyubarsky 2009; Potter & Cotter 2015), we incorporate light-travel time arguments to use the decomposed flares and their emission-frame rise times as proxies for the width of the radio core. We estimate whether the rise times can be related to the jet bulk Lorentz factor, the black hole mass, the accretion disk luminosity, and/or the mass accretion rate.

Our sample includes both the constrained sources of Paper I as well as a larger sample of additional sources for which a constrained timescale was not found. Comparing the results between the decompositions and the PSD timescales, or lack thereof, should help in further understanding the characteristic variability and its dependence on jet physics. The paper is divided as follows: In Section 2, we describe the used data and sample. In Section 3, we describe the methods, followed by the results in Section 4. In Section 5 we discuss the results and their connection to previous studies. Section 6 gives the conclusions obtained from the analysis.

2. Data

The Aalto University Metsähovi Radio Observatory (MRO) hosts a 14-metre telescope that is used to observe AGNs predominantly in 37 GHz. Observations have been conducted regularly for over 40 years, with a few sources having their first observation already in 1979. Because of this long-term monitoring, the MRO sample of light curves allows us to analyse the time variability of the most energetic relativistic jets with high accuracy. More detailed information on MRO and the observed sample is presented in Teräsranta et al. (1998) and in Paper I.

In Paper I, we compiled a sample of 123 sources from the objects monitored at MRO, which had a minimum of 10 years of consecutive observations, a minimum of 100 data points, and at least a 1 Jy maximum flux density. For this work, we constructed a subsample by including the 11 constrained sources with a constrained timescale obtained from the PSD analysis, and a subset of 40 well-sampled unconstrained sources as a further comparison sample. We also included three borderline sources from Paper I, which are sources that had a timescale in the PSD but did not meet all strict criteria to be considered reliable. The additional sources in the analysis were chosen from Paper I with the following criteria: a minimum of 500 data points, no multiple long gaps in the light curve, and prominent flares for the decomposition analysis. Including this larger sample of sources should allow for further confirmation of their timescales and potentially relate the results to the central engine of the sources. The entire sample consists of 54 sources and they are listed in Table A.1.

Because we are interested in comparing characteristic timescales obtained from the PSD to the light-curve decomposition, we used the same cut-off date as in Paper I (1 January 2023), leading to a median light-curve length of ∼39 years. The source types in the sample can be divided based on optical classification into flat-spectrum radio quasars (FSRQs), BL Lac objects (BLOs), and radio galaxies (GALs).

3. Methods

3.1. Light-curve decomposition

To connect our results from the PSD analysis in Paper I to the specific timescales in the light curves, we applied a decomposition algorithm. Decomposing the light curves into individual flares has been shown to be a viable method for estimating radio light curve parameters (e.g. Lähteenmäki & Valtaoja 1999; Hovatta et al. 2009), and it allows us to compare the light curve characteristics to the PSD timescales. Using both the results from the flare decomposition and the PSD analysis may help in understanding whether the light curve is sufficiently long to describe the source behaviour. If no bend is visible in the PSD, it suggests that the amplitude of the long-term variations in the observed light curve is still increasing (e.g. Press 1978), and thus the behaviour of the source cannot be fully characterised.

The Magnetron¹ code (Huppenkothen et al. 2015) decomposes light curves into separate flares, with an exponential rise and decay, with no a priori information on the number of flares or the parameters of the exponential function. The algorithm uses Bayesian hierarchical modelling and diffusive nested sampling (DNest4, Brewer et al. 2011; Brewer & Foreman-Mackey 2018) to determine the number of components and the parameters of individual components independently of each other and to create samples from the posterior. For the list of priors and other assumptions used, we refer the reader to Huppenkothen et al. (2015).

As described in detail in Liodakis et al. (2018), the modification of the Magnetron software to blazar analysis uses the Ornstein-Uhlenbeck (OU) process to model the underlying background variability unrelated to flaring events. The OU process is a mean-reverting stochastic stationary process (Uhlenbeck & Ornstein 1930). If no OU process is modelled for the underlying large-scale variability, the program tries to fit each small variation in the light curve with an exponential flare. The OU background should thus help in avoiding fitting exponential flares to minor variations due to turbulent processes or to the extended diffuse emission from the jet.

Figure 1 shows an example of a single realisation of the decomposition drawn from the posterior sample of source 0415+379. The posterior sample includes several possible decompositions with different numbers of individual flares and varying flare parameters. This way the method accounts also for the uncertainty in the number of flares in any given light curve as often multiple variations of flare superpositions fit the light curve equally well. The number of obtained posteriors varied between ∼250−1250, with a median of 517 samples across all sources. For each observed light curve, the plots of the posterior weights were monitored to decide when a run was successfully finished (see Brewer & Foreman-Mackey 2018). All of the decomposed flares are described by their peak position and amplitude, e-folding timescale (time constant) and skewness. The e-folding timescale is related to the rise of the flare, whereas the decay time is obtained from multiplying the timescale with the skewness factor.

Fig. 1. Example of one decomposition of the 37 GHz light curve of 0415+379 drawn from the posterior sample. The OU process has been subtracted from the decomposition. The flares coloured in blue and cyan have an amplitude higher than the mean amplitudes of all flares in the posterior, and the flares coloured in cyan have amplitudes above 50% of the maximum amplitude in the given posterior. The magenta line shows the fit of the decomposition. In the top right, the inset shows an example of the portion of the flare used to determine the flare duration (cyan) and the portion of the flare used for the rise time (magenta).

Upon visual examination of some of the decompositions, we noticed that some posterior draws included suspiciously large spikes and occasionally extremely low flare amplitudes. The Diffusive Nested Sampling algorithm (Brewer & Foreman-Mackey 2018) used by the Magnetron code, is a Markov chain Monte Carlo method with inherent statistical randomness, so it may be possible that a different run of the code results in a different fit. Therefore, we ran the decomposition on some of the sources twice and found that the results of the test runs were within the standard deviations of the original run. Thus, we concluded that the posterior distributions were sufficiently consistent between different runs. To mitigate the effect of the fast spikes, we required the e-folding timescale of each flare to be longer than the median cadence in the light curve for the flare to be included in the analysis.

3.1.1. Flare durations and separations

We estimate whether the PSD timescale can be associated with something easily measurable in the decomposed AGN light curves. It is possible that multiple timescales can be extracted from the PSD; however, in Paper I, we only probed for a single low-frequency bend in the PSD. To find the potential timescale associated with the PSD timescale, we test two simple parameters obtained from the light curves: the flare durations and the separations of their peaks. The separations obtained from the peaks roughly correspond to the intervals between flare onsets, with the exception of a few sources that only include a small number of exceptionally long flares whose onsets are not included within the observing windows.

Because the red noise nature of AGN PSDs indicates that the amplitude of variations is increasing over longer and longer timescales (e.g. Press 1978), we focus the analysis on the brightest flares. We analyse those flares that have an amplitude above the mean amplitude in a given posterior draw. This decision allows a sufficient number of flares to be included in analysing both the mean durations and the flare peak separations; however, this criterion may underestimate the maximum separations in case there are only a small number of very bright flares. An example of flare decomposition with the mean-amplitude criterion is given in Fig. 1. The figure also shows an example of the flare separation, visualised between two adjacent flares. The durations are determined from when the flare has reached 1% of its maximum amplitude on the rise and decay sides. This limit was also used in Huppenkothen et al. (2015) with count rates. To determine the timescales from the durations and separations of the most dominant flares, we took their means for each posterior, and then similarly the means from all drawn posteriors.

3.1.2. Flare rise times

The brightest flares have been linked to so-called core flares associated with ejections of VLBI components (Savolainen et al. 2002). As a simplified assumption, the rise times of the brightest flares should give the best representation of variability timescales within the radio core. We exclude the fainter flares as they have been associated with knots further downstream of the core (e.g. Savolainen et al. 2002). With the following reasoning and light-travel time arguments, we use the e-folding timescale of the most dominant flares as a proxy for the upper limit of the relative emission region size. The e-folding timescale corresponds to approximately 37% of the flare maximum amplitude, and we quote it as the rise time. This portion of the flare is indicated in Fig. 1.

Here, we narrow down the number of flares included to only consider the ones whose peak amplitudes exceed 50% of the highest-amplitude flare in a given decomposition. This is to ensure that a minimal number of flares being formed further downstream are included in estimating the size of the radio core. From each posterior draw, we take the median rise time of the considered flares and then take the median over the sampled posteriors. To scale the rise times to the emission frame, we used

$$\begin{array}{c}\hfill {\tau }_{\mathrm{rise},\mathrm{em}}=\frac{{\tau }_{\mathrm{rise},\mathrm{obs}}·{\delta }_{\mathrm{var}}}{1+z},\end{array}$$(1)

where τ_rise, obs is the decomposed rise time in the observed frame, δ_var the variability Doppler factor, and z the redshift.

To be able to compare the rise times of the flares and thus the sizes of the emitting regions between sources, we made several simplifications in our assumptions. First, we assume that synchrotron opacity has a negligible effect on the flare timescales. In addition, to compare the obtained results to previous studies, we assume that the emitting region corresponds to the location of the jet transition region and that the upper limit of the emitting region size is a sufficient proxy for the jet width at that location.

3.2. Doppler factors

Lähteenmäki & Valtaoja (1999) discussed different methods for estimating Doppler factors. One approach involves directly measuring the brightness temperatures from VLBI data and comparing them to the theoretical maximum equipartition brightness temperature (Readhead 1994). Another option is to use total flux density flares through light-curve decomposition (e.g. Terasranta & Valtaoja 1994; Lähteenmäki & Valtaoja 1999; Hovatta et al. 2009), where the intrinsic brightness temperature T_b, int is compared to the variability brightness temperature T_b, var.

We estimate the Doppler factors following the prescription in Liodakis et al. (2021), using the cosmological parameters Ω_m = 0.315, Ω_Λ = 0.685, and H₀ = 67.4 km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2020) as

$$\begin{array}{c}\hfill T_{\mathrm{b},\mathrm{var}}=1.47·{10}^{13}\frac{d_{\mathrm{L}}^2\mathrm{\Delta }{S}_{\mathrm{ob}}(\nu )}{{\nu }^2{t}_{\mathrm{var}}^2{(1+z)}^4}K,\end{array}$$(2)

where d_L is the luminosity distance in Mpc, ΔS_ob(v) the flare amplitude in Jy, ν the observing frequency in GHz, t_var the flare rise time in days, and z the redshift (Liodakis et al. 2021). From each posterior draw, the maximum T_b, var is chosen. We use the e-folding timescale as the rise time and consequently the difference between the maximum flux density of a flare and the flux density at the e-folding timescale as ΔS_ob(v). To account for the uncertainty in T_b, int, we follow Liodakis et al. (2018) and use a Gaussian distribution with mean μ = 2.78 ⋅ 10¹¹ K and standard deviation σ = 0.72 ⋅ 10¹¹ K and draw samples from it. The variability Doppler factor δ_var was obtained from

$$\begin{array}{c}\hfill {\delta }_{\mathrm{var}}=(1+z)\sqrt[3]{\frac{T_{\mathrm{b},\mathrm{var}}}{T_{\mathrm{b},\mathrm{int}}}}\end{array}$$(3)

by random sampling from the distribution of T_b, var and T_b, int and constructing the median. The δ_var and the respective 68% confidence intervals are reported in Table A.1. Two sources, 3C84 and ON231, had a Doppler factor less than one, suggesting that their intrinsic brightness temperature is likely lower than indicated in the distribution derived in Liodakis et al. (2018). Thus, we excluded the values from the table. We also note that the Doppler-factor estimate for 0836+710 is high, with a value δ_var = 78.6, and it is possibly an artefact of the decomposition procedure due to some gaps in the beginning of the light curve.

3.3. Lorentz factor

We derived the Lorentz factors for each source through the calculated Doppler factors and the apparent jet speeds (β_app) from the MOJAVE survey (Lister et al. 2018, 2019) using

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{{\beta }_{\mathrm{app}}^2+{\delta }_{\mathrm{var}}^2+1}{2{\delta }_{\mathrm{var}}}·\end{array}$$(4)

As there is no reported redshift for 0716+714 in Lister et al. (2019), we adopt z = 0.2304 from Pichel et al. (2023). For S20109+22, we used z = 0.49 (Koljonen et al. 2024).

Through Eqs. (3) and (4), we were able to estimate the Lorentz factor for 50 sources, reported in Table A.1. The given apparent speeds in units of c have been converted using Ω_m = 0.315, Ω_Λ = 0.685, and H₀ = 67.4 km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2020).

3.4. Central engine parameters

To compare the size of the emission region to the central engine parameters of the sources, we estimated the accretion disk luminosities, L_acc, and black hole masses, M_BH, of the sources mainly using broad-line luminosities (Zamaninasab et al. 2014 and references therein). When available, we use the Hβ line for the accretion disk luminosity, otherwise the MgII line or the OIII line.

3.4.1. Black hole mass

The observed optical continuum emission of radio-loud AGNs includes synchrotron emission from the jet; thus, we needed to use the broad-line luminosities to obtain an estimate for the disk emission. This was done through a power-law relationship between the line luminosities and the continuum, which is constructed through observations of radio-quiet sources. To obtain an estimate for the λL_5000 Å continuum emission, we used the formula from Liu et al. (2006) and for the λL_3000 Å continuum emission from Shen et al. (2011). The estimate for λL_3000 Å was measured from total line luminosities, but according to Shen et al. (2011), the power-law relationship is similar enough that broad emission lines can be used. We only substituted the observed continuum emission with the estimate from the broad emission lines, if the line estimate gives a smaller luminosity.

The formulae for an estimate of M_BH are (Zamaninasab et al. 2014, and references therein)

$$\begin{array}{cc}& M_{\mathrm{BH}}\simeq 4.82{\left[\frac{\lambda L_{5100\mathrm{\AA }}}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right]}^{0.69}{\left[\frac{{\mathrm{FWHM}}_{\mathrm{H}\beta }}{\mathrm{km}{\mathrm{s}}^{-1}}\right]}^2{M}_{\odot },\hfill \end{array}$$(5)

$$\begin{array}{cc}& M_{\mathrm{BH}}\simeq 3.37{\left[\frac{\lambda L_{3000\mathrm{\AA }}}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right]}^{0.47}{\left[\frac{{\mathrm{FWHM}}_{\mathrm{MgII}}}{\mathrm{km}{\mathrm{s}}^{-1}}\right]}^2{M}_{\odot },\hfill \end{array}$$(6)

where FWHM_Hβ is the width of the Hβ, and FWHM_MgII is the width of the MgII emission line.

The majority of the emission lines and full width at half maximums (FWHMs) are obtained from Torrealba et al. (2012), supplemented by values from Shaw et al. (2012), Stickel et al. (1993), and Cohen et al. (1987). We obtained an estimate for the black hole mass for 30 sources; the references for the emission lines are reported in Table A.1. We do not include black hole masses derived through other methods to avoid any additional biases in the analysis.

3.4.2. Accretion disk luminosity and mass accretion rate

The formulae to calculate the accretion disk luminosity L_acc from emission-line luminosities are (Zamaninasab et al. 2014, and references therein)

$$\begin{array}{cc}& {log}_{10}{L}_{\mathrm{acc}}=(12.32\pm 0.32)+(0.78\pm 0.01){log}_{10}{L}_{\mathrm{H}\beta },\hfill \end{array}$$(7)

$$\begin{array}{cc}& {log}_{10}{L}_{\mathrm{acc}}=(16.76\pm 0.26)+(0.68\pm 0.01){log}_{10}{L}_{\mathrm{MgII}},\hfill \end{array}$$(8)

$$\begin{array}{cc}& {log}_{10}{L}_{\mathrm{acc}}=(26.50\pm 0.32)+(0.46\pm 0.01){log}_{10}{L}_{\mathrm{OIII}}.\hfill \end{array}$$(9)

The mass accretion rates can be calculated from the disk luminosities following the relation for thin disks (Shakura & Sunyaev 1973):

$$\begin{array}{c}\hfill \dot{M}=\frac{L_{\mathrm{acc}}}{\eta c^2},\end{array}$$(10)

where Ṁ is the mass accretion rate, L_acc the accretion disk luminosity, and η the radiative efficiency. We adopt the value η = 0.4, similarly to Zamaninasab et al. (2014), albeit it is somewhat unimportant in the upcoming analysis as the efficiency is kept the same for each source. It is noteworthy that this formula is intended for thin disks, and as such it may not work for BLOs. Nevertheless, because we are using the same scaling factor for each source, the mass accretion rate will have the same relationship with the observed values as the accretion disk luminosity. We obtained the accretion disk luminosity for 34 sources.

We calculated the normalised mass accretion rates using the accretion disk luminosity and the black hole mass through the formula

$$\begin{array}{c}\hfill \dot{m}=\frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}=\frac{L_{\mathrm{acc}}}{\eta L_{\mathrm{Edd}}},\end{array}$$(11)

where $L_{\mathrm{Edd}}\approx 1.3·{10}^{47}\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}$ is the Eddington luminosity (e.g. Böttcher et al. 2012).

We note that the accretion disk luminosities obtained through the methods discussed in Zamaninasab et al. (2014) give a larger value than those obtained through spectral energy distribution (SED) modelling (Ghisellini et al. 2011, 2014). This also results in the normalised mass accretion rates being relatively large. However, this should not affect the results of our correlation analysis as the values are calculated in the same way for each source.

4. Results

We decomposed the light curves of 54 sources to estimate the mean flare duration and separation of the brightest flares of each source. These timescales are given in Table A.2.

We used the Kendall’s τ rank correlation to test for monotonic associations between parameters. We opted for a typical significance level of 0.05, but we also report the exact p-values. We considered a correlation weak if 0.1 < τ < 0.3, moderate if 0.3 ≤ τ < 0.6, and strong if τ ≥ 0.6.

4.1. Flare durations and separations

In Paper I, the characteristic timescales were fitted with a bend-timescale grid of 500 days up to 7000 days, with an additional coverage of 100, 200, 400, and 800 days. The coverage was chosen to be sparse because it is difficult to extract exact values from the periodogram, and the uncertainties were correspondingly large. Figure 2 shows the results from the flare decomposition against the obtained PSD timescales. The plots suggest that the flare durations have a better overall one-to-one agreement with the PSD timescale. To quantify this, we performed the Kendall rank correlation test between the PSD timescales and the flare durations and separations. Because most of the PSD timescales of the unconstrained sources are near or at the upper limit of the tested scale (7000 days), we conducted two tests, one with only the constrained sources and one with all of the sources with an estimated PSD timescale. This analysis specifically probes whether the flare durations and/or separations increase with the PSD timescale. We assumed that for the unconstrained sources the PSD timescale is longer than for the constrained sources. The correlation analysis results are reported in Table 1.

Fig. 2. Power spectral density (PSD) timescales against the mean flare durations and separations. The coloured data points are for the constrained sources and the grey ones for the unconstrained sources. The dashed line indicates the one-to-one correspondence between the timescales.

We obtained significant correlations with both the flare durations and separations when using all of the sources. For the 11 constrained sources, the correlation is not significant with the flare durations, likely owing to the two sources 3C120 and 3C454.3 with shorter flare durations compared to the obtained best-fit PSD timescales. For the separations, the correlation with the PSD timescale is strong but does not follow the one-to-one correspondence. Overall, for five of the constrained sources, the mean flare duration matches the PSD timescale better than the flare separation. For four of the sources, both timescales are very close to each other. Mukherjee et al. (2019) show that if both separations between the flares and the electron cooling times (resulting in longer flare decays) are of a similar magnitude, they may both contribute to the bend causing a smoother turnover. The constrained sources are also within the 25 sources with the shortest flare durations. Using the flare separations, the scatter is larger.

The flare durations are, on average, twice the length of the separations. However, for the 11 constrained sources this ratio is smaller with ∼1.3. Figure 3 shows the flare separations plotted against the flare durations for each light curve. The solid line is the equivalence of twice the duration compared to the separation, whereas the dashed line gives the one-to-one equivalence. We obtained a strong positive correlation, with τ ≈ 0.58 and p ≈ 7.9 ⋅ 10⁻¹⁰. Given the longer values obtained for the flare durations, it seems likely that the PSD bend constrained in Paper I is more closely related to the durations. Even if there was a second bend at higher frequencies, we expect the bending power law used in Paper I to rather estimate the low-frequency bend, where the flattening to white noise is expected. However, it is also possible that the flattening observed in some of the sources is not to white noise but to a slope of 1, in which case the flare durations could be probing the high-frequency bend instead. This is potentially the case with 3C120; The best-fit PSD bend of 3C120 is at 2500 days, but its mean flare duration and separation are only ∼500 and ∼1000 days, respectively. If only the very brightest flares are considered, that is, those with amplitudes over 50% of the maximum amplitude in a given posterior draw, the mean separation is ∼2600 days. It is then possible that using a two-bend model for the PSD of 3C120 could reveal a high-frequency bend better associated with the durations. We leave this analysis for future studies.

Fig. 3. Mean flare separations against the mean flare durations obtained from the decomposition of the light curves using the mean-amplitude limit. Source 3C84 is excluded from the plot for visual reasons, due to a much longer estimated mean flare duration and separation compared to the other sources.

Connection with the central engine

Due to the found connection between the light-curve decomposition results and the PSD timescales, we used both the mean flare durations and separations as potential proxies for the PSD timescale. This is because the number of constrained sources (11) is too low for the analysis and is reduced to only eight for the ones with a black hole mass estimate. Motivated by McHardy et al. (2006), who found the PSD timescales of AGNs in the X-ray domain to scale as $\approx M_{\mathrm{BH}}^{1.12}/{\dot{m}}^{0.98}$, we estimated the PSD timescale association with M_BH/ṁ using the black hole masses in units of solar mass (see Fig. 4).

Fig. 4. Black hole mass divided by the normalised mass accretion rate (MBH/ṁ) against the observed mean flare duration (left) and separation (right).

With a sample of 30 sources, we find that both the flare durations and separations have a moderate positive correlation with M_BH/ṁ (Table 2). The correlations may be driven by the black hole mass, as no correlations exist with the mass accretion rate ṁ and either timescale. However, for the flare separations, the correlation with M_BH is marginally weakened (τ = 0.254, p = 0.052) compared to M_BH/ṁ. Thus, we cannot definitively exclude the need for the division by ṁ and a larger sample of sources or more accurate estimates for the black hole mass and accretion rate are required to analyse the connection further.

If the flare durations and separations are corrected to the emission frame, the correlation is weakened for the flare durations and for the flare separations we cannot reject the null hypothesis of no association (Table 2). The implications are discussed in Sect. 5.2.

4.2. Jet and the central engine

We analysed the data for an association between the black hole mass and the width of the transition region approximated from the median rise time, τ_rise, of the most dominant flares. Fig. 5 shows the black hole mass against the rise time in both the observed frame τ_rise, obs and the redshift and Doppler-corrected emission frame τ_rise, em. We applied the Kendall rank correlation analysis to probe for an association between M_BH and τ_rise, em and were not able to reject the null hypothesis (Table 3). However, if the source 0235+164 with weak emission lines and thus uncertain black hole mass is removed from the analysis, a significant weak positive correlation, with τ ≈ 0.29 and p ≈ 0.03, is obtained. A similar correlation, including all of the sources, is obtained for τ_rise, obs and M_BH.

Fig. 5. Black hole mass against the median flare rise time in the observed frame (left) and the emission frame (right).

As in the previous analysis, we used τ_rise, em as a proxy for the emission-region size, and we calculated the accretion disk luminosity for 34 sources with estimated Doppler factors. Figure 6 shows in the top row the accretion disk luminosity L_acc against both the observed and the emission-frame rise times. We found a weak significant positive correlation between L_acc and τ_rise, em (Table 3). The result is the same as for the mass accretion rate Ṁ, when using the same radiative efficiency for all sources. We calculated the normalised accretion rate ṁ using the accretion disk luminosity and the calculated black hole masses for 30 sources. Figure 6 shows ṁ against the observed and emission-frame rise times. Here, we could not reject the null hypothesis of no correlation. Neither parameter has a significant correlation with τ_rise, obs.

Fig. 6. Accretion disk luminosity (top row) and normalised mass accretion rate (bottom row) against the observed median flare rise time (left) and emission-frame rise time (right).

With the 51 sources for which we were able to determine the Lorentz factor, and testing for an association with the emission-frame rise time τ_rise, em, we found a weak positive correlation (Table 3). No association was found with τ_rise, obs. Figure 7 shows the bulk jet Lorentz factor Γ against τ_rise, obs and against τ_rise, em.

Fig. 7. Lorentz factor against the median rise time in the observed frame (left) and the emission frame (right).

5. Discussion

We have analysed the decomposed light curves to understand which timescales could correspond to the PSD timescales found in Paper I and whether they can be connected to the black hole mass and accretion rate. We have also used a simplified assumption that the derived rise times can be used as a proxy for the emission region size of the jet and that it corresponds to flares observed to begin at the transition region. Here, we discuss some of the implications of the potential correlations as well as the caveats in the analysis. Because the number of BLOs and especially GALs is small in the sample, we cannot make strong inferences between the source types.

5.1. Associations with the emission region size

Kovalev et al. (2020) analysed VLBI data to look for jet shape transitions from parabolic to conical. They found the transition region in ten sources, three of which are included in our sample: 0415+379, 3C120, and BL Lac. Comparing the emission-frame rise times of the three sources to the jet widths given in Kovalev et al. (2020), the width of the transition region in parsecs agrees in ordering with the rise times obtained from the light-curve decomposition: In our analysis, BL Lac has the longest emission-frame rise time with 411 days, 0415+379 with 186 days, and 3C120 with 134 days. We note that the estimated rise times in our analysis are medians, and we use the e-folding timescale instead of a double-folding timescale. While the number of sources is too low to make inferences, this result is in agreement with our simplified assumption that the rise times can be used as a proxy for the width of the transition region.

5.1.1. Lorentz factor

Previous studies have suggested that the jet bulk Lorentz factor of AGNs may be higher for sources with longer acceleration regions (e.g. Vlahakis & Königl 2003; Beskin & Nokhrina 2006; Lyubarsky 2009; Kutkin et al. 2019). Consequently, Potter & Cotter (2015) showed that a larger width of the transition region is connected to higher Lorentz factors. A related conclusion was obtained by Malzac et al. (2018), who analysed the behaviour of a BHXRB with a shell-shock model and found that a higher average jet Lorentz factor caused shocks to be formed further downstream of the jet increasing the observed timescales. Their analysis did not consider any specific emission regions in relation to the Lorentz factor, but the conclusion of larger jet Lorentz factors connected to larger jet radii may suggest a similar effect on AGN observations. While we found an association between the Lorentz factor and the rise times, the strength of the correlation is low. This could be in part due to different scaling factors between the length of the accelerating region and the width of the transition region for different sources; Kovalev et al. (2020) also estimated the deprojected break locations from the central black hole and found 0415+379 to have a longer parabolic region compared to BL Lac. This may then indicate a more complex relationship between the length of the accelerating region and the width of the transition region.

5.1.2. Black hole mass

We found a potential weak positive correlation between the rise times and the black hole masses. Assuming a simple relationship with the length of the accelerating region and width at the jet shape transition, we may expect to observe a relationship. This is because some studies have found the location of this transition to coincide with the Bondi radius or the sphere of gravitational influence (SGI) (e.g. Asada & Nakamura 2012; Tseng et al. 2016; Kovalev et al. 2020). However, this relation is likely complicated: Okino et al. (2022) analysed the location of the jet collimation break in 3C273 and compared it to the jet shape transition regions obtained for other sources in relation to the location of the SGI. For 3C273, they found the transition in the jet shape to be downstream of the SGI, and the other sources had a large scatter in the locations of the jet shape transitions around their respective SGIs. Okino et al. (2022) note that this scatter in the locations indicates that the surrounding environment contributes to the location of the transition. Boccardi et al. (2021) suggest that an upstream location of the jet break (from the Bondi radius) could occur especially with low-luminosity sources, and Fariyanto et al. (2025) discuss how this could be connected to a lower accretion rate, which may also affect the collimation profile of the jet.

5.1.3. Accretion disk luminosity

We also found a weak significant positive correlation between the optical accretion disk luminosity (and thus un-normalised mass accretion rate) and flare rise times. In X-ray variability studies of AGNs, it has been shown that a higher X-ray luminosity corresponds to a longer variability timescale (e.g. Barr & Mushotzky 1986; Lawrence & Papadakis 1993; Green et al. 1993). Barr & Mushotzky (1986) compared the doubling times of the intensity of the X-ray emission to the 2−10 keV X-ray luminosity, finding a clear positive correlation. Similarly, Hovatta et al. (2007) found that higher radio luminosities lead to longer variability timescales. Sbarrato et al. (2014) analysed a sample of blazars and GALs and found the BLR luminosity to be closely correlated with radio luminosity. Thus, our results of the connection between accretion disk luminosity and flare rise times are not surprising, given the correlation found in Hovatta et al. (2007) and may point to a common dependency between source luminosity and variability timescales across the electromagnetic spectrum.

5.1.4. Redshift

The found correlations are mostly weak, and their significances in either direction are typically driven by a small number of sources. We examined whether the correlations between the rise time and the jet and central engine parameters could be driven by a common dependency on redshift. The results from the correlation analysis are given in Table 4. The positive correlations between redshift and emission-frame rise time, black hole mass, accretion disk luminosity, and normalised mass accretion rate are significant. Thus, we performed a partial correlation analysis using the ppcor R-package (Kim 2015) to examine whether redshift explains the correlation between the emission-frame rise time and black hole mass and accretion disk luminosity. Indeed, neither correlation is significant when considering redshift as a third variable (τ = 0.126, p = 0.337 and τ = 0.075, p = 0.541, respectively), and we can conclude that the original correlations are driven by redshift. While the Lorentz factor correlation with redshift is not significant, the margin to the significance limit is very small. If a partial correlation accounting for redshift is performed between the Lorentz factor and the rise time, the association weakens to τ = 0.156 and p = 0.110.

5.2. PSD timescales

We found that both the observer-frame flare separations and durations correlate with M_BH/ṁ, which is potentially driven by a correlation with M_BH. The association between the observer-frame rise times and black hole masses is then likely related to the correlation with the flare durations. However, if only FSRQs are considered in the analysis, a statistically significant association persists only between the timescales and M_BH/ṁ. Further analyses are thus required to understand the connection better.

The correlation with the redshift-corrected flare separations and M_BH/ṁ did not persist. This is surprising as the redshifts have a relatively narrow range of values and the resulting corrections to the timescales are moderate; a correlation analysis between M_BH/ṁ and only redshift-corrected flare durations gives a similar result as the analysis between M_BH/ṁ and observer-frame flare durations. Thus, we examined whether a common association with redshift may affect the results. Indeed, the Doppler- and redshift-corrected flare durations have a positive association with redshift whereas the redshift-corrected flare separations have a negative correlation with redshift (Table 5).

While the positive association between M_BH/ṁ and redshift is not significant (τ = 0.230, p = 0.077), it is possible that this association then affects the analysis. Using redshift as a third variable, we find that the correlation between M_BH/ṁ and the Doppler- and redshift-corrected flare durations decreases to τ = 0.204 and p = 0.120, whereas with the redshift-corrected flare separations the correlation is increased to τ = 0.222 and p = 0.090. Nevertheless, we cannot reject the null hypothesis of no correlation.

As McHardy (2009) suggested, a connection with the X-ray PSD bend of a blazar implies that a central-engine source modulates the variations in the jet. This is especially understandable with the flare separations as the rate of ejections of superluminal knots could be related to the central engine properties and matter infall (Marscher et al. 2002; Chatterjee et al. 2009, 2011). On the other hand, the association with the flare durations and the black hole mass could be partially explained by the emission region size being further out from the black hole, as explained in Sect. 5.1.2. The mass accretion rate could then also be related to the location of the emitting region (e.g. Fariyanto et al. 2025).

While we cannot verify the emission-frame timescale association, the results are intriguing as they suggest that the bend seen in the PSDs of radio light curves may have a connection to the central engine, similarly as in the X-ray domain.

5.2.1. Connection between flare durations and separations

Lister (2001) hypothesised that, assuming similarity of sources, higher Doppler-boosting factors would lead to shorter observed flare durations compared to their separations. However, we find no correlation between the Doppler factor and T_dur, obs/T_sep, obs (Fig. 8). One potential explanation is that the flare durations are longer in the emission frame for sources with higher Lorentz factors, related to the increased width of the transition region (Potter & Cotter 2015). This could then naturally counteract the otherwise increasing ratio between the flare durations and separations. Of course, this interpretation assumes that the rate of ejections is not significantly different between sources.

Fig. 8. Ratio of the mean flare durations and separations against the estimated variability Doppler factor. The source 0836+710, with a very high Doppler factor of ∼79, is excluded from the plot.

As mentioned above, the results in Table 5 suggest that the intrinsic flare durations are longer for the more distant sources, and that the intrinsic rate of ejections increases (flare separations decrease) with a higher redshift. The former association is driven by BLOs and GALs, as the association weakens and is no longer significant with FSRQs alone. The negative correlation between the intrinsic rate of ejections and redshift is for both the two source divisions (FSRQs and BLOs/GALs) and the entire sample. The division to BLOs/GALs is due to their generally smaller redshifts in the sample and because the number of sources is too low to isolate them into two groups. The increasing flare durations with redshift may also be a sample bias as we have chosen sources that have discernible flares. Lister (2001) notes that due to severe overalapping of flares in low-beamed sources, they will typically have fairly stable observed flux densities.

Because of this opposite correlation between the emission-frame flare durations and redshift and redshift-corrected flare separations and redshift, the two correlations could counteract each other and explain why we do not see longer observed flare separations with shorter observed flare durations for the more Doppler-boosted sources, which generally are at higher redshifts. On the other hand, if redshift is accounted for through partial correlation, the emission-frame flare durations and redshift-corrected flare separations also correlate weakly (τ ≈ 0.29, p ≈ 3 ⋅ 10⁻³).

5.2.2. VLBI timescales

In Paper I, we found a potential association between the PSD timescale and the longest time a 43 GHz VLBI knot was observed for, and we hypothesised that the timescale may be related to the flare durations. This association is now strengthened with the good correspondence between the flare durations and the PSD timescales. However, in general, the VLBI knot ejections have been associated with the decay times of the flares rather than their durations (e.g. Türler et al. 1999; Lindfors et al. 2006). This is because the flare is expected to peak at the radio core, only after which an associated knot can be observed to separate from the core. We do not have an answer for this apparent conflict; it is possible that these most-dominant flares are low-peaking ones (Valtaoja et al. 1992) in which case the shock becomes optically thin already before the peak amplitude. It is then related to the resolution of the VLBI network and the distance of the source, which determine when the knot is observed. On the other hand, how long the knot is visible in the VLBI image also depends on the sensitivity of the array and the ability to model it as a separate feature, and therefore it may be considered a part of the flare in the 37 GHz single-dish data for a longer time.

Another explanation is that the preliminary correspondence between the PSD and VLBI knot timescale is only at an approximately correct value: In this analysis, we did not isolate and measure the absolute longest high-amplitude flares, but included all of the brightest flares to approximate the most dominant flare durations and separations. Thus, it is possible that the VLBI knot and PSD timescales also correspond to the decay time of a single dominant bright flare and it is approximately of a similar value as the combined durations of a selection of the brightest flares.

6. Conclusions

We have analysed 54 sources to compare the decomposed flare durations and separations either to the constrained PSD timescales or to the unconstrained estimates obtained in Paper I. We found the strongest one-to-one association between the PSD timescale and the mean duration of the brightest flares, whereas the flare separations also correlate with both the PSD timescale and the flare durations. However, as only 11 of the sources in the sample have a constrained PSD timescale, we cannot make definitive statements on whether the durations or separations are more closely associated with the PSD bend.

A connection with the PSD timescale and M_BH/ṁ was found through using the flare durations and separations as proxies for the PSD timescale. While both of the decomposed timescales had a significant association with M_BH/ṁ when using observer-frame values, the association was no longer significant in the emission frame. Further analyses are thus required to understand the connection better and whether the accuracy of the estimated central engine parameters affects the results.

We also found that both the emission-frame flare durations and separations correlate with redshift, with both the flare durations and the rate of flares increasing for higher redshift sources. Curiously, if redshift is accounted for, longer intrinsic flare durations are associated with longer intrinsic flare separations (lower rate of flares). It is unclear, whether this association is then driven by differences in the sources or whether these connections are induced by, for example, a sample bias.

Using the flare rise times as a proxy for the emission region size, we found weak correlations between the emission-frame flare rise times and the jet and central engine parameters. However, the associations appear to be driven by a common dependency on redshift. This does not necessarily mean that no correlation exists, as the analysis may be limited by the sample being both relatively small and including only bright sources at higher redshifts.

Acknowledgments

The authors would like to thank Dr. Yannis Liodakis and Dr. Elina Lindfors for their valuable help when preparing the manuscript. S.K. was supported by Emil Aaltonen Foundation and Jenny and Antti Wihuri Foundation. T.H. acknowledges support from the Research Council of Finland projects 317383, 320085, 345899, and 362571 and from the European Union ERC-2024-COG – PARTICLES – 101169986. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This publication makes use of data obtained at the Metsähovi Radio Observatory, operated by the Aalto University (https://www.metsahovi.fi/opendata/). This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018).

References

Appendix A: Tables

All Tables

All Figures

Fig. 1. Example of one decomposition of the 37 GHz light curve of 0415+379 drawn from the posterior sample. The OU process has been subtracted from the decomposition. The flares coloured in blue and cyan have an amplitude higher than the mean amplitudes of all flares in the posterior, and the flares coloured in cyan have amplitudes above 50% of the maximum amplitude in the given posterior. The magenta line shows the fit of the decomposition. In the top right, the inset shows an example of the portion of the flare used to determine the flare duration (cyan) and the portion of the flare used for the rise time (magenta).

In the text

	Fig. 2. Power spectral density (PSD) timescales against the mean flare durations and separations. The coloured data points are for the constrained sources and the grey ones for the unconstrained sources. The dashed line indicates the one-to-one correspondence between the timescales.
In the text

	Fig. 3. Mean flare separations against the mean flare durations obtained from the decomposition of the light curves using the mean-amplitude limit. Source 3C84 is excluded from the plot for visual reasons, due to a much longer estimated mean flare duration and separation compared to the other sources.
In the text

	Fig. 4. Black hole mass divided by the normalised mass accretion rate (MBH/ṁ) against the observed mean flare duration (left) and separation (right).
In the text

	Fig. 5. Black hole mass against the median flare rise time in the observed frame (left) and the emission frame (right).
In the text

	Fig. 6. Accretion disk luminosity (top row) and normalised mass accretion rate (bottom row) against the observed median flare rise time (left) and emission-frame rise time (right).
In the text

	Fig. 7. Lorentz factor against the median rise time in the observed frame (left) and the emission frame (right).
In the text

	Fig. 8. Ratio of the mean flare durations and separations against the estimated variability Doppler factor. The source 0836+710, with a very high Doppler factor of ∼79, is excluded from the plot.
In the text