Optical monitoring and intra-day variabilities of BL Lac Objects OJ 287

The observations are carried out by the 1.26-m National Astronomical Observatory-Guangzhou University Infrared/Optical Telescope (NGT) at the Xinglong Observatory of the National Astronomical Observatories, Chinese Academy of Sciences. Since 2014, a new optical system has been designed to split the optical beam into optical and near infrared channels. An optical TRIPOL5 instrument and a near infrared PSL camera are attached in two channels, respectively. The TRIPOL5 can observe in three optical bands simultaneously and the PSL camera can observe in near infrared J band. Unfortunately, the near infrared camera has not worked well up to now. The TRIPOL5 use three SBIG STT-8300M cameras with a CCD of 3326 × 2504 pixels and a field of view of 6.0' × 4.5'. The filters adopt standard SDSS g, r, i bands. For OJ 287, the exposure time is 300 s.

The data reduction is carried out by the automated photometry pipeline, which is named RAPP (Huang et al. 2020), which can be introduced as follows. First, the observed images are corrected by bias, flat and dark images. Second, the RAPP will automatically detect the position in each image, and match images according to the information of stars. All images will be used to create an overlaid image. We use the overlaid image to register the position of stars in each observation image, and the positions of stars in different images are obtained. Finally, we carry out the aperture photometry process (the same as APPHOT package in IRAF).

Given K, the number of comparison stars, for each of them (S_i , i = 1,2,...K), we calculate the ith target magnitude (m_i ): m_i = m_i|o + m_i|c – m_i|oc , here m_i|oc is the ith instrumental magnitude of comparison stars, m_i|o is the instrumental magnitude of target magnitude, m_i|c is the standard comparison star magnitude. Considering the whole comparison stars, the target magnitude (m) can be calculated as $m\,=\,\displaystyle \frac{\displaystyle \sum _{i=1}^{K}{m}_{i}}{K}$ with a standard error $\sigma \,=\,\sqrt{\displaystyle \frac{\Sigma {({m}_{i}-m)}^{2}}{K-1}}$. The standard stars of OJ 287 noted as '4, 10, 11', are collected by Fiorucci & Tosti (1996). The comparison stars are listed in Table 1.

For every comparison star, based on VRI magnitudes, we use a least square fitting method to obtain the gri magnitudes, m_ν = alog² ν + blogν + c, here, m_ν is the magnitude at ν-band (ν = V, R and I). The fitting pictures have been shown in Figure 1, in which the black filled circle stand for VRI magnitudes, the red dots stand for gri (g, r, i) magnitudes, and the red lines stand for the least square fitting.

Zoom In Zoom Out Reset image size

3.1.1. Methods

We use the following three methods to analyze the IDV results.

(1) The intra-day optical variabilities (IDVs) can be constrained as the following method. For any two pairs of observations (t_j , m_j , σ_j ), (t_k , m_k ,σ_k ) (j, k = 1,2,...N), we calculated three parameters, time interval: Δ T_jk = |t_j – t_k |, magnitude difference: Δ m_jk = |m_j – m_k |, and the standard deviation: ${\sigma }_{jk}=\sqrt{{\sigma }_{j}{}^{2}+{\sigma }_{k}{}^{2}}$. If Δ m_jk > 3σ_jk , we take Δ m_jk as a real variation and the corresponding time interval Δ T_jk as the timescale. If there are more cases with Δ m_jk > 3σ_jk , we take the shortest Δ T_jk as the timescales as we did in our former papers (Fan et al. 2009a,b,c,2014).

(2) The variability amplitude parameter (Heidt & Wagner 1996) (A_m ),

$$\begin{eqnarray}&&{A}_{m}=\sqrt{{({m}_{\max }-{m}_{\min })}^{2}-{\sigma }_{\max }{}^{2}-{\sigma }_{\min }{}^{2}},\end{eqnarray}$$

here, m_max and m_min are the maximum and minimum magnitudes, and σ_max and σ_min are the corresponding uncertainties. When A_m > 7.5%, the source is variable, and Δ T is the time scale.

