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

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 $\gamma$-ray sources. We collect nearly nine years of $\gamma$-ray data (from MJD 54684 to MJD 57932) with energy range $0.1-300$ GeV from the LAT data server.¹¹1{https://fermi.gsfc.nasa.gov/}. The data was reduced using Fermi Science Tools version v10r0p5, and was analyzed by adopting the unbinned likelihood method. A $15^{\circ}$ region of interest (ROI) centered on the PKS 1502+106 was considered. To count the $\gamma$-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 $\rm{TS}>10$ were extracted from the gtlike result files. In addition, the $\gamma$-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.²²2http://james.as.arizona.edu/~psmith/Fermi 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$\times$2k square pixels (Hu et al., 2014). The field of view is about $11^{\prime}.8\times 11^{\prime}.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,³³3IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy Inc., under contract to the National Science Foundation the apparent magnitudes of the comparison star A were calculated to be $B=16.493\pm 0.029$, ${V}=15.521\pm 0.017$, and $R=14.984\pm 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.⁴⁴4http://ned.ipac.caltech.edu/ 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 measurement errors. The light curves obtained are shown in 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\pi$ 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 $\gamma$-ray, optical, radio, and PD are shown in Figure 2.

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 $\tau\leq\Delta t_{ij}\leq\tau+\delta t$ ($\delta t$ is the bin time), $\bar{a}_{\tau}$ and $\bar{b}_{\tau}$ are the averaged values of the M pair samples, and $\sigma_{a\tau}$ and $\sigma_{b\tau}$ 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 $\gamma$-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 $\gamma$-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 $\beta=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\sigma$ (68.27%), $2\sigma$ (95.45%), and $3\sigma$ (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 $\gamma$-ray data.

Maxmoerbeck et al. (2014a) obtained that the $3\sigma$ significance range in the $\gamma$-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\sigma$ level. First, the main difference between our results and Maxmoerbeck et al. (2014a) stems from the assumed $\beta_{\gamma}=1.6$ in their simulation. Abdo et al. (2010) analyzed the PSD of the $\gamma$-ray light curve with a period of 140 days for PKS 1502+106, and obtained $\beta_{\gamma}=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\sigma$ significance indicates the strong correlation between the $\gamma$-ray and radio light curves. And the variation of $\gamma$-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\sigma$ 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 $\gamma$-ray versus $V$ band (blue diamond) and $\gamma$-ray versus PD (red triangle) in Figure 5. No significant lag is found between $\gamma$-ray and $V$ band flux, while $\gamma$-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 $\gamma$-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 $\gamma$-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\sigma$) of LCCF is located at $-1230$, while the second significant peak (about $2\sigma$) 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\sigma$ 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., $\tau_{p}$ and $\tau_{c}$. $\tau_{p}$ is defined as the lag corresponding to the peak of the LCCF. $\tau_{c}$ is the centroid lag of the LCCF defined as $\tau_{c}=\sum\limits_{i}\tau_{i}C_{i}/\sum\limits_{i}C_{i}$, where $C_{i}$ is the correlation coefficient satisfying $C_{i}>0.8{\rm LCCF}(\tau_{p})$. We repeat $10^{4}$ times to obtain a distribution for $\tau_{p}$ and $\tau_{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\sigma$ range from the distribution (Jiang et al., 2018).

The time delays $\tau_{p}$ and $\tau_{c}$ between different energy bands, together with their average $\langle\tau\rangle$, are shown in Table 1. The most significant time delay is between $\gamma$-ray and radio, i.e., $\tau_{p}=-40_{-8}^{+0}$ or $\tau_{c}=-80_{-11}^{+10}$. Abdo et al. (2010) derived that peak flux of $\gamma$-ray arrives 98 days before that of radio 15 GHz for the giant flare in Epoch I, which is close to our result $\tau_{c}=-80_{-11}^{+10}$. Maxmoerbeck et al. (2014a) obtained $\tau_{p}=-40\pm 13$ for $\gamma$-ray versus radio, which is almost the same as our result $\tau_{p}=-40_{-8}^{+0}$.

We derive $\tau_{p}=-72_{-103}^{+119}$ 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 $\tau_{\rm GP}=64.2$ and $\tau_{\rm DCF}=63_{-51}^{+43}$, respectively. This value is close to our result $\tau_{c}=-64_{-6}^{+7}$ 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 $\gamma$-ray for $\tau_{p}=8_{-43}^{+12}$ or $\tau_{c}=17_{-16}^{+61}$. Abdo et al. (2010) also obtained that $\gamma$-ray leads optical with 4 days using $\tau_{p}$ of DCF, which is roughly in accordance with our result. Within uncertainties, we conclude that $\gamma$-ray and optical photons are most probably emitted from the same region.

We aim to figure out the emission mechanism at optical and $\gamma$-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 $\delta$, are given by⁵⁵5To derive Equations (5)-(7), the convention is that $F_{\nu}\propto\nu^{-\alpha}$. Chatterjee et al. (2012)

where $\alpha_{o}$ is the spectral index at optical band, and $\alpha_{\gamma}$ is that of $\gamma$-ray, and $U^{\prime}_{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+\alpha_{o})\log B+(3+\alpha_{o})\log\delta+C$, where $C$ is a constant. The optical $V$ band radiation is of synchrotron, while $\gamma$-ray photons are produced by SSC or EC in the lepton model. We select pairs of $\gamma$-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$ ($\gamma$-ray) versus $\log\nu F_{\nu}$ ($V$ band) in Figure 6. If variation of $B$ dominates, one has $\Delta\log F_{syn}\sim(1+\alpha_{o})\Delta\log B$, $\Delta\log F_{\rm EC}\sim 0$ and $\Delta\log F_{\rm SSC}\sim(1+\alpha_{o})\Delta\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 $\frac{3+\alpha_{\gamma}}{3+\alpha_{o}}$. All other theoretical slopes are derived and summarized in Table 3.

In Figure 6, it is evident that the logarithm of optical and that of the $\gamma$-ray flux is linearly correlated, the statistic slope is $0.78\pm 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 $\delta$ is the main parameter, one needs to discuss the range of $\frac{3+\alpha_{\gamma}}{3+\alpha_{o}}$ and $\frac{4+2\alpha_{\gamma}}{3+\alpha_{o}}$. The range of $\alpha_{o}$ can be obtained from Figure 7 via $\alpha_{o}=(V-R)/2.5\log(\nu_{V}/\nu_{R})$ , i.e., from about 2.18 to 3.27. In Figure 8, $\alpha_{\gamma}$ ranges from 0.834 to 1.566. Thus, the possible ranges of $\frac{3+\alpha_{\gamma}}{3+\alpha_{o}}$ and $\frac{4+2\alpha_{\gamma}}{3+\alpha_{o}}$ are $[0.61,0.88]$ and $[0.90,1.38]$ respectively. Therefore, both EC and SSC varying with $\delta$ are possible, since the slope $0.78\pm 0.14$ is in these ranges.

In the upper right panel of Figure 3, we obtain a $3\sigma$ correlation between PD and radio light curves. This suggests that variation of PD is correlated with optical flux. So we plot $\log\prod$ versus $\log\nu F_{\nu}$ in the right panel of Figure 6, where the $\log\prod$ is strongly related to $\log\nu F_{\nu}$ with Pearson’s $r=0.77$, and the slope is $0.45\pm 0.03$. In Cawthorne & Cobb (1990), the polarization was shown to be a function of $\theta$, where $\theta$ 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 $\prod\sim A\sin^{n}\theta^{\prime}$, where $n$ is a positive real number, and $\theta^{\prime}$ 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 $\theta^{\prime}=\pi/2$, and fall to minimum at $\theta^{\prime}=0$. In the observer frame, we have

where $\alpha_{\theta}=\log_{\delta}{\sin\theta}$ is the index related to $\delta$. Thus, $\prod$ is maximal for $\theta\sim\frac{1}{\Gamma}$, corresponding to $\theta^{\prime}=90^{\circ}$ in the comoving frame (Nalewajko, 2010). When the blob weaves in our line of sight, $\theta$ will decrease, passing $\frac{1}{\Gamma}$, at which PD achieves its maximum, and then reaches $\theta_{\min}$, at which the flux reaches its maximum. When the blob weaves out, $\theta$ will pass $\frac{1}{\Gamma}$ 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\prod$ and $\log\nu F_{\nu}$, i.e., $\frac{2(1+\alpha_{\theta})}{3+\alpha_{o}}$, can be derived in the case that variations of PD and flux are due to the variation of $\delta$, see Table 3. Since $\delta>1$, $\sin\theta<1$, one has $\alpha_{\theta}<0$. Having $n=2$ and $\alpha_{o}\in[2.18,3.27]$, one has $\frac{2(1+\alpha_{\theta})}{3+\alpha_{o}}<0.39$. For $n=3$, one has $\frac{2(1+\alpha_{\theta})}{3+\alpha_{o}}<0.58$. Thus, the observed slope ($0.45\pm 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 $\theta$. 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^{\circ}$ and $1^{\circ}$, respectively. They also derived that the jet opening angle is $(3.8\pm 0.5)^{\circ}$, and the downstream of jet (away from the core with $1$ mas) bends toward us. Since $\gamma$-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 $\gamma$-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 $\gamma$-ray and optical light curves. This can be understood if $\gamma$-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.

The color index versus radio case shows a peak at the 2$\sigma$ 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\times 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\times 10^{46}{\rm erg\,s}^{-1}(M/10^{8}M_{\bigodot})$, 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_{\nu}>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^{\prime}_{\nu^{\prime}}\propto\nu^{\prime-\alpha_{o}}$ in the jet comoving frame. The observed frequency and flux are relativistically boosted via $\nu\propto\delta(\theta)\nu^{\prime}$ and $F_{\nu}\propto\delta^{3+\alpha_{o}}(\theta)F^{\prime}_{\nu^{\prime}}(\nu)$ (Urry & Padovani, 1995). When $\delta$ 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_{\nu_{V}}/F_{\nu_{R}}=(\nu_{V}/\nu_{R})^{-\alpha_{o}}$. Since $\nu_{V}/\nu_{R}>1$, $F_{\nu_{V}}/F_{\nu_{R}}$ is larger or smaller than unity for $\alpha_{o}<0$ or $\alpha_{o}>0$, respectively. Correspondingly, a BWB or RWB trend will occur for $\alpha_{o}<0$ or $\alpha_{o}>0$ in the optical bands.

Since $\alpha_{o}$ is in the range $[2.18,3.27]$, both $\nu_{V}$ and $\nu_{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 $\sim 0.85\,\rm{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 $-\alpha_{\gamma}$ of $\gamma$-ray are obtained by linearly fitting $\gamma$-ray fluxes of seven energy grids, while the $\nu F_{\nu}$ $\gamma$-ray fluxes (in units of ${\rm erg\,cm^{-2}\,s^{-1}}$) are obtained by integrating over the $0.1-300$ GeV range. The $-\alpha_{\gamma}$ versus $\log\nu F_{\nu}$ is plotted in Figure 8. The linear fit, with slope $-0.480\pm 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 $\gamma$-ray are similar, i.e., turning softer when brighter. Meanwhile, the emitting regions of optical and $\gamma$-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 $\gamma$-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 $\gamma$-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 $\gamma$-ray spectral index may also be due to the curvature effect.