Stochastic Modeling of Multiwavelength Variability of the Classical BL Lac Object OJ 287 on Timescales Ranging from Decades to Hours

Blazars are a major class of active galactic nuclei (AGNs), whose total radiative energy output is dominated by the Doppler-boosted, broadband, and nonthermal emission of relativistic jets launched by accreting supermassive black holes from the centers of massive elliptical galaxies (Begelman et al. 1984; de Young 2002; Meier 2012). The blazar class includes BL Lacertae objects (BL Lac objects) and flat-spectrum radio quasars (FSRQs; see Urry & Padovani 1995). In the framework of the "leptonic" scenario for blazar emission, the radio-to-optical/X-ray segment of the continuum emission is produced by synchrotron radiation of electron–positron pairs (e^±) accelerated up to ∼TeV energies, while the high-frequency X-ray-to-γ-ray segment is widely believed to be caused by the inverse-Comptonization of various circumnuclear photon fields (produced both internally and externally to the outflow) by the jet electrons (e.g., Ghisellini et al. 1998). Alternatively, in the "hadronic" scenario, the high-energy emission continuum could also be generated via protons accelerated to ultrahigh energies (≥EeV) and could produce γ-rays via either direct synchrotron emission or meson decay and synchrotron emission of secondaries from proton–photon interactions (e.g., Böttcher et al. 2013).

Blazars display strong flux variability at all wavelengths from radio to γ-rays, on timescales ranging from decades down to hours, or even minutes. The observed flux changes are often classified broadly into three major categories, namely, "long-term variability" (corresponding to timescales of decades to months), "short-term variability" (weeks to days), and intranight/day variability (timescales <1 day; see, e.g., Wagner & Witzel 1995; Ulrich et al. 1997; Falomo et al. 2014). During the past decade, special attention has been paid to catching and characterizing large-amplitude and extremely rapid (minute-/hour-long) flares in the γ-ray regime, with observed intensity changes of up to even a few orders of magnitude (Aharonian et al. 2007; Aleksić et al. 2011; Foschini et al. 2011; Rani et al. 2013; Saito et al. 2013; Ackermann et al. 2016). However, these are rare, rather exceptional events while, in general, the multiwavelength variability of blazar sources is of a "colored noise" type, meaning larger variability amplitudes on longer variability timescales with only low-level (few percent) flux changes on hourly timescales.

More precisely, the general shapes of the power spectral densities (PSDs) of blazar light curves, which may typically be approximated to first order by a single power law P(f) = A f^−β (where f is the temporal frequency corresponding to the timescale 1/2πf, A is the normalization constant, and β > 0 is the spectral slope), indicate that the flux changes observed in given photon energy ranges are correlated over temporal frequencies (Goyal et al. 2017 and references therein). So far, little or no evidence has been found for flattening of the blazar variability power spectra on the longest timescales covered by blazar monitoring programs (i.e., years and decades), even though such a flattening is expected in order to preserve the finite variance of the underlying (uncorrelated, by assumption) process triggering the variability (but see in this context Kastendieck et al. 2011; Sobolewska et al. 2014). Breaks in the PSD slope (from 0 < β < 2 down to β > 2 at higher temporal frequencies) reported in a few cases may, on the other hand, hint at characteristic timescales related to either a preferred location of the blazar emission zone, the relevant particle cooling timescales, or some global relaxation timescales in the systems (Kataoka et al. 2001; Sobolewska et al. 2014; Finke & Becker 2015). The detection of such break features in blazar periodograms would therefore be of primary importance for constraining the physics of blazar jets. Owing to observing constraints (e.g., weather, visibility), however, the blazar light curves from ground-based observatories are always sampled at limited temporal frequencies. This issue is particularly severe at optical and very high energy (VHE) γ-ray energies, where generally timescales corresponding to ∼12–24 hr can hardly be probed. This difficulty has recently been surmounted in the optical range with the usage of Kepler satellite data, though only for rather limited numbers of blazars/AGNs (Edelson et al. 2013; Revalski et al. 2014).

