Curvature induced polarization and spectral index behavior for PKS 1502+106

A comprehensive study of multifrequency correlations can shed light on the nature of variation for blazars. In this work, we collect the long-term radio, optical and $\gamma$-ray light curves of PKS 1502+106. After performing the localized cross-correlation function analysis, we find that correlations between radio and $\gamma$-ray or $V$ band are beyond the $3\sigma$ significance level. The lag of the $\gamma$-ray relative to 15 GHz is $-60^{+5}_{-10}$ days, translating to a distance $3.18^{+0.50}_{-0.27}$ parsec (pc) between them. Within uncertainties, the locations of the $\gamma$-ray and optical emitting regions are roughly the same, and are away from the jet base within $1.2$ pc. The derived magnetic field in optical and $\gamma$-ray emitting regions is about $0.36$ G. The logarithm of $\gamma$-ray flux is significantly linearly correlated with that of $V$ band fluxes, which can be explained by the synchrotron self-Compton (SSC) process, the external Compton (EC) processes, or the combination of them. We find a significant linear correlation in the plot of $\log\prod$ (polarization degree) versus $\log \nu F_{\nu}$ at $V$ band, and use the empirical relation $\Pi \sim \sin^n \theta'$ ($\theta'$ is the observing angle in the comoving frame blob) to explain it. The behaviors of color index (generally redder when brighter at the active state) and $\gamma$-ray spectral index (softer when brighter) could be well explained by the twisted jet model. These findings suggest that the curvature effect (mainly due to the change of the viewing angle) is dominant in the variation phenomena of fluxes, spectral indices, and polarization degrees for PKS 1502+106.


INTRODUCTION
Blazars are active galactic nuclei (AGNs) with emitting jets pointing to our line of sight, which result in relativistic beaming of radiation (Urry & Padovani 1995). Blazars are famous for their high luminosity, rapid variation, and high polarization. Their broadband spectral energy distributions (SEDs) indicate radiation from radio to γ-ray, characterized by two prominent bumps. It is widely accepted that the first bump is caused by synchrotron radiation, while the second one is due to the inversely Compton scattered photons either from the synchrotron radiation in the jet or from external sources outside the jet. However, the locations of emitting regions for these two bumps have been under intensive debate. The correlation analysis between radio and γ-ray helps shed light on the location of γ-ray emitting regions, based on a series of blazar monitoring programs (Cohen et al. 2014;Fuhrmann et al. 2014;Maxmoerbeck et al. 2014a). However, the variation mechanisms to elucidate the color index and polarization behaviors are still away from convincing. PKS 1502+106 (historically OR 103, S3 1502+10 and 4C 10.39) is a single-sided, core-dominated radio-loud AGN, and is classified as a flat spectrum radio quasar, located at R.A. = 15:04:25.0, decl. = +10:29:39. Its redshift is 1.839 (Smith et al. 1977;Richston et al. 1980). The target underwent a strong high energy outburst in 2008 August, which was reported in ATel #1650, followed by multiwavelength campaigns at radio, visual, ultraviolet, and X-ray bands. The study of the target using Very Long Baseline Interferometry (VLBI) observations indicated the misalignment of radio components (An et al. 2004) and jet bending morphology (Karamanavis et al. 2016a). For this target, the multiple radio bands analysis is used to determine the opacity structure, and helps localize the γ-ray emitting region (Maxmoerbeck et al. 2014a;Fuhrmann et al. 2014;Karamanavis et al. 2016a). Kang et al. (2014) reported that γ-ray emitting regions are most likely beyond the broad line region (BLR) by studying a low frequency peaked blazar sample. They found that SED with seed photons from dust torus is better fitted than those from BLR. However, the SED fitting is not unique in principle. The optical emitting regions cannot be obtained directly from images, since almost all blazars are point sources. The correlation analysis between multiple frequencies may require the use of long-term time series to overcome the difficulty of low spatial resolution.
The intense variation study offers essential insights not only into the emitting regions but also into the intrinsic radiation processes. Abdo et al. (2010) investigated radiation from radio to γ-ray of PKS 1502+106, and found that synchrotron and self-synchrotron Compton (SSC) processes dominate from radio to X-ray in SED, while γ-ray radiation is caused by the external Compton (EC) process. The color index behavior in variation is another aspect to reveal the emission mechanism. Villata et al. (2006) and Sasada et al. (2010) stated that the accretion disk emission contributing to the observed flux leads to the bluer when brighter (BWB) trend in outburst state and the redder when brighter (RWB) trend in active state for 3C 454.3. The similar behavior of CTA 102 studied by Raiteri et al. (2017) also can be explained by the same theory. The variation of polarization is also important and can help us to further constrain the radiation models. In this work, we find a significant correlation between polarization degree (PD) and fluxes for the target. Combining with correlations of multiband light curves, the color index, and γ-ray spectral index behaviors, we propose that the geometrical curvature effect can lead to various variation phenomena for the target in an unified manner.
This paper is organized as follows. In Section 2, the γ-ray, optical V and R band, radio 15 GHz and PD data with periods of nearly nine years are collected. The optical data are calibrated using one meter telescope in the Weihai Observatory (WO). The localized cross-correlation function (LCCF) among different time series are calculated, and time lags between them are derived. In Section 3, we obtain the localization of the γ-ray and optical emitting regions, and discuss their positions relative to the BLR. In Section 4, various variation phenomena and their correlations are discussed. Finally, our conclusion is given in Section 5.

