DETECTION OF POSSIBLE QUASI-PERIODIC OSCILLATIONS IN THE LONG-TERM OPTICAL LIGHT CURVE OF THE BL LAC OBJECT OJ 287

Quasi-periodic oscillations (QPOs) are routinely found in the X-ray light curves of Galactic binary systems (see, e.g., Remillard & McClintock 2006 for a review). In radio-quiet active galactic nuclei (AGNs), the first detection of short (∼1 hr) X-ray QPOs was reported by Gierliński et al. (2008) for the narrow-line Seyfert 1 galaxy REJ 1034+396, with subsequent discoveries (e.g., Alston et al. 2015). In the optical domain, 5.2 year periodicity has recently been found in the light curve of the radio-quiet quasar PG 1302–102 (Graham et al. 2015a), and some other analogous candidates have been identified (Graham et al. 2015b; see also Liu et al. 2015 for the case of the high-redshift radio-loud quasar PSO J334.2028+01.4075, and Zheng et al. 2016 for the radio-quiet quasar SDSS J0159+0105).

Short-timescale quasi-periodic oscillations in optical or X-ray blazar light curves—with characteristic timescales of the orders of hours or even minutes—have been claimed in several cases (e.g., Lachowicz et al. 2009 for PKS 2155–304; Rani et al. 2010 for S5 0716+714); these results have not been confirmed for the other periods in the same sources. On the other hand, analogous findings regarding quasi-periodicity in blazar light curves on longer timescales of months or years seem more robust, as summarized below.

At radio frequencies, harmonics with periods ranging from about one year up to several years or a decade have been reported for AO 0235 (Liu et al. 2006; see also Raiteri et al. 2001), PKS 1510–089 (Xie et al. 2008), NRAO 530 (An et al. 2013), and PKS 1156+295 (Wang et al. 2014, who confirmed the previous analyses by Hovatta et al. 2007, 2008). Also, King et al. (2013) reported persistent ∼150 day periodicity in the 15 GHz light curve of J1359+4011. Similar results have been reported in the optical domain, including Sandrinelli et al. (2014, 2016), who confirmed the ∼315 day period in the light curve of PKS 2155–304 claimed previously by Zhang et al. (2014), or Sandrinelli et al. (2015, 2016) who found marginally significant ($\sim 3\sigma$) periodicity in PKS 0537–441, OJ 287, 3C 279, PKS 1510–089 and PKS 2005–489, on timescales ranging from tens of days up to a few years (often in harmonic relations).

In the X-ray domain, Rani et al. (2009) found ∼17 day and ∼420 day quasi-periodicity in AO 0235+164 and 1ES 2321+419, respectively. In the high-energy γ-ray regime, Sandrinelli et al. (2014) first detected the ∼635 day period in the light curve of PKS 2155–304, most likely in a harmonic relation with the QPO found in the source at optical frequencies. More recently, Ackermann et al. (2015) reported a periodic behavior of the blazar PG 1553+113 with the characteristic timescale of 2.2 yr (similar to the one found at optical frequencies, but not in the radio domain).

BL Lac object OJ 287 (R.A. = 08^h54^m4887, decl = +20^d06^m3064 and z = 0.306) is to date one of the best studied blazars in history. Sparsely monitored for more than a century, and intensely observed during the last 40–50 years, it has become the most famous case of a blazar periodicity, with its double optical outbursts repeating every ∼12 year (Sillanpää et al. 1988; Hudec et al. 2013). Besides this, a number of authors have also claimed the presence of QPOs in the source at various timescales: a $\sim 40\mbox{--}50$ day quasi-periodicity for the source has been claimed at radio frequencies by Wu et al. (2006), and at optical frequencies by Pihajoki et al. (2013). Also Sandrinelli et al. (2016) reported a marginally significant ($\sim 3\sigma$) signal in the optical light curve of OJ 287 with a period of ∼435 days, in addition to a much less significant one of 203 days.

In this paper, we present our analysis and results of a search for quasi-periodical oscillations in the 9.2 year long, and exceptionally densely sampled, optical light curve of the blazar OJ 287, utilizing several widely used statistical methods, including the Lomb–Scargle periodogram (LSP), and weighted wavelet Z-transform (WWZ).

