Multicolor Optical Monitoring of the Quasar 3C 273 from 2005 to 2016

A number of statistical tests have been proposed to prove the validity of IDV reports. Normally, the C-test was often used to analyze the IDV. However, de Diego (2010) has found that the C-test is not a reliable methodology, because it does not have a Gaussian distribution and its criterion is too conservative. At present, the optical IDV is often analyzed by the F-test, χ²-test, modified C-test, and one-way analysis of variance (ANOVA; de Diego 2010; Joshi et al. 2011; Goyal et al. 2012, 2013; Hu et al. 2014; Agarwal & Gupta 2015; Dai et al. 2015). Romero et al. (1999) introduced the variability parameter, C, as the average value between C₁ and C₂:

$$\begin{eqnarray}&&{C}_{1}=\displaystyle \frac{\sigma (\mathrm{BL}-\mathrm{StarA})}{\sigma (\mathrm{StarA}-\mathrm{StarB})},{C}_{2}=\displaystyle \frac{\sigma (\mathrm{BL}-\mathrm{StarB})}{\sigma (\mathrm{StarA}-\mathrm{StarB})},\end{eqnarray} \tag{ 5 }$$

where (BL-StarA), (BL-StarB), and (StarA-StarB) are the differential instrumental magnitudes of the blazar and comparison star A, the blazar and comparison star B, and comparison stars A and B. σ is the standard deviation of the differential instrumental magnitudes. The adopted variability criterion requires C ≥ 2.576, which corresponds to a 99% confidence level. Despite the very common use of C-statistics, de Diego (2010) has pointed out that it has severe problems.

The F-test is thought to be a proper statistics to quantify variability (e.g., de Diego 2010; Joshi et al. 2011; Goyal et al. 2012; Hu et al. 2014; Agarwal & Gupta 2015; Xiong et al. 2016). The F value is calculated as

$$\begin{eqnarray}&&{F}_{1}=\displaystyle \frac{\mathrm{Var}(\mathrm{BL}-\mathrm{StarA})}{\mathrm{Var}(\mathrm{StarA}-\mathrm{StarB})},{F}_{2}=\displaystyle \frac{\mathrm{Var}(\mathrm{BL}-\mathrm{StarB})}{\mathrm{Var}(\mathrm{StarA}-\mathrm{StarB})},\end{eqnarray} \tag{ 6 }$$

where Var(BL-StarA), Var(BL-StarB), and Var(StarA-StarB) are the variances of the differential instrumental magnitudes. The F value from the average of F₁ and F₂ is compared with the critical F value, ${F}_{{\nu }_{{bl}},{\nu }_{* }}^{\alpha }$, where ν_bl and ν_* are the number of degrees of freedom for the blazar and comparison star, respectively ($\nu =N-1$), and α is the significance level set as 0.01 (2.6σ). If the average F value is larger than the critical value, the blazar is variable at a confidence level of 99%. Other alternatives to the standard F-test are ANOVA or the χ² test (e.g., de Diego 2010). de Diego (2010) showed that ANOVA is a powerful and robust estimator for microvariations. We use ANOVA in our analysis because it does not rely on error measurement but derives the expected variance from subsamples of the data. Considering the time of exposure, we bin the data in groups of three or five observations (see Xiong et al. 2016 and de Diego 2010 for detail). This method is used only for light curves with more than 20 observations on a given night; if the measurements in the last group are less than three or five, then it is merged with the previous group. The critical value of ANOVA can be obtained by ${F}_{{\nu }_{1},{\nu }_{2}}^{\alpha }$ in the F-statistics, where ${\nu }_{1}=k-1$ (k is the number of groups), ${\nu }_{2}=N-k$ (N is the number of measurements), and α is the significance level (Hu et al. 2014). When analyzing, we remove the outliers, where S_x is more than 3σ.

The blazar is considered variable (V) if the light curve satisfies the two criteria of the F-test and ANOVA. The blazar is considered probably variable (PV) if only one of the above two criteria is satisfied. The blazar is considered non-variable (N) if none of the criteria are met. We only analyze the data with exposure times of more than 1 minute because much shorter exposure times produce lower quality data and their use makes the analysis of the detection of IDVs much more difficult. The results of our analysis of the IDV are shown in Table 5. From Table 5, it can be seen that there is IDV found on five nights (I band on 2013 April 04, 2015 April 15, 2010 March 19, and 2009 January 17, and V band on 2016 April 25; Figure 1). As an illustration, the specific variabilities on two nights are presented. On 2013 April 04, the largest magnitude change of the I band is ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 055$ in 175 minutes from MJD = 56386.621 to MJD = 56386.742 corresponding to the variability amplitude Amp = 5.39%. On 2015 April 15, 3C 273 brightens by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 038$ in 52 minutes from the beginning of MJD = 57127.656 to MJD = 57127.692, and then fades by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 029$ in 41 minutes from MJD = 57127.692 to MJD = 57127.721. From Table 5 and Figure 1, we can see that the variability amplitudes on these nights are close to 5%. For the R band from 2014 April 25 and the I band from 2015 May 17, the light curves meet the criterion of the F-test, but do not reach the critical value of ANOVA. Though we use different bin sizes (N = 2–6), the two nights still do not reach the critical value of ANOVA. It is can be seen that the light curves of the two nights have large variation with Amp > 30% (Table 5 and Figure 1). The ANOVA test is used to compare the means of a number of samples. For the light curves of the two nights, when 3C 273 flares or darkens, there are only a few data points in the corresponding time interval, which reduces the variations between groups, and causes the ANOVA test to not detect the IDV (see de Diego 2010 for details). Therefore, according to results from the F-test, the two nights are detected as IDVs with a large variability amplitude. On 2014 April 25, 3C 273 brightens by ${\rm{\Delta }}R=0\buildrel{\rm{m}}\over{.} 27$ in 51 minutes from the beginning of MJD = 56772.668 to MJD = 56772.703, and then quickly fades by ${\rm{\Delta }}R=0\buildrel{\rm{m}}\over{.} 27$ in 26 minutes from MJD = 56772.703 to MJD = 56772.721. On 2015 May 17, 3C 273 brightens by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 27$ in 18 minutes from the beginning of MJD = 57159.651 to MJD = 57159.663, and then quickly fades by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 42$ in 5.8 minutes from MJD = 57159.663 to MJD = 57159.667 and again brightens by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 22$ in 5.9 minutes from the beginning of MJD = 57159.667 to MJD = 57159.671. The outliers (S_x > 3σ) only appear in two nights (2015 April 15 and 2016 April 25) and never exceed two data points for per night.