Data Reduction
The Fermi Large Area Telescope (LAT) is a highly sensitive instrument with large viewing fields. Its survey scanning mode views the whole sky every 3 hr. The scanning cadence makes it an ideal piece of equipment to monitor γ-ray sources. We collect nearly nine years of γ-ray data (from MJD 54684 to MJD 57932) with energy range 0.1 − 300 GeV from the LAT data server. 1 The data was reduced using Fermi Science Tools version v10r0p5, and was analyzed by adopting the unbinned likelihood method. A 15 • region of interest (ROI) centered on the PKS 1502+106 was considered. To count the γ-ray background calculation, the galactic diffuse emission model (gll iem v06.fits) as well as the extragalactic isotropic diffuse emission model (iso P8R2 SOURCE V6 v06.txt) were applied in the likelihood analysis. Furthermore, make3FGLxml.py was used to create a source model file. The instrument response function was chosen to be P8R2 SOURCE V6. The flux information of the target under the filtering condition TS > 10 were extracted from the gtlike result files. In addition, the γ-ray spectral index in one time bin was obtained by linearly fitting fluxes of seven energy bins, which are logarithmically divided in the range of 0.1 − 218.7 GeV. This will minimize the correlation between flux and its index in the likelihood analysis.
Steward Observatory (SO) has monitored this target for a long period since 2009 February 24 (MJD 54520), which provides the optical V band, R band, and polarization data. 2 However, as for the photometry data, only the differential magnitudes in V and R bands of the target are available, and the comparison star A in the finding chart was not calibrated (Smith et al. 2009). Using the one-meter Cassegrain telescope at the WO, we performed the photometric observations for the comparison star A and the other six Landolt standard stars on 2019 January 16 (MJD 58499). The telescope was mounted with the Johnson/Cousins set of UBVRI filters and the back-illuminated PIXIS 2048B CCD camera with 2k×2k square pixels (Hu et al. 2014). The field of view is about 11 ′ .8 × 11 ′ .8. The images obtained were corrected with the bias and flat field, and the R band finding chart image for the comparison star A and PKS 1502+106 is presented in Figure 1. Following the standard aperture photometric procedure using IRAF, 3 the apparent magnitudes of the comparison star A were calculated to be B = 16.493 ± 0.029, V = 15.521 ± 0.017, and R = 14.984 ± 0.022, respectively. U band data was abandoned due to low signal-to-noise ratio, while the apparent magnitude in I band was not considered due to its less significant fitting of extinction coefficients. In addition, the galactic extinctions in V and R bands are 0.106 and 0.086 magnitude, respectively. 4 The fluxes of targets are calculated based on the differential photometry data and the calibration results above. The flux errors are inherent only from the differential photometry  Figure 2. The SO used the SPOL to derive the Stokes parameters (q = Q/I and u = U/I). The fractional q and u have been calibrated without considering the interstellar polarization. In this work, we consider only the light curve of polarization degree (PD). The polarization angle (PA) has the nπ ambiguity problem (Marscher et al. 2008;Kiehlmann et al. 2016), and the large gap between observations will reduce its validity of variation especially for the long period.
As for radio data, we collect the calibrated data from the OVRO 40 m monitoring program (Richards et al. 2011). The data period are from 2008 January 8 (MJD 54473) to 2017 November 12 (MJD 58099) with 630 points. The long-term duration will enhance the significance of the correlation analysis. The light curves of γ-ray, optical, radio, and PD are shown in Figure 2.

