Multiwavelength Variability Power Spectrum Analysis of the Blazars 3C 279 and PKS 1510–089 on Multiple Timescales

We analyzed the Fermi-LAT (Atwood et al. 2009, 2013; Bruel et al. 2018) data for 3C 279 and PKS 1510−089 for the 0.1–300 GeV band, for the period 2008 August 8–2020 September 16 with an integration time of seven days using the Fermi Science Tools package Fermitools-conda (version 2.0.18) available from Fermi Science Support Center ²¹ (FSSC) and the latest instrument response function p8r3_source_v3. The standard unbinned likelihood analysis was performed. ²² The procedure started with selecting the events for a time interval (1 week) with standard quality cuts (DATA_QUAL>0)&&(LAT_CONFIG==1) and zenith angle cut of θ < 90° to avoid contamination from Earth's limb using the tasks gtselect and gtmktime. For the target, we selected a 20° × 20° region of interest (ROI) centered at the source catalog position. Next, the exposure map was created from the livetime cube (gtltcube) using the task gtexpmap with a total of 35 logarithmically spaced energy intervals bounded within 0.1–300 GeV. The choice of energy resolution is governed by the FSSC recommendation, as 10 bins per decade was found to be sufficient to accommodate the change in effective area below ∼1 GeV. Finally, in the likelihood analysis (gtlike), all the sources from data release 2 (DR2) of the 10 yr Fermi-LAT catalog (4FGL-DR2; Abdollahi et al. 2020; Ballet et al. 2020) within a 30° radius of the target were included, as were the standard templates for diffuse emission from our Galaxy (gll_iem_v07.fits) and the isotropic γ-ray background (iso_p8r3_source_V3_v1.txt). In the fitting (gtlike), the spectral normalization of the target and sources within 5° of the ROI center were allowed to be free while for the remaining sources the spectral shapes were fixed at their catalog-derived values. The normalization of the Galactic diffuse background was kept frozen while that of the isotropic γ-ray background component was allowed to vary. Specifically, the γ-ray spectra of the targets were modeled with a log-parabolic function given by the 4FGL-DR2 catalog. The light curves were generated including 0.1−300 GeV fluxes for which the test statistic (TS, defined as twice the difference between the log-likelihood, which is maximized by adjusting all the parameters of the model computed with and without including the source; Mattox et al. 1996) was found to be >10. This gives us preliminary light curves for 3C 279 and PKS 1510−089.

Next, we checked for the contamination of the targets' fluxes due to proximity to, or occultation by, the Sun and the Moon as it has been demonstrated that the Sun and the Moon are sources of γ-ray emission (Abdo et al. 2011c, 2012). We extracted the equatorial positions (R.A. and decl.) of the Sun at mission elapsed time (MET) from the spacecraft files and computed the angular distances from the targets in the sky plane. The Sun is known to occult 3C 279 on October 8 each year (Barbiellini et al. 2014). The closest distance between the Sun and PKS 1510−089 is ≈85, occurring on November 11. Similarly, the positions (R.A. and decl.) of the Moon were computed at the METs given in the spacecraft files in the Geocentric Celestial Reference System (GCRS) using the get_moon task from the astropy.coordinates library ²³ and the angular distances between the targets and Moon were computed. The Moon occults 3C 279 while the closest distance between the Moon and PKS 1510−089 is ∼32. For the blazar 3C 279, we extracted all the weekly time bins for which the distances between the targets and the Sun/Moon were found to be smaller than 10° and analyzed them separately. As the Moon moves much faster than the Sun in the sky plane (∼13° day⁻¹ and ∼1° day⁻¹, respectively) and therefore is expected to occult 3C 279 more often than the Sun, we proceeded to account for the contamination due to Sun/Moon to the target fluxes in two steps.

First, we extracted all the weekly time bins for which the angular distances between the Moon and 3C 279 were found to be <10° without having the Sun nearby (i.e., the angular distance between the Sun and 3C 279 was >20°). Of these, we randomly selected a few time bins during which the Moon was passing the target's ROI (separation <10°). We generated the lunar templates for these bins following the FSSC guidance. ²⁴ The description of the Solar System Tools, which model the quiescent emission from the Sun and Moon, is given in Johannesson & Orlando (2013). An all-sky livetime cube was generated for the lunar emission using the task gtltcube2 with 35 logarithmically spaced energy intervals with livetime cube and counts map option from the output of the task gtbin. Next, we generated the livetime cube for the Moon using the task gtltcubesun with body=Moon and instrument angle=1°. As a final step, the template for the Moon emission was generated using the task gtsuntemp using the lunar profile template, ²⁵ and the all-sky livetime cube generated from the task gtltcube2. Finally, the Moon template was included in the model file for the likelihood analysis. The Moon's SED was modeled as a power law for which the normalization was allowed to vary while the index was fixed during the fitting. We compared the target fluxes with and without including the Moon template in the likelihood analysis and found that the change in flux (before Moon—after Moon/before Moon) was comparable to the measurement uncertainties in the weekly time bins with TS > 10. Thus, we conclude that the Moon does not have a significant impact on the weekly binned fluxes of 3C 279.

Next, we proceeded to compute the contamination due to the Sun's proximity to 3C 279. We generated the solar template following the same steps needed to generate the lunar template with appropriate changes in input parameters (such as including body = Sun and instrument angle = 180° in the task gtltcubesun) and the likelihood analysis was performed. The normalization of the power-law SED of the Sun was allowed to vary while the index was fixed during the fitting. We compared the change in flux (before Sun – after Sun/before Sun) and found that for a few time bins, the change in target flux was up to 50% of the uncorrected value when the solar template was included in the model fitting. In two time bins, 3C 279 was detected with TS < 10 (as compared to TS > 10 without including the solar template). In general though we found that the change in flux is lower (<10%) and therefore comparable to measurement uncertainties when the Sun was >2° away from the target. In a few time bins analyzed, the target and the Sun's fluxes were found to be highly correlated (correlation coefficient >0.5), so the flux measurements in these time bins were discarded. The final light curve of the blazar 3C 279 (Figure 1) includes all flux measurements corrected for solar emission and TS > 10.

Zoom In Zoom Out Reset image size

Because the analysis performed for weekly integrated time bins for 3C 279 indicated negligible change when the Sun is >2° away from the target, we did not perform the checks for the Sun's contamination for PKS 1510−089 (closest distance>8°). We also did not correct for the Moon's contamination as weekly integrated fluxes did not show a significant change in flux after the inclusion of the lunar template. Figure 2 shows the final light curve including measurements with TS > 10.

Zoom In Zoom Out Reset image size

