Multi-wavelength Variation Phenomena of PKS 0735+178 on Diverse Timescale

BL Lac object PKS 0735+178 showed some complex multi-wavelength variation phenomena in the previous studies, especially for its color behavior. Bluer-when-brighter, redder-when-brighter and achromatic behavior were all found to be possible long-term trends of PKS 0735+178. In this work, we collected the long-term multi-wavelength data of PKS 0735+178, and also performed a multi-color optical monitoring on intraday timescale. The intraday variability was detected on one night. On the long timescale, a possible 22-day time lag was found between the $R$ and $\gamma$-ray bands. The results of the cross-correlation analysis exhibited strong correlations between various optical bands on both intraday and long timescale. However, only a mild correlation was found between the long-term $\gamma$-ray and $R$-band light curves, which could be interpreted by the different emission mechanisms of the $\gamma$-ray and optical emissions. PKS 0735+178 showed a significant harder-when-brighter in the $\gamma$-ray band, which is consistent with the observed optical bluer-when-brighter trend on both long-term and intraday timescales. We found that the HWB and BWB trends will be enhanced during the active states, especially for the historical low state. Such phenomenon indicates a special activity-dependent color behavior of PKS 0735+178, and it could be well interpreted by the jet emission model.


INTRODUCTION
Blazars are the most violently variable class of active galactic nuclei (AGN) (Urry & Padovani 1995), which have their relativistic jets directed towards the earth. They display violent variability of flux and polarization of non-thermal radiation across the entire electromagnetic spectrum. According to whether there are some broad emission lines, blazars can be further divided into flat-spectrum radio quasars (FSRQs) and BL Lac objects. The broadband spectral energy distributions (SEDs) of blazars exhibit two characteristic bumps. One locates at low frequencies (from radio to the UV or X-ray bands) and is dominated by synchrotron radiation. The other locates at high frequencies (the X-ray and γ-ray bands) and is interpreted as the inverse Compton scattering processes (e.g., Ulrich et al. 1997;Böttcher 2007;Sikora et al. 2009). The color behavior can serve as a tool to investigate the nature of emission of blazars and to probe the physical processes in them. The analysis of color-magnitude trends in the monitoring data of blazars has revealed three common patterns, i.e., bluer-when-brighter (BWB), redder-when-brighter (RWB), and achromatic. Through some comprehensive statistical analysis, BL Lac objects frequently show BWB trend, while FSRQs tend to show RWB trend (e.g., Gu et al. 2006;Wu et al. 2011;Gaur et al. 2012;Isler et al. 2017).
Specially, it seems that PKS 0735+178 didn't follow a simple BWB or RWB trend, and researchers have been arguing about its color behavior for more than two decades. As early as in 2000, Fan & Lin (2000) collected 10-year historical optical data of PKS 0735+178 and found that it showed a BWB trend, and the studies of Rani et al. (2010) and Meng et al. (2018) also showed this trend. However, Sandrinelli et al. (2014) found a positive correlation between the R − H color vs. the H-band magnitude, with a general trend indicating bluer color for decreasing flux, i.e., a RWB trend. In the γ-ray band, a comprehensive study of blazars (including PKS 0735+178) illustrated a significant harder-whenbrighter (HWB) trend for FSRQs, while the BL Lacs showed no propensity toward a HWB or softer-when-brighter trend (SWB) (Williamson et al. 2014). Achromatic behavior or weak correlation in optical bands also appeared in the previous studies on PKS 0735+178 (Gu et al. 2006;Ciprini et al. 2007). Recently, Yuan & Fan (2021) found that PKS 0735+178 showed two opposite spectral behaviors among different bands, i.e., BWB for spectral indices (SIs) vs. g-band flux but RWB for SIs vs. r-and i-band flux. All these results suggested a complex pattern of its color behavior.
In order to explore the specific color/spectral behaviors of PKS 0735+178, we collected its long-term optical and γ-ray data, and we also carried out a multi-color optical monitoring program with the 85cm telescope at Xinglong observatory. The descriptions of observations and data reductions could be found in Section 2. In Section 3, we studied the intraday variability (IDV) of this object and searched for the possible inter-band time lags through the interpolated cross correlation function (ICCF). In Section 4, we analyzed its color/spectral behavior in detail, including the behaviors under different timescales and different activity states, to figure out the primary variation mechanism. Finally, a summary is given in Section 5.