(3) The variability parameter (C_j ), see Romero et al. (1999).

3.1.2. Results

We use the upper three methods to study the intra-day lightcurves, and use the Gaussian function to fit the regions of intra-day lightcurves when showing obvious variabilities. The intra-day variabilities are listed in Table 3. On one day, if there are more than one regions with intra-day variabilities, we list them in Table 4, and noted as, (1:GF), (2:GF), ...

Table 3. The IDV Analysis of OJ 287

Band	Date	ΔT	Δm ± σ	Am	C	Date	ΔT	Δm ± σ	Am	C
(1)	(2)	(3)	(4)	(5)	(6)	(2)	(3)	(4)	(5)	(6)
g	2016/3/26	95.90	0.11 ± 0.01	9.50	2.16	2017/12/10	183.02	0.22 ± 0.06	20.75	2.36
	(1:GF)	95.91	0.11 ± 0.02			(1:GF)	44.08	0.09 ± 0.02
	(2:GF)	67.54	0.07 ± 0.02			(2:GF)	113.21	0.14 ± 0.03
	2016/12/18	248.11	0.08 ± 0.02	7.05	0.95	(3:GF)	64.77	0.17 ± 0.04
	(1:GF)	81.07	0.06 ± 0.02			2017/12/11	142.13	0.14 ± 0.05	13.27	1.65
	(2:GF)	49.52	0.06 ± 0.01			(1:GF)	7.69	0.06 ± 0.02
	2017/12/3	94.46	0.15 ± 0.03	13.45	1.55	(2:GF)	150.22	0.08 ± 0.02
	(1:GF)	52.63	0.11 ± 0.03			(3:GF)	48.77	0.08 ± 0.00
	(2:GF)	74.30	0.11 ± 0.03			(4:GF)	77.36	0.11 ± 0.01
	2017/12/4	129.13	0.16 ± 0.03	14.29	1.71	2019/2/24	344.16	0.16 ± 0.04	15.49	1.27
	2017/12/9	135.94	0.15 ± 0.04	14.10	1.29	(1:GF)	101.23	0.07 ± 0.03
	(1:GF)	28.54	0.06 ± 0.04			(2:GF)	65.30	0.05 ± 0.02
	(2:GF)	99.00	0.11 ± 0.03			(3:GF)	55.92	0.08 ± 0.03
	(3:GF)	49.84	0.09 ± 0.01			2019/3/16	33.70	0.08 ± 0.02	7.34	1.33
	(4:GF)	54.00	0.10 ± 0.02			(1:GF)	26.93	0.08 ± 0.02
						(2:GF)	62.04	0.05 ± 0.01
r	2017/11/26	148.46	0.13 ± 0.04	11.98	10.38	2017/12/11	155.52	0.10 ± 0.03	9.85	1.60
	2017/12/3	114.77	0.11 ± 0.02	10.20	1.89	2019/2/20	371.09	0.26 ± 0.04	12.93	3.28
	2017/12/4	210.24	0.15 ± 0.04	13.87	5.81	(1:GF)	20.17	0.20 ± 0.04
	2017/12/9	108.57	0.10 ± 0.03	9.56	1.56	(2:GF)	115.73	0.07 ± 0.03
	2017/12/10	182.88	0.11 ± 0.02	10.02	1.69	2019/2/24	13.54	0.15 ± 0.04	13.96	2.43
	(1:GF)	33.70	0.06 ± 0.02			2019/3/16	155.23	0.08 ± 0.03	7.88	1.15
	(2:GF)	40.81	0.11 ± 0.02			(1:GF)	27.07	0.07 ± 0.02