In order to extend the gamma-ray PSD temporal range down to hourly timescales we needed to examine periods of high fluxes. So we generated short-duration (3 weeks) LAT light curves using a 3 hr integration binning interval for time ranges when the daily LAT flux from the sources exceeded 1 × 10⁻⁵ ph cm⁻² s⁻¹. We selected three time ranges for the blazar 3C 279, denoted as "Flare 1" (Cutini 2015; Ackermann et al. 2016), "Flare 2" (van Zyl & Ojha 2018), and "Flare 3" (Ojha & van Zyl 2018), respectively. For the blazar PKS 1510−089, we selected one time range, denoted as "Flare 1" (Hungwe et al. 2011; Orienti et al. 2013). In each case, the selected time range was centered around the day of the flare. The analysis was performed in a similar fashion as that of 1 week integration time bins. The final light curves were checked for the goodness of gtlike results by plotting Flux/ΔFlux versus $\sqrt{{TS}}$. For a good fit, one expects a linear trend between these quantities. We removed one flux point from the "Flare 3" and "Flare 1" light curves of 3C 279 and PKS 1510−089, respectively, which did not follow a linear trend. For 3C 279, we checked for the contamination due to the Moon alone as the selected time ranges did not coincide with the passage of the Sun. We note that the change in flux due to inclusion of lunar templates was <30%, i.e., within the measurement uncertainties provided by these measurements. Therefore, we conclude that the occultation of the Moon does not have a significant effect on 3 hr integration bin flux measurements. For PKS 1510−089, the selected time range did not coincide with the passage of the Moon in the ROI (>10°). Therefore, we do not correct for the Moon's effect in the 3 hr integration bin light curves for either of the targets. Figure 3 gives the 3 hr integration bin light curves for "Flare 1", "Flare 2", and "Flare 3" for the blazar 3C 279 (panels (a), (b), and (c)) and "Flare 1" for the blazar PKS 1510−089 (panel (d)).

Zoom In Zoom Out Reset image size

We downloaded the 1998 and 1334 Stdprods (Standard2 data product) for 3C 279 and PKS 1510−089, respectively, from the RXTE quick look analysis from the HEASARC website. ²⁶ These source and background Standard2 spectra were reduced by a pipeline using a standard set of FTOOLS and PCA backgrounds with a general script and employed one single set of filtering criteria for all observations. ²⁷ The Standard2 mode of PCA has 129 energy channels, covering the entire range 2–25 keV of the PCA detector. The spectra were fitted within the energy range 3.0–20.0 keV with a power law of the form $F(E)={F}_{0}{E}^{-{{\rm{\Gamma }}}_{X}}$, where F₀ is the normalization and Γ_X is the photon index, found via xspec with χ² statistics. In the spectral analyses, the Galactic absorption corresponding to neutral hydrogen column density in the direction of the source was fixed (Kalberla et al. 2005). The unabsorbed photon fluxes were obtained within the 3–20 keV energy range by fitting wabs*cpflux*powerlaw. We averaged the fluxes using a binning interval of 1 day. Figures 1(b) and 2(b) show these X-ray photon fluxes for 3C 279 and PKS 1510−089, respectively. A few earlier studies report part of the RXTE-PCA data sets examined by us for 3C 279 (Hartman et al. 2001; Larionov et al. 2008; Böttcher et al. 2007; Chatterjee et al. 2008; Abdo et al. 2010b) and PKS 1510−089 (Marscher et al. 2010; Castignani et al. 2017).