Optical: Intraday
The intraday monitoring was performed with the 85cm telescope at Xinglong Station of the National Astronomical Observatories Chinese Academy of Sciences (NAOC). This telescope uses the prime focus optical design with a focal ratio f /3.3. The CCD was a 2048 × 2048 chip with a field of view of ∼ 32.8 × 32.8 arcmin 2 . We observed PKS 0735+178 for 3 nights in Dec. 2020, and ∼400 data points were collected. We have listed the relevant observation information in Table 1.
As shown in Fig. 1, PKS 0735+178 and six field stars are labeled. The data reduction followed the standard process including bias-subtraction, flat-fielding and instrument magnitudes extraction with the Image Reduction and Analysis Facility (IRAF) 2 software. Through aperture photometry, we extracted the instrument magnitudes of PKS 0735+178 and six field stars. In order to obtain the best aperture, the photometry was carried out with ten different aperture radii range from 1 to 3 times of the full width at half-maximum (FWHM), and the inner and outer radii of the sky annulus were accordingly set as 6 and 8 times of the FWHMs, respectively. We finally selected 1.6×FWHM, which gave minimum standard deviations of the differential magnitudes between the check star and the comparison star. In order to minimize the error of calibration, star 1 was chosen as the comparison star, which is the brightest one of all six. Moreover, the closest star 2 was chosen as the check star. The fluxes of PKS 0735+178 and check star were calibrated with respect to a comparison star. The calibrated magnitudes of our target and star 2 were presented in Fig. 2. For clarity, the B-and V -band light curves are shifted by −0.3 mag and −0.4 mag, respectively. Similar shifts of −1.1 mag and −0.5 mag were also applied respectively to the B-and R-band light curves in the small panels.
where N b is the number of data points in the blazar's light curve), (6) the degree of freedom between groups in the enhanced F -test (the number of values between all the groups that are free to vary, and it is calculated by ν2 = ( k j=1 Nj) − k, where Nj is the number of data points in the jth comparison star's light curve and k is the number of comparison stars), (7-8) F and the critical value Fc in the enhanced F -test, (9-10) two degrees of freedom in the ANOVA test (their meanings are the same as those in the enhanced F -test, but calculated by ν1 = g − 1 and ν2 = N b − g, respectively, where g is the number of groups), (11-13) F of PKS 0735+178, F of star C and Fc in the ANOVA test, (14) variable or not, (15) variability amplitude and associated uncertainty, respectively.

Optical: Long-term
We collected the historical optical data from various archives, including those of the Steward Observatory (Smith et al. 2009) 3 , the Small and Moderate Aperture Research Telescope System (SMARTS) 4 , the Catalina Real-Time Transient Survey (CRTS) (Drake et al. 2009) 5 , the Katzman Automatic Imaging Telescope (KAIT) 6 and the All-Sky Automated Survey for Supernovae (Jayasinghe et al. 2019) 7 . Notice that the data from SMARTS are actually raw images after bias subtracted and flat field calibrated, which can be queried from NOIRLab Astro Data Archive 8 , and the subsequent processing is the same as in Sec. 2.1. The collected data are in the V -, R-, and I-bands, but we will not consider the I-band due to the low sampling. If more than one data point is present in a night, the average magnitudes are taken.
2.3. γ-ray energy bins on each time bin. The variation behavior of spectra will be studied by pairing the simultaneous flux and SI Γ γ . The long-term light curves in all three bands were presented in Fig. 3. Moreover, the measurements with TS < 10 were labeled by "upper limits" (red triangles).