We began optical monitoring of the blazar OJ 287 in 2006 September at three sites: Mt. Suhora Observatory of the Pedagogical University in Kraków (SUH), the Astronomical Observatory of the Jagiellonian University (KRK) and the Astronomical Observatory of the University of Athens (GR).

Here we report the data taken until the end of the 2014/2015 observing season. Most of the observations were gathered by taking a series of about 10 frames each night from which nightly averages were calculated. We obtained 562 points at SUH, 164 at GR, and 123 at KRK. All three sites are equipped with CCD cameras and sets of wide-band filters (Bessell 1990). We took most of the data in the R filter, occasionally also measuring the target colors at SUH and KRK. After calibrating images for bias, dark and flat-field (usually taken on the dusk or dawn sky) with IRAF, we extracted magnitudes with the C-Munipack (Hroch 1998; Motl 2004) software. All scientific images were visually inspected. We used GSC 1400–222 ($R=13\buildrel{\rm{m}}\over{.} 74$) as the comparison star, for which the constancy was checked against the nearby star GSC 1400–444.

To fill some of the gaps in our observations, optical observations in the R filter from the Small and Moderate Aperture Research Telescope System (SMARTS) were also included; the details of data acquisition and reduction procedures are discussed in Bonning et al. (2012). To match our observations with those from SMARTS, a small offset of 0.02 mag, comparable to instrumental errors, was subtracted from the SMARTS data. Furthermore, the magnitudes were converted into flux units using the zero points for the filter R given in Bessell et al. (1998).

3.1. The Light Curve

The analyzed optical light curve of the source spans about 9.2 yr, from 2006 to 2015 (53997.590–57352.578 MJD), with a total of 1338 high-quality observations. The mean brightness and the standard deviation during the observation period are 14.47 and 0.48 mag, respectively. Excluding nine summer gaps (for a total of 799 days) due to the source visibility, the average data sampling rate is ∼1.92 day⁻¹. The data acquisition rate indicated by the cumulative number of observations and the locations of the summer gaps are shown in Figure 1. The observed nearly linear increase in the cumulative number of observations with time, except for a few shorter time intervals manifesting in the graph as flat plateaus, is clearly an indicator of the highly regularly sampled data collected for the whole observation epoch.

Zoom In Zoom Out Reset image size

As revealed by the resulting source light curve presented in Figure 2, the blazar exhibits considerable variability during the period: the amplitude of peak-to-peak oscillations (Heidt & Wagner 1996) is calculated to be 3.00 mag, while the fractional variability (Edelson et al. 2002; Vaughan et al. 2003) is 44.25 ± 0.04%.

Zoom In Zoom Out Reset image size

3.2. Periodicity Search

3.2.1. Lomb–Scargle Periodogram

The LSP, one of the most widely used methods in QPO analysis, modifies the periodogram such that the least-squares fitting of sine waves of the form ${X}_{f}(t)=A\cos \omega t+B\sin \omega t$ is minimized (Lomb 1976; Scargle 1982). The modified periodogram is expressed as

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{2}\left\{\displaystyle \frac{{\left[{\displaystyle \sum }_{i}{x}_{i}\cos \omega ({t}_{i}-\tau )\right]}^{2}}{{\displaystyle \sum }_{i}{\cos }^{2}\omega ({t}_{i}-\tau )}+\displaystyle \frac{{\left[{\displaystyle \sum }_{i}{x}_{i}\sin \omega ({t}_{i}-\tau )\right]}^{2}}{{\displaystyle \sum }_{i}{\sin }^{2}\omega ({t}_{i}-\tau )}\right\},\end{eqnarray} \tag{ 1 }$$

where τ is given by

$$\begin{eqnarray}&&\tan (2\omega \tau )=\displaystyle \frac{{\displaystyle \sum }_{i}\sin \omega {t}_{i}}{{\displaystyle \sum }_{i}\cos \omega {t}_{i}}.\end{eqnarray} \tag{ 2 }$$