The archival data from the Swift-XRT (Gehrels et al. 2004), consisting of 415 and 262 pointed observations made between 2005 September 17 and 2016 July 17, for the blazars 3C 279 and PKS 1510−089, respectively, were also analyzed by us. We used the latest version of the calibration database (CALDB) and version 6.19 of the heasoft package. ²⁸ For each data set, we used the level 2 cleaned event files of the "photon counting" (PC) and "window timing" (WT) data acquisition modes generated using the standard xrtpipeline tool. The details of the XRT data analysis are given in Goyal et al. (2018, 2020), and we briefly outline the procedure here. For the PC mode data, the source and background spectra were generated using a circular aperture with appropriate region sizes and grade filtering using the xselect tool. An aperture radius of 47'' around the source position and four source-free regions of 118'' radius were used to estimate the source and background spectra, respectively. For the WT mode data, the source region of 47'' circular aperture and an annular region with inner and outer radii of 187'' and 281'', respectively, for background estimation were selected to estimate the source and background spectra. Following the XRT tutorial, ²⁹ appropriate corrections were made in the "BACKSCAL" keyword in the extracted source and background spectra to account for the different sizes of the source and background estimation regions. We also checked for the recommended pileup limit for the PC (∼0.5 count s⁻¹) and WT (∼100 count s⁻¹) modes. The ancillary response matrix was generated using the task xrtmkarf for the exposure map generated by xrtexpomap. All the source spectra were then binned over 20 points and corrected for the background using the task grppha. For each exposure, the spectral analysis between 0.3 and 10 keV was performed in an identical manner to that of the RXTE-PCA data. The unabsorbed 0.3–10 keV photon fluxes were averaged using a 1 day binning interval, and these light curves are shown in Figure 1(c) and Figure 2(c) for the targets 3C 279 and PKS 1510−089, respectively. In addition, we show the fitted Γ_X for the RXTE-PCA and Swift-XRT data sets analyzed in Figures 1(c) and 2(c) for the targets 3C 279 and PKS 1510−089, respectively. A few earlier studies report part of the Swift-XRT data sets analyzed by us for 3C 279 (Abdo et al. 2010b; Hayashida et al. 2015; Larionov et al. 2020) and PKS 1510−089 (Abdo et al. 2010c; D'Ammando et al. 2011; Nalewajko et al. 2012).

The blazars 3C 279 and PKS 1510−089 have been the targets of regular monitoring on daily timescales at multiple optical (B, V, R, and I) and infrared (J, H, and K) bands and here we combine many of these observations. We begin by using publicly available B, V, R, J, and K-band light curves from the Small and Moderate Aperture Research Telescope System (SMARTS; Bonning et al. 2012) program coordinated by Yale University ³⁰ for the period ∼2008 until 2017. We also take V-, R-, I-, J-, H-, and K-band light curves, covering the period 2005–2012, from the Rapid Eye Mounting telescope (REM; Zerbi et al. 2001) which were published in Sandrinelli et al. (2016, S16 hereafter). Many R-band optical photometric data were also obtained from the Tuorla Blazar monitoring program for the period 2004–2016 ³¹ (Takalo et al. 2008) and those published in (Chatterjee et al. 2008, C08 hereafter) for the period 1996–2004. Additionally, for the blazar PKS 1510−089 R-band photometric data were obtained from the VLBA-BU-BLAZAR blazar monitoring program for the period 2005–2018 ³² (Jorstad & Marscher 2016). For the blazar 3C 279, additional R-band photometric data for the period 2018–2019 were obtained from the Skynet Robotic Telescope Network program (SKYNET). ³³

There is a great deal of overlap between the SMARTS and the REM programs' data in the V, R, J, and K bands, and they complement each other well at the V, R, and J bands, with the pairs of magnitudes within 0.1 mag on the same nights, consistent with the calibration uncertainties provided by these programs. Moreover, the R-band photometric data obtained from the VLBA-BU and the Tuorla monitoring programs are comparable to other data sets within the measurement uncertainties (<0.1 mag) on the same nights. As a majority of these programs have mean sampling intervals >1 day, we have combined the data sets at the V, R, and J bands and average magnitudes were computed with 1 day binning intervals. This is done to avoid obtaining poor sampling on intranight timescales as well as to obtain homogeneous sampling on daily timescales. On the other hand, the K-band magnitudes from the SMARTS program has a magnitude offset ≥0.1 mag from that provided by S16 on the same nights. On visual inspection, the K-band SMARTS light curve appears more variable than the S16 light curve on similar timescales. Therefore, in this study, we use only the K-band data provided by S16. On a few occasions, fluxes at B, I, H, and K bands were also obtained multiple times on a given night, and we have averaged those observations over a 1 day binning interval.

Next, the optical and infrared magnitudes, m, were converted to flux densities using m_zero × 10^{−0.4× m} , where m_zero (=4063 Jy, 3636 Jy, 3064 Jy, 2635 Jy, 1590 Jy, 1020 Jy, and 600 Jy for the B, V, R, I, J, H, and K bands, respectively) refer to zero-point magnitude fluxes of the photometric system (Glass 1999). The errors in fluxes were derived using standard error propagation (Bevington & Robinson 2003). The resulting optical- and infrared-band light curves for the targets 3C 279 and PKS 1510−089 are presented in Figures 1(d), (e) and 2(d), (e), respectively.

GHz-band radio light curves of the targets were gathered from the Aalto University Metsähovi Radio Observatory (MRO) at 36.9 and 22.0 GHz (Teraesranta et al. 1998), the UMRAO at 14.5, 8.0, and 4.8 GHz (Aller et al. 1999), and the OVRO monitoring programs at 15 GHz (Richards et al. 2011). All the light curves gathered were binned using 1 day binning intervals as the typical sampling intervals of these programs are much larger (>3 days to 2 weeks). The radio light curves of the blazars 3C 279 and PKS 1510−089 for periods (1995–2017) that overlap monitoring with other frequencies are presented in Figures 1(f) and 2(f), respectively. The radio light curves for the full durations of the monitoring are shown in panels (a) and (b) of Figure 4.

Zoom In Zoom Out Reset image size

The rms-normalized periodogram is given as the squared modulus of its DFT for the evenly sampled light curve f(t_i ), observed at discrete times t_i and consisting of N data points,

$$\begin{eqnarray}\begin{array}{rcl}P({\nu }_{k}) & = & \displaystyle \frac{2\,T}{{\mu }^{2}\,{N}^{2}}\left\{{\left[\sum _{i=1}^{N}f({t}_{i})\,\cos (2\pi {\nu }_{k}{t}_{i})\right]}^{2}\right.\\ & & \left.+{\left[\sum _{i=1}^{N}f({t}_{i})\,\sin (2\pi {\nu }_{k}{t}_{i})\right]}^{2}\,\right\},\end{array}\end{eqnarray} \tag{ 1 }$$

where μ is the mean of the light curve and is subtracted from the flux values, f(t_i ). The DFT is computed for evenly spaced frequencies ranging from the total duration of the light curve, T, down to the Nyquist sampling frequency (ν_Nyq). Specifically, in our analysis, the frequencies corresponding to ν_k = k/T with k = 1,...,N/2, ν_Nyq = N/2T, and T = t_N − t₁ are considered. The constant noise floor level from measurement uncertainties is computed using (Isobe et al. 2015)

$$\begin{eqnarray}&&{P}_{\mathrm{stat}}=\displaystyle \frac{2\,T}{{\mu }^{2}\,N}\,{\sigma }_{\mathrm{stat}}^{2},\end{eqnarray} \tag{ 2 }$$

where, ${\sigma }_{\mathrm{stat}}^{2}={\sum }_{j=1}^{j=N}{\rm{\Delta }}f{({t}_{j})}^{2}/N$ is the mean variance of the measurement uncertainties for the flux values, Δf(t_j ), in the observed light curve at times t_j , with N denoting the number of data points in the observed light curve. As expected, the noise floor level due to measurement uncertainties has a white-noise character as the variability power contributed by measurement uncertainties alone is equal at all variability frequencies.

The unevenly sampled time series is rendered evenly sampled through linear interpolation between two consecutive data points with the interpolation interval, T_in, roughly 5–10 times smaller than the mean observed sampling interval; otherwise, the "spectral window function" corresponding to the sampling times gives a nonzero response in the Fourier domain, resulting in spurious powers in the periodograms (see Appendix A of Goyal 2020). In this work, the PSDs have been generated for the actual duration of the light curves down to the mean (observed) sampling intervals, T_mean. We note that interpolation is a smoothing operation, so it alters the high-frequency PSD of a data set, and this would be very obvious if we reported the spectrum up to 1/(2*T_in). The effects are reduced by cutting it off at 1/(2*T_mean), but it does not necessarily solve the problem. It depends on how much the actual sampling varies, or equivalently what the typical distance between the actual data point and the interpolate point used in the spectrum at 1/(2*T_mean) is. Because the light curves are relatively densely sampled, the PSDs have adequate coverage of many different frequencies. We expect a little extra noise because we do not have as many data pairs per frequency as would be the case for even sampling. These ramifications are not problematic in our case as we have plenty of coverage of the frequencies, and gridding the data is simply expedient. We demonstrate this qualitatively by calculating the typical separation of each sampling point from the even (i.e., T_mean) grid and computing its variance. For our light curves, we computed the square root of the variance of the sampling pattern normalized by T_mean, Var_sp = variance^1/2/T_mean (Column 7 of Table 2). This value should be zero for precisely evenly sampled data. As can be seen from Table 2, these values are close to zero for weekly binned and <2 for 3 hr binned Fermi-LAT light curves for both targets, as they are better sampled. For the optical and infrared light curves, the values range between ∼3.0 and ∼4.0 because of the relatively large (4–6 months) systematic gaps in monitoring due to the target's proximity to the Sun. For the X-ray- and GHz-band radio light curves, these values range between ∼1.4 and ∼3.0 for the majority of the gathered light curves, except for the 14.5 GHz UMRAO light curve for the target PKS 1510−089, for which it is 3.4. As only relatively small (weeks) systematic gaps are expected in the RXTE, Swift, and the radio data because of pointing constraints due to the target's proximity to the Sun, these light curves can be considered better sampled than the optical ones. Because for the great majority of our data these Var_sp values are relatively small numbers, our cutting of the spectrum down to 1/(2*T_mean) does not alter its shape.

Table 2. Summary of the Observation and PSD Analysis

Light Curve (Energy Range)	Epoch of Monitoring	Nobs	Ttotal	Tmean	Tin	Varsp	${\mathrm{log}}_{10}({P}_{\mathrm{stat}})$	${\mathrm{log}}_{10}({\nu }_{k})$ Range	β ± err	pβ b
	(Start-End)		(yr)	(day)	(day)		$\left(\tfrac{\mathrm{rms}}{\mathrm{mean}}\right)$ 2day	(day−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
3C 279
Fermi-LAT (0.1–300 GeV)	2008 Aug 8–2020 Sep 16	601	12.1	7.3	1.0	0.25	−0.84	−3.64 to −1.40	1.0 ± 0.7	0.595
RXTE-PCA (3–20 keV)	1996 Jan 22–2011 Dec 30	1751	15.9	3.3	0.5	3.16	−0.96	−3.76 to −1.11	1.2 ± 0.3	0.496
Swift-XRT (0.3–10 keV)	2006 Jan 13–2017 Jun 28	230	11.5	18.2	0.5	2.57	−0.43	−3.62 to −1.78	1.4 ± 0.7	0.999
Optical (B)	2008 May 17–2017 Jul 19	624	9.2	5.4	0.5	3.05	−1.41	−3.52 to −1.28	1.8 ± 1.3	0.567
Optical (V)	2005 Apr 7–2017 Jul 19	775	12.3	5.8	0.5	3.26	−0.21	−3.65 to −1.20	2.1 ± 1.4	0.651
Optical (R)	1995 Nov 25–2019 Jun 4	1829	23.5	4.7	0.5	3.17	−1.24	−3.93 to −1.07	1.4 ± 0.5	0.635
Optical (I)	2005 Apr 2–2011 Jan 30	183	6.2	12.5	0.5	2.82	−0.15	−3.35 to −1.53	1.6 ± 1.8	0.924
Infrared (J)	2005 Apr 9–2017 Jul 18	727	12.3	6.2	0.5	3.34	−1.37	−3.65 to −1.20	2.1 ± 1.5	0.832
Infrared (H)	2005 Apr 9–2011 Jun 23	197	6.2	11.5	0.5	3.37	−0.83	−3.35 to −1.53	1.8 ± 1.8	0.994
Infrared (K)	2005 Apr 9–2011 Jun 23	170	6.2	13.3	0.5	3.58	−0.31	−3.35 to −1.53	2.9 ± 1.6	0.992
MRO (36.9 GHz)	1979 Dec 21–2017 Jun 19	1740	37.5	7.9	0.5	2.97	−1.68	−4.13 to −1.48	1.7 ± 0.3	0.968
MRO (22.0 GHz)	1980 Jun 19–2004 Jun 21	751	24.0	11.7	0.5	2.55	−1.38	−3.94 to −1.49	2.9 ± 1.1	0.225
OVRO (15 GHz)	2008 Jan 8–2020 Dec 26	535	13.0	8.9	0.5	1.35	−2.11	−3.67 to −1.43	2.5 ± 1.1	0.813
UMRAO (14.5 GHz)	1974 May 8–2012 Jun 23	1567	38.1	8.9	0.5	2.00	−2.39	−4.14 to −1.49	2.6 ± 1.0	0.977
UMRAO (8.0 GHz)	1965 Aug 13–2012 May 16	1653	46.8	10.3	0.5	1.59	−2.40	−4.23 to −1.52	2.4 ± 0.9	0.871
UMRAO (4.8 GHz)	1978 Apr 15–2012 Jun 14	1101	34.2	11.3	0.5	1.64	−2.46	−4.09 to −1.64	2.1 ± 0.6	0.309
Flare 1 a (0.1–300 GeV)	2015 Jun 5–2015 Jun 26	80	0.057	0.24	0.042	1.84	−1.45	−1.29 to +0.11	1.4 ± 3.9	0.500
Flare 2 a (0.1–300 GeV)	2018 Jan 7–2018 Jan 28	146	0.057	0.14	0.042	0.38	−1.88	−1.31 to +0.30	2.0 ± 2.8	0.808
Flare 3 a (0.1–300 GeV)	2018 Apr 13–2018 May 4	157	0.057	0.14	0.042	0.27	−2.35	−1.31 to +0.30	1.5 ± 2.9	0.680
PKS 1510−089
Fermi-LAT (0.1–300 GeV)	2008 Aug 8–2020 Sep 16	603	12.1	7.3	1.0	0.30	−0.76	−3.64 to −1.40	0.9 ± 0.6	0.300
RXTE-PCA (3–20 keV)	1996 Dec 13–2011 Dec 30	1252	15.1	4.4	0.5	1.47	−0.79	−3.73 to −1.08	1.4 ± 0.3	0.482
Swift-XRT (0.3–10 keV)	2006 Aug 4–2017 Jun 26	181	10.9	22.0	0.5	3.04	−0.25	−3.59 to −1.76	1.5 ± 1.0	0.994
Optical (B)	2008 May 17–2017 Jun 14	542	9.1	6.1	0.5	3.29	−2.26	−3.52 to −1.28	1.5 ± 0.5	0.185
Optical (V)	2006 Jan 10–2017 Jun 14	611	11.4	6.8	0.5	3.48	−1.27	−3.62 to −1.38	1.2 ± 0.3	0.003
Optical (R)	2005 Mar 25–2018 Jul 10	1207	13.3	4.0	0.5	3.12	−2.15	−3.68 to −1.03	0.7 ± 0.4	0.532
Optical (I)	2006 Jan 11–2012 May 28	225	6.4	10.3	0.5	3.75	−1.09	−3.36 to −1.12	0.7 ± 0.8	0.245
Infrared (J)	2006 Feb 17–2017 Jun 11	675	11.3	6.1	0.5	3.86	−1.75	−3.61 to −1.37	0.2 ± 1.9	0.317
Infrared (H)	2006 Feb 17–2012 Jun 1	304	6.3	7.5	0.5	4.06	−1.44	−3.36 to −1.32	1.3 ± 4.4	0.003
Infrared (K)	2006 Feb 22–2012 Jun 1	223	6.3	10.3	0.5	3.59	−1.27	−3.36 to −1.52	0.2 ± 2.1	0.011
MRO (36.9 GHz)	1983 Apr 16–2017 Jun 19	885	34.2	14.1	0.5	2.51	−0.97	−4.09 to −1.65	1.8 ± 0.7	0.559
MRO (22.0 GHz)	1984 Jul 12–2004 Jun 21	239	19.9	30.5	0.5	1.44	−0.52	−3.86 to −2.03	1.4 ± 0.5	0.386
OVRO (15.0 GHz)	2009 Apr 2–2020 Dec 27	427	11.7	10.0	1.0	1.18	−2.17	−3.63 to −1.39	2.1 ± 1.0	0.505
UMRAO (14.5 GHz)	1974 Aug 14–2012 Jun 27	1345	37.9	10.3	0.5	3.41	−1.70	−4.14 to −1.48	1.6 ± 0.4	0.145
UMRAO (8.0 GHz)	1974 Sep 4–2012 May 18	1147	37.7	12.0	0.5	1.68	−1.55	−4.13 to −1.48	1.6 ± 0.3	0.164
UMRAO (4.8 GHz)	1979 Mar 23–2012 Jun 16	672	33.2	18.1	0.5	1.56	−1.19	−4.08 to −1.64	1.7 ± 0.3	0.199
Flare 1 a (0.1–300 GeV)	2011 Oct 14–2011 Nov 4	74	0.057	0.27	0.042	0.95	−1.34	−1.31 to +0.09	1.5 ± 3.6	0.978

Notes. (1) The instrument/band (energy range/frequency band used for observations in parentheses), (2) the beginning and ending dates of the light curve used, (3) the number of data points in the light curve, (4) the total duration of the observed light curve, (5) the mean sampling interval of the light curve (= Col. 4/Col. 3), (6) the interpolation interval used in the PSD analysis, (7) square root of the variance of the sampling pattern for the data normalized by T_mean, (8) the noise level in the PSD due to the measurement uncertainty, (9) the temporal frequency range covered by the binned logarithmic power spectrum, (10) the best-fit power-law slope of the PSD derived using the PSRESP method along with the corresponding errors representing the 98% confidence limit (see Section 3.2), and (11) the corresponding p_β .

^a Refers to Fermi-LAT light curves made with 3 hr integration bins (Figure 3), ^b A power-law model is considered a poor fit for p_β < 0.1 as the corresponding rejection confidence for the model is ≥90% (Section 3.2).

Download table as:  ASCIITypeset image

The sampling distribution and the corresponding spectral window functions of the observed light curves are shown in Appendices A1 and A2 for the blazars 3C 279 and PKS 1510−089, respectively. In our analysis, to avoid the distortions introduced due to the finite duration of the light curve, known as "red-noise leakage," the PSDs are generated using the "Hanning" window function (e.g., Press et al. 1992; Max-Moerbeck et al. 2014). The distortions in the PSD due to discrete sampling of the time series, i.e., aliasing, are not serious for red-noise-dominated time series and usually contribute a small amount of power at all frequencies (Uttley et al. 2002). The "raw" periodograms, obtained using Equation (1), provide a noisy estimate of the spectral power (as it consists of independently distributed χ² variables with two degrees of freedom (DOF); Timmer & Koenig 1995; Papadakis & Lawrence 1993; Vaughan et al. 2003); therefore, we average a number of them to obtain a reliable estimate. A binning factor of 1.6 is used, with the representative frequency taken as the geometric mean of each bin in our analysis (except for the R-band and 8.0 GHz PSDs of 3C 279 for which a binning factor of 1.4 is used). Next, the "true" power spectrum in the log–log space is obtained by offsetting the observed power spectrum by the expectation value of a χ² distribution with two DOF, which is equal to −0.25068 (Vaughan 2005).

As the PSDs are subjected to various distortions due to Fourier transformation, the PSRESP method is an efficient approach to obtain a best-fit spectral model (e.g., Uttley et al. 2002; Chatterjee et al. 2008; Max-Moerbeck et al. 2014; Meyer et al. 2019). We have previously employed this method in Goyal (2020, 2021); therefore, we only briefly discuss its main features here. In this method, an (input) PSD model is tested against the observed PSD, and the estimation of best-fit model parameters and their uncertainties is performed by varying the model parameters. A large number of mock light curves are generated with a known underlying power-spectral shape via Monte Carlo (MC) simulations. These simulated light curves are rebinned to mimic the observed data and subjected to the Fourier transformation in exactly the same manner as that of the observed data. The averaging of a large number of such PSDs gives the mean of the distorted model (input) power spectrum while the standard deviation around the mean periodogram gives the errors on the modeled (input) power spectrum. The goodness of fit of the model is estimated by computing two functions, similar to χ², defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{obs}}^{2}=\sum _{{\nu }_{k}={\nu }_{\min }}^{{\nu }_{k}={\nu }_{\max }}\displaystyle \frac{{\left[\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}({\nu }_{k})-{\mathrm{log}}_{10}{P}_{\mathrm{obs}}({\nu }_{k})\right]}^{2}}{{\rm{\Delta }}\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}{\left({\nu }_{k}\right)}^{2}}\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{\chi }_{\mathrm{dist},\ {\rm{i}}}^{2}=\sum _{{\nu }_{k}={\nu }_{\min }}^{{\nu }_{k}={\nu }_{\max }}\displaystyle \frac{{\left[\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}({\nu }_{k})-{\mathrm{log}}_{10}{P}_{\mathrm{sim},{\rm{i}}}({\nu }_{k})\right]}^{2}}{{\rm{\Delta }}\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}{\left({\nu }_{k}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

Here, ${\mathrm{log}}_{10}{P}_{\mathrm{obs}}$ and ${\mathrm{log}}_{10}{P}_{\mathrm{sim},{\rm{i}}}$ are the observed and the simulated log-binned periodograms, respectively; $\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}$ and ${\rm{\Delta }}\overline{{\mathrm{log}}_{10}{P}_{\mathrm{sim}}}$ are the mean and the standard deviation obtained by averaging large number of PSDs; k represents the number of frequencies in the log-binned power spectrum (ranging from ${\nu }_{k,\min }$ to ${\nu }_{k,\max }$), while i runs over the number of simulated light curves for a given β. The ${\chi }_{\mathrm{obs}}^{2}$ determines the minimum χ² for the model compared to the data, and the ${\chi }_{\mathrm{dist}}^{2}$ values determine the goodness of the fit corresponding to the ${\chi }_{\mathrm{obs}}^{2}$, as neither ${\chi }_{\mathrm{dist},{\rm{i}}}^{2}$ nor ${\chi }_{\mathrm{obs}}^{2}$ follow a standard χ² distribution. For this, the ${\chi }_{\mathrm{dist}}^{2}$ are sorted in ascending order. The probability, or p_β , that a given model can be rejected is then given by the percentile of ${\chi }_{\mathrm{dist}}^{2}$ distribution above which ${\chi }_{\mathrm{dist}}^{2}$ is found to be greater than ${\chi }_{\mathrm{obs}}^{2}$ for a given β (Uttley et al. 2002; see also Chatterjee et al. 2008 who call p_β as a success fraction). A large value of p_β or the success fraction represents a good fit in the sense that a large number of random realizations represent the shape and slope of the intrinsic (input) PSD. Subsequently, the PSRESP method allows for an estimation of best-fit model as well as the goodness of fit. It essentially uses the MC approach toward a frequentist estimation of the quality of the model compared to the data and so is a good approach when the fit statistics are not well understood (Press et al. 1992).

We use the method of Emmanoulopoulos et al. (2013) for light-curve simulations as it preserves both the probability density function (PDF) of the flux distribution as well as the underlying power spectral shape and not just the power spectral shape as is the case with the method of Timmer & Koenig (1995). In our light-curve simulations, we assume a single power-law spectral shape with a given β and a log-normal flux distribution. For this purpose, the mean and the standard deviation were computed by fitting a Gaussian function to the logarithmically transformed flux distributions, and the measurement errors were incorporated by adding a Gaussian random variable with mean 0 and standard deviation equal to the mean error of the measurement uncertainties on the observed flux values. In such a manner, 1000 light curves are simulated in the β range 0.1–3.0, with a step of 0.1 (and in the range 0.1–4.0 for the 3 hr integration bin LAT light curves). The best-fit PSD slope for the observed PSD is given by the one with the highest p_β value. The error on the best-fit PSD slope is obtained by fitting a Gaussian function to the p_β distribution curve (e.g., Chatterjee et al. 2008; Bhatta et al. 2016a; Bhattacharyya et al. 2020) and is given as 2.354σ, where σ is the standard deviation of the Gaussian function. Therefore, the errors provide a 98% confidence interval on the best-fit PSD slope (Goyal 2021).

We note that a common approach to derive PSDs of an unevenly sampled time series is the Lomb–Scargle periodogram (LSP; Lomb 1976; Scargle 1982), which is also frequently used in the AGN variability literature (Rani et al. 2009, 2010; Bhatta et al. 2016b; Sandrinelli et al. 2016; Castignani et al. 2017; Sandrinelli et al. 2017; Gaur et al. 2018; Gupta et al. 2019; Żywucka et al. 2020; Bhatta & Dhital 2020; Tarnopolski et al. 2020; Li et al. 2021). However, with numerical simulations, we demonstrate in Appendix C that the LSP method does not usually reproduce the shape of the underlying power spectrum correctly up to the highest (Nyquist) frequency, whereas the DFT method does. Therefore, the DFT method (with linear interpolation) is a more appropriate choice as far as the characterization of the PSD shape is concerned when the available sampling is not very close to uniform.

Details on the light curves used for the analysis and the PSDs derived for the given light curve, along with the best-fit PSDs and the maximum p_β , are summarized in Table 2. Figures 5 and 6 present the corresponding best-fit PSDs for the analyzed light curves. The corresponding p_β distribution curves are given in Appendix B for the blazars 3C 279 and PKS 1510−089, respectively. In our analysis, we have not subtracted the constant noise floor level. The PSRESP method also allows us to compute the rejection confidence for the input PSD shape; the maximum probability lower than x% means that the rejection confidence (1 − p_β ) is higher than (100 − x)% for the (input) PSD model. In our analysis, we set the rejection confidence at 90% for the input PSD model. Moreover, a direct comparison of the average fractional variability power for a given timescale is given by the ν_k P(ν_k ) versus ν_k curves for the acceptable PSD fits in Figure 7 for blazars 3C 279 (panels (a), (b), and (c)) and PKS 1510−089 (panels (d), (e), and (f)). The quantity ν_k P(ν_k ) is equivalent to the square fractional variability (see Goyal 2020 and references therein). In addition, we note that Goyal (2021) derived the R-band intranight variability PSDs for the blazar 3C 279 (monitored on 2006 January 26, 2006 February 28, and 2009 April 20) and PKS 1510−089 (monitored on 2009 May 1) using the PSRESP method. On two occasions (2006 January 26 and 2006 February 28), the intranight PSDs for the blazar 3C 279 showed a good fit to a power-law model (p_β > 0.1) while the intranight PSD obtained on 2009 April 20 gave a bad fit (p_β < 0.1). Similarly, the intranight PSD obtained for the blazar PKS 1510−089 gave a bad fit. Therefore, we extend the frequency range of the optical power spectrum to minute timescales for the blazar 3C 279 using the 2006 January 26 and 2006 February 28 intranight PSDs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be seen from Table 2 that the derived 98% confidence limits on the best-fit PSD slopes using the PSRESP method are large in some cases. Therefore, it is fortunate that we can obtain representative PSD slopes (and the corresponding errors) on emission frequencies by averaging several PSD slopes at nearby frequencies. Such an approach is possible solely due to the availability of the large number of light curves analyzed by us in this study. For this, we average all the PSD slopes at different radio and optical frequencies; this gives the mean PSD slope at two different synchrotron frequencies while the combination of two X-ray PSD slopes and one GeV PSD slope gives an estimate of the mean PSD slope at IC frequencies of the emission spectrum (see Section 1). As the UMRAO 14.5 GHz and OVRO 15.0 GHz data sets are essentially the same (Aller et al. 1985; Richards et al. 2011), for the calculation of the mean PSD slope at radio-synchrotron frequencies we did not use the OVRO estimates because the UMRAO 14.5 GHz PSD covers a wider range of variability frequencies. For averaging, we used the PSD slopes that were successfully represented by a single power-law model (i.e., p_β > 0.1). The mean PSD slope, β_mean, was computed in a straightforward manner while the error was computed using the MC bootstrap method as follows. For a sample of PSD slopes, a slope is drawn from a Gaussian distribution of mean and standard deviation equal to β and σ of the best-fit PSD estimate, respectively (note that Table 2 reports errors equal to 2.354σ of the fitted Gaussian function). The β_mean is computed. These two steps are repeated 1000 times. The error is computed as the standard deviation of the distribution of β_mean values. Table 3 gives the summary of the mean PSD slopes and the corresponding standard deviation (1σ error) for the synchrotron and IC portions of the emission continuum for these two blazars.

The γ-ray PSD slopes computed with the PSRESP technique used here for the blazars 3C 279 and PKS 1510−089 employing the decade-long light curves are consistent with those given in previous studies (β ∼0.6 and ∼1.0 for 3C 279 and PKS 1510−089, respectively; Meyer et al. 2019). Nilsson et al. (2018) reported the optical R-band PSD slopes ∼1.5 and ∼0.9 for the blazars 3C 279 and PKS 1510−089, respectively, which are also consistent within errors with our estimates. Chatterjee et al. (2008) reported the PSD slopes for the blazar 3C 279 at X-ray, R-band optical, and radio frequencies as ∼2.3, ∼1.7, and ∼2.3, respectively, using the decade-long RXTE-PCA data, R-band, and 14.5 GHz UMRAO data sets for the period 1996–2007. Note that our longer-duration X-ray, optical R-band, and 14.5 GHz light curves include data sets reported in Chatterjee et al. (2008). The X-ray PSD slope obtained in their analysis of RXTE data is significantly steeper, β = 2.3 ± 0.3, as compared to our result of 1.2 ± 0.3. Our X-ray PSD slope is corroborated by the Swift-XRT data set over a similar timescale (β = 1.4 ± 0.7); therefore, this discrepancy cannot be due to the different lengths of light curves used in Chatterjee et al. (2008) and in our analyses for the PCA data. Good matches (within uncertainty), however, are found between the optical and radio PSD slopes between our analysis and that of Chatterjee et al. (2008). Using the full-duration UMRAO light curves for both blazars, Park & Trippe (2017) reported the PSD slopes, β ∼2.5, ∼2.2, and ∼1.6 at 14.5 GHz, 8.0 GHz, and 4.8 GHz, respectively, for 3C 279 and β ∼ 1.5, ∼1.3, and ∼1.4 at 14.5 GHz, 8.0 GHz, and 4.8 GHz, respectively, for PKS 1510−089; these are consistent within errors with our results.

We note that Ryan et al. (2019) performed direct modeling of 9.5 yr long, weekly binned γ-ray light curves for a sample of blazars, including the targets studied here, using the CARMA process (Kelly et al. 2014). Their analysis revealed a break on a year-like timescale in the PSD of the blazar PKS 1510−089 while the PSD of 3C 279 showed no hint of a break (meaning that the relaxation timescale could be longer than the length of the time series for this source). Our PSD analysis for the blazar PKS 1510−089 contradicts the results of Ryan et al. (2019) in the sense that our γ-ray PSD is best fitted with a single power-law model with a reasonably good probability (p_β = 0.3; Table 2). For the blazar 3C 279, Ryan et al. (2019) report no hint of a break, indicating a good match between our PSD modeling and that reported from the CARMA modeling.

We now mention the reported QPOs on year-like timescales in γ-ray and optical/infrared light curves for the period 2005–2014 claimed by Sandrinelli et al. (2016) for the blazars 3C 279 (V, R, and K bands) and PKS 1510−089 (V, R, J, and K bands). We also use the same data set presented in Sandrinelli et al. (2016) and supplement it with other data sets to increase the length of the analyzed light curves. Our γ-ray PSDs are well represented by the single power-law spectral shapes with high confidence (p_β ≥0.76). Moreover, the optical/infrared PSDs of the blazar 3C 279 are represented by the single power-law spectral form adopted in our analysis with high confidence (p_β > 0.6). For the blazar PKS 1510−089, the PSDs obtained at the B, R, I, and J bands are also well represented by a single power-law spectral form with decent confidence (p_β > 0.18, meaning that the rejection confidence for the single power-law spectral shape is <82%). The PSDs obtained at the V, H, and K bands, where the data are more limited, give poor fits, with rejection confidence ≥95%. The single power-law PSDs obtained in our analysis argue against the presence of reported QPOs at the V, R, and K bands for the blazar 3C 279 and at the R and J bands for the blazar PKS 1510−089. This discrepancy between our results and those of Sandrinelli et al. (2016) could arise from the use of different methods of PSD estimation (LSP in their analysis versus DFT in our analysis) as well as the longer lengths of the light curves used in our study. If a QPO were always present, a longer-duration light curve should reveal a higher significance for it, but if a QPO were transitory, then additional data would weaken its significance. The lack of reproducibility is not uncommon in blazar studies as many previously claimed QPOs in the literature could not be confirmed when different methods or longer-duration light curves were used for periodicity analyses (for recent discussions, see Goyal et al. 2018; Covino et al. 2019, 2020; Tarnopolski et al. 2020; Yang et al. 2021). Further investigation of this point is beyond the scope of this paper.

As the single-zone models do not explain the statistical properties of long-term multiwavelength variability adequately (Section 1), based on the different statistical characters of multiwavelength variability exhibited by the BL Lac object PKS 0735+178 (β ∼ 2 for radio/optical emissions but ∼1 at HE γ-ray energies), we hypothesized in Goyal et al. (2017) that the entire broadband emission is produced in an extended and turbulent jet. We noted there that these results would arise if a single stochastic process is responsible for the variability at synchrotron frequencies while a linear superposition of two or more stochastic processes with widely different relaxation timescales produce variability at IC frequencies. We speculated that the drivers behind the synchrotron variability can be stochastic fluctuations in the jet plasma conditions (such as variations in the magnetic field or bulk Lorentz factor), leading to dissipation of energy over wide spatial scales. The radiative response of the accelerated particles will be delayed with respect to the input perturbations, thereby producing the red-noise variability at synchrotron frequencies on timescales smaller than the relaxation timescales of the process. Because we did not recover the flattening of PSDs on longer timescales (≥1000 days), we noted that the driver should be related to some global magnetohydrodynamic timescale. At IC frequencies, however, due to inhomogeneities in the population of photons available for upscattering, an additional stochastic process with relaxation timescale equal to a light-crossing time ∼1 day (a valid approximation for a boosted jet with bulk Lorentz factor 30) was envisaged. IC variability is expected to show red-noise character on timescales shorter than a ∼day in our model. Our hypothesis was supported in Goyal (2020) where the VHE γ-rays also exhibited β ∼ 1 for the blazars Mrk 421 and PKS 2155−304 whose TeV emission is often concluded to be due to SSC processes (Abdo et al. 2011a; Abramowski et al. 2012; Aleksić et al. 2012; Madejski et al. 2016; Abdalla et al. 2020; Dmytriiev et al. 2021).

The best-fit PSD slopes for all of the individual emission frequencies that we report here range between β ∼ 0.2–2.9; however, the statistically significant slopes span a narrower range, 0.7–2.9, which correspond to the flicker and red-noise statistical characters of variability. We note that in many cases, the reported 98% confidence limit on the PSD estimate is large (Table 2), which renders a strict statistical classification of the variability unfeasible. The mean PSD slopes, however, obtained at radio-synchrotron and optical-synchrotron frequencies correspond to statistical characters consistent with strict red-noise β ∼2 (except for PKS 1510−089 for which the mean β at optical/infrared frequency is ∼0.8) and to a strict flicker-noise (β ∼ 1) statistical character at IC emission frequency on timescales ranging from decades/years to weeks/days. As such, their broadband variability behaviors basically can be understood within our model, as summarized in the previous paragraph and given in more detail in Goyal et al. (2017).

For PKS 1510−089 the optical/infrared variability seems to be better represented by flicker noise whenever a single power-law form was found to be an adequate fit. This is not expected within our model (Section 5.3) as the emission at optical/infrared frequencies is presumed to be dominated by synchrotron processes for such FSRQ type blazars, although the radio variability is well represented by red-noise processes for the same source and, of course, that is synchrotron radiation. Therefore, our model cannot be considered inadequate based on the flicker-noise character at synchrotron frequencies because in such a case we expect the radio frequencies to mimic the infrared/optical frequencies.

However, the FSRQ PKS 1510−089 is well known to have a prominent thermal accretion disk component in its SED at UV/infrared frequencies, not only during quiescent states but even when flaring (e.g., D'Ammando et al. 2011; Nalewajko et al. 2012; Aleksić et al. 2014; Ahnen et al. 2017); therefore, one expects different variability properties if the jet emission (as the source shows significant optical polarization variability at times; Itoh et al. 2016; Beaklini et al. 2018) is diluted by the accretion disk component. Assuming that the majority of the disk emission originates within 100 Schwarzschild radii of the accretion disk (Shakura & Sunyaev 1976) with viscosity parameter α = 0.01, the characteristic light-crossing, gas orbital, and thermal/viscous timescales for this blazar turn out to be ≤1 day, ∼31–73 days, and ∼1.4–3.2 yr, respectively (Kelly et al. 2009) for a black hole of mass range (3–7)×10⁷ M_⊙ (Table 1). Because the optical/infrared PSDs cover 6–12 yr down to weeks timescale, such instabilities operating on these timescales seem to be valid candidates for the production of fluctuations in the disk emission in this source (e.g., Wiita 2006). Therefore, the flatter optical/infrared PSD slopes for this source could be reconciled if the emission at optical/infrared frequencies is driven by a linear combination of stochastic process(es) driving the disk emission at shorter timescales (days to years) and the jet emission, which could include longer relaxation timescales (≈1000 days) (see, in this regard, Kelly et al. 2011; Sobolewska et al. 2014).

The combination of γ-ray measurements made over the long-term (obtained with 1 week integrations) and short-term (obtained with 3 hr integrations for short durations only when source fluxes were high) allowed us to characterize the γ-ray variability PSDs across four decades. The mean PSD slope at short timescales is consistent with those at longer timescales, within the large uncertainties of the former, and as such does not require the presence of a break as apparently seen in some blazars around ∼days timescales (Sobolewska et al. 2014; Ryan et al. 2019). The optical PSD for the blazar 3C 279 could be characterized across six decades due to the inclusion of intranight PSDs on two occasions (2006 January 26 and 2006 February 28; Goyal 2021). The continuity of normalized power spectra across six decades indicates that a continuous (single) stochastic process appears to drive the variability at optical frequencies.

In none of the analyzed light curves, however, do we recover the relaxation timescale of the process driving the variability, which should entail seeing β changing from ≥1 to ∼0 at frequencies below the inverse of that timescale. The γ-ray and optical/infrared PSDs are in particular well represented by single power-law spectral shapes with high confidence up to timescales of the order of a decade or more. As the majority of HE jet emission is believed to originate at ∼1 pc distance from the supermassive black hole (see Sikora et al. 2009; Harvey et al. 2020, and references therein), the corresponding relaxation timescales for the γ-ray variability should turn out to be a few months in the observer's frame (assuming a typical bulk Lorentz factor of 10; Lister et al. 2016). The absence of PSD breaks in these blazars implies that the γ-ray emission originates, like the optical/infrared emission, from a rather extended region, in accordance with the expectation of our hypothesis discussed in Section 5.3.

The noise-like appearance of blazar light curves initiated efforts to model the jet emission due to many cells behind a shock (Marscher 2014; Calafut & Wiita 2015; Pollack et al. 2016). The detailed model developed by Marscher (2014) provides plausible flux light curves at both synchrotron and SSC emission frequencies but did not provide the PSDs, making a direct comparison difficult with these multiband variability observations. In contrast, the modeling done by Calafut & Wiita (2015) and Pollack et al. (2016) computes light curves as well as PSDs, but only at synchrotron emission frequencies. For relativistic hydrodynamic simulations of 2D jets where the emission and its variability are shaped by the combination of bulk Lorentz factor fluctuations and the turbulence within the jets, PSD slopes in the range 2.1 to 2.9 due to bulk Lorentz factor changes and in the range 1.7 to 2.3 due to turbulence are expected (Pollack et al. 2016). In this scenario, the majority of radio-optical PSD slopes (synchrotron emission frequencies) obtained in our study can be reconciled with those attributed primarily to changes in bulk Lorentz factors for the blazar 3C 279 and to turbulence for the blazar PKS 1510−089. However, such models would need to be extended to show that they naturally can produce significantly flatter PSD slopes at IC energies, before considering them actually supported by our analysis.

In summary, we have explored the statistical properties of blazar light curves by means of variability PSDs using good photometric quality (∼1%–15% flux accuracy), extensive multiwavelength data sets (≥15 different frequencies covering 13 decades of the emission spectrum) available for the blazars 3C 279 and PKS 1510−089 over timescales exceeding a decade. The majority of PSDs show a good fit to single power-law spectral shapes over temporal frequencies ranging between ∼10⁻⁴ and ∼10⁻¹ (day⁻¹). However, we note that PSDs on temporal frequencies longer than 10^−2.7 (day⁻¹) are poorly sampled, providing only a few up to 10 cycles in the PSDs.

The mean PSD slopes at radio and optical-synchrotron emission frequencies are steeper from those at IC emission frequencies at the ≥2σ level on timescales ranging from decades to weeks/days for the blazar 3C 279. On similar broad timescales, the mean PSD slopes at radio and optical-synchrotron emission frequencies and at IC emission frequencies turn out to be statistically indistinguishable (≤2σ level) from one another for the blazar PKS 1510−089. The γ-ray variability PSDs covers ∼4 decades of the variability spectrum for these two blazars where the short-term variability PSDs connect smoothly to the long-term variability PSDs. The optical variability power spectrum covers ∼6 decades for the blazar 3C 279, and the normalization of intranight PSDs appears to be a simple extrapolation from longer timescales.