Light Curves
From visual inspection of Figure 2 and 3, we can find that the optical light curves are well correlated with each other on both intraday and long timescales, and a significant IDV is exhibited on MJD 59188. However, on the long timescale, the correlation between the optical band and the γ-ray band is clearly weaker than that between two optical bands. A prominent flare was observed during MJD 56500-56800 in the V and R bands. A historical low state (17.23 mags) was observed in the V band during MJD 57350-57450, which was ∼0.2 mag fainter than the previous lowest state (∼17 mags during MJD 50400-50600, see Ciprini et al. (2007)). During MJD 59188-59191, the source was the faintest on MJD 59190 with 15.95 mags in the R band and reached the brightest state on MJD 59188 with 15.78 mags. During our long-term optical observations, PKS 0735+178 reached the brightest state on MJD 56649 in the R band, which was ∼ 0.3 magnitude fainter than its historical brightest state on MJD 59551, i.e., six days before the neutrino event IC211208A. We noticed that all the optical flares over the past decades have far lower intensity than the most recent flare possibly associated with the neutrino event IC211208A. It brightened by ∼1 mag in 50 days since MJD 59500, and then dimmed at nearly the same rate (GCN 31529), which was about three times the rate of that flare during MJD 56500-56800. In the γ-ray band, the object is 44 times brighter in the highest state than in the lowest state.
In order to determine the variability on intraday timescale, we used two robust statistical tests, namely, the powerenhanced version of F -test (de Diego 2014) and analysis of variance (ANOVA) (de Diego et al. 1998), which were widely used in the past IDV studies (e.g. Gaur et al. 2015;Polednikova et al. 2016;Meng et al. 2017;Zhang et al. 2018;Pandey et al. 2020). The detailed algorithms can be found in the corresponding references. If the F exceeds the critical value F c at the 99% confidence level, the null hypothesis will be rejected and the existence of variability (V) will be confirmed, otherwise it is non-variable (N). To avoid detecting spurious variability, an additional ANOVA Figure 1. Finding chart of PKS 0735+178 in the R band. The labels "T" and "1-6" represent PKS 0735+178, the comparison star, the check star, and four field stars, respectively. test was performed to the check star. Only both tests of PKS 0735+178 show "V" and the test of check star is "N", PKS 0735+178 is genuinely variable ("V"), otherwise it is non-variable ("N"). The final results were listed in Table. 1, including two corresponding degrees of freedom, the F value and the critical value F c . PKS 0735+178 only showed IDV on MJD 59188. We found that the observing duration on MJD 59188 was the longest among three nights, and most of the variabilities on this night occurred in the second half night. Longer sampling would naturally have higher probability to detect IDV. We also calculated the variation amplitudes on this night, by where A max and A min are the maximum and minimum magnitudes, respectively, and σ is the measurement error (Heidt & Wagner 1996). The propagating uncertainty of Amp is given by σ Amp = ( Amax−Amin Amp ) · (σ 2 Amax + σ 2 Amin ), where the σ Amax and σ Amin are the measurement errors corresponding to A max and A min , respectively. The amplitudes and associated uncertainties were also given in Table 1. We can find that the variation amplitude tends to decrease with decreasing frequency, which is consistent with the common tendency, i.e., the higher the energy band, the greater the IDV amplitude (Webb et al. 1998;Meng et al. 2017).
On the long-term trend, three distinct γ-ray flares were observed during MJD 55950-56420, 56420-56780 and 57840-58040, respectively, and the second flare was observed quasi-simultaneously in the V and R bands. According to the activities, we divided the γ-ray light curve into five phases, namely, phase A (the 1st quiescent, and from MJD 54839  to 55720), phase B (the 1st active state, and from MJD 55720 to 57280), phase C (the 2nd quiescent, and from MJD 57280 to 57850), phase D (the 2nd active state, and from MJD 57850 to 58230) and phase E (the 3rd quiescent, and from MJD 58230 to the end), which were distinguished by vertical dashed lines in Fig. 3. The optical light curves were segmented according to the γ-ray phase.   Table 2). The solid and dotted lines indicate the measured lags and uncertainties, respectively. The corresponding lag calculation results are listed in each panel. The dashed lines indicate the smoothing of the weighted CCCDs and the histograms are the unweighted CCCDs. The gray regions highlight the ranges of lags included in the final calculation, and the removal efficiency fout were also given.