OJ 287 (J2000.0, α = 08^h54^m48875, δ = +20°06'3664; redshift z = 0.3056; Nilsson et al. 2010) is a typical example of a "low-frequency-peaked" BL Lac object with positive detections in the GeV and TeV photon energy ranges (Abdo et al. 2010; O'Brien 2017). It is highly polarized in the optical band (PD_opt > 3%; Wills et al. 2011) and exhibits a flat-spectrum radio core with a superluminal parsec-scale radio jet, both being characteristic of blazars (Wills et al. 1992; Lister et al. 2016). A supermassive black hole binary was claimed in the system, based on the evidence for a ∼12 year periodicity in its optical and radio light curves (Sillanpaa et al. 1996; Valtaoja et al. 2000; Valtonen et al. 2016); in addition, hints for a quasi-periodicity, with a characteristic timescale of ∼400–800 days, have been reported for the blazar based on the decade-long optical/near-infrared and high-energy γ-ray light curves. (See, in the multiwavelength context, Sandrinelli et al. 2016 and Bhatta et al. 2016 for the most updated list of claims of quasi-periodic oscillation (QPO) detections in blazars, in general.) OJ 287 is, in fact, one of the few blazars for which good-quality, long-duration optical monitoring data are available, dating back to around 1896 (Hudec et al. 2013). It is also one of the few blazars that have been observed by the Kepler satellite, for a continuous monitoring duration of 72 days with cadence becoming as small as 1 minute. Hence, OJ 287 is an outstanding candidate for characterizing statistical properties of optical flux changes on timescales ranging from ∼100 years to minutes.

Here, for the first time, we present the optical PSD of OJ 287 covering—with no gaps—about six decades in temporal frequency, by combining the 117 year long optical light curve of the source (using archival as well as newly acquired observations with daily sampling intervals) with the Kepler satellite data. The source has also been monitored in the radio (GHz) domain with a number of single-dish telescopes, in X-rays by the space-borne Swift's X-Ray Telescope (XRT), and in the high-energy γ-ray range with the Large Area Telescope (LAT) on board the Fermi satellite. Here, we utilize these massive data sets to derive the radio, X-ray, and γ-ray PSDs of OJ 287 and compare them with the optical PSD.

In Section 2, we describe in more detail all the gathered data as well as the data-reduction procedures. The data analysis and the results are given in Sections 3 and 4, respectively, while a discussion and our main conclusions are presented in Section 5.

2.1. High-energy $\gamma$-rays: Fermi-LAT

2.2. X-Rays: Swift-XRT

2.3. Optical: Ground-based Telescopes and Kepler

The long-term optical data presented in this work have been gathered from several sources and monitoring programs listed in Table 1, including newly acquired measurements, together resulting in a very long optical light curve ranging from 1900 to 2017 February. We note that, starting from 2015 September 2, the blazar OJ 287 has been the target of the dense multiwavelength optical monitoring campaign "OJ287-15/16 Collaboration" led by S. Zola, which was undertaken because of the predicted giant outburst in the system related to the ∼12-year-long periodicity of a putative supermassive black hole binary (see Sillanpaa et al. 1996; Valtonen et al. 2016 and references therein). Details of the "OJ287-15/16 Collaboration" monitoring program, including the list of observatories, are given in Valtonen et al. (2016). All of the data were taken from the start of 2016 through 2017 February using the skynet robotic telescope network⁷⁵ and the Mt. Suhora telescopes, with external observers from Greece, Ukraine, and Spain (see Zola et al. 2016). For these newly acquired optical data, including also the Kraków quasar monitoring program,⁷⁶ data reduction was carried out using the standard procedure in the Image Reduction and Analysis Facility (iraf)⁷⁷ software package.

Table 1.  Optical Observations of OJ 287 Made from Ground-based Observatories Included in the Present Work

Database	Monitoring Epoch (UT)	Filter	N	References
(1)	(2)	(3)	(4)	(5)
Harvard College Observatorya	1900 Oct 4–1988 Nov 16	B	272	Hudec et al. (2013)
Sonneberg Observatory	1930 Dec 20–1971 May 25	V	91	Valtonen & Sillanpää (2011)
Partly historicalb	1971 Mar 25–2001 Dec 29	R	3717	Takalo (1994), this work
Catalina sky surveyc	2005 Sep 4–2013 Mar 16	V	606	this work
Perugia and Rome database	1994 Jun 3–2001 Nov 5	R	802	Massaro et al. (2003)
Shanghai Astronomical Observatory	1995 Apr 19–2001 Dec 29	R	71	Qian & Tao (2003)
Tuorla monitoringd	2002 Dec 7–2011 Apr 14	R	1525	Villforth et al. (2010), this work
Krakow quasar monitoringe	2006 Sep 19–2017 Feb 20	R	1155	Bhatta et al. (2016), this work

Notes. Columns: (1) name of the observatory/university/monitoring program; (2) period covered by the monitoring program (start–end); (3) observing filter; (4) number of the collected data points; and (5) references for the data (either full or partial data sets).