i	2017/12/4	236.28	0.15 ± 0.05	15.76	1.88	2017/12/11	243.36	0.13 ± 0.02	12.13	1.30
	2017/12/9	149.48	0.20 ± 0.05			(1:GF)	101.20	0.12 ± 0.03
	2017/12/10	203.18	0.20 ± 0.04	19.13	1.62	(2:GF)	20.35	0.11 ± 0.03
	(1:GF)	74.97	0.09 ± 0.03			2017/12/26	256.32	0.11 ± 0.05	51.78	1.59
	(2:GF)	47.23	0.12 ± 0.02			(1:GF)	115.20	0.10 ± 0.02

Column (1): Band; Col.(2): Date; Col.(3): the variable timescale, ΔT, in units of minutes; Col.(4): the variable value and the corresponding error, Δm ± σ, in units of mag; Col.(5): Am; Col.(6): C.

Table 4. Intra-day Periods of OJ 287

Date	Band	PS	JUR	DCF	Po
		(min)	(min)	(min)	(min)
(1)	(2)	(3)	(4)	(5)	(6)
2017/12/10	g	94.97 ± 14.82	92.69 ± 9.65	93.69 ± 19.48,	93.78 ± 3.66
				212.21 ± 27.16
	r	103.5 ± 32.07	92.69 ± 7.43,	88.2 ± 19.48,	94.80 ± 12.32
			185.52 ± 8.46	189.16 ± 21.4

Column (1): Date; Col.(2): Band; Col.(3): PS, power spectrum results, in units of min; Col.(4): JUR, Jurkevich results, in units of min; Col.(5): DCF, DCF results, in units of min; Col.(6): P^o , in units of min.

At g band, there are 9 days showing the IDVs, which are 2016/03/26, 2016/12/18, 2017/12/03, 2017/12/04, 2017/12/09, 2017/12/10, 2017/12/11, 2019/02/24 and 2019/03/16, see Figure 3 (the 1st, 2nd lines), with the variable timescales (Δ T) in the range from 7.69 min (Δ m = 0.06 ± 0.02 mag) to 344.16 min (Δ m = 0.16 ± 0.04 mag).

Zoom In Zoom Out Reset image size

At r band, there are 9 days showing the IDVs, which are 2017/11/26, 2017/12/03, 2017/12/04, 2017/12/09, 2017/12/10, 2017/12/11, 2019/02/20, 2019/02/24 and 2019/03/16 (the 3rd, 4th lines), see Figure 3, with the variable timescales (Δ T) in the range from 13.54 min (Δ m = 0.15 ± 0.04 mag) to 371.09 min (Δ m = 0.26 ± 0.04 mag).

At i band, there are 5 days showing the IDVs, which are 2017/12/04, 2017/12/09, 2017/12/10, 2017/12/11 and 2017/12/26 (the 5th lines), see Figure 3, with the variable timescales (Δ T) in the range from 20.35 min (Δ m = 0.11 ± 0.03 mag) to 243.36 min (Δ m = 0.13 ± 0.02 mag).

The relations between the gri variable timescales (ΔT) and variable values (Δm) are shown in Figure 4. At g band, Δ m = (3.05 ± 0.96) × 10⁻⁴Δ T + (0.07 ± 0.01), r = 0.51, p = 0.004, see the top-left panel in Figure 4; at r band, Δ m = (2.87 ± 0.91) × 10⁻⁴Δ T + (0.08 ± 0.01), r = 0.54, p = 0.04, see the top-right panel in Figure 4; at i band, Δ m = (1.42 ± 1.16) × 10⁻⁴Δ T + (0.10 ± 0.01), r = 0.40, p = 0.25, see the bottom-left panel in Figure 4. The correlation analysis imply that at g and r bands, Δ T and Δ m lie strong correlations; at i band, Δ T and Δ m lie no correlation.

Zoom In Zoom Out Reset image size

For all the data, Δ m = (2.91 ± 0.66) × 10⁻⁴Δ T + (0.08 ± 0.009), r = 0.52, p = 5.22 ± 10⁻⁵, see the bottom-right panel in Figure 4, which show that Δ m and Δ T have strong correlation.

For blazars, it is difficult to judge the periodicity because of the uneven distribution of observations. To solve this problem, we can have two solutions. (1) By selecting appropriate fitting and interpolation methods, the non-uniform observations can be transformed into uniform data distribution. (2) We choose the appropriate method to deal with the non-uniform observations.