This weights the data points and hence accounts for the problem with unevenly spaced data set. The calculated LSP of the source light curve in the R band is presented in Figure 3, which clearly shows two prominent features around the periods of 410 ± 38 and 789 ± 58 days. The associated uncertainties are the half-widths at half-maxima (HWHMs) of Gaussian fits to the corresponding periodogram peaks after subtracting the mean periodogram due the colored noise determined from Monte Carlo (MC) simulations, as explained in detail in Section 4.

Zoom In Zoom Out Reset image size

3.2.2. Weighted Wavelet Z-transform

Although the LSP method accounts for irregular spacing in a time series, the method does not take into account time fluctuations in the periodic signal information. In real astronomical systems, QPOs may develop and evolve in frequency and amplitude over time. In such cases, the wavelet transform method proves to be a more useful tool, and indeed it is frequently employed in the analysis of blazar sources (e.g., Hovatta et al. 2008; Bhatta et al. 2013).

The method attempts to fit sinusoidal waves localized in both time and frequency domains, which can be scaled in frequency and shifted in time (Torrence & Compo 1998). However, for unevenly spaced data the method still may not yield reliable results due to the fluctuations of the local number density of the data points, an important parameter in the transform. As a result, a false time evolution of the characteristic frequency may be recovered. In this context, Foster (1996) suggested that the problem can be handled efficiently by rescaling of the wavelet functions such that they can be viewed as weighted projections on the trial functions ${\phi }_{1}(t)=1(t)$, ${\phi }_{2}=\cos [\omega (t-\tau )]$, and ${\phi }_{2}=\sin [\omega (t-\tau )]$ by the weight given as ${w}_{i}={e}^{-c{\omega }^{2}{({t}_{i}-\tau )}^{2}}$, where c acts as a fine-tuning parameter usually chosen around 0.0125. This modified approach is known as the WWZ. With V_x and V_y as weighted variations of the data and the model function, respectively, the WWZ power can be written as

$$\begin{eqnarray}&&\mathrm{WWZ}=\displaystyle \frac{({N}_{\mathrm{eff}}-3){V}_{y}}{2({V}_{x}-{V}_{y})},\end{eqnarray} \tag{ 3 }$$

where N_eff is the effective number of the data points (for further details see Foster 1996).

We computed the WWZ power of the R-band flux of OJ 287 as a function of time and characteristic period, with the conventional value c = 0.0125, using the WWZ software available at the American Association of Variable Stars Observers (AAVSO).⁶ The top panel in Figure 4 presents the color-scaled WWZ power of the source in the time-period plane. The figure indicates that the power for the characteristic periods centered at 419 ± 52 days and 817 ± 130 days is highly significant, confirming the result obtained using the LSP. The uncertainties given here for the periods are evaluated as the means of the HWHMs of the Gaussian fits centered around a maximum power along the time (τ) axis; these uncertainties are relatively large because the maximum powers "meander" along τ.

Zoom In Zoom Out Reset image size

The QPO corresponding to the period of 817 ± 130 days starts from the beginning of the observing epoch and gradually dies out toward the end of the observations, although we note that the high power at the beginning could have been affected to some extent by edge effects. On the other hand, a QPO behavior around the period of 419 ± 52 days shows up a little further from the start of the campaign, and decays well before the end of the observations. The bottom panel of Figure 4 presents the average oscillation power at a given period. As shown, while in the lower frequency domain the signal seems weaker and transitory, at higher frequencies it becomes more coherent and persistent.

Although the above-described LSP and WWZ methods are widely used in time series analysis in astrophysics, the effects related to uneven sampling of a light curve—which can sometimes produce spurious peaks in the periodogram that can be mistaken for a real signal—should be analyzed in detail, especially when searching for QPO oscillations. In addition, there is a pitfall that blazar variability is, in general, of a colored noise type, with larger amplitude fluctuations at longer variability timescales; as a result, high-amplitude peaks in the lower-frequency region of the periodograms can sometimes appear to be QPO features (see in this context Press 1978). Hence, in order to establish a robust significance of the potential QPO oscillations revealed by both the LSP and WWZ methods in OJ 287, both the uneven sampling of the source light curve and the colored noise-type behavior of the blazar have to be addressed carefully. This can be done by means of MC simulations of a large number of red-noise light curves of the source, obtained by randomizing both the amplitude and the phase of the Fourier components, following Timmer & Koenig (1995).