Data Analysis
One popular approach to calculate the correlation between unevenly sampled light curves is the discrete correlation function (DCF; Edelson & Krolick (1988)). The absolute value of DCF can be larger than unity. Another normalized correlation function is invented by Welsh (1999), which is named as the LCCF, and is given by where a i and b j are time series, M denotes the number of (a i , b j ) pairs which satisfy the condition τ ≤ ∆t ij ≤ τ + δt (δt is the bin time),ā τ andb τ are the averaged values of the M pair samples, and σ aτ and σ bτ are the corresponding standard deviations. Following the definition of DCF uncertainty, the error of the LCCF coefficient is taken to be the standard deviation of the local M samples, i.e., The values of LCCF are in the range [−1, 1], which is a good property to be used in significance estimation. Compared with DCF, the LCCF is more efficient to pick up physical signals (Max-Moerbeck et al. 2014b). Thus, LCCF will be used to analyze correlations between various time series of PKS 1502+106. The sampling of the optical observation is extremely uneven, while radio observation has relatively even samplings. Thus, we bin the optical and radio light curves with 4 and 7 day intervals, respectively. The 7 day bin for the optical curve has also been tested, and no obvious change in the LCCF result is evident. In the reduction procedure, a time step of 7 day has already been considered to produce the γ-ray light curve. The light curve of PD is not binned as that of fluxes, since the binning process for PD is of nonlinear and could lead to spurious correlations. The rebinning procedure can smooth the LCCF profile and the significance levels. Although interpolation can reduce the red-noise leakage problem, it can also bring spurious signals especially when there are large observation gaps. Thus, no interpolation has been applied to all observed data. The LCCFs of γ-ray, optical, PD, V − R (magnitude) versus radio are plotted in Figure 3. The range of lag time is taken to be [−2000, 2000] with an 8 day bin. The Monte Carlo (MC) simulation to produce the significance levels is essential to interpret the cross-correlation result. In our recipe, we simulate 10,000 artificial light curves using the TK 95 algorithm with β = 2.3 (Timmer & Koenig 1995), which is an ensemble of the radio light curve. Each light curve contains 3000 points separated by equal bins of 2 days. Then LCCFs between artificial light curves and the observed one are calculated to produce a distribution at each lag bin. The 1σ (68.27%), 2σ (95.45%), and 3σ (99.73%) significance levels are marked with olive dashsed dotted, red dotted and royal blue dashed lines in Figure 3, respectively. This step is different from that of Cohen et al. (2014) and Maxmoerbeck et al. (2014a), who used the completely observed and completely artificial pairs, respectively. The aim of significance estimation is to find the chance probability for a correlated physical signal in a random sample. In some sense, our procedure can reduce the impact of sampling affects from observation, and can escape from the red-noise leakage problem in PSD for optical and γ-ray data. Maxmoerbeck et al. (2014a) obtained that the 3σ significance range in the γ-ray versus radio MC approaches 0.9, which is larger than our result. They concluded that the peak of correlation coefficients is below the 3σ level. First, the main difference between our results and Maxmoerbeck et al. (2014a) stems from the assumed β γ = 1.6 in their simulation. Abdo et al. (2010) analyzed the PSD of the γ-ray light curve with a period of 140 days for PKS 1502+106, and obtained β γ = 1.3. A flatter PSD will decrease LCCF coefficients, and further reduce the significance range. This reason is evident in the comparison between PKS 1502+106 and other two sources (AO 0235+1164 and B2 2308+34). Second, our LCCF peak is a little higher than Maxmoerbeck et al. (2014a), since we use nearly 9 yr of data to perform LCCF. The long duration of the observation will enhance coefficients of correlation, when two curves have physical coherence. The TK95 algorithm considers only the PSD to produce the artificial light curve, whose fluxes have a Gaussian distribution. However, the fluxes of observed light curves usually are non-Gaussian distributed. Such property is characterized by the probability density function (PDF), which can also affect the confidence level. Emmanoulopoulos et al. (2013) indicated that the significance for the peak of DCCF is more conservative, if both PSD and PDF are considered in the simulation. Thus, significance levels in our recipe are still underestimated to some extent.
In the upper left panel of Figure 3, the sharp peak with beyond 3σ significance indicates the strong correlation between the γ-ray and radio light curves. And the variation of γ-ray leads a variation of radio of about dozens of days. In both the upper right panel and lower left panel, there are two peaks beyond the 3σ level, and the most significant peak shows a plateau, which spans less than 300 days. To elucidate the appearance of the plateau, we calculate LCCFs between optical and radio in Epochs I and II, respectively. The plateau occurs in Epoch I, and disappears in Epoch II, see Figure 4. We also plot LCCFs of γ-ray versus V band (blue diamond) and γ-ray versus PD (red triangle) in Figure 5. No significant lag is found between γ-ray and V band flux, while γ-ray leads PD for about 50 days, but the peak coefficient of LCCF (about 0.5) is less significant. In this plot, there is no plateau. Notice that there are several complete flares in both the γ-ray and V band light curves in Epoch III. Hence, the appearance of the plateau is most probably due to the missing rising phase of the giant flare at V band in epoch I, because the correlation between a monotonically decreasing curve and a complete flare curve is invariant in the shift of lag time. The plateau brings trouble when the lag time is estimated. It is most likely that variations of V band and γ-ray are simultaneous, and they both lead to variation of radio. Variation of PD is significantly correlated with variation of radio flux, which helps us to understand the variation mechanism. Detailed lag times will be analyzed in the following.
In the lower right panel of Figure 3, the most significant peak (beyond 3σ) of LCCF is located at −1230, while the second significant peak (about 2σ) is located around zero. Connected with the significance analysis above, the peak around −1230 is a spurious signal. It is impossible that the variation color index leads that of the flux for nearly a thousand days. The magnitude of variation for the color index has no linear relation with that for flux generally. It is shown that V − R has larger magnitude also in quiescent state in Figure 7, leading to a curved shape of color index light curve in Epoch II. The mismatched flares in two time series will lead to a spurious signal in LCCF. Even so, the second peak tells us that the color index varies in a similar cadence with the flux. For this target, the sparse samplings for optical observation hinder us from performing a time delay analysis between them.
The time lag estimation, together with its 1σ error range, is based on the model-independent Monte Carlo method proposed by Peterson et al. (1998). The procedure considers the flux randomization (FR) and random subset selection (RSS; Peterson et al. (1998); Lasson (2012)). Two kinds of time delays are considered, i.e., τ p and τ c . τ p is defined as the lag corresponding to the peak of the LCCF. τ c is the centroid lag of the LCCF defined as τ c = where C i is the correlation coefficient satisfying C i > 0.8LCCF(τ p ). We repeat 10 4 times to obtain a distribution for   Figure 2. From top to bottom, the light curves of γ-ray, optical V band and R band, radio 15GHz, and PD are plotted, respectively. The two vertical dashed lines (MJD 55444 and 57015) divide light curves into three periods, namely Epoches I, II and III, which correspond to the giant flare, the quiescent, and the active state, respectively. τ p and τ c . To better locate time delays, we set that the lag range of LCCF to [−600, 600] and the lag bin to 4 days. The errors of time delays are of the 1σ range from the distribution (Jiang et al. 2018).
The time delays τ p and τ c between different energy bands, together with their average τ , are shown in Table 1. The most significant time delay is between γ-ray and radio, i.e., τ p = −40 +0 −8 or τ c = −80 +10 −11 . Abdo et al. (2010) derived that peak flux of γ-ray arrives 98 days before that of radio 15 GHz for the giant flare in Epoch I, which is close to our result τ c = −80 +10 −11 . Maxmoerbeck et al. (2014a) obtained τ p = −40 ± 13 for γ-ray versus radio, which is almost the same as our result τ p = −40 +0 −8 . We derive τ p = −72 +119 −103 for V band versus radio, which has large uncertainties due to the plateau. Karamanavis et al. (2016a) measured the time lags among different radio frequencies using both DCF and Gaussian process regression (GP) for the target. The time delays between 15 and 142.33 GHz are τ GP = 64.2 and τ DCF = 63 +43 −51 , respectively. This value is close to our result τ c = −64 +7 −6 for optical versus radio. So the radiation of 142.33 GHz and optical band may originate from the same optically thin regime. We obtain that optical V band lags behind γ-ray for τ p = 8 +12 −43 or τ c = 17 +61 −16 . Abdo et al. (2010) also obtained that γ-ray leads optical with 4 days using τ p of DCF, which is roughly in accordance with our result. Within uncertainties, we conclude that γ-ray and optical photons are most probably emitted from the same region.
[b!]  Note-τp and τc are all in units of days, τ is the mean of τp and τc. A negative lag indicates that the former leads the latter.   Note-D is the average of Dp and Dc in units of parsec. The positive value indicates that the emitting region of the former is in the upstream of the latter.