So we use the following process to deal with the periodicity of blazars. (1) Choose the suitable method to analyze the light curves with uneven observations and obtain the possible period (P^o ). (2) Based on the regress analysis and interpolation to analyze the light curve, we obtain the uniform data distribution. (3) We obtain the periodicity (P^u ) of uniform data distributions. (4) Finally, we compare P^o with P^u , and if P^o is consistent with P^u , we choose P^o as the period, otherwise, P^o might be unreasonable.

3.2.1. Method

3.2.2. Results

1. Intra-day periods

We use the three methods to calculate every intra-day lightcurves, with the main results in Table 4.

For OJ 287, there are possible periods on 2017 December 10 at g and r band. We use the Locfit Regression (Feigelson & Babu 2012) to fit the intra-day light curves and obtain the uniform data distributions, which can be analyzed by the upper three methods, and obtain the period of uniform data distributions (P^u ). The upper results are shown in Figure 5.

Zoom In Zoom Out Reset image size

For the results calculated by the three different methods, we choose the common parts which are consistent with each other within the error range as the period. From Figure 5, we can see that P^o is consistent with P^u , which are noted by the red dotted lines, so we choose P^o as the period, at g band, P_g = 91.99 ± 5.66 min, at r band, P_r = 93.30 ± 3.15 min. The periods from two different bands are consistent with each other, which imply that the physical origins of g and r emissions are correlated.

2. Long-term periods

To analyze the long-term variability more comprehensively, we combine our observations with the available data from the literature (Pollock et al. 1971; Schaefer et al. 1980; Lloyd 1984; Sillanpää et al. 1988; Jia et al. 1995; Soltan et al. 1996; Jia et al. 1998; Xie et al. 1988; Bai et al. 1999; Ghosh et al. 2000; Wu et al. 2006; Zheng et al. 2008; Villforth et al. 2010; Dai et al. 2011; Pihajoki et al. 2013; Rakshit et al. 2017), and SMARTs (the Small and Moderate Aperture Research Telescope System).

We combine the whole available data in Figure 6 (the upper panel), in which, the black dots stand for the data from the literature, and the red dots stand for our observations.

Zoom In Zoom Out Reset image size

Based on power spectrum, the results are P_p = 0.48 ± 0.05, 1.91 ± 0.05, 12.02 ± 0.41 and 18.35 ± 0.48 yr; based on Jurkevich method, the results are P_J = 1.31 ± 0.27, 3.93 ± 0.13, 5.90 ± 0.53, 11.93 ± 0.4 and 22.96 ± 0.55 yr; based on DCF method, the results are P_D = 12.13 ± 1.82, 24.25 ± 2.28 and 37.04 ± 1.8 yr. We choose the common parts, and obtain the suspected period P^o = 12.02 ± 0.41 yr.

We use the Locfit Regression (Feigelson & Babu 2012) to fit the long-term light curves and show the results in Figure 6 (red line in the upper panel). After calculation, the possible period is P^u = 12.80 ± 0.95 r, which is consistent with P^o , so period of OJ 287 might be P = 12.02 ± 0.41 yr.

3.3.1. Methods

3.3.2. Results

After calculation, there are 4639 spectral indices (α), which are in the range from 0.39 ± 0.02 to 1.24 ± 0.02, with the averaged value $\bar{\alpha }=0.75\pm 0.11$.

At g band, α = –(1.91 ± 0.05) × 10⁻² F_g + (0.87 ± 0.003), with r = –0.50, p = 1.49 × 10⁻⁴, see Figure 7 (the left panel).

Zoom In Zoom Out Reset image size

At r band, α = –(1.67 ± 0.05) × 10⁻² F_r + (0.87 ± 0.003), with r = –0.47, p = 7.55 × 10⁻⁵, see Figure 7 (the middle panel).