First, in order to characterize the exact form of the underlying colored noise of the source (needed as an input for simulating the source light curves), we followed the power response method (PSRESP; see Uttley et al. 2002), which is widely used in modeling AGN power spectra (e.g., Chatterjee et al. 2008; Edelson et al. 2014; Bhatta et al. 2016; Kapanadze et al. 2016; Stone et al. 2016). In short, the method attempts to fit the binned periodogram with a model power spectral density (PSD) which maximizes the probability that the given model PSD can be accepted.

The discrete Fourier periodogram (DFP) of the light curve and the periodogran binned in logarithm frequency to reduce the scatter, are both presented in the left panel of Figure 5; the horizontal orange dashed line marks the Poisson noise level. We attempted to model the periodogram with a broken power-law of the form $P{(f)\propto {[1+(f/{f}_{b})}^{2}]}^{-\beta /2}$, where β and f_b correspond to the spectral slope and the break frequency, respectively. This particular "knee" model is designed to suppress the excess power at long variability timescales, given the fact that the power in the lower-frequency domain of the periodogram cannot increase indefinitely, in accord with the analysis of the long-term optical light curve of PKS 2155+304 presented by Kastendieck et al. (2011), who found the characteristic break timescale in the blazar of the order of a few years.

Zoom In Zoom Out Reset image size

In order to determine the best-fit PSD model, the observed binned DFP (blue triangles in the left panel of Figure 5) was compared with the mean periodogram of 1000 simulated light curves for each pair of trial values of β and f_b. In particular, as a measure of the goodness of fit, we calculated the probabilities of a model being accepted for 10 trial β ranging from 0.8 to 1.7 with an interval of 0.2, and 7 trial f_b ranging from 900 to 1500 days with an interval of 100 days (a total of $10\times 7\times 1000=70\mathrm{,000}$ simulated light curves). The resulting distribution of the probabilities for the corresponding trial break periods and PSD slopes are presented in the the color-coded contour plot in the right panel of Figure 5. The resulting best-fit model (the highest probability of 0.73), corresponding to $1/{f}_{b}=1400\pm 55$ days and $\beta =1.3\pm 0.2$, is represented by the magenta points connected by the dashed line in the left panel of Figure 5. The uncertainties associated with the slope and the break period were calculated here by taking half the vertical and horizontal ranges of the highest probability ($\gt 0.7$) region denoted in the right panel of Figure 5 in red.

The best-fit PSD model described above was next used to simulate numerous light curves of the source, which were re-sampled according to the observed light curve. The distribution of 1000 such light curves was studied to establish the significance of the peaks seen in the observed LSP; in particular the 99% significance contour was determined by the 99th percentile of a distribution of the power in the periodogram of the simulated light curves. As shown in Figure 3 (top-left panel), this significance turns out to be 99.8% and 99.3% for the peak around 410 and 789 days, respectively. In order to validate further the QPO detection, sinusoidal oscillations with 419 day and 780 day periods and a small amplitude ($\sim 0.7{\sigma }_{\mathrm{flux}}$) were superimposed on the previously generated light curves. The corresponding average periodograms resulting from the $1000$ simulations are shown as green curves in Figure 3: 410 day period only in the top-right panel, 780 day period only in the bottom-left panel, and both periods together in the bottom-right panel.

Similarly, in order to establish the significance of the observed WWZ power of the light curve, WWZ analysis was carried out on the 1000 simulated light curves from the above-described reference PSD model, and their distribution of WWZ power averaged over τ (time; day) in the considered frequency-space (alternatively period-space) were studied for the significance estimates. In particular, the observed average WWZ power was compared against the 90% and 99% confidence contour (90th and 99th percentile at a given frequency) from the distribution of the τ-averaged WWZ power for the simulated light curves represented in the lower panel of Figure 4 by the orange and red dashed-curves, respectively. The significance of the of the peaks centered around ∼400 and ∼800 days were in this way evaluated as 99.7% and 99.0%, respectively.