Zoom In Zoom Out Reset image size

To sum up, IDV is detected on seven nights. For the seven nights, we also check the color variations on intraday timescales. However, the results from three statistical tests do not show the corresponding color variations on intraday timescales. There are 14 nights detected as PV, excluding the seven nights. When estimating the variability amplitude (Amp), we only consider the night detected as "V" and "PV". The distributions of the variability amplitude from different bands are given in Figure 2, which shows that the variability amplitudes for most nights are less than 10%, and for four nights, more than 20% (also see Table 5). The correlations between the variability amplitudes and the source average brightness are shown in Figure 3. The results from the analysis of the Spearman rank indicate that there are no significant correlations between variability amplitude and brightness (significance level P > 0.05, N = 11, 7, 9 for the I, R, and V bands). However, there are not enough data to support the results. Making use of Equation (4), we calculate the DC of the IDV. The value of the DC is 10.84% for the V case and 53.65% for the PV+V cases. When considering the nights with time spans >4 hr, the value of the DC is 14.17% for the V case and 55.24% for the PV+V cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Long-term light curves and color index variations are given in Figure 4. From Figure 4, we can see that on the whole, the quasar becomes dark from 2005 to 2015 but also bright between different years; in 2015, the quasar reaches a low flux state; compared with 2015, the quasar in 2016 has a tendency to brighten; there is a color index variability for long timescales. Over the 12 years, the overall magnitude and color index variabilities are ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 67$, ${\rm{\Delta }}R=0\buildrel{\rm{m}}\over{.} 72$, ${\rm{\Delta }}V=0\buildrel{\rm{m}}\over{.} 68$, and ${\rm{\Delta }}(V-R)=0\buildrel{\rm{m}}\over{.} 25$, respectively. The magnitude distributions in the V, R, and I bands are 1310 < V < 1242, 1309 < R < 1237, and 1259 < I < 1193, respectively. The average values of the magnitude and color index are $\langle I\rangle =12\buildrel{\rm{m}}\over{.} 178\pm 0\buildrel{\rm{m}}\over{.} 118$, $\langle R\rangle =12\buildrel{\rm{m}}\over{.} 645\,\pm 0\buildrel{\rm{m}}\over{.} 107$, $\langle V\rangle \,=12\buildrel{\rm{m}}\over{.} 791\pm 0\buildrel{\rm{m}}\over{.} 106$, and $\langle V-R\rangle =0\buildrel{\rm{m}}\over{.} 126\pm 0\buildrel{\rm{m}}\over{.} 023$, respectively.

Zoom In Zoom Out Reset image size

Using the z-transformed discrete correlation function (ZDCF; Alexander 1997), we perform the interband correlation analysis and search for possible interband time delay. The ZDCF estimates the cross-correlation function (CCF) in the case of non-uniformly sampled light curves, and deals with sparsely and unequally sampled light curves better than both the interpolated cross-correlation function (ICCF) and the discrete correlation function (DCF; Edelson et al. 1996; Alexander 1997; Giveon et al. 1999; Roy et al. 2000). The ZDCF attempts to correct the biases that affect the original DCF by using equal-population binning. The ZDCF involves three steps (Alexander 1997; Oscoz et al. 2001). (i) All possible pairs of observations are sorted according to their time lag, and binned into equal-population bins of at least 11 pairs. Multiple occurrences of the same point in a bin are discarded so that each point appears only once per bin. (ii) Each bin is assigned its mean time lag and the intervals above and below the mean that contain 1σ of the points each. (iii) The correlation coefficients of the bins are calculated and z-transformed. The error is calculated in z-space and transformed back to r-space. The time lag corresponding to the maximum value of the ZDCF is assumed as the time delay between both components. The ZDCF code of Alexander (1997) can automatically set how many bins are given and used to calculate the interband correlation and the ACF. For each correlation, a Gaussian fitting is made to find the central ZDCF points. The time where the Gaussian profile peaks denotes the lag between the correlated light curves (Wu et al. 2012). The results that can determine the time delay are displayed in Figure 5. From Figure 5 and the results of the Gaussian fitting, the corresponding time lags are −0.002 ± 0.001, 0.002 ± 0.002, −0.007 ± 0.003, 0.003 ± 0.002, 0.011 ± 0.004, 0.014 ±0.003, 0.026 ± 0.006, and −0.005 ± 0.004, respectively. Considering the time resolutions, we do not find significant time lags between the V-band magnitude and I-band magnitude. For the rest of the nights, there are not enough data points or good Gaussian profiles to determine time lags.