^aB-band measurements listed by Hudec et al. (2013), after applying quality cuts and removing the upper limits. ^bPartly displayed in Figures 3(a)–(c) of Takalo (1994), converted to the R band by using a constant color difference. ^c http://www.lpl.arizona.edu/css/ ^d http://users.utu.fi/kani/1m/ ^e http://www.as.up.krakow.pl/sz/oj287.html

Download table as:  ASCIITypeset image

The procedure starts with pre-processing of the images through bias subtraction, flat-fielding, and cosmic-ray removal. The instrumental magnitudes of OJ 287 and the standard calibration stars listed by Fiorucci & Tosti (1996) in the image frames were determined by aperture photometry using apphot. This calibration was then used to transform the instrumental magnitude of OJ 287 to a standard photometric system. Our data have quoted photometric uncertainties of ∼2%–5%, arising mainly from large calibration errors in the estimated magnitudes of the stars in the field (Fiorucci & Tosti 1996). A typical 0.2 mag calibration uncertainty is assumed for B-band photographic magnitudes listed by Hudec et al. (2013). In the cases when during a given night the flux has been measured multiple times with different telescopes, we have averaged the measurements over the daily intervals. For the flux measurements obtained in B and V, fixed color differences of B − R = 0.87 mag and V − R = 0.47 mag were used to convert to photometric R-band magnitudes in the standard Landolt photometric system (Takalo et al. 1994). Finally, for a given R-band magnitude m_R, the R-band flux was derived as $3064\times {10}^{-0.4{m}_{R}}\,\mathrm{Jy}$, where 3064 Jy is the zero-point flux of the photometric system (Glass 1999). The uncertainties in the R fluxes were derived using standard error propagation (Bevington & Robinson 2003). The resulting long-term R-band historical light curve of OJ 287 is presented in Figure 1.

Zoom In Zoom Out Reset image size

OJ287 was also observed during Campaign 5 of the Kepler's ecliptic second-life (K2) mission. The Kepler spacecraft contains a 0.95 m Schmidt telescope with a 110 deg² field of view imager having a pixel size of 4''. It is in a heliocentric orbit currently about 0.5 au from Earth and provides high-cadence, very high precision (1 part in 10⁵) photometry for rather bright stellar targets (Howell et al. 2014). Campaign 5 lasted for 72 days, starting on 2015 April 27 and ending on 2015 July 10, data accumulating with both long (29.4 minutes) and short (58.85 s) cadences. The data are publicly available from the Barbara A. Mikulski Archive for Space Telescopes (MAST) and stored in a fits table format.⁷⁸ Timestamps are provided in Barycentric Julian Date (BJD). The long- and short-cadence data are not independent; that is why we used only the latter, which provide higher temporal resolution. In the short-cadence mode, for each timestamp a target mask is provided, while an optimal aperture is not, meaning that the light curves are not extracted. To estimate the fluxes and their uncertainties, we applied our customized scripts based on the tasks daofind and phot in IRAF. We applied an aperture with a radius of 4 pixels. The background has been already subtracted during the in-house processing. The extracted light curve was subject to 4σ clipping. The uncertainties were estimated following the recipe given in the phot manual.

K2 data analysis must struggle with onboard systematics, such as the thruster firings, reaction-wheel jitter, temperature-change-induced CCD module sensitivity, and solar-pressure-induced drift (Howell et al. 2014). Most significant are the "thruster firings," which cause targets to drift across the detector and usually occur at intervals of 6 hr. The light curves are also affected at the few percent level by differential velocity aberration (DVA), which originates in the motion of the spacecraft in its annual orbit around the Sun. However, the DVA and thruster firings need only be of concern if they mask the intrinsic variability present in the light curve. In our case, the thruster firings add to the nominal variations at ∼6 hr intervals and the DVA has relevance on multi-week timescales. However, these rarely seem to exceed the intrinsic variations on those timescales, and we have decided not to attempt to correct for them. Unfortunately, the analysis approaches used as countermeasures for the drifts that are appropriate for detecting exoplanets also remove real astrophysical brightness variations on the same timescales (e.g., Revalski et al. 2014). Hence, we assume that the K2 light curve found in this way is dominated by the intrinsic variability. Although some influence from temperature changes can still be present in the light curve, this systematic could not be removed owing to the lack of calibration files provided by the archive. The K2 short-cadence light curve of OJ 287 produced in this fashion is shown in Figure 2.

Zoom In Zoom Out Reset image size

2.4. Radio Frequencies: UMRAO and OVRO