Cross-Correlation Analysis
We performed the correlation analysis to search for the possible inter-band time lags by using ICCF 9 (Peterson et al. 1998a,b). Through shift and linear interpolation, the Pearson coefficient r between the two light curves will be calculated at each given time shift τ . The optimal time lag and its uncertainties were obtained via the flux randomization/random subset sampling (FR/RSS) method, using Monte Carlo (MC) realizations, as discussed in Peterson et al. (2004). Those realisation with a peak correlation coefficient r peak < 0.5 will be labeled as a failed measurement, and the failure rate f f ail will be recorded. The centroid lag will be estimated using points surrounding the r peak , and we can obtain the cross correlation centroid distribution (CCCD) after 10000 MC realizations. The search range of time lag was empirically fixed as [−200, 200] days and [−60, 60] minutes for long-term and intraday timescales, respectively.
However, the alias of multiple significant peaks has been reported in the previous studies (Grier et al. 2017;Zhang et al. 2018;Homayouni et al. 2019;Li et al. 2019). Thus we followed Grier et al. (2017) and introduced the weighted lag identification. The weight of each data point in CCCD were defined as P = [N (τ )/N (0)] 2 , where N (τ ) corresponds to the number of overlapping epochs at a time shift τ . After a Gaussian kernel estimation (KDE), we can obtain the peak with the largest area in CCCD, and the data points beyond the range of this primary peak will be removed. The final lag was estimated as the median of the unweighted & unsmoothed CCCD within the range of the primary peak, and the 16 th and 84 th percentiles were chosen as the lower and upper uncertainties, respectively. Empirically, a reliable result should satisfy f f ail ≤ 30%, f out ≤ 50% (the ratio of the removed data points to the total data points) and r max ≥ 0.5. All the measurement results that satisfied these criteria were presented in Fig. 4, and corresponding f f ail , f out , the maximum correlation coefficient r max and time lags were listed in Table 2. The negative lags indicate that the former band leads the latter one.
On the intraday timescale, all the time lags were near-zero within the error range. A similar result was obtained between the V and R bands on the long timescale. The non-detection of the lag between variations in different optical bands may be because the emitting regions of these optical bands are too close to each other to result in the time lag (Carini et al. 2011;Wu et al. 2012;Meng et al. 2017). Moreover, according to the data compiled from Sandrinelli et al. (2014), a near-zero lag (−3.6 +4.2 −4.2 days, see Appendix B) was also found between the R and J bands. The only one non-zero result was the 22-day lag between the R and γ-ray bands. On the long timescale, the r max between the γ-ray and R bands is only 0.5, which is much less than the r max = 0.82 between two optical bands. Several studies also found that the correlation between the optical and γ-ray bands was much weaker than that between various optical bands (e.g., Chatterjee et al. 2013;Cohen et al. 2014;Liodakis et al. 2019;Rajput et al. 2020). This is mainly due to the different origins of the optical and γ-ray emissions, the former is dominated by the synchrotron radiation while the latter is by the inverse Compton process. Alternatively, Rajput et al. (2020) demonstrated that the changes of the magnetic field strength could lead to this phenomenon.