LOCATIONS OF THE OPTICAL AND γ-RAY EMITTING REGIONS [b!]
Assuming a canonical jet, a perturbation propagates along the jet, emitting γ-ray (for instance) at the upstream and radio emission at the downstream. The observed time delay between γ-ray and radio depends on the red-shift z, the viewing angle θ, the velocity of perturbation β = v/c, and the distance D between their emitting regions (see Figure  3 in (Maxmoerbeck et al. 2014a)). Analytically, the distances between different emission regions at different energy bands in the rest frame of the quasar are derived as (Kudryavtseva et al. 2011;Maxmoerbeck et al. 2014a) where β app is the apparent velocity in the observer frame, c is the light speed, T is the time delays (τ p or τ c ) between different bands in observer frame, and z is the redshift, and θ is the viewing angle between jet axis and observing line of sight. For PKS 1502+106, Hovatta et al. (2009) obtained β app = 14.7 using the fastest knot feature. Pushkarev et al. (2009) derived θ = 4 • .7 by the VLBA observation. The target has a redshift z = 1.839 (Smith et al. 1977;Abdo et al. 2010). The jet parameters are estimated by variation of radio fluxes and knot features in images, and can vary in a range for different knot observations. Corresponding to τ p and τ c , we define two kinds of relative distance D p and D c via Equation 3. The derived relative distances, i.e., D p , D c and their average D , are summarized in Table 2. We obtain that the γ-ray emitting region is at upstream of the jet, separating from the radio emitting region with distance 3.18 +0.50 −0.27 parsec (pc) (corresponding to τ ). Meanwhile, relative distance between γ-ray and V band is 0.66 +1.54 −1.94 pc, which indicates that the emitting regions of optical and γ-ray are very close in jet.
To derive distances between emitting regions and the jet base, we need to refer to r core , which denotes the distance between the 15 GHz core region and the jet base. Pushkarev et al. (2012) presented r core = 8.19 pc for PKS 1502+106. Karamanavis et al. (2016b) presented r core = 3.8 ± 0.7pc for 15GHz emissions, using the time lag based core shift measurement, which combines the proper motion and time lags to derive core position offset. The mass of the black hole in PKS 1502+106 is about 10 9 M ⊙ , and the BLR region size is about 0.1 pc. If one takes r core = 8.19 pc, then the γ-ray emitting region is located about 5 pc away from the jet base, which is far away from the BLR region. If one takes r core = 3.8 ± 0.7pc, then the distance between the γ-ray (possibly and optical) emitting region and the jet base is less than 1.2 pc. The possibility that γ-rays are emitted inside the BLR region cannot be ruled out.
Additionally, the magnetic field of emitting zones is derived from the formula B = B 1 r −1 , where B 1 is given by which represents the magnetic filed at 1 pc away from the jet base. For this target, Pushkarev et al. (2012) presented B 1 = 0.69G via core shift measurement, while Karamanavis et al. (2016b) presented B 1 = 0.18G via time lag core position offset. Referring to parameters given by Karamanavis et al. (2016b), the magnetic field in γ-ray and optical emitting regions is about 0.36G.