The radio data were obtained from the University of Michigan Radio Astronomy Observatory (UMRAO) 26 m dish at 4.8, 8.0, and 14.5 GHz and the 40 m telescope at the Owens Valley Radio Observatory (OVRO) at 15 GHz. The UMRAO fluxes at 4.8, 8.0, and 14.5 GHz were typically measured weekly (Aller et al. 1985), from 1979 March 23 to 2012 June 15, from 1971 January 27 to 2012 May 16, and from 1974 June 20 to 2012 June 23, respectively. The OVRO light curve at 15 GHz was sampled twice a week (Richards et al. 2011), during the period from 2008 January 8 to 2016 November 11. Discussions of the corresponding observing strategies and calibration procedures can be found in Aller et al. (1985) for the UMRAO data and in Richards et al. (2011) for the OVRO data. Figure 3 (bottom panel) shows the resulting long-term radio light curves of OJ 287, compared with the high-energy γ-ray, X-ray, and optical R-band light curves for the overlapping monitoring epochs. The zoomed-in version of the plot is shown in Figure 4 to highlight the Swift-XRT and Fermi-LAT data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 presents a summary of these observational data.

Table 2.  Observational Summary for the Light Curves of OJ 287

Data Set	N	Tobs	ΔTminutes	ΔTmax	ΔTmed	ΔTmean	${\sigma }_{\mathrm{stat}}^{2}$	log(Pmed)	log(Pmean)	log(Frequency Range)
		(years)	(day)	(day)	(day)	(day)	(rms2)	(rms2 day)	(rms2 day)	(day−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Fermi-LAT	199	8.5	14	70	14	15.5	3.0 10−1	+0.78	+0.81	−3.5 to −1.9
Swift-XRT	239	11	1	528	2.9	18.2	2.3 10−2	−0.86	−0.07	−3.6 to −2.1
Optical (all)	3490	117	1	1728	1.6	12.2	8.9 10−2	−0.54	+0.33	−4.6 to −1.4
Optical (trun.)	3238	46	1	342	1.5	5.2	1.5 10−2	−1.37	−0.81	−4.2 to −0.9
Kepler	109408	0.2	0.00068	0.062	0.00067	0.00068	1.2 10−5	−7.78	−7.78	−1.9 to +1.2
OVRO	529	8.7	1	87	3.2	6.1	5.1 10−4	−2.48	−2.20	−3.5 to −1.5
UMRAO1	1364	33	2	468	6.0	10.2	8.8 10−4	−1.97	−1.74	−4.1 to −1.5
UMRAO2	1300	41	2	218	7.0	11.6	9.0 10−4	−1.89	−1.67	−4.2 to −1.6
UMRAO3	978	38	2	185	7.9	12.4	1.3 10−3	−1.69	−1.49	−4.1 to −1.9

Note. Columns: (1) data set; subscripts 1, 2, and 3 (respectively) refer to the 14.5, 8.0, and 4.8 GHz observing frequencies for the UMRAO data sets; (2) number of data points in the observed light curve; (3) total duration of the light curve; (4) minimum sampling interval of the light curve; (5) maximum sampling interval of the light curve; (6) median sampling interval of the light curve; (7) mean sampling interval of the observed light curve (duration of the monitoring divided by the number of data points); (8) mean variance of the light curve due to measurement uncertainties; (9) median noise floor level in the PSD due to measurement uncertainty (Equation (3)); (10) mean noise floor level in the PSD due to measurement uncertainty (Equation (3)); and (11) temporal frequency range covered in PSD analysis, above the median noise floor level.

Download table as:  ASCIITypeset image

Power spectral analysis of astrophysical sources typically invokes Fourier decomposition methods, where a source light curve is represented by a sum of a set of sinusoidal signals with random phases and amplitudes, which correspond to various timescales of a source's variability in the time series (e.g., Timmer & Koenig 1995). As such, the constructed PSD is a Fourier transform without the phase information. However, PSDs generated using Fourier decomposition methods can be distorted owing to aliasing and red-noise leak. Aliasing arises from the discrete sampling of a time series, while the red-noise leak appears because of the finite length of a light curve. This problem is particularly severe if a time series is not evenly sampled, as the response of a spectral window (i.e., the discrete Fourier transform of the sampling times) is in such a case unknown in the Fourier domain (e.g., Deeming 1975). Therefore, in order to derive reliable PSDs, an evenly sampled time series has to be obtained through a linear interpolation from an unevenly sampled data. Even though this procedure introduces interpolated data in a time series, the underlying PSD parameters can then be recovered up to a typical (mean) sampling interval of an unevenly sampled time series (see the discussion in Goyal et al. 2017 and, in particular, the Appendix therein).

We emphasize that the most common approaches in the literature to deal with unevenly spaced data in order to obtain PSDs using Fourier transform methods, such as the Lomb–Scargle periodogram (LSP; Lomb 1976, Scargle 1982) and the Fourier transform (FT) of the autocorrelation function (ACF; Edelson & Krolik 1988), do not reproduce correctly the slope of the underlying power spectrum. We demonstrate this in Appendix A by simulating blazar light curves using the method of Emmanoulopoulos et al. (2013). We note here that for the LSP "least squares fitting" method (see in this context the discussion in VanderPlas 2018), a variability power is added at the highest frequencies—limited by the mean sampling interval—due to a small number of the available frequencies in the periodogram, or in other words a small number of degrees of freedom (DOF); in the periodograms derived using FT of the ACF, on the other hand, the effect of the spectral window function becomes progressively more prominent for progressively less evenly spaced data (see Appendix A for details).

Since our aim is to characterize the variability properties of OJ 287 over an extremely broad range of temporal frequencies (∼6 decades), using the 117 year long optical light curve, albeit with extremely uneven sampling, instead of standard Fourier decomposition methods here we use a certain statistical model to fit the light curve in the time domain, and thus to derive the source power spectrum.⁷⁹ Specifically, we use the publicly available Continuous-time Auto-Regressive Moving Average (CARMA) model⁸⁰ by Kelly et al. (2014), which is a generalized version of the first-order Continuous-time Auto-Regressive (CAR(1)) model (also known as an Ornstein-Uhlenbeck process). In the CAR(1) model, the source variability is essentially described as a damped random walk—that is, a stochastic process with an exponential covariance function $S({\rm{\Delta }}t)={\sigma }^{2}\,\exp (-| {\rm{\Delta }}t/\tau | )$ defined by the amplitude σ and the characteristic (relaxation) timescale τ (Kelly et al. 2009). Meanwhile, in the CARMA model, the measured time series y(t) is approximated as a process defined to be the solution to the stochastic differential equation

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{p}y(t)}{{{dt}}^{p}}+{\alpha }_{p-1}\displaystyle \frac{{d}^{p-1}y(t)}{{{dt}}^{p-1}}+....+{\alpha }_{0}y(t)\\ &&=\,{\beta }_{q}\displaystyle \frac{{d}^{q}\epsilon (t)}{{{dt}}^{q}}+{\beta }_{q-1}\displaystyle \frac{{d}^{q-1}\epsilon (t)}{{{dt}}^{q-1}}+....+\epsilon (t),\end{eqnarray} \tag{ 1 }$$