At i band, α = –(1.25 ± 0.04) × 10⁻² F_i + (0.86 ± 0.003), with r = –0.44, p = 1.37 × 10⁻⁵, see Figure 7 (the right panel).

For OJ 287, Fan et al. (2009a) obtained the IDVs with the timescales in the range from 10 min to 2 hour with the variations from 0.11 mag to 0.75 mag. Gaur et al. (2012) found evident intra-day variabilities. In this work, we obtain that, at g band, the variable timescales are in the range of 7.69∼344.16 min; at r band, the variable timescales are in the range of 13.54∼210.24 min; at i band, the variable timescales are in the range of 6.62∼276.77 min.

Some works claim to find the intra-day periods, but have not been confirmed (e.g., Frohlich 1973; Visvanathan & Elliot 1973; De Diego & Kidger 1990). For example, Visvanathan & Elliot (1973) obtained a period of ∼40min at optical band for OJ 287.

In this work, for OJ 287, we obtain the consistent periods at g and i bands on 2017 September 10. OJ 287 is famous for the discovery of the convincing period about 12 yr (Sillanpää et al. 1988). Many models have been proposed to explain this period, which can be classed into two types: dynamical models and geometrical models (Liu & Wu 2002; Kushwaha 2020). The dynamical models show that the period is caused by the accretion dynamics in the supermassive binary black holes (SBBHs) (Sillanpää et al. 1988; Lehto & Valtonen 1996; Valtonen et al. 2008), while the geometrical models consider that the period caused by the Doppler boosted jet emission from the jet procession (Katz 1997; Villata et al. 1998; Britzen et al. 2018; Qian 2018).

Goyal et al. (2018) analyzed the long-term optical light-curves (∼117 yr), but did not find the period of 12 yr. In this work, based on our observations and the data from the available literature, we analyze the long-term lightcurves and obtain the period P = 12.02 ± 0.41 yr, which is consistent with that of Lehto & Valtonen (1996)'s.

For this source, we find that there are strong correlations between the IDV timescales (Δ T) and the corresponding variable values (Δ m). Based on the results from Fan et al. (2009b), we can obtain Δm = (–0.002 ± 0.001) Δ T + (0.44 ± 0.07), with r = –0.40, p = 0.11, which is different from the result in this work, seeing Figure 8.

Zoom In Zoom Out Reset image size

The origin of IDV timescales might come from the inner part of the blazar, and they are close to the black hole. In this sense, the IDV timescale can be used to estimate the mass (Abramowicz & Nobili 1982). For this source, we obtain the intra-day period P ≈ 95 min.

Based on the expression from Fan et al. (2014), the emission regions can be calculated as follows.

(1) $r=\displaystyle \frac{6GM}{{c}^{2}}$ (thin accretion disk with Schwarzschild black hole);

(2) $r=\displaystyle \frac{4GM}{{c}^{2}}$ (thick accretion disk with Schwarzschild black hole);

(3) $r=1.48\times {10}^{5}(1+\sqrt{1-{a}^{2}})\displaystyle \frac{M}{{M}_{\odot }}$ (Kerr black hole, with the radius of the event horizon r, angular momentum parameter a).

If we take the period, P as the time for the light travel in the innermost stable orbit, then 2 π r = c P/(1 + z), so the black hole mass should be

(1) $M=3.18\times {10}^{5}\displaystyle \frac{P}{1+z}{M}_{\odot }$ (thin accretion disk with Schwarzschild black hole);

(2) $M=4.77\times {10}^{5}\displaystyle \frac{P}{1+z}{M}_{\odot }$ (thick accretion disk with Schwarzschild black hole);

(3) $M=1.93\times {10}^{6}\displaystyle \frac{P}{1+z}{M}_{\odot }$ (Kerr black hole, with the radius of the event horizon r, angular momentum parameter a).