Zoom In Zoom Out Reset image size

The autocorrelation function (ACF) is defined by

$$\begin{eqnarray}&&{ACF}(\tau )\equiv \langle (m(t)-\langle m\rangle )\cdot (m(t+\tau )-\langle m\rangle )\rangle ,\end{eqnarray} \tag{ 7 }$$

where the brackets denote a time average. The ACF measures the correlation of the light curve with itself, shifted in time, as a function of the time lag τ (Giveon et al. 1999). The zero-crossing time is the shortest time it takes the ACF to fall to zero (Alexander 1997). If there is an underlying signal in the light curve, with a typical timescale, then the width of the ACF peak near zero time lag will be proportional to this timescale (Giveon et al. 1999; Liu et al. 2008). The zero-crossing time is a well-defined quantity and used as a characteristic variability timescale (e.g., Netzer et al. 1996; Alexander 1997; Giveon et al. 1999; Liu et al. 2008; Liao et al. 2014). The width of the ACF may be related to the characteristic size scale of the corresponding emission region (Abdo et al. 2010a; Chatterjee et al. 2012). Another function used in variability studies to estimate the variability timescales is the first-order structure function (SF; e.g., Trevese et al. 1949) defined by

$$\begin{eqnarray}&&\mathrm{SF}(\tau )\equiv \sqrt{\langle {(m(t+\tau )-m(t))}^{2}\rangle }.\end{eqnarray} \tag{ 8 }$$

There is a simple relation between the ACF and the SF,

$$\begin{eqnarray}&&{\mathrm{SF}}^{2}(\tau )=2(V-\mathrm{ACF}(\tau )),\end{eqnarray} \tag{ 9 }$$

where V is the variance of the light curve (Giveon et al. 1999). We therefore perform only an ACF analysis on our light curves. The ACF was estimated by ZDCF. We only analyze the nights detected as IDV. The results used to estimate the characteristic variability timescales are given in Figure 6. Following Giveon et al. (1999), Liu et al. (2008), and Liao et al. (2014), we use a least-squares procedure to fit a fifth-order polynomial to the ACF, with the constraint that the ACF(τ = 0) = 1. The fitting results show that the detected variability timescales are 0.107, 0.023, 0.005, and 0.028 days for the I band on 2013 April 04, R band on 2014 April 25, and I bands on 2015 May 17 and 2015 April 15. In Section 3.1, we have given the timescales corresponding to the change of brightness on these nights (also see Figure 1). By comparison, we obtain that the variability timescales from the ACF analysis are consistent with the timescales corresponding to the change of brightness.

Zoom In Zoom Out Reset image size