where (t) is the Gaussian (by assumption) "input" white noise with zero mean and variance σ², the parameters α₀...α_p−1 are the autoregressive coefficients, the parameters β₁...β_q are the moving average coefficients, and finally α_p = 1 and β₀ = 1. The case with p = 1 and q = 0 corresponds to the CAR(1) process (with the relaxation timescale τ = 2π/α₀); hence, a CARMA(p, q) model describes a higher-order process when compared with CAR(1).

For an in-depth discussion of the CARMA model, the reader is referred to Kelly et al. (2014). Here, we only note that in this approach, for a given light curve y(t), one derives the probability distribution of the (stationary) CARMA(p, q) process via Bayesian inference, and in this way one calculates the corresponding power spectrum

$$\begin{eqnarray}&&P(f)={\sigma }^{2}{\left|\displaystyle \sum _{j=0}^{q}{\beta }_{j}{(2\pi {if})}^{j}\right|}^{2}{\left|\displaystyle \sum _{k=0}^{p}{\alpha }_{k}{(2\pi {if})}^{k}\right|}^{-2},\end{eqnarray} \tag{ 2 }$$

along with the uncertainties. Kelly et al. (2014) provide the adaptive Metropolis MCMC sampler routine to obtain the maximum-likelihood estimates. The quality of the fit is assessed by standardized residuals: if the Gaussian CARMA model is correct, the residuals should form a Gaussian white-noise sequence, for which the ACF is normally distributed with mean zero and variance 1/N, where N is the number of data points in the measured time series. Importantly, since here a light curve is directly modeled in the time domain—unlike in a PSD analysis using the standard FT methods—the resulting PSD shape should, in principle, be free from "uneven sampling" effects.

Here, we employ the generalized "corrected" Akaike information criterion (AICc; Hurvich & Tsai 1989), which is based on the maximum-likelihood estimate of the model parameters, including penalizing against overfitting due to the model complexity for finite sample sizes, to chose the order of the CARMA(p, q) process, although other criteria could be applied in this context as well. As such, models having AICc values <2 can be considered as sufficiently close to the null hypothesis, while models having AICc values >10 are not supported (Burnham & Anderson 2004). The CARMA software package from Kelly et al. (2014) that we employed finds the maximum-likelihood estimates of the model parameters by running 100 optimizers with random initial sets of model parameters, and then selects the order (p, q) that minimizes the AICc.