We use the upper method to obtain the lower limit of central black hole mass. For OJ 287, z = 0.3056, and P ≈ 95 min, so the corresponding masses of the black holes should be more than 2.31 ∼ 14.04 (× 10⁷ M_⊙). Cao (2003) pointed out the black hole mass of OJ 287 less than 6.16 × 10⁸ M_⊙, which is consistent with our result. Gupta et al. (2017) used the variable timescale to obtain the mass 6.5 ∼ 16.7(× 10⁷ M_⊙), which is very close to ours.

Considering the binary black hole system, with the semi-major axes being b₁ and b₂, the value of b₁ + b₂ can be calculated by the Kepler's law,

$$\begin{eqnarray}&&{P}^{2}=\displaystyle \frac{4{\pi }^{2}{({b}_{1}+{b}_{2})}^{3}}{G(M+m)},\end{eqnarray} \tag{ 1 }$$

which can be transformed as:

$$\begin{eqnarray}&&P \sim 1.72(1+z){M}_{8}{}^{-1/2}{r}_{16}^{3/2}{(1+\displaystyle \frac{{m}_{8}}{{M}_{8}})}^{-1/2}\,{\rm{yr}},\end{eqnarray} \tag{ 2 }$$

here, M₈ and m₈ are masses of the primary and secondary black hole, in units of 10⁸ M_⊙, P is period, G is the gravitational constant, r₁₆ = a₁ + a₂, in units of 10¹⁶ cm.

Based on the work of Sillanpää et al. (1988), r ∼ 0.1 pc. In this work, P = 12.02 yr, the secondary black hole mass m₈ = 0.231 ∼ 1.404. Based on Equation (2), we can obtain the primary black hole mass M₈ ≈ 74, that is, M = 7.4 × 10⁹ M_⊙, which is close to the value 1.8 × 10¹⁰ M_⊙ calculated from Valtonen et al. (2012). Sillanpää et al. (1988) give the ratio $(\displaystyle \frac{m}{M})$ is 0.004, our value is about 0.003 ∼ 0.019. The SMBH estimated by the IDV timescales gives a lower limit, so that mass is smaller than the mass estimated here.

For OJ 287, our results show the BWB behaviors at g, r and i bands. Dai et al. (2011) used the observations during an optical outburst to find a BWB chromatism. Bonning et al. (2012) found the BWB behaviors, as well as the redder-when-brighter (RWB) behaviors. Siejkowski & Wierzcholska (2017) did not reveal any BWB or RWB behaviors whenever the source was in any state of activity. Soltan et al. (1996) found that the color indices were stable during the outburst in 1993–1994.

To study the difference of the relations between F_ν and α when the source is in the different state of brightness, we divide the F_ν -α distributions into some sub-distributions according to the fixed flux span (F_ν|0), which can be calculated by the Jurkevich method from the F_ν -α distributions. Then, we calculate the correlations between F_ν and α from every sub-distribution.

For OJ 287, F_g|0 = 1.16 mJy, F_r|0 = 1.38 mJy, F_i|0 = 1.98 mJy, which have been placed in the top right corner of Figure 9. The correlations between F_ν and α from every sub-distribution have been noted in Figure 9, in which the red lines stand for the linear fitting. Most of the sub-distributions show BWB, especially when the source is in the state of brightness.

Zoom In Zoom Out Reset image size

We use the DCF method to analyze the time delays among the gri lightcurves and we plot the results in Figure 10. For the gross lightcurves, we cannot find time delay among gri bands, see Figure 10 (the first line), and find sub-structures, which are noted by 'a' and 'b', and might come from the different stages of the lightcurves.

Zoom In Zoom Out Reset image size

The lightcurves of OJ287 can be divided into six denser stages, which are noted as '1', '2', '3', '4', '5' and '6', see the green dotted rectangles in Figure 2. We analyze the time delays from different stages, and plot the results in the lower four lines of Figure 10. In some stages ('2', '3', '4' and '5'), there are time delays, but the delay time is different, for example, τ_ri|2 = 12.14 ± 1.88 min, τ_ri|3 = 8.03 ± 3.77 min, τ_ri|45 = –11.46 ± 6.66 min; in stage '6', no time delay.