In this section, we analyze the relationships between the magnitude and color index for intraday timescales, short-term timescales, and whole time-series data sets. The whole time-series data sets include data from 2005 to 2016, but exclude the data from 2006 to 2008 due to non-quasi-simultaneous measurements. For the color index, we correct for the Galactic extinction. The correction factors of the Galactic extinction are from Schlafly & Finkbeiner (2011). We concentrate on the V − R index and V magnitude because the V − R index versus V magnitude is frequently studied. When exploring the relationships, we only analyze the data with more than nine color indices obtained on intraday timescales and quasi-simultaneous measurements. The results of correlations between the V − R index and V magnitude are given in Table 6. As an example, Figure 7 shows the correlations between the V − R index and V magnitude on intraday timescales. The results from Table 6 show that most of the nights have strong correlations between the V − R index and V magnitude on intraday timescales, and the rest of the nights have moderate or weak correlations. Therefore, a bluer-when-brighter (BWB) chromatic trend is dominant for 3C 273 on intraday timescales. The BWB trend exists for short-term timescales and intermediate-term timescales but different timescales have different correlations between the V − R index and V magnitude (Table 6). The correlation between the V − R index and V magnitude in the whole time-series data sets is shown in Figure 8. The result of the correlation analysis from Table 6 shows that there is no BWB trend in the whole time-series data sets. We have the opportunity to explore the relationship between the BWB trend and the length of timescales in the V band because of long-term optical monitoring. We use the coefficient of correlation from error-weighted linear regression analysis to indicate the intensity of the BWB trend and the sum of the monitoring time spans on intraday timescales to indicate length of timescales. Analysis of the Spearman rank shows that there is significant anticorrelation between the BWB trend and length of timescales (r = −0.491, P = 0.002; see Figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to search for periodicity, well-covered long-term light curves are required. Soldi et al. (2008) presented an update of 3C 273's database hosted by the ISDC.¹² 3C 273's database, which was first published by Turler et al. (1999), is one of the most complete multiwavelength databases currently available for the quasar. The time spans of the optical V-band database are from 1968 to 2005. However, there are no V-band data from 1998 to 1999. We also compiled V-band data in 1998–1999 from Fan et al. (2014). Our data from 2005 to 2016 are updates and supplementary to the above data. For the data from Fan et al. (2014) and the present work, we convert magnitudes to fluxes in the Geneva photometric system as described in Turler et al. (1999), without any additional correction (Vcorr = 0) because the data from Fan et al. (2014) and the present work used standard Johnson photometric system while Soldi et al. (2008) and Turler et al. (1999) used the Geneva photometric system. For data from Soldi et al. (2008), we consider data of Flag = 0 in order to use only the best quality data and exclude the data with the influence of strong synchrotron flares (Flag = 1). Combining all data, we search the periodic signals of 3C 273 in 48 year time spans. The long-term light curves (bin = 1 day) are given in Figure 10. When there is more than one data on intraday timescales, we use the standard deviation as the error.

Zoom In Zoom Out Reset image size

For the periodicity analysis of AGNs, unevenly sampled data, frequency-dependent red noise, flares of high activity, and total monitoring time need to be considered (e.g., Schulz & Statteger 1997; Schulz & Mudelsee 2002; Vaughan 2005, 2010; Fan et al. 2014; Sandrinelli et al. 2016). For the data, the strong synchrotron flares are excluded. The time spans from the data are 48 years. Therefore, when analyzing periodicity, we mainly consider unevenly sampled data and red noise. In addition, in order to reduce false periodicity produced by short timescale variability, we pay more attention to periodicity in year timescales. In order to increase the reliability of the periodicity analysis, we use three methods (REDFIT: Schulz & Mudelsee 2002; Lomb–Scargle method: Lomb 1976 and Scargle 1982; Jurkevich method: Jurkevich 1971) to explore periodicity. The periodicity can be considered valid only if three methods have consistent results. Schulz & Mudelsee (2002) presented a computer program (REDFIT3.8e¹³ ) estimating red-noise spectra directly from unevenly spaced time series by fitting a first-order autoregressive (AR1) process. The program can be used to test if peaks in the spectrum of a time series are significant against the red-noise background from an AR1 process, and removes the bias of this Fourier transform for unevenly spaced data by correcting for the effect of correlation between Lomb–Scargle Fourier components. However, this program has two underlying assumptions: (i) the noise background recorded in a time series can indeed be approximated by an AR1 process; (ii) the distribution of data points along the time axis is not too clustered. For our data, the results of the non-parametric runs test indicate that the spectrum is consistent with an AR1 model (rtest = 139 falls inside the 98% acceptance region [129, 170]; see Schulz & Mudelsee 2002 and usage38e for further details). Figure 10 also shows that the distribution of data points along the time axis is not too clustered. So, it is appropriate to use REDFIT3.8e for periodicity analysis. The Lomb–Scargle method is commonly used to detect periodicity in unevenly sampled time series. The periodogram is a function of circular frequency ω, and is defined by the formula (Li et al. 2016)

$$\begin{eqnarray}{P}_{X}(\omega ) & \equiv & \displaystyle \frac{1}{2}\left\{\displaystyle \frac{{[{{\rm{\Sigma }}}_{i}X({t}_{i})\cos \omega ({t}_{i}-\tau )]}^{2}}{{{\rm{\Sigma }}}_{i}{\cos }^{2}\omega ({t}_{i}-\tau )}\right.\\ & & +\left.\displaystyle \frac{{[{{\rm{\Sigma }}}_{i}X({t}_{i})\sin \omega ({t}_{i}-\tau )]}^{2}}{{{\rm{\Sigma }}}_{i}{\sin }^{2}\omega ({t}_{i}-\tau )}\right\},\end{eqnarray} \tag{ 10 }$$

where X(t_i) (i = 0, 1..., N₀) is a time series. τ is calculated by the equation

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{2\omega }{\tan }^{-1}\left[\displaystyle \frac{{{\rm{\Sigma }}}_{i}\sin 2\omega {t}_{i}}{{{\rm{\Sigma }}}_{i}\cos 2\omega {t}_{i}}\right],\end{eqnarray} \tag{ 11 }$$

where ω = 2πν. Thus, the periodogram is a function of frequency ν. For a true signal X(t_i), the power in P_X(ω) would present a peak, or the power of a purely noise signal would be an exponential distribution. For a power level z, the False Alarm Probability (FAP) is calculated by (Scargle 1982; Press et al. 1994, p. 569; Li et al. 2016)

$$\begin{eqnarray}&&p(\gt z)\approx N\cdot \exp (-z),\end{eqnarray} \tag{ 12 }$$

where N is the number of data points. The Jurkevich method is an alternative and powerful method that is insensitive to the mean shape of the periodicity, less sensitive than Fourier analysis to uneven sampling. The Jurkevich method involves testing a run of trial periods around which the data is folded. The data is binned in phases around each trial period, and the total variance for all the points in the bin computed. A good period will give a much reduced variance (${V}_{{\rm{m}}}^{2}$). No firm rule exists for assessing the significance of a minimum in the ${V}_{{\rm{m}}}^{2}$. But a good guide is the fractional reduction of the variance $\left(f=\tfrac{1-{V}_{{\rm{m}}}^{2}}{{V}_{{\rm{m}}}^{2}}\right)$. Generally, a value f ≥ 0.5 (${V}_{{\rm{m}}}^{2}\leqslant 0.667$) implies that there is a very strong periodicity in the data (Jurkevich 1971; Kidger et al. 1992). We use half of the FWHM of the peak as errors of the periodicity. For REDFIT, when estimating errors of the periodicity, we choose the minimum of the peak height above the average power level of the red noise in the vicinity of the investigated period. The results of the periodicity analysis are given in Figure 11. The results of REDIFIT show that there is a broad peak (the center frequency V_peak = 2.55 × 10⁻⁴ day⁻¹) with 90% significance levels and above the red-noise power level. The corresponding periodicity is 3918 ± 1112 days in year timescales. The results of the Lomb–Scargle method show that there are three significant peaks, P₁ (4961 ± 1072 days), P₂ (3848 ± 471 days), and P₃ (2793 ± 325 days) above 0.01 FAP levels. The results of the Jurkevich method show that from 2749 (±36) days to 5594 (±44) days, there are some strong periodicity signals. The strong periodicity signal of 3950 ± 14 days also appears in the results of the Jurkevich method. We do not consider periods >5840 days because the periods (considered) are limited to three times periods less than the total time spans (48 years). From the above results, we should note that the results from the REDFIT and Lomb–Scargle methods have large errors, and that for the Jurkevich method, there are some possible periods with range from 2749 days to 5594 days. Therefore, within the errors of the periods, the periods estimated by the three methods have consistent results. A possible quasi-periodicity of P = 3918 ± 1112 days is found.

Zoom In Zoom Out Reset image size

Long-term multiband observations of the quasar 3C 273 were performed from 2005 to 2016 in the V, R, and I bands. The magnitude distributions in the V, R, and I bands are 1310 < V < 1242, 1309 < R < 1237, and 1259 < I < 1193, respectively. The average values of the magnitude are $\langle I\rangle =12\buildrel{\rm{m}}\over{.} 178\pm 0\buildrel{\rm{m}}\over{.} 118$, $\langle R\rangle =12\buildrel{\rm{m}}\over{.} 645\pm 0\buildrel{\rm{m}}\over{.} 107$, and $\langle V\rangle =12\buildrel{\rm{m}}\over{.} 791\pm 0\buildrel{\rm{m}}\over{.} 106$, respectively. Toone (2004) observed this source during the period of 1980–2004 and presented the yearly averaged light curve. The magnitude distributions in the V band from Toone (2004) are 1312–125 with $\langle V\rangle =12\buildrel{\rm{m}}\over{.} 86\pm 0\buildrel{\rm{m}}\over{.} 13$. The magnitude distributions in 2000–2008 from Fan et al. (2009) are 13567–12204, 13313–12014, and 12669–11628 for the V, R, and I bands respectively. The magnitude distributions in 1998–2008 from Fan et al. (2014) are 129–125, 128–123 and 122–1185 for V, R and I bands respectively. The magnitude distributions in 2003–2005 from Dai et al. (2009) are 12824–12663, 12714–12338 and 12318–12055 in V, R and I bands with the mean magnitudes 12698, 12441, and 12139, respectively. The magnitude distributions in 1996 from Xie et al. (1999) are 139–1286, 1298–126, and 1247–1214 in the V, R, and I bands, respectively. Ikejiri et al. (2011) showed that the V-band magnitude distributions between 2008 and 2010 are 1275–1251. From the above references, we obtain that the maximum and minimum magnitudes for this source data are 139–122, 1331–1201, and 1267–1163 for the V, R, and I bands respectively. Through comparisons, we can find that our distributive ranges of magnitude in the V, R, and I bands are within the ranges of the maximum and minimum magnitudes from previous observations. From Figure 10, we also get the same conclusion. Moreover, our long-term multiband observations are updates and supplementary to the quasar data. During our observation, on the whole the quasar becomes dark from 2005 to 2015 but also bright between different years. In 2015, the quasar reaches a low flux state. Compared with 2015, the quasar in 2016 has a tendency to brighten. Over the 12 years, the overall magnitude variabilities are ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 67$, ${\rm{\Delta }}R=0\buildrel{\rm{m}}\over{.} 72$, and ${\rm{\Delta }}V=0\buildrel{\rm{m}}\over{.} 68$, respectively. These results indicate that unlike other blazars, 3C 273 shows moderate magnitude variations on a timescale of many years.

In order to accurately determine the optical IDV, we use two statistical tests to analyze the observation data. The blazar is considered variable only if the light curve satisfies the two criteria of the F-test and ANOVA. During the analysis, we remove the outliers, that is, S_x is more than 3σ. Through these tests and processing, we attain the reliable result that IDV is found on five nights. On 2013 April 04, the largest magnitude change of the I band is ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 055$ in 175 minutes, corresponding to the variability amplitude Amp = 5.39%. The variability amplitudes for the rest of the four nights are close to 5%. Rani et al. (2011) presented two nights detected as an IDV (Amp ∼ 5%) using the C- and F-tests. For the R band from 2014 April 25 and the I bands from 2015 May 17, the light curves meet the criteria of the F-test, but do not reach the critical value of ANOVA. We still consider the two nights detected as do IDV because the light curves of the two nights have large variations with Amp > 30% and ANOVA easily misses this kind of IDV since there are only a few data points in the time interval of flaring or darkening (see Section 3.1). On 2015 May 17, 3C 273 quickly faded by ${\rm{\Delta }}I=0\buildrel{\rm{m}}\over{.} 42$ in 5.8 minutes. Fan et al. (2009) noted that the quasar showed IDVs of 0.55 mag in 13 minutes and 0.81 mag in 116 minutes. The results we also obtain show that the variability amplitudes for most nights are less than 10% and for four nights, more than 20%. The value of the DC is 10.84% for the V case and 53.65% for the PV+V cases. Extensive IDV studies of different subclasses of AGNs revealed that the occurrence of IDV in a blazar observed on a timescale of <6 hr is ∼60%–65%; if the blazar is observed more than 6 hr, then the possibility of IDV detection is 80%–85% (Gupta & Joshi 2005; Rani et al. 2011). Then it is quite likely that the value of the DC is related to the time spans of observation. When considering the nights with time spans >4 hr, the value of the DC is 14.17% for the V case and 55.24% for the PV+V cases. Therefore, the quasar has low variability amplitudes and DC, which are different from other FSRQs. However, we still need more nights with longer time spans to confirm the results further. Although our results show that there are no significant correlations between the variability amplitude and brightness, we do not consider this result reliable because the result is not supported by enough data.

We notice that the variability amplitudes for the non-variable nights are non-zero and greater than 5% or even 10% (Table 5). The variability amplitudes from three nights with exposure time ∼1 minute (V bands on 2007 March 27, 2014 April 21, and 2016 April 24) are more than 5%. If we choose the bin = 3 minutes, the variability amplitudes decrease to 2.4%, 6.17%, and 3.48% respectively. Gopal-Krishna et al. (2003) and Goyal et al. (2013) have found that IRAF and DAOPHOT produce nominal errors that are too small by factors of ∼1.5. It is possible that the non-zero variability amplitudes for the non-variable nights are due to very small values of σ in Equation (3). Assuming variability amplitudes below 1% for the non-variable nights, the underestimated factors range from 1 to 4.7, with average values ∼2. If the original variability amplitudes from Table 5 for the non-variable nights are below 1%, the underestimated factors are set to 1, i.e., the value of σ is not underestimated. For the nights detected as the V case, using the underestimated factors ∼2, we recalculate the variability amplitudes. The results show that compared with the original variability amplitudes from Table 5, the recalculated variability amplitudes do not change significantly. Therefore, for some nights, the uncertainties in Tables 2–4 are underestimated by factors of ∼2. On the other hand, the ANOVA (such as the I band on 2015 May 17 or nights of low sampling rate) and standard F-test could fail to detect IDV (Joshi et al. 2011). This may also partly explain the problem.

From the results of the ACF, the detected characteristic variability timescales are 0.107, 0.023, 0.005, and 0.028 days for the I band on 2013 April 04, the R band on 2014 April 25, and the I bands on 2015 May 17 and 2015 April 15, respectively. Given the minimum values of the lags that can be measured, we cannot claim 0.023, 0.005, and 0.028 days as the significance of the timescales. Since the flares are essentially random and only one or two are seen on 2013 April 04, the variability timescale of 0.107 days (154 minutes) is considered as a possible variability timescale. We also use a least-squares procedure to fit fourth-order and sixth-order polynomials to the ACF on 2013 April 04. The fitting results of different order polynomials show that the differences in the variability timescales are within 2 minutes.

For blazars, the intrinsic origin of IDV is relativistic jet activities or accretion disk instabilities (Marscher & Gear 1985; Chakrabarti & Wiita 1993; Mangalam & Wiita 1993; Marscher et al. 2008). Extrinsic mechanisms include interstellar scintillation and gravitational microlensing. Interstellar scintillation causes radio variations at low frequencies (Heeschen et al. 1987) and gravitational microlensing is important on timescales of weeks to months (Agarwal & Gupta 2015). Therefore, we do not consider extrinsic mechanisms to explain the origin of the IDV. For 3C 273 in the optical band, there is the blue-bump emission. Paltani et al. (1998) presented two possible variable components (B and R) to explain the blue-bump emission. The R component could be synchrotron emission not directly related to the radio–millimeter jet (Soldi et al. 2008). The B component could be explained by reprocessing on an accretion disk (Paltani et al. 1998). The results from Paltani et al. (1998) showed that the B component has a maximum timescale of variability of about 2 years (a mean of 0.46 years) and the R component has a longer timescale of variability. Thus, the two variable components cannot explain the origin of the IDV. In the outburst state, the origin of the IDV can be attributed to the shock-in-jet model—relativistic shocks propagating through a relativistic jet of plasma (Marscher & Gear 1985; Marscher et al. 2008). The observed radiation is predominantly non-thermal Doppler-boosted jet emission. The flare arises as disturbances in the flow of a jet cause a shock wave to propagate along the jet (Marscher & Gear 1985). Models based on the instabilities on the accretion disk could give rise to IDV when the blazar is in the low state, because in the low state, any contribution from the jets is weak (Rani et al. 2011). For the quasar, the instabilities on the accretion disk could yield IDV in the modest state because it has an obvious blue-bump emission. From Figure 10, we see that on 2013 April 04 and 2014 April 25, the quasar has the mean flux level, but on 2015 May 17 it  has a low flux level. On 2015 May 17, a large variability amplitude (41.93%) in 5.8 minutes is found. Although the quasar is in a low flux level on this night, it is impossible to use disk-based models to explain the type of variability, because for luminous quasars with black hole masses in excess of 10⁸ M_⊙, even the fastest disk-based timescale can be a few hours (Joshi et al. 2011). This type of variability is likely to have an association with jet activities. The shocks propagate down relativistic jets, sweeping emitting regions. If the emitting regions have large intrinsic changes (magnetic field, particle velocity/distribution, a large number of new particles injected), then we could see a large flare on a very short variability timescale. Another possible explanation is the Doppler-factor change of the emitting region. The helical jet structure (Gopal-Krishna & Wiita 1992) may cause the Doppler-factor change on a very short variability timescale. On 2014 April 25, the large intrinsic changes in emitting regions are likely to explain the IDV. On 2013 April 04, the small variability amplitude (∼5%) in 154 minutes can be explained by turbulence behind an outgoing shock along the jet (Agarwal & Gupta 2015) or  hot spots or disturbances in or above accretion disks (e.g., Chakrabarti & Wiita 1993; Mangalam & Wiita 1993; Gaur et al. 2012). The interpretations of IDV for the other nights are similar to that on 2013 April 04. If the observed 154 minutes indicate an innermost stable orbital period from the accretion disk, an upper limit can be obtained for the mass of the central black hole (e.g., Xie et al. 2004b; Fan et al. 2014; Dai et al. 2015). From Fan (2005) and Fan et al. (2014), the innermost stable orbit depends on the black hole and the accretion disk, and $r\leqslant c{\rm{\Delta }}{t}_{\min }/(1+z)$, and then the upper limits of the black hole masses are (i) ${M}_{8}\leqslant 1.2\times \tfrac{{\rm{\Delta }}{t}_{\min }(\mathrm{hr})}{1+z}$ for a thin accretion disk surrounding a Schwarzschild black hole; (ii) ${M}_{8}\leqslant 1.8\times \tfrac{{\rm{\Delta }}{t}_{\min }(\mathrm{hr})}{1+z}$ for a thick accretion disk surrounding a Schwarzschild black hole;  and (iii) ${M}_{8}\leqslant \left(\tfrac{7.3}{1+\sqrt{1-{a}^{2}}}\right)\left(\tfrac{{\rm{\Delta }}{t}_{\min }(\mathrm{hr})}{1+z}\right)$ for the case of a Kerr black hole, where M₈ is the black hole mass in units of 10⁸ M_⊙ and a is an angular momentum parameter. Then we obtain that the upper limits of black hole mass are M₈ ≤ 2.66 and M₈ ≤ 3.99 for thin and thick accretion disks, and M₈ ≤ 16.18 for the extreme Kerr black hole case (a = 1). The previous reverberation mapping results obtained M₈ = 2.3–5.5 (Kaspi et al. 2000) and M₈ = 65.9 (Paltani & Turler 2005). Therefore, based on the assumption that the observed 154 minutes indicate an innermost stable orbital period, our upper limits on the black hole masses are consistent with reverberation mapping results.

Fan et al. (2009) found that the quasar shows a time delay between the V and I bands. However, the results of our ZDCF analysis indicate that there are no significant time lags between the V-band magnitude and the I-band magnitude. The optical interband time delay supports the shock-in-jet model. For our results, these factors (e.g., sampling rate, time length of observation, flat light curve) can cause the significant time delay to not be found.

Over the 12 years, the overall color index variability is ${\rm{\Delta }}(V-R)=0\buildrel{\rm{m}}\over{.} 25$. The average value of the color index is $\langle V-R\rangle =0\buildrel{\rm{m}}\over{.} 126\pm 0\buildrel{\rm{m}}\over{.} 023$. From our results, the BWB chromatic trend is dominant for 3C 273 and appears at different flux levels for intraday timescales. The BWB trend exists for short-term timescales and intermediate-term timescales but different timescales have different correlations between the V − R index and the V magnitude. There is no BWB trend for our whole time-series data sets (12 years). The results from Dai et al. (2009) implied that the quasar has the BWB trend on both intraday and long-term timescales (three years). Ikejiri et al. (2011) found that 3C 273 exhibits the BWB trend in their whole time-series data sets (three years). Fan et al. (2014) found that the spectrum becomes flat when the source becomes bright in the overall trend (about five years). However, there is a steeper-when-brighter trend in the details. When considering the V-band magnitude range from Dai et al. (2009), we still find the BWB chromatic trend in the long-term timescales for our data. When considering the V-band magnitude range from Ikejiri et al. (2011), we find a weak BWB chromatic trend for our data. The difference may be due to the different samples and different color indices—Ikejiri et al. (2011) used V − J while we used V − R. From Figure 7 in Fan et al. (2014), there is a gap close to the V-band flux F_V = 28 mJy. From our Figure 10, we find that there are some data close to the V-band flux F_V = 28 mJy. Then if the new data are included in Figure 7 of Fan et al. (2014), the results of Fan et al. (2014) could change. In addition, we notice that compared with previous results, the time spans of our results are longer. The above discussions can explain the differences between our results and previous results in the long-term timescales.

The BWB behavior is most likely to support the shock-in-jet model. According to the shock-in-jet model, as the shock propagating down the jet strikes a region with a large electron population, radiations at different visible colors are produced at different distances behind the shocks. High-energy photons from the synchrotron mechanism typically emerge sooner and closer to the shock front than the lower frequency radiation, thus causing color variations (Agarwal & Gupta 2015). When we consider two distinct synchrotron components, the BWB trend could be explained if the flare component has a higher V_peak than the underlying component (Ikejiri et al. 2011). For 3C 273, there is a weak host galaxy contribution, which is a redder non-variable component. Also, the Doppler-boosted jet emission almost invariably swamps the light from the host galaxy. Gravitational microlensing is important on weeks to months timescales and is achromatic. The optical emission may be contaminated by thermal emission from the accretion disk and the surrounding regions, especially for quasars with a blue bump. The quasar 3C 273 has a blue bump, which flattens the spectral slope in the optical region. Since the thermal contribution is larger in the blue region, the composite spectrum would be flatter than the non-thermal component. Then, when the object is brightening, the non-thermal component has a more dominant contribution to the total flux, and the composite spectrum steepens (Gu et al. 2006). Consequently, in the low flux level, the quasar could have a redder-when-brighter (RWB) trend. For a very low flux state where the thermal radiation of the disk dominates the total flux, the accretion disk model can explain the BWB trend (Li & Cao 2008; Gu & Li 2013; Liu et al. 2016). So the BWB trend on intraday timescales is most likely explained by the shock-in-jet model, and is also possibly due to two distinct synchrotron components or the accretion disk model. In addition, we find that there is a significant anticorrelation between the BWB trend and the length of timescales (Figure 9). Sun et al. (2014) found that the color variability of quasars is prominent on timescales as short as ∼10 days, but gradually reduces toward timescales of up to years, i.e., the variable emission on shorter timescales is bluer than that on longer timescales. They proposed the thermal accretion disk fluctuation model to explain the anticorrelation that fluctuations in the inner, hotter region of the disk are responsible for short-term variations, while longer-term and stronger variations are expected from the larger and cooler disk region. The BWB trend can also be interpreted in terms of two components: a "mildly chromatic" longer-term component and a "strongly chromatic" shorter-term one, which can be likely due to Doppler-factor variations on a convex spectrum and intrinsic phenomena, respectively (Villata et al. 2002, 2004; Gu et al. 2006). The long-term achromatic trend could be due to the superposition of different components (jet components, accretion disk components, gravitational microlensing effect; Fan et al. 2008; Ikejiri et al. 2011; Bonning et al. 2012; Agarwal & Gupta 2015; Gaur et al. 2015; Xiong et al. 2016).

In order to increase the reliability of the periodicity analysis, we appropriately dealt with all of these factors: the best quality data excluding the influence of strong synchrotron flares, three methods for periodicity analysis, unevenly sampled data, red noise, periodicity in year timescales, and periods three times less than total time spans. Converting magnitudes to fluxes in different photometric systems could bring uncertainties that can cause a small (false) variability in the light curve. In this case, if we search for periodicity on short or intraday timescales, the small variability could result in false periodicity. However, the small variability will not cause a significant impact on periodicity in year timescales because the periodicity mainly represents the total change in year timescales. So, analyzing the periodicity in year timescales can reduce the effect from the uncertainties. When analyzing V-band data, a possible quasi-periodicity of P = 3918 ± 1112 days is found. However, it is very possible that a non-periodic model could provide a better fit. Fan et al. (2014) used data over 100 years from the B band (from 1882 to 2012) to analyze the periods of the quasar. Their results obtained possible periods of P = 21.10 ± 0.14, 10.90 ± 0.14, 13.20 ± 0.09, 7.30 ± 0.10, 2.10 ± 0.06, and 0.68 ± 0.05 years for 3C 273. Smith & Hoffleit (1963) suggested a period of P = 12.7–15.2 years for 3C 273. Babadzhanyants & Belokon (1993) analyzed the B-band data of 1887–1991 and found a period of P = 13.4 years. Vol'vach et al. (2013) obtained periods of 11.7 ± 2.5, 7.2 ± 0.8, 4.9 ± 0.3, and 2.8 ± 0.3 years based on the optical B-band light curve of 1963–2011. Therefore, within the errors, the quasi-period of P = 3918 ± 1112 days is consistent with the previous works. For the long-term period of 3C 273, some possible interpretations are as follows: the orbit of a perturbing object, precession of jet, a rotating helical jet structure, and a binary black hole model (Lehto & Valtonen 1996; Villata & Raiteri 1999; Rani et al. 2009; Vol'vach et al. 2013; Fan et al. 2014; Marscher 2014; Sandrinelli et al. 2016).