Finally, we note that the noise floor level in the derived PSD (Equation (2)) resulting from statistical fluctuations caused by measurement errors, is calculated as

$$\begin{eqnarray}&&{P}_{\mathrm{stat}}=2\,{\rm{\Delta }}t\,{\sigma }_{\mathrm{stat}}^{2},\end{eqnarray} \tag{ 3 }$$

where Δt is the sampling interval and ${\sigma }_{\mathrm{stat}}^{2}={\sum }_{j=1}^{j=N}{\rm{\Delta }}y{({t}_{j})}^{2}/N$ is the mean variance of the measurement uncertainties in the flux values y(t_j) in the observed light curve at times t_j.

The flux distributions of blazar light curves can be modeled nonlinearly, in the sense that they often can be represented as $y(t)=\exp [l(t)]$, where l(t) is a linear Gaussian time series (see, e.g., Edelson et al. 2013; Kushwaha et al. 2016; Abdalla et al. 2017; Liodakis et al. 2017). Hence, we have logarithmically transformed the light curves of OJ 287 analyzed here (Figures 1–3), and then modeled them as Gaussian CARMA(p, q) processes. For each light curve, the minimum (p, q) order was selected by minimizing the AICc values on the grid p = 1, .., 7 and q = 0, .., p − 1. As such, we ran the MCMC sampler for ${ \mathcal N }$ iterations with the first ${ \mathcal N }/2$ iterations discarded as a burn-in. Next, we employed the Gelman & Rubin (1992) method as a diagnostic to analyze the chain convergence using the multiple-chain approach, allowing us to compare the "within-chain" and "between-chain" variances. The number of iterations was chosen to derive the potential scale reduction factor for all of the model parameters to be <1.001. We select as the best-fit model the one produced by the pair of (p, q) values having the lowest order within the range in which models are statistically indistinguishable from each other (i.e., minimum AICc <10; see Section 3). Figure 15 in Appendix B shows, for comparison, the power spectra obtained for different (p, q) parameters consistent with the null hypothesis of the CARMA process for the analyzed Fermi-LAT light curve (see Figure 3(a)).

The results of the CARMA model fitting are presented in Figures 5 (high-energy γ-rays), 6 (X-rays), 7–9 (optical), and 10–13 (radio). We plot the measured time series along with the modeled values based on the best-fit CARMA process (panel (a)), the standardized residuals and their distribution compared with the expected normal distributions (panel (b)), the corresponding ACFs compared with the 95% confidence regions for a white-noise process (panel (c)), the squared ACFs compared with the 95% confidence regions for a white-noise process (panel (d)), the AICc values for different (p, q) pairs (panel (e)), and the resulting PSDs with 2σ confidence regions, as well as noise floor levels P_stat marked by horizontal lines (panels (f)). As shown, all of the analyzed light curves are well represented by Gaussian CARMA processes, as the residuals from the model fitting follow the expected normal distributions with the ACFs and the squared ACFs lying within 2σ intervals for most of the temporal lags. Note that because some of the light curves are sparsely sampled (in particular, the historical optical and the Swift-XRT light curves), we estimated the noise floor level (Equation (3)) with either "mean" or "median" sampling intervals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 presents the PSDs for the entire available long-term (1900–2017) monitoring optical data set, with typical daily sampling (hereafter "optical(all)"). The majority of the optical data obtained before 1970 (∼7% of the data points) have measurement uncertainties of the order of 20%, owing to large calibration errors resulting from observations recorded on photographic plates (see Hudec et al. 2013 for a discussion of error estimation). The overall noise floor level in the derived long-term PSD is relatively high, first because of larger measurement uncertainties, and second because the mean sampling interval is large, >12 days. Therefore, we also derived the long-term PSD for the data obtained using only the good-quality photomultiplier tubes and CCD photometric measurements for the period 1970–2017, with typical measurement uncertainties ∼2%–5% and a typical sampling of 5 days (hereafter "optical(trun.)"), as shown in Figure 8. Finally, the PSD corresponding to the continuous 72 day long monitoring Kepler data, with sampling down to ∼1 minute (see Table 2), is presented in Figure 9. The optical PSD of the blazar obtained by combining the PSDs generated with optical(all), optical(trun.), and the Kepler data covers an unprecedented frequency range of nearly 6 dex (from 117 years down to hour-long timescales), without any gaps. The lower-frequency segment of the Kepler PSD is found to be consistent with a simple extrapolation of the optical PSD following from the long-term monitoring of the source.

