Optical Variability Modeling of Newly Identified Blazar Candidates behind Magellanic Clouds

The LSP (Lomb 1976; Scargle 1982; Press & Rybicki 1989; VanderPlas 2018) is a way for constructing a power spectral density (PSD) for arbitrarily spaced data (unevenly sampled time series). For an LC with N observations x_k at times t_k, it is computed as

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{LS}}(\omega ) & = & \displaystyle \frac{1}{2{\sigma }^{2}}\left[\displaystyle \frac{{\left({\displaystyle \sum }_{k=1}^{N}({x}_{k}-\bar{x})\cos [\omega ({t}_{k}-\tau )]\right)}^{2}}{{\displaystyle \sum }_{k=1}^{N}{\cos }^{2}[\omega ({t}_{k}-\tau )]}\right.\\ & & +\left.\displaystyle \frac{{\left({\displaystyle \sum }_{k=1}^{N}({x}_{k}-\bar{x})\sin [\omega ({t}_{k}-\tau )]\right)}^{2}}{{\displaystyle \sum }_{k=1}^{N}{\sin }^{2}[\omega ({t}_{k}-\tau )]}\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where ω = 2π f is the angular frequency, τ ≡ τ(ω) is defined as

$$\begin{eqnarray}&&\tau (\omega )=\displaystyle \frac{1}{2\omega }\arctan \left[\displaystyle \frac{{\sum }_{k=1}^{N}\sin (2\omega {t}_{k})}{{\sum }_{k=1}^{N}\cos (2\omega {t}_{k})}\right],\end{eqnarray} \tag{ 2 }$$

$\bar{x}$ and σ² are the sample mean and variance.

The lower limit for the sampled frequencies is ${f}_{\min }=1/({t}_{\max }-{t}_{\min })$, corresponding to the length of the time series. The upper limit, f_max, would be the Nyquist frequency, the same as in the case of the Fourier spectrum, if the data were uniformly sampled. For unevenly spaced data, the common choices for a pseudo-Nyquist frequency are somewhat arbitrary (VanderPlas 2018). The upper limit herein is set to correspond to variations on timescales of one day, i.e., the maximal frequency is f_max = 1 day⁻¹. This is motivated by the fact that the most common time step between consecutive observations is scattered around one day. The sampling during each observing night was very nonuniform, and the data uncertainties obscure any short-term variability, hence the study of intranight variability is beyond the scope of this work. The total number of sampling frequencies is set to be ${N}_{P}={n}_{0}\tfrac{{f}_{\max }}{{f}_{\min }}$, with n₀ = 10 employed hereinafter. The Poisson noise level, coming from the statistical noise due to uncertainties in the measurements, Δx_k, is given by

$$\begin{eqnarray}&&{P}_{\mathrm{Poisson}}=\displaystyle \frac{1}{2{\sigma }^{2}}\sum _{k=1}^{N}{\rm{\Delta }}{x}_{k}^{2}.\end{eqnarray} \tag{ 3 }$$

To fit a PSD model in log–log space, one needs to take into account that the evenly spaced frequencies f are no longer uniformly spaced when logarithmized, i.e., their density is greatly increased at higher values, where the Poisson noise can be expected to be significant. A straightforward least squares fitting would then rely mostly on points clustered in one region of the log f values, i.e., at high frequencies. To circumvent this problem, binning is applied. The values log f of the raw LSP are binned into equal-width bins, with at least six points in a bin, and the representative frequencies are computed as the geometric mean in each bin. The PSD value in a bin is taken as the arithmetic mean of the logarithms of the respective PSD values in that bin (Papadakis & Lawrence 1993; Isobe et al. 2015).

Fits of different models are compared using the small sample Akaike Information Criterion (AIC_c) given by

$$\begin{eqnarray}&&{\mathrm{AIC}}_{c}=2p-2{ \mathcal L }+\displaystyle \frac{2(p+1)(p+2)}{N-p-2},\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal L }$ is the log-likelihood, p is the number of parameters, and N is the number of fitted points (Akaike 1974; Hurvich & Tsai 1989; Burnham & Anderson 2004). For a regression problem,

$$\begin{eqnarray}&&{ \mathcal L }=-\displaystyle \frac{1}{2}N\mathrm{ln}\displaystyle \frac{\mathrm{RSS}}{N},\end{eqnarray} \tag{ 5 }$$

where RSS is the residual sum of squares, p is an implicit variable in the log-likelihood. A preferred model is one that minimizes AIC_c.

What is essential in assessing the relative goodness of a fit in the AIC_c method is the difference, ${{\rm{\Delta }}}_{i}={\mathrm{AIC}}_{c,i}-{\mathrm{AIC}}_{c,\min }$, between the AIC_c of the i-th model and the one with the minimal value, AIC_c,min. If Δ_i < 2, then there is substantial support for the i-th model (or the evidence against it is worth only a bare mention), and the proposition that it is an adequate description is highly probable. In other words, both models are equally good, and one cannot decide which one is better based only on the information criterion. A possible decision might be to choose the simpler model. Subsequently, if 2 < Δ_i < 4, then there is strong support for the i-th model. When 4 < Δ_i < 7, there is considerably less support, and models with Δ_i > 10 essentially exhibit no support.