Optical V Band and γ-Ray
We aim to figure out the emission mechanism at optical and γ-ray bands. The observed fluxes of synchrotron, SSC and EC, which mainly depend on three parameters, i.e., the particle number density N e , the magnetic field strength B and the Doppler factor δ, are given by 5 Chatterjee et al. (2012) where α o is the spectral index at optical band, and α γ is that of γ-ray, and U ′ ext is the energy density of external photons in the jet comoving frame. Taking the logarithm of these fluxes, three parameters can be disentangled. For instance, log F syn = log N e + (1 + α o ) log B + (3 + α o ) log δ + C, where C is a constant. The optical V band radiation is of synchrotron, while γ-ray photons are produced by SSC or EC in the lepton model. We select pairs of γ-ray and V band fluxes observed on the same date (the time difference is less than 2 days), and plot log E 2 dN/dE (γ-ray) versus log νF ν (V band) in Figure 6. If variation of B dominates, one has ∆ log F syn ∼ (1 + α o )∆ log B, ∆ log F EC ∼ 0 and ∆ log F SSC ∼ (1 + α o )∆ log B. Thus, one can expect that the slope in the plot is 1 for the SSC process, and no linear correlation should be found for the EC process. If the Doppler factor is the dominant variable, then one can derive the slope for the SSC versus synchrotron case to be 3+αγ 3+αo . All other theoretical slopes are derived and summarized in Table 3.
Note-The predicted correlations in the plot of γ-ray versus optical fluxes in logarithm for different processes are plotted.  Figure 6. Left panel is the plot of log E 2 dN/dE (γ-ray) vs. log νFν (V band). We subtract a base level flux 3.09 × 10 −13 erg cm −2 s −1 from the νFν at V band, because the nonvariable component (including contributions from host galaxy etc.) should not be included in the variation analysis. And the right panel is the plot of log Π ∼ log νFν . The orange squares, gray circles and pink triangles represent the data in Epochs I, II, and III, respectively.
In Figure 6, it is evident that the logarithm of optical and that of the γ-ray flux is linearly correlated, the statistic slope is 0.78 ± 0.14 with Pearson's r = 0.71. The upper limit of the slope approaches 1. Referring to Table 3, if the variation is caused by N e , the predicted slope is 1 for EC and 2 for SSC. Thus, we can conclude that the variation of N e in the SSC process cannot explain the slope. One can also rule out the case in which the variation of fluxes is due to the change of B in EC process. If the Doppler factor δ is the main parameter, one needs to discuss the range of .38] respectively. Therefore, both EC and SSC varying with δ are possible, since the slope 0.78 ± 0.14 is in these ranges.