Figure 14 presents the composite multiwavelength PSDs of OJ 287, truncated below the median noise floor level for each data set. As seen, the PSDs derived using GHz-band radio light curves follow each other closely, on timescales ranging from decades to weeks/months. Also, there is a remarkable similarity between the radio and optical PSDs on variability timescales longer than ∼1 year. On the other hand, on timescales shorter than one year, the X-ray and radio power spectra seem to resemble each other, with the variability amplitudes smaller than those observed in the optical. Importantly, the overall shape of the γ-ray PSD is clearly distinct from the (generally speaking, colored-noise type) PSDs derived at lower-frequency bands, at least over timescales probed by the length of the LAT light curve, down to the median LAT sampling interval. In particular, in the high-energy γ-ray range we see rather uncorrelated flux changes on timescales longer than $\simeq {157}_{-91}^{+350}\,\mathrm{days}$, manifesting as a flat ("white-noise") segment in the derived variability power spectrum. We note in this context that some low-frequency flattening can be noticed in basically all of the PSDs modeled with CARMA; this is due to the fact that, by definition (Equation (1)), the model applied includes the uncorrelated Gaussian ("driving") term (t). However, only in the case of the γ-ray power spectrum could such a flattening be considered significant (note the 2σ confidence regions marked as blue areas in Figures 5–13, especially in the X-ray power spectrum shown in Figure 13).

Zoom In Zoom Out Reset image size

The main findings from our analysis of the multiwavelength and multi-epoch measurements of OJ 287 can be summarized as follows.

• (i)  
  On timescales ranging from tens of years down to months, the power spectra derived at radio frequencies based on single-dish observations indicate a colored-noise character of the source variability.
• (ii)  
  Owing to the inclusion of the Kepler data, we were able to construct for the first time the optical variability power spectrum of a blazar without any gaps across ∼6 dex in temporal frequencies. The modeled power spectrum is well represented by a higher-order CARMA process, on timescales ranging from 117 years down to <1 hr, but with no significant narrow features that could be identified with QPOs (see the next section below).
• (iii)  
  The power spectrum at X-ray photon energies, based on relatively sparsely sampled Swift-XRT data, could be characterized essentially only on timescales shorter than a year (down to months), where it resembles the radio power spectrum of the target.
• (iv)  
  The high-energy γ-ray power spectrum of OJ 287, modeled using the Fermi-LAT data, is noticeably different from the radio, optical, and X-ray power spectra of the source. In particular, in the framework of the CARMA modeling, this power spectrum is consistent with the CAR(1) process, according to the minimum AICc criterion adopted in this study, with an uncorrelated Gaussian (white) noise above the relaxation timescale of ∼150 days and a correlated (colored) noise on timescales shorter than 150 days.

5.1. Periodicity in the Long-term Optical and Fermi-LAT Light Curves

5.2. Characteristic Variability Timescales

5.3. Multiwavelength Power Spectra

We thank the anonymous referees for constructive comments on the manuscript. The authors also thank B. C. Kelly and J. Ballet for useful discussions during the revision of the paper.

A.G. and M.O. acknowledge support from the Polish National Science Centre (NCN) through the grant 2012/04/A/ST9/00083. D.K.W. and A.G. acknowledge the support from 2013/09/B/ST9/00026. Ł.S. and V.M. are supported by Polish NSC grant UMO-2016/22/E/ST9/00061. M.S. acknowledges the support of 2012/07/B/ST9/04404. S.Z. acknowledges the support of 2013/09/B/ST9/00599 and 2017/27/B/ST9/01855. R.H. acknowledges GA CR grant 13-33324S. A.S. and M.So. were supported by National Aeronautics and Space Administration (NASA) contract NAS8-03060 (Chandra X-ray Center). M.So. also acknowledges Polish NCN grant OPUS 2014/13/B/ST9/00570. A.V.F. is grateful for support from National Science Foundation (NSF) grant AST-1211916, NASA grant NNX12AF12G, the TABASGO Foundation, the Christopher R. Redlich Fund, and the Miller Institute for Basic Research in Science (U.C. Berkeley). Research at Lick Observatory is partially supported by a generous gift from Google.