Periodicity in the light curves of blazar sources (and, more generally, of other types of AGN), may be related to the period of a perturbing object in a binary black hole system (eg Lehto & Valtonen 1996; Graham et al. 2015a), or a jet precession caused by either closely orbiting binary black holes or warped accretion disks (see the discussion in, e.g., Graham et al. 2015b; Sandrinelli et al. 2016). Moreover, coherent helical/non-ballistic motions of relativistic blobs within blazar jets could be responsible for quasi-periodicity of blazar sources (e.g., Camenzind & Krockenberger 1992; Mohan & Mangalam 2015). Finally, jet modulation by various instabilities developing within the innermost parts of accretion disks could, in principle, result in quasi-periodic modulation of the observed jet emission (see the discussion in Liu et al. 2006; Wang et al. 2014).

In the specific context of OJ 287, previous claims of the detection of (quasi-)periodic oscillations on year-like timescales—besides the 60 year modulation advocated by Valtonen et al. (2006) and the 12 year cycle discussed by Sillanpää et al. (1988)—have been made by a number of authors. For example, Hughes et al. (1998) identified a persistent oscillation in the radio light curve of the source, with a characteristic timescale of about 410 days. This result was supported by the wavelet analysis attempted by Hovatta et al. (2008). More recently, using the data for the partly overlapping observing epoch considered in this paper, Sandrinelli et al. (2016) claimed to have detected $\simeq 435\,\mathrm{day}$ and $\simeq 412\,\mathrm{day}$ period in the NIR-optical and γ-ray light curve of the blazar, respectively. The detection of the same signal in more than one band—radio, optical, and γ-rays—enforces the physical relevance of the findings. It is therefore interesting to indicate at this point a possible analogy between the year-like periodicity we see in OJ 287 and the low-frequency "C-type" QPOs commonly found in Galactic X-ray binaries. These QPOs come in harmonic peaks sticking out of a power spectrum at frequencies just above/around the transition from the "white noise" to the "pink/red noise" segments of the PSDs, and are typically linked to the Lense–Thirring precession of the innermost parts of accretion disks (see, e.g., Stella & Vietri 1998; Motta et al. 2011).

Considering a possible explanation for the observed QPOs in OJ 287, we note that since binary black holes have been claimed to shape the 12 year periodicity of the optical light curve in the source (Valtonen et al. 2008; but see also Villforth et al. 2010 for a critical review), the same scenario could not account for the ∼400 day period found in our analysis. Jet precession due to a warped accretion disk, on the other hand, as well as a "grand-design" helical magnetic field in the OJ 287 jet, seem both in conflict with the erratic jet wobbling observed on parsec scales by Agudo et al. (2012) in the high-frequency radio domain at the time of our monitoring program. Hence we conclude that the likely explanation for the ∼400 day period in the blazar (with the possible accompanying ∼800 day harmonic) could be a jet modulation by the innermost parts of the accretion disk. Interestingly, in the case of magnetically arrested disks, the characteristic timescale of QPOs in a jet with production efficiency set by the rotating and unstable ("chocking") magnetic field accumulated at the saturation level around the horizon of a spinning black hole, as seen in recent MHD simulations (Tchekhovskoy et al. 2011; McKinney et al. 2012), corresponds to tens/hundreds of the gravitational radius light-crossing times, $\sim 100\,{r}_{g}/c\sim {10}^{6}\,({M}_{\mathrm{BH}}/{10}^{9}\,{M}_{\odot })$ s. This is of the correct order of magnitude assuming ${M}_{\mathrm{BH}}\simeq 2\times {10}^{10}\,{M}_{\odot }$ as suggested by Valtonen et al. (2008) for OJ 287. Even so, in such a case the year-like periodicity claimed for radio and optical light curves of other blazars (most likely hosting smaller black holes) would remain a puzzle, unless lower black hole spins are considered, allowing the accommodation of longer timescales for smaller black hole masses.

We acknowledge the support by the Polish National Science Centre through the grants: DEC-2012/04/A/ST9/00083 and 2013/09/B/ST9/00599. We are grateful to C. Porowski, T. Szymanski, D. Jableka, E. Kuligowska, J. Krzesinski, A. Baran and A. Kuzmicz for carrying out some of the observations. We gratefully acknowledge Marian Soida for his valuable assistance in the computational part of the work. We thank the anonymous referee whose comments and suggestions greatly improved the analysis. This paper has made use of up-to-date SMARTS optical/near-infrared light curves that are available online.⁷