Polarization and Optical V Band
In the upper right panel of Figure 3, we obtain a 3σ correlation between PD and radio light curves. This suggests that variation of PD is correlated with optical flux. So we plot log versus log νF ν in the right panel of Figure 6, where the log is strongly related to log νF ν with Pearson's r = 0.77, and the slope is 0.45 ± 0.03. In Cawthorne & Cobb (1990), the polarization was shown to be a function of θ, where θ is the angle between the observer's line of sight and the moving direction of the radiative blob. So we take the empirical relation that ∼ A sin n θ ′ , where n is a positive real number, and θ ′ is the viewing angle in comoving frame of the blob. Raiteri et al. (2013) considered n = 2, based on the study of Lyutikov et al. (2005), which tells that PD reaches its maximum when θ ′ = π/2, and fall to minimum at θ ′ = 0. In the observer frame, we have ∼ δ n sin n θ ≡ δ n(1+α θ ) , where α θ = log δ sin θ is the index related to δ. Thus, is maximal for θ ∼ 1 Γ , corresponding to θ ′ = 90 • in the comoving frame (Nalewajko 2010). When the blob weaves in our line of sight, θ will decrease, passing 1 Γ , at which PD achieves its maximum, and then reaches θ min , at which the flux reaches its maximum. When the blob weaves out, θ will pass 1 Γ again, leading to another peak in the light curve of PD. Thus, one peak of flux is accompanied by two peaks of PD for one outburst for the line of sight inside the beaming cone of the blob. However, the sampling of PD is sparse, so that we cannot verify this correspondence from the two light curves directly. A better sampling example, i.e., the intra-day variation study of S5 0715+714, seems to support such correspondence, see Figure 1 in Chandra et al. (2015). Based on Equation 8, the linear correlation between log and log νF ν , i.e., 2(1+α θ ) 3+αo , can be derived in the case that variations of PD and flux are due to the variation of δ, see Table 3. Since δ > 1, sin θ < 1, one has α θ < 0. Having n = 2 and α o ∈ [2.18, 3.27], one has 2(1+α θ ) 3+αo < 0.39. For n = 3, one has 2(1+α θ ) 3+αo < 0.58. Thus, the observed slope (0.45 ± 0.03) can be explained from Equation 8. The index n is model dependent, the correlation analysis here can constrain the polarization model.
The variations of N e and B cannot lead to the variation of PD in models. The significant correlation suggests that Doppler factor is the key parameter that leads to the variations of fluxes and PD. Finally, the variation of the Doppler factor is mainly due to the variation of observing angle θ. For this target, it was found that the jet component position angles are nonlinearly distributed, and the jet position angles depend on frequencies at radio bands (An et al. 2004). Karamanavis et al. (2016a) studied the dynamical structures based on VLBI images at different radio frequencies, and found that the viewing angle of the inner and outer jet are 3 • and 1 • , respectively. They also derived that the jet opening angle is (3.8 ± 0.5) • , and the downstream of jet (away from the core with 1 mas) bends toward us. Since γ-ray and optical bands are significantly correlated with radio 15GHz, and they are located in the upstream of the radio core, one can expect that radiative blobs of the γ-ray and optical follow a curved trajectory. Note that the radio light curve has a long-term increasing trend in Epoch III, which does not occur in the γ-ray and optical light curves. This can be understood if γ-ray and optical emitting regions have the same viewing angles, which are different from that of the downstream radio emission regions. The curvature of the jet leads to the different orientation of emitting zones for different frequencies (Raiteri et al. 2017). The moving direction of the radiative blob changes when the blob propagates along the jet, which can be realized in the jet precession models (Abraham 2000;Sobacchi et al. 2017). It is noted that if the trail of emitting regions is helical, the distances from the jet base obtained from time lags would be the upper limit. The method to measure the helical trajectory will be an interesting future work.