The Hurst exponent H (Hurst 1951; Mandelbrot & van Ness 1968; Katsev & L'Heureux 2003; Tarnopolski 2016; Knight et al. 2017) measures the statistical self-similarity of a time series x(t). It is said that x(t) is self-similar (or self-affine) if it satisfies

$$\begin{eqnarray}&&x(t)\doteq {\lambda }^{-H}x(\lambda t),\end{eqnarray} \tag{ 6 }$$

where λ > 0 and $\doteq$ denotes equality in distribution. Self-similarity is connected with long-range dependence (memory) of a process via the autocorrelation function for lag k

$$\begin{eqnarray}&&\rho (k)=\displaystyle \frac{1}{2}\left[{\left(k+1\right)}^{2H}-2{k}^{2H}+{\left(k-1\right)}^{2H}\right].\end{eqnarray} \tag{ 7 }$$

When H > 0.5, ρ(k) decays to zero as $k\to \infty$ so slowly that its accumulated sum does not converge (i.e., $\rho (k)\propto | k{| }^{-\delta }$, 0 < δ < 1), and x(t) is then called a persistent process, i.e., the autocorrelations persist for a prolonged time. The persistency here relates to the higher probability of an increase (decrease) to be followed by another increase (decrease) in the short term, rather than alternating. Such a process is characterized by positive correlations at all lags. A process with H < 0.5 is antipersistent, also referred to as a mean-reverting series, i.e., having a tendency to quickly return to its long-term mean. Its autocorrelation, ρ(k), is summable. Both persistent and antipersistent cases can be stationary and nonstationary (we consider weak stationarity herein).

The archetypal processes with long-range dependence are the fractional Gaussian noise (fGn, a stationary process) and fractional Brownian motion (fBm, a nonstationary process with variance growing ∝t^2H with time; the increments of an fBm constitute an fGn with the same H). There is a discontinuity of H at the border between the two, where an fGn with H ≲ 1 is very similar to an fBm with H ≳ 0, as illustrated in Figure 2, and the two cases can be easily misidentified. For H = 0.5, fGn and fBm reduce to white noise and Brownian motion, respectively. The properties of H can be summarized as:

• 1.  
  0 < H < 1,
• 2.  
  H = 0.5 for an uncorrelated process,
• 3.  
  H > 0.5 for a persistent (long-term memory, correlated) process,
• 4.  
  H < 0.5 for an antipersistent (short-term memory, anticorrelated) process.

Zoom In Zoom Out Reset image size

For irregularly sampled data, the Hurst exponent can be obtained without any interpolation of the examined time series (Knight et al. 2017) with the use of the lifting wavelet transform algorithm called lifting one coefficient at a time (LOCAAT). The algorithm aims at producing a set of wavelet-like coefficients, ${\{{d}_{{j}_{r}}\}}_{r}$, whose variance obeys the relation

$$\begin{eqnarray}&&{\mathrm{log}}_{2}\mathrm{var}({d}_{{j}_{r}})=\alpha \cdot {j}^{* }+\mathrm{const}.,\end{eqnarray} \tag{ 8 }$$

where j* is an equivalent of the wavelet scale, constructed for a set of j_r coefficients. Eventually, the value of α is obtained via a linear regression of Equation (8), and is linearly related with H via $H=\tfrac{\alpha -1}{2}$ when α ∈ (1, 3), and $H=\tfrac{\alpha +1}{2}$ when α ∈ (−1, 1) (the two most common instances, related to fBm- and fGn-like signals, respectively).⁵ Note the discontinuity of H for the two ranges of α (see Figure 2). We used the package liftLRD⁶ implemented in R to estimate H. The standard errors, obtained via bootstrapping, are returned by the package as well.

The Abbe value (von Neumann 1941b, 1941a; Williams 1941; Kendall 1971; Mowlavi 2014; Tarnopolski 2016) is defined as

$$\begin{eqnarray}&&{ \mathcal A }=\displaystyle \frac{\tfrac{1}{N-1}{\sum }_{k=1}^{N-1}{\left({x}_{k+1}-{x}_{k}\right)}^{2}}{\tfrac{2}{N}{\sum }_{k=1}^{N}{\left({x}_{k}-\bar{x}\right)}^{2}}.\end{eqnarray} \tag{ 9 }$$

It quantifies the smoothness of a time series by comparing the sum of the squared differences between two successive measurements with the standard deviation of the time series. It decreases to zero for time series displaying a high degree of smoothness, while the normalization factor ensures that ${ \mathcal A }$ tends to unity for a purely noisy time series (more precisely, for a white noise process).

Three consecutive data points, ${x}_{k-1},{x}_{k},{x}_{k+1}$, can be arranged in six ways; in four of them, they will create a peak or a valley, i.e., a turning point (Kendall & Stuart 1973; Brockwell & Davis 1996). The probability of finding a turning point in such a subset is therefore 2/3, and the expected value for a random data set is ${\mu }_{T}=\tfrac{2}{3}(N-2)$—the first and last ones cannot form turning points. Let T denote the number of turning points in a time series, and ${ \mathcal T }=T/N$ be their frequency relative to the number of observations. ${ \mathcal T }$ is asymptotically equal to 2/3 for a purely random time series (Gaussian noise). A time series with ${ \mathcal T }\gt 2/3$ (i.e., with raggedness exceeding that of a white noise) will be deemed more noisy than white noise. Similarly, a time series with ${ \mathcal T }\lt 2/3$ will be less ragged than Gaussian noise. The maximal value of ${ \mathcal T }$ asymptotically approaches unity for a strictly alternating time series.

The ${A-T}$ plane (Tarnopolski 2016) was initially introduced to provide a fast and simple estimate of the Hurst exponent. It is able to differentiate between different types of colored noise, P(f) ∝ 1/f^β, characterized by different values of β. In Figure 3 we show the locations in the ${A-T}$ plane of colored noise plus Poisson noise spectra of the form P(f) ∝ 1/f^β+C, where the term C is introduced to account for the uncertainties in the measurements, manifesting through the Poisson noise level from Equation (3).

Zoom In Zoom Out Reset image size

The following models are fitted to the LSP of each object:

• 1.  
  Model A: a power law (PL) plus Poisson noise

  $$\begin{eqnarray}&&P(f)=\displaystyle \frac{{P}_{\mathrm{norm}}}{{f}^{\beta }}+C;\end{eqnarray} \tag{ 10 }$$
• 2.  
  Model B: a smoothly broken PL (SBPL) plus Poisson noise (McHardy et al. 2004; Alston et al. 2019; Alston 2019):

  $$\begin{eqnarray}&&P(f)=\displaystyle \frac{{P}_{\mathrm{norm}}{f}^{-{\beta }_{1}}}{1+{\left(\tfrac{f}{{f}_{\mathrm{break}}}\right)}^{{\beta }_{2}-{\beta }_{1}}}+C,\end{eqnarray} \tag{ 11 }$$

where the parameter C is an estimate of the Poisson noise level from Equation (3), f_break is the break frequency, and β₁, β₂ are the low- and high-frequency indices, respectively.

Parameters of the fits are gathered in Table 1, where ${T}_{\mathrm{break}}=1/{f}_{\mathrm{break}}$. The uncertainties are the standard errors of the fit, and ΔT_break comes from the law of error propagation. The best model is chosen based on AIC_c values.⁷ When Δ_i < 2, both models are equally good. Only 10 FSRQs clearly yield model A as the better one (with J0512−7105 consistent with a nearly flat PSD⁸ ); for 13 FSRQs, model B is preferred. In the case of BL Lacs, 14 are best described by model A, and only in one instance (J0538−7225) was model B pointed at.

The distributions of the exponents β from model A of all 44 objects are displayed in Figure 4. For FSRQs, the exponent β mostly lies in the range (1,2), with the mode at 1.5. BL Lacs are slightly flatter, spanning mostly the range (1, 1.8), with the mode at about 1.2; one object has a flat PSD.⁹ On the other hand, three BL Lacs have steeper PSDs, with β ∼ 3–4. Fits of the models are shown in Figures 11 and 12 in Appendix A.

Zoom In Zoom Out Reset image size

To check the reliability of the fits, we performed Monte Carlo simulations described in Appendix B. In short, one can expect the β values from model A to be spread over a wide range for input β ≲ 2, and slightly overestimated for β ≳ 2. Fitting to model B yields statistically reliable estimates, although with a nonnegligible spread and greater uncertainties. Finally, in Appendix C we tested for spurious detections of model B due to irregularly sampled LCs by imposing the same sampling as in the LCs of the examined FSRQs. We found that owing primarily to statistical fluctuations, model B should appear only nine times, half of the actual case of 18 instances of model B being a plausible result. We can therefore conclude that a subpopulation of FSRQs, with PSDs given by model B, is very likely present in our sample.

Our BL Lac–type objects are on average dimmer than FSRQs: $\bar{I}=20.09\pm 0.18\,\mathrm{mag}$ and $\bar{I}=19.00\pm 0.20\,\mathrm{mag}$, respectively (Żywucka et al. 2018). This fact may be the reason, for a majority of BL Lacs, model A is more adequate, while FSRQs require a more complex model B, i.e., the variability properties of BL Lacs could not be constrained on shorter timescales due to domination of Poisson noise at such scales. To test this, we calculated the critical frequency f₀ at which the Poisson noise has the same power as the PL component, ${P}_{\mathrm{norm}}/{f}_{0}^{\beta }=C\Rightarrow {f}_{0}={({P}_{\mathrm{norm}}/C)}^{1/\beta }$. For all FSRQs, $\mathrm{log}{f}_{0}\gtrsim -2.4$, and the corresponding β ∈ (1, 2). In case of BL Lacs, however, for the three objects with β > 2, their log f₀ < −2.8, and I ≳ 20.5 mag (see Figure 5 and Żywucka et al. 2018). This implies that their high β values are artifacts of fitting model A to PSDs that are significantly dominated by Poisson noise over a wide range of timescales, including the longest ones, where only one to two points contribute to the PL part while fitting model A; hence, they are just statistical fluctuations. Note, however, that there are still three other BL Lacs that are even dimmer (J0039−7356, J0441−6945, and J0545−6846), but appear to be well described by model A. Overall, in all instances where model B was chosen as the better one, β₂ ≳ 3 (except for J0552−6850, which yields an unusually high T_break = 1201 ± 417, and J0602−6830, whose T_break = 538 ± 579 is not constrained well due to a huge uncertainty).

Zoom In Zoom Out Reset image size

Figure 6 presents the relations between T_break and the exponents β₁ and β₂ from model B. There are 18 FSRQs and 4 BL Lacs (including J0444−6729) that yielded reasonable fits. The Pearson correlation coefficient for the β₂ − T_break relation is r = −0.49; after discarding the FSRQs with T_break > 500 days and large errors, it reduces to r = −0.01. For the β₁ − T_break relation, the respective correlation coefficients are −0.78 and −0.33. This indicates no significant correlation between the high-frequency PSD index β₂ and T_break, and moderate correlation between the low-frequency index β₁ and T_break. According to this, the relative power in the long timescale variability increases with decreasing T_break.

Zoom In Zoom Out Reset image size

The PL form of a PSD is indicative of a self-affine stochastic process, characterized by the Hurst exponent H, underlying the observed variability. The H values are listed in Table 2, and displayed graphically in Figure 7. We find that most objects have H ≤ 0.5 within errors, indicating short-term memory. Four BL Lacs (whose PSDs are best described by model A) and two FSRQs (J0512−7105 with a flat PSD, and J0512−6732, a secure blazar candidate with PSD given by model B, with an exceptionally short T_break = 67 ± 10) yield H > 0.5, implying long-term memory. Very few objects are characterized by H ≈ 0.5, so the modeled stochastic process is not necessarily uncorrelated. There is also a number of FSRQs and a few BL Lacs with H ≲ 0.2, i.e., close to the discontinuity value on the border between fGn- and fBm-like processes (see Figure 2), hence their H estimates are uncertain. In general, the autocorrelation functions drop to zero after timescales comparable to T_break, above which the system becomes decorrelated (Caplar & Tacchella 2019). This strongly suggests that models admitting long-range dependence (Tsai & Chan 2005; Tsai 2009; Feigelson et al. 2018) should be considered candidates for the underlying stochastic processes governing the observed variability of blazar LCs.

Zoom In Zoom Out Reset image size

Table 2.  Hurst Exponents and $({ \mathcal A },{ \mathcal T })$ Locations of the Newly Identified FSRQ- and BL Lac–Type Blazar Candidates

Number	Object	H	${ \mathcal A }$	${ \mathcal T }$
(1)	(2)	(3)	(4)	(5)
FSRQ-Type Blazar Candidates
1	J0054−7248	0.42 ± 0.02	0.83 ± 0.02	0.670 ± 0.011
2	J0114−7320a	0.06 ± 0.04	0.033 ± 0.002	0.653 ± 0.015
3	J0120−7334a	0.04 ± 0.03	0.027 ± 0.002	0.653 ± 0.012
4	J0122−7152	0.24 ± 0.04	0.24 ± 0.02	0.643 ± 0.012
5	J0442−6818a	0.04 ± 0.03	0.022 ± 0.002	0.664 ± 0.014
6	J0445−6859	0.31 ± 0.05	0.59 ± 0.03	0.667 ± 0.017
7	J0446−6758	0.45 ± 0.05	0.30 ± 0.02	0.652 ± 0.011
8	J0455−6933	0.25 ± 0.05	0.22 ± 0.02	0.671 ± 0.019
9	J0459−6756	0.21 ± 0.03	0.17 ± 0.01	0.668 ± 0.012
10	J0510−6941	0.04 ± 0.03	0.20 ± 0.01	0.637 ± 0.013
11	J0512−7105	0.63 ± 0.05	0.90 ± 0.05	0.717 ± 0.019
12	J0512−6732a	0.85 ± 0.03	0.26 ± 0.02	0.640 ± 0.014
13	J0515−6756	0.38 ± 0.02	0.82 ± 0.03	0.670 ± 0.012
14	J0517−6759	0.29 ± 0.04	0.53 ± 0.03	0.663 ± 0.013
15	J0527−7036	0.04 ± 0.03	0.07 ± 0.01	0.655 ± 0.013
16	J0528−6836	0.26 ± 0.05	0.19 ± 0.01	0.651 ± 0.011
17	J0532−6931	0.07 ± 0.03	0.013 ± 0.001	0.626 ± 0.009
18	J0535−7037	0.42 ± 0.03	0.89 ± 0.03	0.630 ± 0.011
19	J0541−6800	0.08 ± 0.05	0.30 ± 0.02	0.667 ± 0.012
20	J0541−6815	0.21 ± 0.03	0.27 ± 0.02	0.643 ± 0.012
21	J0547−7207	0.03 ± 0.02	0.09 ± 0.01	0.655 ± 0.014
22	J0551−6916a	0.06 ± 0.04	0.046 ± 0.004	0.607 ± 0.016
23	J0551−6843a	0.11 ± 0.04	0.065 ± 0.005	0.623 ± 0.014
24	J0552−6850	0.22 ± 0.05	0.25 ± 0.02	0.706 ± 0.013
25	J0557−6944	0.49 ± 0.05	0.26 ± 0.03	0.684 ± 0.014
26	J0559−6920	0.22 ± 0.03	0.03 ± 0.02	0.647 ± 0.014
27	J0602−6830	0.07 ± 0.04	0.10 ± 0.01	0.634 ± 0.014
BL Lac–Type Blazar Candidates
1	J0039−7356	0.44 ± 0.03	0.96 ± 0.02	0.660 ± 0.010
2	J0111−7302a	0.35 ± 0.04	0.73 ± 0.03	0.680 ± 0.012
3	J0123−7236	⋯	0.95 ± 0.03	0.681 ± 0.012
4	J0439−6832	0.60 ± 0.06	0.83 ± 0.03	0.658 ± 0.014
5	J0441−6945	0.45 ± 0.03	0.80 ± 0.04	0.660 ± 0.014
6	J0444−6729	0.21 ± 0.05	0.51 ± 0.04	0.641 ± 0.019
7	J0446−6718	0.58 ± 0.05	0.94 ± 0.03	0.655 ± 0.015
8	J0453−6949	0.48 ± 0.04	0.93 ± 0.03	0.672 ± 0.012
9	J0457−6920	0.26 ± 0.03	0.66 ± 0.03	0.651 ± 0.012
10	J0501−6653a	0.29 ± 0.04	0.39 ± 0.02	0.643 ± 0.014
11	J0516−6803	0.62 ± 0.05	0.99 ± 0.03	0.670 ± 0.014
12	J0518−6755a	0.18 ± 0.03	0.26 ± 0.02	0.665 ± 0.015
13	J0521−6959	0.23 ± 0.04	0.48 ± 0.03	0.659 ± 0.016
14	J0522−7135	0.39 ± 0.04	0.88 ± 0.03	0.657 ± 0.014
15	J0538−7225	0.28 ± 0.04	0.32 ± 0.02	0.659 ± 0.014
16	J0545−6846	0.58 ± 0.05	0.87 ± 0.04	0.658 ± 0.019
17	J0553−6845	0.36 ± 0.04	0.58 ± 0.02	0.634 ± 0.013

Notes. Columns: (1) number of the source; (2) source designation; (3) Hurst exponent; (4) Abbe value; (5) ratio of turning points.

^aSources considered secure blazar candidates by Żywucka et al. (2018).

Download table as:  ASCIITypeset image

We note an interesting correlation between H and the bolometric luminosities for FSRQ-type candidates (see Table 3 and Section 5.4), displayed in Figure 8. The Pearson coefficient is r = −0.41; after discarding J0512−6732, a bright outlier with H = 0.85 ± 0.03, we obtain r = −0.70. While potentially this could link the persistence properties of an LC (via H) and physical processes governing the radiative output (via L_bol), we propose a more straightforward explanation: in our FSRQ sample, it turns out the dimmer the object, the lower the L_bol (r = −0.89), and so, as we argued in Section 4.1, the more the LC is consistent with white noise (Poisson noise level dominates the PSD). Indeed, most objects with PSDs given by model A lie roughly around H ≈ 0.5. On the other hand, J0512−6732 does not follow this scheme, because it is one of the brightest FSRQ candidates in our sample and yields a remarkably high value of H. Recall that this source's PSD is better described by model B, with an exceptionally short T_break = 67 ± 10 days. Finally, cases with H ≲ 0.2 are dubious due to the discontinuity of H (see Figure 2), so they might as well yield high values of H. It is therefore unclear whether this one outlier is a statistical fluctuation, or a hint of a subpopulation of bright, long-term memory blazars.

Zoom In Zoom Out Reset image size

Table 3.  Bolometric Luminosities and BH Mass Estimates of FSRQ-Type Blazar Candidates

Number	Object	$\mathrm{log}{M}_{\mathrm{BH}}$ a	$\mathrm{log}{L}_{\mathrm{bol}}$	$\mathrm{log}{M}_{\mathrm{BH}}$ b
		(M⊙)	($\mathrm{erg}\,{{\rm{s}}}^{-1}$)	$({M}_{\odot })$
(1)	(2)	(3)	(4)	(5)
FSRQ-Type Blazar Candidates
1	J0054−7248	⋯	44.54 ± 0.13	⋯
2	J0114−7320c	(9.32,10.45)	46.15 ± 0.05	9.31 ± 0.31
3	J0120−7334c	(9.06,10.14)	46.80 ± 0.07d	9.53 ± 0.34
4	J0122−7152	(8.97,10.11)	45.52 ± 0.07d	8.86 ± 0.26
5	J0442−6818c	(9.27,10.24)	46.19 ± 0.05	9.26 ± 0.30
6	J0445−6859	⋯	46.37 ± 0.03	⋯
7	J0446−6758	⋯	45.75 ± 0.07	⋯
8	J0455−6933	(9.20,10.19)	46.24 ± 0.04	9.29 ± 0.30
9	J0459−6756	(9.01,10.18)	46.19 ± 0.07d	9.27 ± 0.31
10	J0510−6941	(9.22,10.16)	46.21 ± 0.03	9.25 ± 0.30
11	J0512−7105	⋯	44.05 ± 0.12	⋯
12	J0512−6732c a	(8.49,9.40)	46.90 ± 0.03	9.33 ± 0.33
13	J0515−6756	⋯	44.58 ± 0.06	⋯
14	J0517−6759	(9.16,10.25)	44.50 ± 0.14	8.39 ± 0.22
15	J0527−7036	(8.18,9.73)	45.90 ± 0.04	8.79 ± 0.30
16	J0528−6836	⋯	47.15 ± 0.02	⋯
17	J0532−6931	(8.76,9.90)	47.15 ± 0.02	9.56 ± 0.36
18	J0535−7037	⋯	45.06 ± 0.09	⋯
19	J0541−6800	(9.16,10.47)	45.59 ± 0.07d	9.04 ± 0.29
20	J0541−6815	(9.03,9.98)	46.43 ± 0.03	9.31 ± 0.31
21	J0547−7207	(9.28,10.40)	45.76 ± 0.07d	9.09 ± 0.28
22	J0551−6916c	(9.01,10.00)	46.95 ± 0.03	9.60 ± 0.35
23	J0551−6843c	⋯	46.52 ± 0.04	⋯
24	J0552−6850	(9.74,10.84)	46.19 ± 0.04	9.59 ± 0.32
25	J0557−6944	⋯	44.73 ± 0.10	⋯
26	J0559−6920	(9.02,10.15)	46.15 ± 0.07	9.24 ± 0.31
27	J0602−6830	< 10.79	46.23 ± 0.03	9.44 ± 0.38

Notes.

^aAssuming α = 0.1; the lower limit is for T_break − ΔT_break and a_⋆ = 0; the upper limit is for T_break+ΔT_break and a_⋆ = 0.998. See text for details. ^bEstimates from the fundamental plane of AGN variability. Columns: (1) number of the source; (2) source designation; (3) range of BH mass based on Equation (12); (4) bolometric luminosity; (5) mass of BH based on Equation (15). ^cSources considered secure blazar candidates by Żywucka et al. (2018). ^dMissing uncertainty of I magnitude; the error of log L_bol is estimated as the mean error of other objects.

Download table as:  ASCIITypeset image

The locations of FSRQs and BL Lacs in the ${ \mathcal A }\mbox{--}{ \mathcal T }$ plane are gathered in Table 2, and displayed in Figure 9. The uncertainties are computed by bootstrapping. Most objects fall in the region occupied by PL plus Poisson noise processes, with various Poisson noise levels. Three FSRQs are interestingly placed: J0535−7037 is marginally below the pure PL line, while J0512−7105 and J0552−6850 are above the limiting ${ \mathcal T }=2/3$ line. The LSP implies that J0535−7037 has a PL PSD with β ≈ 1, i.e., pink noise. J0512−7105 has a flat PSD, but yields H > 0.5. The PSD of J0552−6850, however, shows a clear flattening on timescales greater than 1200 days, and has H < 0.5. They are also distant in the ${ \mathcal A }\mbox{--}{ \mathcal T }$ plane, with ${ \mathcal A }=0.90$ and ${ \mathcal A }=0.25$, respectively.

Zoom In Zoom Out Reset image size

One thing to bear in mind is that both ${ \mathcal A }$ and ${ \mathcal T }$ are order statistics, i.e., they are insensitive to the spacing between consecutive data points. One way of justifying the usage of the ${ \mathcal A }\mbox{--}{ \mathcal T }$ plane in the case of irregularly sampled time series is by noting that any LC comes from sampling a continuous process; hence, a continuum of data between any two observations is missing. In spite of this fact, it can be generally expected that the true characteristics of an analyzed system are properly captured by the observations, and the interrelation between available data points catches the overall behavior of the system.

Among our 44 blazar candidates, 41 are located in the region of the ${ \mathcal A }\mbox{--}{ \mathcal T }$ plane occupied by PL plus Poisson noise processes; one FSRQ is located marginally below (less noisy than white noise), and two FSRQ candidates are above the line ${ \mathcal T }=2/3$ (more noisy than white noise). While not of direct interpretation herein, these two objects clearly stand out in this context. Moreover, BL Lac candidates are clearly characterized by higher ${ \mathcal A }$ values than FSRQ ones (means of 0.71 ± 0.06 and 0.29 ± 0.05, respectively). This, as discussed in Sections 4.1 and 4.2, is consistent with dimmer objects being more noise polluted. Recent developments of the ${ \mathcal A }\mbox{--}{ \mathcal T }$ methodology (Zunino et al. 2017; Zhao & Morales 2018) prove it to be a useful tool in time series analysis.

The optical emission of blazars is expected to be predominantly due to the synchrotron radiation of highly relativistic electrons and positrons accelerated in situ within the blazar emission zone. In the majority of theoretical models and scenarios proposed to account for the PL shape of blazars' PSDs, it is assumed that the energy density of the synchrotron-emitting electrons fluctuates about a mean value, with the distribution, which is also a PL in the frequency domain (e.g., Marscher 2014). Such fluctuations could be related to the dominant electron acceleration process at work, e.g., fluctuations in the bulk Lorentz factor at the base of a jet, which determine properties of the internal shocks developing further along the outflow and dissipating the jet bulk kinetic energy into the internal energy of the jet electrons (see, e.g., Malzac 2013, 2014). That is, one assumes a given PL form of the variability power spectrum of electron fluctuations, to match the observed PSD in a given range of the electromagnetic spectrum. As discussed by Finke & Becker (2014, 2015), for the electron variability power spectrum ∼1/f^β, the index of the corresponding PSD of the synchrotron emission is also β on timescales longer than the characteristic cooling timescale of the emitting electrons, but β+2 on the timescales shorter than that.

This situation is analogous to the case of disk-dominated systems, where perturbations in the local disk parameters (e.g., magnetic field), leading to the enhanced energy dissipation, shape the observed PSD of the thermal disk emission. In particular, Kelly et al. (2009, 2011) discussed how uncorrelated (Gaussian) fluctuations within the disk may result in the observed PSD of the ∼1/f⁰ form on timescales longer than the characteristic relaxation (thermal) timescale in the system, and ∼1/f² on the timescales shorter than this.

Based on the above reasoning and scenarios, for the sources most likely dominated by the disk emission, i.e., FSRQs and FSRQ candidates from our sample, we propose a possible interpretation for the breaks and high β₂ in the obtained PSDs by connecting them to the dynamical change of the accreting matters orbits at the inner edge of the disk, related with the thermal timescale t_th. We explore the implications of such association in the subsequent sections.

Chatterjee et al. (2008) analyzed an optical (R filter) LC of the FSRQ 3C 279, spanning 11 yr. They found its PSD to be well described by a PL with β = 1.7. Four Kepler AGNs observed over 2–4 quarters exhibited steep PSDs, with β ≈ 3, within timescales of 1–100 days (Mushotzky et al. 2011). On the other hand, four other radio-loud AGNs, including three FSRQs, observed over 8–11 quarters, were characterized by β ∈ (1, 2) (Wehrle et al. 2013; Revalski et al. 2014). Simm et al. (2016) examined ≈90 AGNs from the Pan-STARRS1 XMM–COSMOS survey, and obtained in most cases good fits with a broken PL, with break timescales of 100–300 days, a low-frequency PL index β₁ ∈ (0, 2), and a high-frequency index β₂ ∈ (2, 4). Recently, Aranzana et al. (2018) examined short-term variability of 252 Kepler AGNs, resulting in β ∈ (1, 3.5), with a mode at 2.4, although the displayed PSDs clearly exhibit flattening at timescales ≳10^5–5.5 s. Smith et al. (2018) examined 21 Kepler AGNs, spanning 3–14 quarters, and in six cases found breaks within ∼10–50 days in their PSDs. In particular, a mild correlation between M_BH and T_break was observed, contrary to Simm et al. (2016). Finally, Caplar & Tacchella (2019) investigated the PSDs of ≈2200 AGNs from the Palomar Transient Factory survey, and obtained β ∈ (1.5, 4).

Kastendieck et al. (2011) used three methods to investigate the long-term LC of a BL Lac object, PKS 2155−304, spanning the years 1934–2010, i.e., LSP, the SF, and the multiple fragments variance function (MFVF). Interestingly, they found $\beta ={5.0}_{-3.0}^{+1.7}$ using LSP, $\beta ={1.6}_{-0.2}^{+0.4}$ with SF, and $\beta ={1.8}_{-0.2}^{+0.1}$ with MFVF. All three methods give comparable results within errors with a break to an assumed white noise at timescales ≳1000 days. For six blazars (five FSRQs and one BL Lac), Chatterjee et al. (2012) obtained β ≈ 2, with an exception of the FSRQ PKS 1510−089, which yielded a much flatter PSD with β = 0.6. An R filter LC of the BL Lac PKS 0735+178 exhibits a purely red noise PSD, i.e., with β = 2 (Goyal et al. 2017). Similarly, R filter LCs of 29 BL Lacs and two FSRQs, spanning ≈10 yr, were analyzed by Nilsson et al. (2018), who obtained β ∈ (1, 2).

Overall, it appears that AGNs either are characterized by PSDs in the form of model A, with β ∈ (1, 3), occasionally with a flatter PSD, or exhibit breaks (model B) on timescales of 100–1000 days, with a steeper PL component at high frequencies. The blazar candidates examined herein fall into those two categories: objects well described by model A yield β ∈ (1, 2), while those better characterized by model B exhibit break timescales within roughly 100–400 days, low-frequency PL index β₁ ≲ 1, and a steep PL component at higher frequencies, β₂ ∈ (3, 7), reaching as high as nine (see Table 1). Such steep PSDs effectively imply no variability on the associated timescales, because the power drops drastically from the conventional PL at lower frequencies to the Poisson noise level. This means there is a sharp cutoff at T_break below which variability on short timescales is wiped out (excluding the region of Poisson noise domination). Therefore, these faint, distant sources might constitute a different, peculiar class of blazars.

In sources for which the observed optical emission is dominated by the radiative output of accretion disks rather than jets, as is in fact expected for FSRQ-type objects, the break timescale T_B can, in principle, be connected to the BH mass and spin. Such a break might appear when the matter inspiraling in the disk transitions from bound orbits to a free fall occurring for distances smaller than the inner radius of the disk, and can explain the high β₂. Therefore, by assuming that the disk extends sharply to the innermost stable circular orbit (ISCO) located at radius r = r_ISCO (Mohan & Mangalam 2014):

$$\begin{eqnarray}\begin{array}{rcl}{T}_{B} & = & 2\pi ({r}_{\mathrm{ISCO}}^{3/2}+{a}_{\star })(1+z)\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{{c}^{3}}\\ & = & 0.359\,{m}_{9}(1+z)f({a}_{\star })\mathrm{day},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${m}_{9}={M}_{\mathrm{BH}}/({10}^{9}{M}_{\odot })$, ${a}_{\star }={Jc}/{{GM}}_{\mathrm{BH}}^{2}$ is the dimensionless spin, and J the angular momentum of the BH. For a prograde rotation, $f({a}_{\star })$ is given by Bardeen et al. (1972):

$$\begin{eqnarray}&&f({a}_{\star })={r}_{\mathrm{ISCO}}^{3/2}+{a}_{\star },\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&{r}_{\mathrm{ISCO}}=3+{Z}_{2}-\sqrt{\left(3-{Z}_{1}\right)\left(3+{Z}_{1}+2{Z}_{2}\right)},\end{eqnarray} \tag{ 13b }$$

$$\begin{eqnarray}&&{Z}_{1}=1+{\left(1-{a}_{\star }^{2}\right)}^{1/3}\left[{\left(1+{a}_{\star }\right)}^{1/3}+{\left(1-{a}_{\star }\right)}^{1/3}\right],\end{eqnarray} \tag{ 13c }$$

$$\begin{eqnarray}&&{Z}_{2}=\sqrt{3{a}_{\star }^{2}+{Z}_{1}^{2}}.\end{eqnarray} \tag{ 13d }$$

We consider an accretion disk with a viscosity parameter α (Novikov & Thorne 1973; Shakura & Sunyaev 1973; Page & Thorne 1974); then by associating the break timescale T_break from model B with the thermal timescale t_th, we get (Czerny 2006; Lasota 2016):

$$\begin{eqnarray}&&{T}_{\mathrm{break}}\sim {t}_{\mathrm{th}}\sim {\alpha }^{-1}{t}_{{\rm{K}}},\end{eqnarray} \tag{ 14 }$$

where ${t}_{K}=2\pi /{\rm{\Omega }}$ is the orbital periodicity of a Keplerian motion on a circular orbit with radius r, with angular frequency ${\rm{\Omega }}={({r}^{3/2}+{a}_{\star })}^{-1}$ around the BH. Therefore, in Equation (12), ${T}_{B}\sim {t}_{{\rm{K}}}\sim \alpha {T}_{\mathrm{break}}$.

With the values and uncertainties of T_break from Table 1 and redshifts from Żywucka et al. (2018), we obtain a set of BH masses and spins satisfying Equation (12). In Table 3, the 68% confidence intervals¹⁰ for log M_BH are given assuming α = 0.1 (Liu et al. 2008). The lower limits of the intervals are obtained for T_break − ΔT_break and a_⋆ = 0; the upper limits are obtained for T_break+ΔT_break and a_⋆ = 0.998, the maximal spin of an accreting BH (Thorne 1974). The value of α is accurate to a multiplicative factor of ∼3 (Grzędzielski et al. 2017), which can change log M_BH by ±log 3 ≈ ±0.48. Moreover, the emission need not necessarily come from the ISCO, as, e.g., in a truncated disk model. By assuming, e.g., r = 2r_ISCO, the log M_BH estimates are lowered by 0.45 for a_⋆ = 0, and by 0.3 for a_⋆ = 0.998. For r = 3r_ISCO, the respective factors are 0.7 and 0.5. Hints at the existence of truncated disks around supermassive BHs (SMBHs) were obtained for the radio galaxies: 3C 120 with a relatively low log M_BH = 7.74 (Cowperthwaite & Reynolds 2012; Lohfink et al. 2013), and 4C+74.26 with log M_BH = 9.6 (Gofford et al. 2015; Bhatta et al. 2018), suggesting that the mass estimates from Table 3 may be lowered by this account (but can be simultaneously increased if α > 0.1). On the other hand, our estimates are consistent with masses of bright Fermi blazars that are within 8 ≲ log M_BH ≲ 10 (Ghisellini et al. 2010b), with FSRQs on average more massive than BL Lacs. A more recent sample of bright FSRQs only has a similar log M_BH distribution (Castignani et al. 2013). Moreover, some BH masses of FSRQs are known to attain high values, even up to log M_BH = 10.6 (Ghisellini et al. 2010a). Therefore, our upper limits on the BH masses seem reasonable, and are consistent with the upper limits derived theoretically by Inayoshi & Haiman (2016) and King (2016). They could be further constrained with dedicated observations.

Most SMBHs inhabiting radio-quiet galaxies appear to have spins a_⋆ ≳ 0.5 (McClintock et al. 2011; Reynolds 2013, 2014); hence, the masses are inclined to lie in the upper half of the presented intervals. Quasars are expected, on average, to exhibit accretion efficiency η > 0.1 (Soltan 1982), corresponding to a_⋆ > 0.67 (Sądowski 2011). Elvis et al. (2002) argued that SMBHs should yield η > 0.15, i.e., a_⋆ > 0.88. Even though it is not clear what spins' range should be expected for radio-loud galaxies, blazars in particular, it seems reasonable to expect the rotation rates to be high.

In Figure 10, the effect of viscosity is presented for α = 0.03, 0.1, 0.3, for some representative values of z and T_break. Overall, if the viscosity is not constrained tightly, log M_BH is uncertain to an additive factor ≳2. On the other hand, if the BH mass and spin can be obtained with other methods, e.g., continuum-fitting, Fe K_α line, or X-ray reflection for the spin (McClintock et al. 2011; Middleton 2016; Kammoun et al. 2018), and reverberation mapping (Peterson 2014) for the mass, Equation (12) can be used to estimate the viscosity parameter α. In Figure 10, such a hypothetical scenario is presented for log M_BH and a_⋆ with errors within the gray bands, with the black dot denoting the exact values, fixed in the simulation. It can be found that for α = 0.05, one gets agreement between all relevant parameters.

Zoom In Zoom Out Reset image size

Possibilities of retrograde rotation in AGNs were considered before (Garofalo et al. 2010). We find, however, that then the dependence of log M_BH on a_⋆ is weakly negative, obviously coinciding for a_⋆ = 0 with predictions from the prograde scenario. Whether the rotation is prograde or retrograde is another factor to take into account, although known BH spins suggest prograde rotation is more common.

McHardy et al. (2006) discovered a relation between the break timescale, BH mass, and bolometric luminosity to be:

$$\begin{eqnarray}&&\mathrm{log}\displaystyle \frac{{T}_{\mathrm{break}}}{1\,\mathrm{day}}=A\mathrm{log}\displaystyle \frac{{M}_{\mathrm{BH}}}{{10}^{6}{M}_{\odot }}-B\mathrm{log}\displaystyle \frac{{L}_{\mathrm{bol}}}{{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}+C,\end{eqnarray} \tag{ 15 }$$

where A = 2.1 ± 0.15, B = 0.98 ± 0.15, C = −2.32 ± 0.2.

Marshall et al. (2009) subsequently compared T_break from Equation (15) with the estimate from the PSD, and obtained good agreement. However, Equation (15) was obtained using X-ray data, hence the T_break therein refers to the X-ray PSD. Indeed, Carini & Ryle (2012) found that there is a discrepancy when the T_break is derived from the optical PSD. It is, however, still unclear whether ultimately the optical and X-ray characteristic timescales are the same or not (Smith et al. 2018).

Nevertheless, we employed Equation (15) to estimate the BH masses. The bolometric luminosities were calculated according to Kozłowski (2015),¹¹ and are gathered in Table 3. For computing L_bol, the latest cosmological parameters within a flat ΛCDM model (Planck Collaboration et al. 2018) are employed: ${H}_{0}\,=67.4\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, Ω_m = 0.315, Ω_Λ = 0.685. The M_BH estimates are mostly consistent with values obtained in Section 5.3, except for J0517−6759, for which we obtain a discrepancy of ≳1 order of magnitude. This can imply, assuming that Equation (15) gives a proper mass, that the accretion disk is truncated, and/or is characterized by a low viscosity parameter.