This research has made use of data from the University of Michigan Radio Astronomy Observatory, which has been supported by the University of Michigan and by a series of grants from the NSF, most recently AST-0607523, and NASA Fermi grants NNX09AU16G, NNX10AP16G, and NNX11AO13G. The OVRO 40 m Telescope Fermi Blazar Monitoring Program is supported by NASA under awards NNX08AW31G and NNX11A043G, and by the NSF under awards AST-0808050 and AST-1109911. Based on observations obtained with telescopes of the University Observatory Jena, which is operated by the Astrophysical Institute of the Friedrich-Schiller University.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include NASA and the Department of Energy (DOE) in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan, and the KA Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE Contract DE-AC02-76SF00515.

In our test simulations, an artificial light curve is generated using the method of Emmanoulopoulos et al. (2013), assuming a pure red-noise (a = 2) power spectrum. The data points are evenly sampled with a sampling interval of 1 day, and the total duration of the simulated time series is 1000 days. We then kept 30% of the data selected at random times to mimic an unevenly sampled data set. The simulated light curves are presented in Figure 15(a). Figure 15(b) shows the corresponding power spectra derived using the LSP method and the FT of the ACF, while Figure 15(c) shows the response of the spectral window for the FT of the ACF. We refer the reader to Scargle (1982) and Edelson & Krolik (1988) or the Appendix of Goyal et al. (2017) for the mathematical details of these methods. The derived power spectrum is fitted with simple power-law form $P({\nu }_{k})\propto {\nu }_{k}^{-a}$ where a is the slope of the power spectrum over the frequencies corresponding to 1000 days down to 2 days. As shown, the derived PSDs (slopes a = 0.3 and −0.1 for the LSP and the FT of the ACF, respectively) are very different from the true PSD (a = 2). The flattening of the derived power spectra in the high-variability frequency range in the case of the LSP is due to the DOF problem, which introduces artificial power in the high-frequency range of the spectrum. Note that this is not caused by the addition of the Gaussian noise (white noise; a = 0), as the simulated light curve is free of such an effect. In the case of the FT of the ACF, which is free of such a "missing data points" problem, the flattening is instead caused by a substantial amount of variability power provided by the spectral window function. For FT of the ACF method, the resulting power spectrum is convolved with the spectral window function, leading to flatter power-law slopes. In principal, one can de-convolve the power spectrum knowing the spectral window function (Figure 15(c)). As shown, the spectral window function becomes more and more noisy toward higher frequencies and therefore the de-convolved power spectrum will be more and more uncertain toward higher and higher frequencies.

Zoom In Zoom Out Reset image size

The order of the CARMA(p, q) process is chosen based on how close the model is to the data using the minimum AICc criterion adopted in the present study. It has been argued that the models for various pairs of (p, q) values for which the minimum AICc is within 10 are not statistically indistinguishable from each other. Figure 16 shows power spectra corresponding to a few sets of (p, q) parameters of the analyzed Fermi-LAT light curve. As these overlap substantially, we choose the lowest-order model—CARMA(1, 0)—as the best-fitting model to describe the high-energy γ–ray variability in OJ 287.

Zoom In Zoom Out Reset image size

Here, we discuss in more detail the lack of the claimed ∼12 year quasi-periodicity in the analyzed optical light curve. To demonstrate the robustness of the CARMA modeling in detecting QPO features against the background of red noise in the power spectrum for a finite time series covering only a few periods, an artificial light curve is generated with three components: (1) a pure red-noise (β = 2) power spectrum using the method of Emmanoulopoulos et al. (2013), (2) a sinusoidal component with a 12 year period, and (3) a Gaussian white noise with mean 0 and standard deviation 0.1 representing measurement uncertainty. The data points are evenly sampled with a sampling interval of 1 day, and the total duration of the simulated time series is 14,600 days (40 years, corresponding to our relatively better sampled optical light curve starting from 1970; see Table 2). Next, we kept 20% of the data selected at random times to mimic an unevenly sampled data set. Figures 17(a)–(e) present the results of CARMA modeling on our simulated light curve, while Figure 17(f) shows the computed power spectrum. Since the simulated light curve covers only ∼3 periods, we do not detect a clear peak on the ∼12 year timescale against the background of a red-noise power spectrum and Gaussian white noise, in accordance with the discussion by Vaughan et al. (2016).

Zoom In Zoom Out Reset image size

Using the historical optical(all) light curve (1900–2017) analyzed here, we calculated the power spectrum using the DFT method (see Goyal et al. 2017 for mathematical details on computing PSD using the DFT method). Figure 18 shows the resulting power spectrum obtained using CARMA modeling and the DFT method. Within a 3σ confidence interval, the two methods give compatible results. A mild peak ∼12 years from the DFT method is within the 3σ confidence intervals estimated from the CARMA modeling.

Zoom In Zoom Out Reset image size