Observational Findings
In principle, the color index (CI) should be obtained by pairing the simultaneous data points. Hence, our quasisimultaneous intraday light curves were linearly interpolated to get the intraday CIs. The intraday color behaviors on MJD 59188 were shown in Fig. 5. In order to comprehensively consider the errors in both magnitudes and CIs, we introduced the BCES estimator (Akritas & Bershady 1996). Only when the absolute value of the correlation coefficient exceeded 0.2 and the statistical significance of the correlation exceeded 99% (i.e., p < 0.01), the target was considered to show a confident color-magnitude correlation. It is clear that all the diagrams in Fig. 5 show significant BWB trends.
However, since our long-term optical data were collected from multiple datasets, we can only pair two data points on the same night. Even so, all the obtained CIs clustered in the period of MJD 55000-55400. Thus, we interpolated the long-term light curves with the damped random walk (DRW) method, which was proven to model quasar light curve behavior quite well on time scales of months to years (Gaskell & Peterson 1987;MacLeod et al. 2010;Zu et al. 2013). More details can be found in Appendix A. We interpolated the R-band light curve to obtain the simultaneous data points with those of the V band and calculated the CIs. For comparison, the original CIs without interpolation  were also calculated. The variation of CIs was presented in Fig. 3. The long-term color trends of PKS 0735+178 were presented in Fig. 6, including the overall trend and four individual trends 10 . The BWB behaviors were exhibited on the overall phase, phases A, B and C. Such behaviors were also proven by the color-magnitude correlation without interpolation (red dots and dashed line in the top-left panel). Our results are consistent with the common color trend of BL Lac objects (Fan & Lin 2000). We noticed that the BWB trend exhibited in phase B, which contained a significant flare component, was stronger than that of phase A. Several studies have demonstrated that the BWB trend of some blazars will be enhanced during active states, e.g., S5 0716+714 (Wu et al. 2007), PG 1553+113 (Agarwal et al. 2021) and AO 0235+164 (Wang & Jiang 2020). As a comparison, the phase A that contained only one incomplete flare presented a mild BWB trend, and an achromatic trend was exhibited in phase D that lacked a flare component.
The most complex color behavior was presented in phase C. From visual inspection, two distinct components were found in the color-magnitude correlation, i.e., a strong BWB trend and an achromatic one. To further investigate the color behavior, we divided the variations in phase C into low and high states by the mean magnitude (∼ 16.23 mags) of this period, as shown in Fig. 7. Particularly, the low state corresponds to the historical low state mentioned above. PKS 0735+178 showed no trend in the high state and a strong BWB trend in the low state. One can find that the correlation coefficient of the low state reached 0.664, which means the most significant BWB trend was observed in the historical low state.  Figure 6. Long-term optical color-magnitude correlations of PKS 0735+178. The panels from top to bottom and left to right represent the overall phase and four individual phases, respectively, and the corresponding formula, correlation coefficient r and p value were given in the top-left corner of each panel. We also fitted those data points without interpolation, which are labeled by the red dots and dashed line in the top-left panel, and the specific parameters were given by the red legends.
The γ-ray spectral behaviors were presented in Fig. 8. In the past decade, PKS 0735+178 showed a significant HWB trend, which is consistent to the BWB trend in optical bands. Moreover, the overall HWB trend tends to saturate toward higher flux (see the red dashed lines in the first panel). This is a common phenomenon in the X-ray and γ-ray bands (e.g., Xue & Cui 2005;Giebels et al. 2007;Weng et al. 2020;Acciari et al. 2021), which suggests a more efficient particle acceleration mechanism in the high state (Giebels et al. 2007). As for the individual trends, two significant HWB behaviors appeared in phases A and B, and two low significant HWB behavior appeared in phases C and D. We noticed that the HWB trend of phase B with flare components was obviously stronger than those of phases A, C and D. Such phenomenon is consistent with the situations in the optical bands. The flare in phase D was the shortest one and most of the data points were concentrated in the highest state, so the exhibited HWB trend was not so significant. The spectral hardening during the flare state indicates the increasing of detected high energy photons, which is expected if the inverse Compton peak shifts to higher energy (Shah et al. 2019(Shah et al. , 2021. Moreover, a SWB trend was found with a low significance level in phase E. Since the lack of a radiative efficiently accretion disk, the color behavior of BL Lac objects were always interpreted by the jet emission model (Marscher & Gear 1985). In the frame of the leptonic model of blazars, the electrons are accelerated near the root of the relativistic jet and the magnetic field is compressed therein. This process leads to the observed flux and spectral variabilities (Marscher & Gear 1985;Qian et al. 1991). Additionally, the weak intraday BWB behavior also could be explained by the superposition of many distinct new variable components (Gaur et al. 2015).

Interpretation and Discussions on Color Behavior
Furthermore, the significant BWB trend in the historical low state suggests that the acceleration process in flare and low activity states may have different origins. In principle, we should consider carefully the contamination from the host galaxy in such a low state. However, the brightness of the host galaxy (> 20.8 mag, Scarpa et al. (2000)) was too faint with respect to the nucleus, and therefore it was unlikely to affect the color behavior. Yan et al. (2013) found that the shock acceleration is dominant in the low state, while stochastic turbulence acceleration is dominant in the flare state.
According to the scenario proposed by Virtanen & Vainio (2005), the electrons are accelerated at shock front and move along the jet. The shock acceleration will continuously affect electrons, both in the low and high states. However, the turbulence appears rather strong in the downstream jet, and the stochastic turbulence can have considerable contributions to the acceleration process of electrons (Virtanen & Vainio 2005;Marscher 2013). The shock-accelerated electrons are efficiently re-accelerated in the downstream jet through a stochastic turbulence, and the increasing radiation contributes to the total emission and the blazar will turn to a flare state (Yan et al. 2013;Feng et al. 2020). Moreover, the high-energy γ-ray photons were also thought to be originated in the downstream region (Jorstad et al. 2001). Both shock acceleration and stochastic turbulence acceleration lead to variations in the magnetic field inside the jet, which further change the brightness of the blazar (Kirk et al. 1998). If the inter-band amplitude difference in the low state is larger than that of the flare state, the corresponding color behavior will be more significant (Dai et al. 2011).
Based on a review of previous studies, Ciprini et al. (2007) proposed an achromatic trend of PKS 0735+178. Unfortunately, they lacked the variation data during the 2001 outburst. So they weren't able to check the spectral behavior. Thus, their data suggested a rather constant SI and presented an achromatic trend. This is not in conflict with our conclusions. A similar achromatic trend caused by the incomplete data was also reported by Gu et al. (2006). The RWB trend proposed by Sandrinelli et al. (2014) needs an in-depth discussion. We found that their data also did not contain any complete flare. In our analysis, such properties should be the features of an achromatic or weak BWB trend rather than a clear RWB trend. We noticed that their variation analysis involved the infrared J-band. Isler et al. (2017) mentioned that the color behavior depended on the ratio of the contributions from the accretion disk and relativistic jet, i.e., the ratio of thermal emission to non-thermal emission. In the infrared band, the thermal emission is thought to partly come from the dusty molecular torus that reprocesses the photons from the broad line region and the accretion disk into the infrared region (Antonucci 1993;Perlman et al. 2008). This torus component increases the proportion of the thermal emission and causes the source to exhibit a RWB trend. Moreover, Safna et al. (2020) also showed that the RWB seems a common trend among blazars when an infrared band was involved.
The two opposite color trends reported by Yuan & Fan (2021) seem to interfere by the clustering of the data points (see Fig. 9 therein). We noticed that they didn't take a nightly average for those data points, which would have affected the final result by superimposing the intraday trend on the long-term trend. To verify our idea, we reconstructed the color-magnitude correlations of their data with our method, as shown in Fig. 9. Unfortunately, all the three diagrams didn't show any significant correlations, which were mainly caused by the small sample of data points and associated large uncertainties. So we couldn't conclude which color-magnitude correlation exhibited during their observations.

SUMMARY
In this work, we monitored PKS 0735+178 with the 85cm telescope at Xinglong observatory for three nights, and we also collected 11-years multi-wavelength light curves of this source. The inter-band time lags were calculated through ICCF and simply discussed the correlations between various bands. We analyzed the color and spectral behavior to investigate the variation mechanism. Our conclusions are summarized as follows: • IDV was observed on MJD 59188. Based on the results of ICCF, the only one non-zero lag of 22.1d was detected between the R and γ-ray bands. Too close emitting regions prevent a non-zero lag been detected between two optical bands.
• The correlation between the γ-ray and R bands is much weaker than that between the V and R bands. This is mainly due to the different origins of the optical and γ-ray emission.
• We found that PKS 0735+178 showed a BWB trend in the optical bands on both long-term and intraday timescales, and a HWB trend was found in the γ-ray band. These color or spectral behaviors will be enhanced during the active states. Such color/spectral behaviors could be naturally explained by the jet emission model.
• The most significant BWB trend was found during the historical low state. It could be explained by the larger difference of the inter-band variation amplitude of the trough state compared to that of the flare state, which may be caused by the acceleration process with different origins. Moreover, the saturation of HWB trend in the high state suggests an efficient particle acceleration mechanism.
Either as a blazar with special color behaviors or as a potential or even promising neutrino emitter, PKS 0735+178 is an interesting target and can be monitored intensively at multiple electromagnetic wavelengths and, especially, with IceCube or future neutrino telescope (e.g., Cubic Kilometre Neutrino Telescope), in order to disclose the mechanisms and regions responsible for the neutrino and electromagnetic emissions. Improved localization and sensitivity on neutrino detection is required, and a joint electromagnetic-neutrino observations will be a key in identifying flaring neutrino sources.
The covariance function of a DRW has an exponential form where ∆t is the time interval between two epochs, σ DRW is the DRW amplitude, and τ DRW is the DRW timescale. The DRW interpolation was implemented through JAVELIN 11 algorithm (Zu et al. 2011). First we fitted the R-band light curve to constrain the σ DRW and τ DRW . The V -band light curve was considered as a shifted, smoothed and scaled version of the R-band light curve and shared the same τ DRW and σ DRW . In this process, four additional parameters are added, i.e., the time lag, the kernel width, the transfer function amplitude, and the ratio between the two bands. Through the Markov chain Monte Carlo (MCMC) method, we can obtain the best fitting parameters and use them to predict the two light curves. Furthermore, we restricted the time lags between V and R bands in [−10, 10] days, and the DRW timescale have been restricted in [100, 300] days, which is found for a larger sample of quasars from the Sloan Digital Sky Survey (MacLeod et al. 2010). The final predictions of the two light curves were presented in Fig. A1.
To check whether the interpolated light curve preserves the properties of the original light curve, its time-domain power density spectrum (PSD), i.e., structure function (SF), was calculated. The SF describes the growth of variability with time and is widely used in the analysis of AGN variabilities (e.g., Ulrich et al. 1997;Abdo et al. 2010;Ackermann et al. 2011). Compared to the PSD in the frequency domain, the SF is less subject to sampling problems in the presence of very irregular time series (Simonetti et al. 1985;Paltani et al. 1997). One of the most common form of SF is defined as where τ is the time difference between the data pair, m(t i ) and σ err (t i ) are the magnitude and associated uncertainty at time t i , respectively (Bauer et al. 2009). The uncertainty of SF is defined as where N i is the number of data point pairs in the bin. In our settings, a linear bin τ = 4 days was adopted, and the time scales τ > 200 days was discarded. We calculated the SFs of the R-band interpolated light curve (SF i ) and original light curve (SF o ), and the comparison between the two SFs were presented in Fig. A2. Notice that our CIs were obtained through the interpolated R-band light curve and the original V -band light curve, and thus the comparison between the two SFs in the V band was not given. Except for the deviation in the low τ case (τ < 32 days), the two SFs are well consistent with each other, which proves that the interpolation preserves the properties of the original light curve. One possible reason for the low-τ deviation is the poor fit of the DRW model to the short timescale light curves.   Figure A2. The SFs of the R-band interpolated light curve and original light curve, which are denoted by the black crosses and red dots, respectively.
As stated in Sec. 4.1, we obtained the CIs without interpolation during MJD 55000-55400, and a comparison between the two CIs would also be a good visualization to the impacts of interpolation. As shown in Fig. A3, the CIs obtained with and without interpolation are presented. The time baseline of the two CIs was restricted to the same period, i.e., MJD 55000-55400. It's clear that the color-magnitude correlation with interpolation is only sightly weaker than that without interpolation. It suggests that DRW interpolation does not produce a spurious correlation. Unfortunately, we don't have more data to extend the time baseline to assess its performance on longer timescale.  Figure A3. The comparison between the CIs obtained with and without interpolation, which were denoted by the black and red dots, respectively. The corresponding formula, correlation coefficient r and p value are given.

B. ICCF OF R-J BANDS
We used the raw data from Sandrinelli et al. (2014) to implement the ICCF. We restricted the search range of time lags to [−200, 200] days and 10000 MC realizations were applied. The final result showed a 0% failure rate, 0% removal efficiency and a maximum correlation coefficient ∼ 0.8, which indicate a reliable measurement, and a −3.6 +4.2 −4.2 days lag was detected (see Fig. B4).