Variations of Color Index and γ-Ray Spectral Index
The color index versus radio case shows a peak at the 2σ significance level. Thereby, it is necessary to explore the relationship between the color index and fluxes at different frequencies. As for the quiescent state (Epoch II), the large range of spectral index and tiny range of flux below 0.25 mJy make it too difficult to determine the trend of the color index variation. The sparse distribution of the spectral index is probably caused by extra emissions from disk, BLR, torus, or the combination of them. Abdo et al. (2010) derived that the bolometric luminosity of BLR is about 3.7 × 10 45 erg s −1 , equivalent to 0.03 mJy flux at V band. Considering the Eddington limit of the quasar luminosity (Shapiro & Teukolsky 1983), i.e., L Edd = 1.3 × 10 46 erg s −1 (M/10 8 M ), this predicts an extreme 1 mJy flux, if all energy is released at V band. The broadband emission will significantly reduce the flux at V band. Besides, photons from the accretion disk have very small PD, which cannot explain the PD variation at optical band. Due to z = 1.839, the radiation of dust torus will mainly shift to far infrared. Thus, contributions from these components can be ignored when the target is in the active state.
The diagram of V − R versus fluxes at V band is shown in Figure 7. Points in Epoch II (the quiescent state) and Epoch III (the active state) are marked with gray and pink color, respectively. The color index in the quiescent state is scatter, while it has a weak RWB trend (with Pearson's r = 0.32) in the active state. For F ν > 0.9 mJy, both a saturation and BWB trend are likely. In theory, there are several models that can explain the RWB trend. First, Villata et al. (2006) found an RWB trend for 3C 454.3, and the color index is saturated when fluxes approach the maximum (Villata et al. 2006;Sasada et al. 2010). The explanation for the RWB trend is that both accretion disk and jet contribute to the fluxes. The accretion disk contributes a bluer component to the broadband SED, while the synchrotron emission of the jet contributes a redder component to the continuum. Secondly, the RWB phenomenon could also be interpreted by the shock in the jet model (Kirk et al. 1998). When the cooling time scale is roughly the same as the accelerating time scale, the simulation of light curves shows a RWB trend. The BWB trend is due to the fact that the cooling time scale is larger than the accelerating time scale for electrons. However, the accelerating time scale is determined by the strength of the shock. The twisted jet model is the third model to explain the color index behavior. Suppose the spectrum of radiation is of power law F ′ ν ′ ∝ ν ′−αo in the jet comoving frame. The observed frequency and flux are relativistically boosted via ν ∝ δ(θ)ν ′ and F ν ∝ δ 3+αo (θ)F ′ ν ′ (ν) (Urry & Padovani 1995). When δ increases, peak frequency will move from lower frequency to the higher one, the spectral index at a fixed wavelength will undergo variation. For the observing V and R bands, the amplitude of the Doppler factor is the same for both V and R, so one has F νV /F νR = (ν V /ν R ) −αo . Since ν V /ν R > 1, F νV /F νR is larger or smaller than unity for α o < 0 or α o > 0, respectively. Correspondingly, a BWB or RWB trend will occur for α o < 0 or α o > 0 in the optical bands.
Since α o is in the range [2.18, 3.27], both ν V and ν R are higher than peak frequency in SED, which is evident in the broadband SED plot given by Abdo et al. (2010). Gupta et al. (2016) showed that the peak frequency is more variable than other frequencies, which can be interpreted in this twisted jet model. After the peak frequency passes through the observed frequency, the object shows the BWB trend. For PKS 1502+106, there is an RWB trend below ∼ 0.85 mJy, and a less significant BWB trend beyond that, see Figure 7. Combined with correlation analysis, the curvature effect is a better choice to explain the color index behavior. Other models cannot be ruled out by the color index analysis alone. The spectral indices −α γ of γ-ray are obtained by linearly fitting γ-ray fluxes of seven energy grids, while the νF ν γ-ray fluxes (in units of erg cm −2 s −1 ) are obtained by integrating over the 0.1 − 300 GeV range. The −α γ versus log νF ν is plotted in Figure 8. The linear fit, with slope −0.480 ± 0.036 and Pearson's r = −0.816, indicates that the spectral index is anticorrelated with the flux. This is a softer when brighter (SWB) trend. Thus, the variations of spectra at optical and γ-ray are similar, i.e., turning softer when brighter. Meanwhile, the emitting regions of optical and γ-ray bands are close, which is presented in the previous section. It is likely that the SWB trend is due to the intrinsic property, such as the evolution of injected particle distribution. By studying the continuous equation of injected particles, it was shown that the spectral slope of most energetic particles is steeper than that of the less energetic particles, regardless of whether the radiation is in the fast cooling or slow cooling phases (Ghisellini et al. 2002;Ghisellini & Tavecchio 2008).
Such intrinsic spectral evolution can also explain the RWB trend, but the amplification of fluxes needs more injected electrons. However, if the variation is mainly due to N e , the slope in Figure  6 is predicted to be 1 for EC and 2 for SSC. Another possible reason for the SWB behavior is the curvature effect. The analysis of flux behavior also applies to the γ-ray case, the RWB and SWB trend can both be derived in the modulation of Doppler effects. Abdo et al. (2010) (see Figure 11) showed that the peak frequencies of broadband SED for the low and high bumps are lower than the observing optical and γ-ray bands, respectively. Based on the SED, both RWB and SWB trends are natural results of the curvature effect. An et al. (2004) and Karamanavis et al. (2016a) presented that the jet of PKS 1502+106 is twisted. Thus, the observed variation of γ-ray spectral index may also be due to the curvature effect.

CONCLUSION
We gather the multifrequency data of PKS 1502+106, including nearly nine years of data of γ-ray, optical, and radio. From LCCF calculations, we find that the γ-ray, optical V band, and PD light curves are correlated with the radio 15 GHz light curve with significance larger than 3σ. Based on the FR/RSS MC procedure, the γ-ray leads the radio with 60 +5 −10 days, and leads V band with 13 +37 −30 days. According to LCCFs of both whole period and separated period data, we learn that the sparse sampling is responsible for the plateau which appears in the LCCF of V band versus radio. The distance between γ-ray and radio core regions is 3.18 +0.5 −0.27 pc, which is consistent with the result of Karamanavis et al. (2016b). The optical and γ-ray emitting regions are almost the same. Referring to the distance from the 15 GHz core region to the jet base, the γ-ray emitting region is located less than 1.2 pc away from the jet base. NNX09AU10G, NNX12AO93G, and NNX15AU81G. This research has made use of data from the OVRO 40 m monitoring program (Richards et al. (2011)) which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.