Exploring the origin of multiwavelength activities of high-redshift FSRQ PKS 1502+106 during 2014-2018

TL;DR

This study investigates the origins of multiwavelength flares in high-redshift FSRQ PKS 1502+106, revealing magnetic reconnection and shock processes as key triggers through detailed analysis of gamma-ray and optical data.

Contribution

It provides the first detailed multi-band analysis of PKS 1502+106's flares, constraining physical parameters and proposing magnetic reconnection and shock as mechanisms.

Findings

Discovered a 1-hour gamma-ray flare with spectral hardening.

Linked optical polarization to shock compression during 2017 outburst.

Constrained emission region size and magnetic field strength.

Abstract

The origin of the multi-band activities (outbursts/flares) of blazars is still a heavily debated topic. Shock and magnetic reconnection have long been considered as possible triggers for the multi-band activities. In this paper, we present an exploration of the origin of multi-band activities for a high-redshift (z =1.8385) FSRQ PKS 1502+106. Utilizing multi-band data from radio to $\gamma$-ray and optical polarization observations, we investigate two dramatic activities in detail: a $\gamma$-ray dominated outburst in 2015 and an optical dominated outburst in 2017. Our main results are as follows. (I) A fast $\gamma$-ray flare with a flux-doubling time-scale as short as 1-hr in 2015 is discovered. Based on the variability time-scale, the physical parameters of the flaring region (e.g, minimum Doppler factor, emission region size, etc.) are constrained. At the peak of the flare, the $\gamma$-ray spectrum hardens to $\Gamma_{\gamma} = 1.82\pm0.04$ and exhibits an obvious curvature/break characteristic that is caused by the typical "cooling break". Modelings of multi-band SEDs reveal a very hard electronic energy spectrum with the electronic spectral index of $1.07\pm0.53$. This result suggests that this fast $\gamma$-ray flare may be triggered by magnetic reconnection. (II) During the outburst in 2017, the optical polarization degree and optical fluxes show a very tight correlation. By analyzing Stokes parameters of polarization observations, our results show that this outburst could be triggered by a transverse shock with a compression ratio of $\eta> 2.2$, and the magnetic field intensity of the shock emission region is about $0.032$ G.

Tables

Table 1. Table 1: Summary of Swift-XRT Observations of PKS 1502+106

Obs-ID	Date (MJD)	Net Exp.	Flux	${\Gamma }_X$	C-stat./dof
		(s)	(${10}^{-12}$ erg cm-2 s-1)
(1)	(2)	(3)	(4)	(5)	(6)
00036388020a	2015-07-08 (57,211.63)	4722.40	1.29 $\pm$ 0.16	${1.85}_{-0.13}^{+0.13}$	119.4 (137)
00036388021b	2015-07-12 (57,215.14)	1982.85	1.18 $\pm$ 0.27	${1.70}_{-0.25}^{+0.25}$	44.3 (67)
00081590001c	2015-07-14 (57,217.74)	1376.01	1.86 $\pm$ 0.51	${1.08}_{-0.37}^{+0.36}$	47.8 (46)
00036388022d	2015-07-15 (57,218.34)	2976.76	0.41 $\pm$ 0.11	${1.80}_{-0.27}^{+0.30}$	47.6 (45)
00036388023	2015-07-18 (57,221.27)	2904.34	1.36 $\pm$ 0.25	${1.40}_{-0.19}^{+0.20}$	69.2 (82)
00036388024	2015-07-21 (57,224.13)	2881.86	1.73 $\pm$ 0.33	${1.23}_{-0.17}^{+0.17}$	66.3 (92)
00036388025	2015-07-24 (57,227.59)	1882.97	1.54 $\pm$ 0.32	${1.54}_{-0.20}^{+0.20}$	61.8 (70)
00036388027	2016-05-23 (57,531.82)	2976.76	1.42 $\pm$ 0.24	${1.44}_{-0.18}^{+0.17}$	79.2 (86)
00092241001	2016-06-18 (57,557.09)	794.15	1.04 $\pm$ 0.41	${1.80}_{-0.41}^{+0.42}$	12.3 (26)
00092241002	2016-06-26 (57,565.45)	1038.87	0.58 $\pm$ 0.18	${1.95}_{-0.49}^{+0.63}$	19.4 (25)
00092241003	2016-07-04 (57,573.42)	891.53	0.86 $\pm$ 0.65	${1.81}_{-0.60}^{+0.60}$	15.5 (24)
00092241004	2016-07-12 (57,581.27)	1023.90	0.71 $\pm$ 0.36	${1.43}_{-0.57}^{+0.57}$	23.1 (25)
00092241005	2016-07-20 (57,589.24)	976.44	1.25 $\pm$ 0.51	${1.40}_{-0.43}^{+0.42}$	9.9 (24)
00092241006	2016-07-28 (57,597.35)	976.44	0.92 $\pm$ 0.42	${1.52}_{-0.51}^{+0.52}$	13.9 (24)
00092241007	2016-12-19 (57,741.20)	978.94	0.81 $\pm$ 0.29	${1.45}_{-0.53}^{+0.61}$	35.7 (28)
00092241008	2016-12-27 (57,749.90)	933.99	1.56 $\pm$ 0.55	${1.15}_{-0.38}^{+0.37}$	36.5 (32)
00092241009	2017-01-04 (57,757.55)	966.45	0.96 $\pm$ 0.34	${1.46}_{-0.42}^{+0.41}$	22.8 (29)
00092241010	2017-01-12 (57,765.18)	961.46	1.55 $\pm$ 0.55	${1.50}_{-0.30}^{+0.31}$	32.8 (41)
00092241011	2017-01-19 (57,772.56)	1091.31	1.31 $\pm$ 0.40	${1.50}_{-0.33}^{+0.34}$	31.2 (36)
00092241012	2017-01-28 (57,781.13)	1068.84	0.79 $\pm$ 0.21	${1.82}_{-0.39}^{+0.40}$	18.6 (30)
00092241014	2017-02-08 (57,792.64)	566.89	1.94 $\pm$ 0.83	${1.19}_{-0.38}^{+0.38}$	21.2 (26)
00093162001	2017-06-23 (57,927.76)	1083.82	1.24 $\pm$ 0.31	${1.18}_{-0.50}^{+0.59}$	20.3 (26)
00093162002	2017-06-30 (57,934.86)	1078.83	1.28 $\pm$ 0.35	${1.68}_{-0.40}^{+0.45}$	22.8 (35)
00093162003	2017-07-07 (57,941.38)	1098.81	1.96 $\pm$ 0.63	${1.06}_{-0.34}^{+0.33}$	43.8 (44)
00093162004	2017-07-14 (57,948.48)	1003.91	1.72 $\pm$ 0.70	${1.45}_{-0.41}^{+0.39}$	35.3 (36)
00093162005	2017-07-21 (57,955.33)	1023.89	1.30 $\pm$ 0.38	${1.84}_{-0.33}^{+0.37}$	29.2 (39)
00093162006	2017-07-28 (57,962.04)	1048.86	1.00 $\pm$ 0.32	${1.75}_{-0.54}^{+0.60}$	21.1 (29)
00093162007	2017-12-22 (58,109.64)	966.45	1.92 $\pm$ 0.46	${1.91}_{-0.28}^{+0.27}$	34.3 (53)
00036388028	2017-12-27 (58,114.29)	2647.14	1.67 $\pm$ 0.23	${1.67}_{-0.17}^{+0.18}$	93.3 (99)
00093162008	2017-12-29 (58,116.62)	1008.91	1.59 $\pm$ 0.32	${1.69}_{-0.28}^{+0.30}$	37.4 (46)
00036388029	2017-12-30 (58,117.02)	4729.87	2.17 $\pm$ 0.19	${1.71}_{-0.10}^{+0.10}$	130.5 (180)
00036388030	2018-01-02 (58,120.14)	3141.60	1.60 $\pm$ 0.21	${1.64}_{-0.15}^{+0.15}$	95.5 (109)
00093094001	2018-01-04 (58,122.98)	1493.37	2.71 $\pm$ 0.68	${1.04}_{-0.22}^{+0.22}$	57.4 (67)
00093162009	2018-01-05 (58,123.00)	704.26	2.32 $\pm$ 0.58	${1.55}_{-0.31}^{+0.39}$	44.8 (42)
00093162012	2018-01-19 (58,137.26)	1086.32	2.04 $\pm$ 0.53	${1.26}_{-0.28}^{+0.27}$	30.7 (47)
00093162013	2018-01-21 (58,139.39)	1303.59	2.00 $\pm$ 0.49	${1.33}_{-0.24}^{+0.24}$	42.1 (56)
00093162014	2018-01-26 (58,144.76)	1066.34	0.82 $\pm$ 0.24	${2.02}_{-0.40}^{+0.39}$	39.7 (33)
00093162015	2018-02-02 (58,151.81)	973.94	2.35 $\pm$ 0.73	${1.03}_{-0.34}^{+0.34}$	31.2 (41)
00094003001	2018-06-20 (58,289.90)	1133.77	1.00 $\pm$ 0.37	${1.50}_{-0.42}^{+0.56}$	24.9 (32)
00094003002	2018-06-27 (58,296.33)	1463.41	1.18 $\pm$ 0.38	${1.30}_{-0.35}^{+0.35}$	34.0 (33)
00094003004	2018-07-11 (58,310.61)	978.94	1.16 $\pm$ 0.85	${1.19}_{-0.61}^{+0.57}$	30.0 (26)
00094003005	2018-07-18 (58,317.39)	1136.27	1.47 $\pm$ 0.45	${1.21}_{-0.35}^{+0.36}$	36.9 (37)
00094003006	2018-07-25 (58,324.23)	1068.84	0.97 $\pm$ 0.35	${1.58}_{-0.49}^{+0.50}$	34.4 (29)
00094003007	2018-08-01 (58,331.13)	936.00	1.25 $\pm$ 0.43	${1.61}_{-0.36}^{+0.36}$	27.9 (32)
00094003008	2018-08-08 (58,338.45)	1088.82	0.83 $\pm$ 0.24	${2.11}_{-0.36}^{+0.41}$	15.8 (26)

Table 2. Table 2: Observation of XMM-Newton EPIC-pn for PKS 1502+106

Obs-ID	Date (MJD)	Net Exp.	Flux	${\Gamma }_X$	${\chi }_{\mathrm{red}}^2$/dof
		(ks)	(${10}^{-12}$ erg cm-2 s-1)
(1)	(2)	(3)	(4)	(5)	(6)
0763790501	2015-07-14 (57,217.83)	13.7	$1.61\pm 0.04$	${1.67}_{-0.02}^{+0.02}$	0.91 (159)

Table 3. Table 3: Summary of the Likelihood Fitting Results of γ 𝛾 \gamma -ray Spectra at Different Phases of the Flare in 2015.

Power Law (PL)
Phase	$F_{0.1-300\mathrm{GeV}}$	$\Gamma$			$-\log \left(\mathrm{Likelihood}\right)$	${\chi }_{\mathrm{red}}^2$
	(${10}^{-6}$ ph cm-2 s-1)
Pre-flare	1.51 $\pm$ 0.11	1.99 $\pm$ 0.05			8761.88	1.36
Flare	2.21 $\pm$ 0.16	1.80 $\pm$ 0.04			6289.04	1.56
Post-flare	1.17 $\pm$ 0.11	2.01 $\pm$ 0.07			7738.03	1.19
Log Parabola (LP)
Phase	$F_{0.1-300\mathrm{GeV}}$	$\alpha$	$\beta$		$-\log \left(\mathrm{Likelihood}\right)$	${\chi }_{\mathrm{red}}^2$	${\mathrm{TS}}_{\mathrm{curve}}$
	(${10}^{-6}$ ph cm-2 s-1)
Pre-flare	1.38 $\pm$ 0.12	1.78 $\pm$ 0.09	0.11 $\pm$ 0.04		8757.59	2.18	8.58
Flare	2.00 $\pm$ 0.17	1.56 $\pm$ 0.11	0.08 $\pm$ 0.03		6285.18	4.32	7.72
Post-flare	1.05 $\pm$ 0.12	1.76 $\pm$ 0.13	0.12 $\pm$ 0.06		7735.03	1.41	6.00
PLExpCutoff (PLEC)
Phase	$F_{0.1-300\mathrm{GeV}}$	${\Gamma }_{\mathrm{PLEC}}$	$E_{\mathrm{cutoff}}$		$-\log \left(\mathrm{Likelihood}\right)$	${\chi }_{\mathrm{red}}^2$	${\mathrm{TS}}_{\mathrm{curve}}$
	(${10}^{-6}$ ph cm-2 s-1)		(GeV)
Pre-flare	1.43 $\pm$ 0.12	1.84 $\pm$ 0.09	12.6 $\pm$ 7.1		8758.44	0.95	6.88
Flare	2.05 $\pm$ 0.16	1.61 $\pm$ 0.08	16.3 $\pm$ 7.2		6283.91	0.92	10.26
Post-flare	1.08 $\pm$ 0.12	1.78 $\pm$ 0.12	7.6 $\pm$ 4.3		7734.97	0.80	6.12
Broken PowerLaw (BPL)
Phase	$F_{0.1-300\mathrm{GeV}}$	${\Gamma }_1$	${\Gamma }_2$	$E_{\mathrm{break}}$	$-\log \left(\mathrm{Likelihood}\right)$	${\chi }_{\mathrm{red}}^2$	${\mathrm{TS}}_{\mathrm{curve}}$
	(${10}^{-6}$ ph cm-2 s-1)			(GeV)
Pre-flare	1.40 $\pm$ 0.15	1.76 $\pm$ 0.14	2.26 $\pm$ 0.18	0.9 $\pm$ 0.4	8758.38	0.56	7.00
Flare	2.02 $\pm$ 0.17	1.57 $\pm$ 0.12	2.04 $\pm$ 0.14	1.3 $\pm$ 0.6	6285.55	1.15	6.98
Post-flare	1.06 $\pm$ 0.12	1.72 $\pm$ 0.17	2.38 $\pm$ 0.22	1.0 $\pm$ 0.5	7734.79	0.59	6.48

Table 4. Table 4: Optimal Parameters of the Multiwavelength SED Modelings of PKS 1502+106 at Different Epochs of the γ 𝛾 \gamma -ray Flare in 2015.

Phase	$B$	$R_{\mathrm{size}}$	$p_1$	${\gamma }_{\mathrm{br}}$	$P_r$	$P_B$	$P_e$	$P_p$
	(G)	(${10}^{15}$ cm)		(103)	(${10}^{45}$ erg s-1)	(${10}^{44}$ erg s-1)	(${10}^{44}$ erg s-1)	(${10}^{45}$ erg s-1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
Pre-flare	1.7 $\pm$ 0.6	6.3 $\pm$ 0.7	2.23 $\pm$ 0.24	1.2 $\pm$ 0.9	7.4	5.8	11.1	27.1
Flare	1.0 $\pm$ 0.4	6.7 $\pm$ 0.8	1.78 $\pm$ 0.28	1.4 $\pm$ 0.8	13.3	2.1	15.1	18.5
Post-flare I	1.7 $\pm$ 0.5	3.6 $\pm$ 0.5	1.07 $\pm$ 0.53	0.4 $\pm$ 0.1	5.0	2.0	9.3	7.9
Post-flare II	1.7 $\pm$ 0.2	4.0 $\pm$ 0.4	1.74 $\pm$ 0.19	0.4 $\pm$ 0.1	5.5	2.4	10.6	16.3
Post-flare III	1.6 $\pm$ 0.5	8.3 $\pm$ 1.5	2.17 $\pm$ 0.29	1.2 $\pm$ 1.0	5.9	9.8	6.5	14.7

Table 5. Table 5: Polarization Properties of the Variable Component

Stage	$R_{Q-I}$	$r_{Q-I}$	$R_{U-I}$	$r_{U-I}$	$P_{\mathrm{var}}$	${\Theta }_{\mathrm{var}}$
					($\%$)	($\mathrm{deg}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
Rising			0.90	$0.13\pm 0.02$	$49.7\pm 1.4$	$82.4\pm 1.3$
Peak	-0.99	$-0.48\pm 0.01$	-0.89	$-0.22\pm 0.02$	$52.8\pm 1.7$	$102.3\pm 1.2$
Declining			0.59	$0.08\pm 0.05$	$48.6\pm 1.8$	$85.3\pm 3.0$

Equations

$$F(t)=F_c+F_0\times (e^{\frac{t_0-t}{T_r}}+e^{\frac{t-t_0}{T_d}}{)}^{-1},$$

$${\delta }_{\min }=(\frac{{\sigma }_{\mathrm{T}}{d}_l^2(1+z{)}^2{f}_{ϵ}{E}_{\max }}{4t_{var}^{obs}{m}_e{c}^4}{)}^{1/6},$$

$$R\le c{t}_{var}^{obs}\frac{{\delta }_{\min }}{1+z}=1.4\times 10^{15}\mathrm{cm}.$$

$$N(\gamma)=N_{0}n(\gamma)=\left\{\begin{array}[]{rcl}N_{0}\gamma^{-p_{1}}&&\gamma_{\min}<\gamma\leq\gamma_{\mathrm{br}}\\ N_{0}\gamma^{-p_{2}}\gamma_{\mathrm{br}}^{p_{2}-p_{1}}&&\gamma_{\mathrm{br}}<\gamma\leq\gamma_{\max},\\ \end{array}\right.$$

$$P_i=\pi R_{\mathrm{size}}^2{\Gamma }^{2}c{U}_i\approx \pi R_{\mathrm{size}}^2{\delta }^{2}c{U}_i,$$

$$p\approx \frac{\alpha +1}{\alpha +5/3}\frac{(1-{\eta }^{-2}){\sin }^2\Psi }{2-(1-{\eta }^{-2}){\sin }^2\Psi }.$$

$$\Psi ={\tan }^{-1}\{\sin \theta /\Gamma (\cos \theta -\sqrt{1-{\Gamma }^{-2}})\}.$$

$$t_{\mathrm{cooling}}\approx 4.75\times 10^2\frac{(1+z)}{\delta {\nu }_{\mathrm{GHz}}{B}^3}\mathrm{days},$$

Full text

galaxies: active — galaxies: jets — radiation mechanisms: non-thermal — quasars: individual (PKS 1502+106)

††software: Fermi Science Tools (Asercion & Science Support Center, 2014), HEAsoft (Heasarc, 2014), SAS (Gabriel et al., 2004), XSPEC (Arnaud, 1996), LMFIT (Newville et al., 2016), Astropy (Astropy Collaboration et al., 2013), SciPy (Jones et al., 2001-), NumPy (Walt et al., 2011), Matplotlib (Barrett et al., 2005)

1 INTRODUCTION

Blazars, including BL Lac objects and flat-spectrum radio quasars (FSRQs), are the most extreme type of radio-loud active galactic nuclei (AGN); their relativistic jets point close to our line of sight, for which they possess distinctive observational characteristics, such as large amplitude and rapid variability, high and variable polarization, and compact radio emission (see, e.g., Urry & Padovani, 1995; Massaro et al., 2015; Padovani et al., 2017). Their radiation spans the entire electromagnetic wave, and as non-catastrophic sources, they can be monitored for a long time, thus becoming one of the unique objectives in the era of multi-messenger astronomy (see, e.g., Rani et al., 2019; Burns et al., 2019).

Multi-band variability in blazars is an important window to understand the non-thermal radiation mechanisms and the energy dissipation of relativistic jets. Generally, the variability of blazars exhibits random behavior (e.g., Abdo et al., 2010b; Chatterjee et al., 2012; H.E.S.S. Collaboration et al., 2017), similar to the variability of radio-quiet AGNs (e.g., Markowitz et al., 2003). However, in some activities (outbursts/flares), the multi-band variability of blazars presents recognizable patterns (e.g., Marscher et al., 2008). These activities possess a variety of characteristics, with variability time-scales ranging from minutes to years (see, e.g., Raiteri et al., 2013; Ackermann et al., 2016; Shukla et al., 2018). Using the shortest variability time-scale of these activities, we can constrain the physical properties of the emission region, such as the size, the magnetic field intensity, and the location of the high-energy emission region (e.g., Böttcher & Reimer, 2004; Yan et al., 2018; Zhu et al., 2018). However, the physical origin of these activities is still not well understood. A variety of scenarios have been proposed (see Böttcher 2019 for a review), and these scenarios can be roughly classified into models associated with shock (i.e., shock-in-jet models, e.g., Marscher & Gear 1985; Böttcher & Dermer 2010), models related to turbulence/magnetic reconnection (e.g., Giannios et al., 2009; Giannios, 2013; Sironi et al., 2015), models in which jets interact with external environments (e.g., jet-star collision model; Barkov et al. 2012; Araudo et al. 2013), and geometric models (e.g., Larionov et al., 2013; Raiteri et al., 2017). It can be noted that shock and magnetic reconnection are the two main internal causes triggering these multi-band activities.

PKS 1502+106 is a powerful FSRQ at a high redshift of $z=1.8385$ (Adelman-McCarthy et al., 2008). High-resolution very-long-baseline interferometry (VLBI) observations at mm wavelengths reveal that it possesses a compact, core-dominated morphology with a one-sided, bent, magnetically dominated parsec-scale jet; the jet has a fast superluminal motion with the estimated Doppler factors ranging from $\sim 7$ up to $\sim 50$ along the jet, showing obvious acceleration characteristic (Karamanavis et al., 2016a). At the beginning of the Fermi satellite operation (August 2008), PKS 1502+106 exhibited a sudden high-energy flare, becoming the second brightest extragalactic source in the $\gamma$-ray sky at that time. This was followed by high and variable fluxes over the subsequent months, triggering an intensive multi-band campaign covering the radio, optical, ultraviolet (UV), and X-ray bands. Abdo et al. (2010a) utilized these synergistic observation data to analyze this dramatic outburst. Their results suggest that PKS 1502+106 is a sub-GeV peaked blazar, and its energy dissipation probably occurs within the broad-line region (BLR).

A minor renewed $\gamma$-ray activity observed by the Fermi-LAT was announced again in Jan 2009 through ATel #1905. After this, PKS 1502+106 entered a quiescent stage of up to six years. Until mid-2015, a highest-level renewed $\gamma$-ray activity was detected by the Fermi-LAT (see ATel #7801), and associated outbursts from near-infrared to UV bands were also observed (see ATel #7804 and ATel #7783). Subsequently, in December 2017, Steward Observatory observed extremely high optical polarization from PKS 1502+106. The optical polarization degree is as high as $(47.4\pm 0.1)\%$, which is among the highest levels of polarization ever observed for blazars (see ATel #11047). Long-term optical photometric monitoring indicates that PKS 1502+106 was undergoing a prominent optical outburst. Motivated by these dramatic multi-band activities, in this paper, we systematically investigate the multiwavelength activities of the source during 2014–2018 using multi-band data from radio to $\gamma$-ray bands as well as optical polarization observations. Especially, two dramatic outbursts, a $\gamma$-ray dominated outburst in 2015 and an optical dominated outburst in 2017, are analyzed in detail to explore the physical origin of their multiwavelength activities and the physical properties of emission regions.

This paper is organized as follows. Section 2 describes multi-band observations and data reduction. Section 3 presents the global characteristics of the multi-band variability of the source during 2014–2018. The $\gamma$-ray dominated outburst in 2015 is investigated and discussed in Section 4, and the optical dominated outburst in 2017 is analyzed and discussed in Section 5. In Section 6, we summarize the main results of this work. Cosmological parameters of $H_{0}=70~{}$km s*-1Mpc-1*, $\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$ (Planck Collaboration et al., 2016) are adopted in this work.

2 Observations and Data Reduction

In this section, we present the multiwavelength observations (from radio to $\gamma$-rays bands) of PKS 1502+106 from October 2014 to October 2018 and the processes of data reduction.

2.1 Gamma-ray Observations: Fermi-LAT

Using the Fermi data server111https://fermi.gsfc.nasa.gov/ssc/data/, the newest Pass 8 data of PKS 1502+106 from October 10, 2014 to October 10, 2018, are acquired. Following the standard procedure222https://fermi.gsfc.nasa.gov/ssc/data/analysis/, we use the Fermi Science Tools v10r0p5 with the P8R2_SOURCE_V6 instrument response function to analyze the acquired data. The 0.1–300 GeV events (evclass=128 and evtype=3) are extracted within a $10^{\circ}$ region of interest (ROI) centered on the location of PKS 1502+106 (R.A.=226.104, decl.=10.494, J2000) with gtselect. In order to eliminate the Earth’s limb events, the recommended quality cuts, (DATA_QUAL$==$1)$\&\&$(LAT_CONFIG$==$1), and a zenith angle cut at $90^{\circ}$ are applied. To prepare for variability and spectral analyses at different time-scales, a binned likelihood analysis to the total acquired data is first performed with gtlike to obtain an initial spectral model. In this analysis, the input model file describing ROI is created using the Fermi-LAT third source catalog (3FGL; Acero et al. 2015), together with the latest isotropic background model, iso_P8R2_SOURCE_V6_v06 and the galactic diffuse emission model, gll_iem_v06.333https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html The model file contains sources within ROI+$10^{\circ}$ from the target. For sources lying within 10*∘* from the center of ROI, their photon indexes and normalized parameters are left free to vary during the model fitting. For sources lying within $10^{\circ}$–$20^{\circ}$ from the center of ROI, their parameters are kept fixed to the 3FGL catalog values. In addition, the normalization of diffuse background components is kept free during the model fitting.

Unbinned likelihood analyses are applied to extract fluxes and spectra under different time-scale requirements. In unbinned likelihood analyses, the initial spectral model obtained in the above binned likelihood analysis is used as the input model file. Considering the need for the convergence of likelihood fitting under short time-scale, for sources with the maximum-likelihood test statistic (TS; Mattox et al. 1996) less than 25, their photon indexes and normalized parameters are set to be fixed regardless of their location. The fixed/free settings of parameters for other sources are consistent with that used in the binned likelihood analysis. Unless otherwise stated, the $\gamma$-ray spectrum of PKS 1502+106 is considered as a simple power-law form throughout the analyses.

2.2 X-ray Observations: Swift-XRT and XMM-Newton

Using the astronomical archival data retrieval service provided by HEASARC,444https://heasarc.gsfc.nasa.gov/docs/archive.html we search for the available Swift observations of PKS 1502+106 from October 10, 2014 to October 10, 2018. Fifty Swift observations are retrieved. Among them, 11 observations pointed to PKS 1502+106, and 39 observations pointed to Mkn 841. Mkn 841 is a Seyfert type-1 galaxy located about 7*′* SW of PKS 1502+106, which allows the observations of Mkn 841 to provide the serendipitous X-ray data for PKS 1502+106. In order to obtain reliable photometric results, among the 50 observations, we only analyze the XRT data of 45 observations with net exposure time greater than 0.5 ks. The 45 observations are all made in the photon counting mode. The XRT data of these observations are processed by the task xrtpipeline (v0.13.4). The calibration files (CALDB v20180710) and the standard filtering and screening criteria are used in this process.

We use the xselect tool to extract source photometry and spectrum, and we use the xrtmkarf task to produce ancillary response file for each observation. A circular region with a radius of $47^{\prime\prime}$ centered at the object is used as the extracted area of the source, while a source-free annular region with a radius of $70^{\prime\prime}$ slightly away from the object is chosen as the extracted area of background.

We use XSPEC v12.9 to fit the extracted 0.3–10 keV spectrum of each observation. Because of the low number of events from the source, events are not grouped, and $C$ statistics (Cash, 1979) is employed to determine the best fit in the fittings. A power-law model modified by Galactic absorption ($N_{\mathrm{H}}=2.20\times 10^{20}$ cm*-2*) is used in the fittings. The XRT observation logs of PKS 1502+106 are summarized in Table 1, and the 0.3–10 keV model fluxes555The energy detection range of the XRT is 0.2–10 keV. However, in most observations, only 0.5–8 keV events (photons) are detected (due to the reduction of the effective area at high energy and low energy ends), which means that the model is mainly constrained by the 0.5–8 keV events., photon indices, and their $1\sigma$ uncertainties obtained from the spectral fittings are also listed in Table 1.

In the 45 observations, the count rate of PKS 1502+106 (see Figure 3) is lower than the typical lower limit of the count rate ($\sim 0.5$ counts s*-1*) that will produce pile-up, suggesting that pile-up is unlikely to occur in these observations. To further confirm whether pile-up exists, we analyze each observation following the standard pile-up analysis process666http://www.swift.ac.uk/analysis/xrt/pileup.php, and we do not find any pile-up effect in these 45 observations. This result excludes the possibility that the hard photon indices ($\sim 1.0-1.3$) in Table 1 are due to the pile-up effect.

XMM-Newton performed an observation in the Imaging mode for Mkn 841 on July 14, 2015. At this time, PKS 1502+106 happened to be in the brightest stage in the $\gamma$-ray band (see Section 4 below). This observation, therefore, provides an important serendipitous X-ray data with a net exposure time of up to $\sim 13$ ks for PKS 1502+106. The European Photon Imaging Camera (EPIC) on board the XMM-Newton is composed of three co-aligned X-ray telescopes, which simultaneously observe a source by accumulating photons in the three CCD-based instruments: the twins MOS1 and MOS2 and the pn. Here we only consider the EPIC-pn data as it is most sensitive. Following "The XMM-Newton ABC Guide" (v4.7), we use the XMM-Newton Science Analysis System (SAS v16.1.0) to extract the photometry and spectrum of PKS 1502+106. We extract the $0.3$–$10$ keV source events from a circular region with a radius of $45^{\prime\prime}$ centered on the object, and we extract the background events from a source-free circular region, with a radius of $70^{\prime\prime}$, slightly away from the object. Using the epatplot task, we exclude the presence of pile-up. The standard background count rate threshold (RATE$<=0.4$) for screening good time intervals and the standard event pattern and flag filtering criteria (FLAG $==0$ && PATTERN $<=4$) are applied in the extraction process. The tasks rmfgen and arfgen are used to create redistribution matrix file and ancillary response file, respectively. The extracted X-ray spectrum is grouped into at least 25 counts per energy bin to ensure the validity of $\chi^{2}$ statistics. Because the X-ray spectrum that obtained by the XMM-Newton has more photons, we first attempt to use a power-law model with both intrinsic absorption and Galactic absorption to fit it. The fitting result shows that the intrinsic absorption of PKS 1502+106 is very weak (the $90\%$ upper limit of $N_{\mathrm{H}}$ is $0.9\times 10^{20}$ cm*-2*), so we finally consider the power-law model only modified by Galactic absorption to fit its X-ray spectrum. The observation log and the X-ray spectrum fitting result are listed in Table 2.

2.3 UV/Optical Observations: Swift-UVOT and Steward Observatory

Another instrument on the Swift, UVOT, provides ultraviolet/optical data. In the 50 Swift observations carried out between October 10, 2014 and October 10, 2018, except for those observations where PKS 1502+106 falls outside the field of view, the UVOT data of 36 observations are available. The Swift-UVOT has three optical-band filters (V, B, U) and three UV-band filters (W1, M2, W2). Of the 36 available observations, 12 observations obtain six band data, and the remaining 24 observations only obtain the W2-band data. For each observation, we sum multiple images in the same filter with the task uvotimsum and then perform aperture photometry with uvotsource. Source counts are extracted from a circle with $5^{\prime\prime}$ radius centered on the object, while background counts are extracted from a circle with $20^{\prime\prime}$ radius located in a source-free region close to the object. The observed magnitudes are corrected for Galactic extinction using the reddening coefficient of $E_{(B-V)}=0.0275$ mag from Schlafly & Finkbeiner (2011) and Fitzpatrick (1999) reddening law with $R_{V}=3.1$. The corrected observed magnitudes are then converted into fluxes using the zero points (Breeveld et al., 2011).

As a part of the Fermi multiwavelength support programme, the Steward Optical Observatory of the University of Arizona perform long-term optical monitoring for PKS 1502+106 (Smith et al. 2009). The optical R-band and V-band photometric data and polarimetric data (optical linear polarization degree and polarization angle) of PKS 1502+106 can be easily acquired from the public archival database777http://james.as.arizona.edu/%7epsmith/Fermi/. 153 R-band and V-band photometric data and 188 polarimetric data are obtained from the optical monitoring between October 10, 2014 and October 10, 2018. Since the comparison star (SDSS J150418.48+102757.6) of PKS 1502+106 is not calibrated, only the differential magnitudes of PKS 1502+106 relative to the comparison star are provided in the database. In order to recover the apparent magnitudes of PKS 1502+106, we calculate the apparent magnitudes of the comparison star by utilizing its SDSS photometric data.888SDSS J150418.48+102757.6 is a star and its SDSS magnitudes in the u, g, r, and i bands are $17.932\pm 0.011$ mag, $15.984\pm 0.004$ mag, $15.225\pm 0.004$ mag, and $14.947\pm 0.004$ mag, respectively (Alam et al., 2015). Following the recipe of Lupton et al. (2005), its apparent magnitudes in the R and V bands can be converted from the SDSS magnitudes. The estimated apparent magnitudes of the comparison star in the R and V bands are $15.0$ mag and $15.5$ mag, respectively. Based on the estimated apparent magnitudes, we obtain the apparent magnitudes of PKS 1502+106, while we consider a large error ($0.1$ mag) on the recovered magnitudes account for the uncertainties introduced in the apparent magnitude estimation of the comparison star. The recovered magnitudes of PKS 1502+106 are corrected for Galactic extinction and then converted into fluxes.

2.4 Radio Observations: OVRO 40m Telescope

PKS 1502+106 is observed in radio band by the Owens Valley Radio Observatory (OVRO) 40m Telescope as a part of Fermi multiwavelength support programme (Richards et al., 2011). Through the database999http://www.astro.caltech.edu/ovroblazars/index.php?page=home provided by the OVRO, we obtain the radio data of PKS 1502+106 detected at 15 GHz from October 10, 2014 to October 10, 2018.

3 Global characteristics of multi-band variability

In this section, we first summarize and analyze the global characteristics of the multiwavelength variability of PKS 1502+106 from October 2014 to October 2018. Then, two significant outbursts (see Figure 1) in this period are investigated in detail in Sections 4 and 5 below.

3.1 Multiwavelength light curves

Figure 1 presents the multiwavelength light curves of PKS 1502+106 from October 2014 to October 2018. The first panel shows the 7-day binning light curve in the $\gamma$-ray band. Before April 30, 2015, PKS 1502+106 were in a quiescent stage with an average $\gamma$-ray flux of $1.1\times 10^{-7}$ ph cm*-2* s*-2*.

Between May 2015 and September 2015, a significant $\gamma$-ray outburst ($\gamma$-ray dominated outburst) occurred in PKS 1502+106. The peak flux in the $\gamma$-ray band was $1.54\times 10^{-6}$ ph cm*-2* s*-2* at July 13, 2015 (for 7-day binning data). During this period, the X-ray and optical/ultraviolet fluxes showed associated changes. There was a gap period in the OVRO radio observations. Nevertheless, it can still be seen that the radio fluxes did not show significant changes, with a fractional variability101010The fractional variability is usually used to reflect variability amplitude. It and its $1\sigma$ uncertainty are calculated using Eqs. 10 and 11 in Vaughan et al. (2003). of only ($6.8\pm 0.2$)%. During the outburst, the average optical linear polarization degree of PKS 1502+106 was ($15.1\pm 0.5$)%, and the fractional variability of the polarization degree was ($29.6\pm 2.7$)%; the polarization angle (PA) did not seem to have a significant swing of $>180^{\circ}$, and the maximum rotation amplitude was $\sim 44.9^{\circ}$. Here, to solve the $\pm 180^{\circ}$ ambiguity, following the recipe of Larionov et al. (2013), we added/subtracted $180^{\circ}$ each time that the subsequent PA value is $>90^{\circ}$ less/more than the preceding one. However, it should be noted that due to sparse polarization observations, the $\pm 180^{\circ}$ ambiguity of PA may not be completely eliminated. There is still a possibility that the PA rotates an angle of $>90^{\circ}$ within a few days (see, e.g., Marscher et al. 2008; Larionov et al. 2013; Kiehlmann et al. 2016), and such rotation events will be handled improperly by the current processing method, resulting in an underestimation of the variability of the PA.

During the period from September 20, 2015 to April 30, 2017, the $\gamma$-ray emission of PKS 1502+106 was in a plateau with average flux slightly higher than that of the quiescent stage (about four times higher); the $\gamma$-ray fluxes showed some minor fluctuations, with a fractional variability of ($38.0\pm 2.7$)%. In the X-ray and UV/optical bands, there were also only some minor fluctuations and no significant outburst. There were large random variations in the optical polarization degree, with a fractional variability of $(45.4\pm 2.2)$%.

The second significant outburst occurred between May 2017 and March 2018. It can be clearly seen that the optical emission of PKS 1502+106 showed a clear outburst profile, and the maximum flux in the R band reached $(3.1\pm 0.3)$ mJy, about six times that of the quiescent stage. During this period, X-ray emission and polarization degree showed associated changes, but $\gamma$-ray fluxes did not show significant variations, with a fractional variability of ($28.1\pm 1.6$)%, so we called this outburst an optical dominated outburst.

3.2 Multiwavelength correlations

By investigating the correlation between different energy bands, we can obtain some information about the radiation mechanism and the location of the emission region in different bands (see, e.g., Max-Moerbeck et al., 2014a; Ramakrishnan et al., 2016; Liodakis et al., 2018). Therefore, we employ discrete correlation function (DCF, Edelson & Krolik 1988), which is one of the most common methods for studying cross-correlation between two time series with uneven spacing, to perform cross-correlation analyses on the multi-band data of PKS 1502+106. The significance of DCF is estimated by a Monte Carlo (MC) method introduced by Max-Moerbeck et al. (2014b). In MC simulations, we generate 1000 random light curves for each band with a simple power-law power spectral density (PSD) ($P(\nu)\propto\nu^{-\beta}$) following the algorithm described in Timmer & Koenig (1995); the spectral index $\beta$ for each band is determined by fitting the PSD of the actual observed light curve.

Figure 2 shows the DCF between the $\gamma$-ray emission and the emissions in the R, V, radio bands and the optical polarization degree over the whole observation period, respectively. The 7-day binning interval is adopted in the analyses. The X-ray and UV bands are not included because they have fewer observation data. We do not find any significant ($>3\sigma$) correlations between the $\gamma$-ray band and the optical bands (R and V bands). This may be due to the fact that PKS 1502+106 is in a steady state for most of the whole observation period. In addition, the thermal radiation from the accretion disk may also result in the weakening of the correlation between the optical emission and non-thermal $\gamma$-ray emission. Max-Moerbeck et al. (2014a) performed a cross-correlation analysis for the $\gamma$-ray data and 15 GHz radio data of PKS 1502+106 observed between 2008 to 2011 (MJD 54600–55800). They found a $2.25\sigma$ confidence-level positive correlation between the leading $\gamma$-ray emission and the lagging radio emission with a time delay of $(40\pm 13)$ days. However, our results show that there is no significant ($>2\sigma$) correlation between the $\gamma$-ray emission and 15 GHz radio emission for the data observed from 2014 to 2018. This result implies that even for the same source, the correlation between different bands will behave differently at different periods (Böttcher, 2019). During 2008–2011, the delay of radio emission and its associated changes with $\gamma$-ray emission could be explained by a moving disturbance propagating along the jet from inner region to downstream extended radio emission region (Max-Moerbeck et al., 2014a). During 2014–2018, the increase of radio fluxes have no obvious correlation with the $\gamma$-ray or other bands (see Figure 1). At this time, the variability behavior of the radio emission may only be related to the local variation of downstream extended radio emission region. Recent studies have found that in some blazars, their $\gamma$-ray emission is associated with their optical polarization degree and the rotation of PA (see, e.g., Itoh et al., 2016; Blinov et al., 2018). Since the ambiguity of the PA can not be completely eliminated (see Section 3.1 above), here we only analyze the cross-correlation between the polarization degree and the $\gamma$-ray emission. We do not find any significant ($>3\sigma$) correlation between the two physical quantities.

4 $\gamma$-ray dominated outburst in 2015

During the period from May to September 2015, the Fermi-LAT detected a highest-level renewed $\gamma$-ray activity in PKS 1502+106 after the $\gamma$-ray activities in 2008–2009 (see ATel #7592, ATel #7801 as well as Figure 1). The 3-day binning light curves in the 0.1–1 GeV and 1–300 GeV band during this period are shown in Figure 3. It can be clearly seen that a rapid GeV flare occurred in PKS 1502+106 from July 8 to July 17, 2015. The highest flux in the 1–300 GeV band is about six times higher than that in the average state, reaching ($3.2\pm 0.3)\times 10^{-7}$ ph cm*-2* s*-1* (for 3-day binning data). Figure 3 also presents the evolution of the $\gamma$-ray spectral index (0.1–300 GeV) with time. Similar to the GeV flares of other blazars (e.g., 3C 454.3 (Ackermann et al., 2010), OT 081 (Kim et al., 2018), CTA 102 (Gasparyan et al., 2018)), the $\gamma$-ray spectrum hardens when the source becomes brighter. At the peak of the flare, the $\gamma$-ray spectral index reached the hardest $\Gamma_{\gamma}=1.82\pm 0.04$, which is rarely seen in FSRQs. The distribution of high energy events ($>$ 10 GeV) detected from PKS 1502+106 is also presented in Figure 3. These high energy events are extracted using the gtsrcprob tool within the $0.5^{\circ}$ ROI. Most of the high-energy photons were observed during the flare, and their energy mainly ranges from 10 to 20 GeV. Below we will combine the multi-band data to analyze this fast GeV flare to explore the possible origin mechanism of the flare and the physical properties of flaring region.

4.1 The time profile of the flare

High time-resolution light curve can eliminate the degeneracy caused by the superposition of multiple short time-scale variations and reveals more basic physical information. Therefore, on the premise of sufficient photon statistics, we obtain a 3-h binning $\gamma$-ray light curve during the flare (see Figure 4). Simultaneous optical/ultraviolet/X-ray observations and the changes of the X-ray photon index (0.3–10 keV) with time are also shown in this figure. The $\gamma$-ray flare is very fast, and there seems to be a related change in the X-ray band after the $\gamma$-ray flare. We divide the light curve into three phases (Pre-flare, Flare, and Post-flare) as shown in Figure 4.

To quantitatively determine the rising and decaying time-scales of the flare, we fit the time profile of the flare with the following double exponential form function (Abdo et al., 2010b):

[TABLE]

where $F_{c}$ is the quiescent flux; $F_{0}$ is the flux at time $t_{0}$, representing the amplitude of flare; $T_{r}$ and $T_{d}$ determine the rising and decaying time-scales of flare, respectively. The fitting result is shown in Figure 5. The rising and decaying time-scales of the best fitting and their $1\sigma$ uncertainties are $T_{r}=(6.21\pm 1.7)$ hr and $T_{d}=(1.46\pm 0.5)$ hr, respectively. The flare possesses an obviously asymmetric shape (a relatively slow rise followed by a very rapid decay), with the symmetry parameter of $\xi=(T_{d}-T_{r})/(T_{d}+T_{r})=-0.61$. This implies that the injected/accelerated high-energy electrons are quickly cooled or escaped.

Based on the $\gamma$-$\gamma$ opacity argument, the minimum Doppler factor can be estimated numerically from the maximum photon energy $E_{\max}$ detected during the flare (Dondi & Ghisellini, 1995; Poutanen & Stern, 2010). Assuming that the $\gamma$-$\gamma$ interaction optical depth at $E_{\max}$ is one, the minimum Doppler factor can be calculated as

[TABLE]

where $\sigma_{\mathrm{T}}$ and $d_{l}$ are the Thomson scattering cross-section and the luminosity distance, respectively. $t_{var}^{obs}$ is the shortest observed doubling/halving time-scale, and it is approximately equal to $\ln{2}\times T_{d}$(see, e.g, Gasparyan et al., 2018). During the flare, the maximum photon energy $E_{\max}$ is detected to be $\sim 20$ GeV, and the contemporaneous X-ray flux (i.e., $f_{\epsilon}$) is $\sim 1.5\times 10^{-12}$ erg cm*-2* s*-1* (see Figure 4), so that $\delta_{\min}\approx 37$. For blazars, the bulk Lorentz factor $\Gamma$ of radiating blob generally have $\Gamma\sim\delta$, so an upper limit on the viewing angle of the jet can also be estimated: $\theta\lesssim 1/\delta_{\min}=1.5^{\circ}$. Karamanavis et al. (2016a) performed ultra-high angular resolution mm-VLBI observations for PKS 1502+106. The observations show that PKS 1502+106 possesses a parsec-scale jet with the estimated Doppler factors ranging from $\sim 7$ to $\sim 50$ along the jet. The viewing angle differs between the inner and outer jet, with the former at $\theta\sim 3^{\circ}$ and the latter at $\theta\sim 1^{\circ}$. The parameters constrained by the doubling/halving time-scales are compatible with that estimated from the mm-VLBI observations.

Based on the doubling/halving time-scale, the upper limit on the physical size of the flaring region can be given by

[TABLE]

This result indicates that the emission region of this flare is very compact, similar to findings in some low-redshift FSRQs (e.g., PKS 1510-089 (Tavecchio et al., 2010; Prince et al., 2017), CTA 102 (Shukla et al., 2018)). For a conical jet, if the emission is produced across the entire jet area, such a compact emission region implies that the flaring region is very close to the central engine; the upper limit on the distance of the flaring region from the central engine is $d\sim 2R/\theta\sim 2c\delta_{\min}^{2}t_{var}^{obs}/(1+z)=0.03~{}\mathrm{pc}$ (Abdo et al., 2011).

4.2 The evolution of $\gamma$-ray spectrum

To quantify the curvature of $\gamma$-ray spectra at different phases of the flare and search for possible spectral evolution, we perform unbinned likelihood spectral fittings for the $\gamma$-ray spectra at the Pre-flare, Flare, and Post-flare phases with power law (PL), log parabola (LP), broken PL (BPL), and power law with an exponential cutoff (PLEC) function forms. These functions are commonly used to analyze the $\gamma$-ray spectra of blazars (see, e.g., Ackermann et al., 2010; Abdo et al., 2011). The fitting results of the $\gamma$-ray spectra at different phases are presented in Figure 6, and the corresponding fitting parameters are listed in Table 3. As mentioned in Abdo et al. (2011), although the likelihood analysis can return the log-likelihood value along with the fitting parameters, the log-likelihood value does not provide an absolute goodness-of-fit evaluation. Therefore, the reduced $\chi^{2}$ value of each fitting is also calculated and listed in Table 3. Following Nolan et al. (2012), we use the $\mathrm{TS}_{\mathrm{curve}}$ value, which is calculated by $\mathrm{TS}_{\mathrm{curve}}=2(\log\mathcal{L}(\mathrm{LP/BPL/PLEC})-\log\mathcal{L}(\mathrm{PL}))$, to evaluate the significance of spectral curvature.

As shown in Figure 6, at the Flare phase, the $\gamma$-ray spectrum is rather hard, and a curvature/truncation feature emerges. The significance level of the curvature/break is $\mathrm{TS}_{\mathrm{curve}}=10.26$ ($\sim 3.2\sigma$) under the PLEC function form. In all cases, the fitting of the LP function is the worst than that of other functions, and the PLEC function is the best fitting form. There is no obvious evolution in the break/cutoff energy between the three phases and only a slight increase in the break/cutoff energy at the Flare phase. This feature is similar to that found in low-redshift FSRQs (e.g., 3C454.3 (Abdo et al., 2011), PKS 1510-089 (Prince et al., 2017)).

There are four scenarios that are usually used to explain the observed spectral curvature/break (see, e.g., Georganopoulos et al., 2001; Abdo et al., 2009; Poutanen & Stern, 2010; Ackermann et al., 2010; Stern & Poutanen, 2014): 1) the "cooling break" resulting from the radiative losses of high-energy electrons, 2) the attenuation of high-energy photons by the extragalactic background light (EBL), 3) the transition of Compton scattering cross-section from the Thomson regime to the Klein-Nishina (K-N) regime, and 4) the photon-photon pair absorption associated with the He ii Lyman continuum (LyC) or H LyC. For PKS 1502+106, although it has a relatively high redshift (z = 1.838), the EBL at this redshift will only affect the spectra at energies $\gtrsim 20$ GeV in the observation frame (Domínguez et al., 2011). In addition, the typical break energy of spectrum produced by the photon-photon pair absorption associated with the He ii LyC is $\sim 1.8$ GeV in the observation frame, and that associated with the H LyC is $\sim 7$ GeV in the observation frame (Stern & Poutanen, 2014). The spectral break energies observed in the three phases are $\sim 1$ GeV (see BPL fittings in Table 3), which is less than the typical break energies in the above two cases. This implies that the observed spectral curvature/break here is unlikely to be caused by the external photon-photon pair absorption. It can be found that the difference between the high-energy spectral index and the low-energy spectral index obtained by the BPL fitting in the three phases is $\sim 0.5$ (see Table 3), which is consistent with the expectation of the typical "cooling break" (Georganopoulos et al., 2001; Abdo et al., 2009). It is suggested that the observed spectral curvature/break is most likely caused by the radiative losses of high-energy electrons.

4.3 The modeling and evolution of multi-band SED

To further investigate the underlying trigger mechanism of the rapid $\gamma$-ray flare, we model the simultaneous multi-band SEDs at different epochs. The multi-band SEDs at different epochs are shown in Figure 7. Here the SEDs of Post-flare I, Post-flare II and Post-flare III correspond to three different X-ray observations at the Post-flare phase, respectively (see Figure 4); they have the same data in the UV/optical and $\gamma$-ray bands. In each SED, the X-ray data points are the actual observed X-ray data after grouping, and each energy bin has no less than 15 photons. These observed X-ray data points are used to perform SED modeling, while a butterfly plot obtained from the X-ray spectral fitting is displayed in the X-ray band as a reference. The grey data points are archive data collected by the ASDC SED Builder Tool, an online service developed at the ASI Science Data Center (Stratta et al., 2011).111111The archive data in the $\gamma$-ray band are the integral data of Fermi-LAT three-year observations (i.e. three-year binning data), so the source exhibits the $\gamma$-ray detection in the $>50$ GeV band even if the source is in a low state. On the contrary, since the integral time of the $\gamma$-ray data obtained during the flare is very short (2-day bin), even if the source is in a high state, only the upper limits are obtained in the $>50$ GeV band. These archive data are not involved in the SED modeling and are only for reference. At high redshift, the Lyman-alpha forest produced by neutral hydrogen absorption shifts to the low-frequency band. For a neutral hydrogen cloud with $z\sim 1.8$, their Lyman-alpha absorption starts at about 343 nm in the observation frame, which means that the UV photometric data of the Swift-UVOT may suffer from the influence of Lyman-alpha absorption (see, e.g., Bottacini et al. 2010; Rau et al. 2012). Therefore, the UV data in here are considered as the lower limits.

Below we model the SEDs in the framework of a conventional one-zone homogeneous leptonic model (see, e.g., Ghisellini & Tavecchio 2009; Yan et al. 2016; Ding et al. 2017, for detailed model description). It should be noted that lepto-hadronic models can also reproduce the multi-band emission of blazars (e.g., Böttcher, 2012; Böttcher et al., 2013). In the lepto-hadronic models, the emission of X-ray to high-energy $\gamma$-ray is generated by proton-photon interactions (Mannheim, & Biermann, 1992) or relativistic proton synchrotron radiation (Aharonian, 2000). Typically, however, the lepto-hadronic models are difficult to account for fast $\gamma$-ray flare owing to the long cooling time scales of protons (see Böttcher, 2012). In addition, for FSRQs, the lepto-hadronic models usually require an extremely high (super-Eddington) proton power (see, e.g., Sikora et al. 2009; Zdziarski, & Böttcher 2015). In view of these, we will not explore the lepto-hadronic scenario in detail. In the leptonic model, we assume that multi-band emission is produced by relativistic electrons in a spherical blob with radius $R_{\mathrm{size}}$ and filled with a uniform magnetic field $B$. The radiating blob moves in the direction close to the line of sight, and the bulk Lorentz factor has $\Gamma\sim\delta$. We assume that the electron energy distribution is a commonly used broken power-law form as follows:

[TABLE]

where $N_{0}$ is the number of emitting particles per unit volume; $\gamma_{\min}$, $\gamma_{\mathrm{br}}$, and $\gamma_{\max}$ are the electron Lorentz factors for minimum, break, and maximum, respectively; $p_{1}$ and $p_{2}$ are indices of the power law below and above the break energy. In the leptonic model, the observed high-energy emission is produced by the inverse Compton (IC) scattering of relativistic electrons. The seed photons for IC process could be from the local synchrotron emission (i.e., synchrotron self-Compton, SSC) or from external fields (EC), such as broad-line region (BLR) and dusty torus (DT). As shown in Figure 7, in most cases, the shapes of the energy spectra in the X-ray band and the $\gamma$-ray band are similar. Such SED can not be reproduced well with SSC model only and needs the involvement of an EC component. The type of soft photons in the EC process depends on the location of the emission region. For PKS 1502+104, its monochromatic luminosity at 1350 Å is $L_{1350}=7.8\times 10^{46}$ erg s*-1* cm*-2* (Shaw et al., 2012). According to the C IV radius–$L_{1350}$ relation obtained from the reverberation mapping of luminous quasars at high redshift (Lira et al., 2018), the BLR radius of PKS 1502+106 is estimated to be $R_{\mathrm{BLR}}\approx 0.11$ pc. This size is larger than the upper limit on the distance of the emission region from the central engine (see Section 4.1), suggesting that the photons coming from the BLR are dominant in the EC process. The central black hole mass ($M_{\mathrm{BH}}=7.9\times 10^{8}~{}M_{\odot}$) and the BLR luminosity ($L_{\mathrm{BLR}}=1.47\times 10^{45}$ erg s*-1* cm*-2*) of PKS 1502+106 are estimated by Xiong & Zhang (2014) based on the C IV emission line. Following Ghisellini et al. (2010), the luminosity of accretion disk can be estimated as $L_{\mathrm{disk}}\approx 10\times L_{\mathrm{BLR}}=1.47\times 10^{46}$ erg s*-1* cm*-2*. Based on these values, the radiation spectrum from a standard thin disk is calculated and shown in Figure 7 (a) Shakura & Sunyaev (1973). It can be seen that even in the UV/optical band, the contribution of the emission from the disk is much smaller than the observed SED, indicating the observed SEDs are dominated by non-thermal radiation. For these reasons, we finally adopt SSC+EC (target photons coming from BLR) model to fit observed SEDs. The BLR photon density in the comoving/jet frame can be calculated as $U_{\mathrm{BLR}}=\frac{\Gamma^{2}L_{\mathrm{BLR}}}{4\pi cR_{\mathrm{BLR}}^{2}}\approx 0.03\delta^{2}~{}\mathrm{erg~{}cm}^{-3}$, and a effective BLR temperature $T_{\mathrm{BLR}}\sim 5\times 10^{4}$ K is adopted in the model. In addition, based on the EBL model in Domínguez et al. (2011), we consider the effect of EBL absorption in the GeV band in the model calculation.

There are nine free parameters in the model. Six of them specify the electron energy distribution ($N_{0}$, $\gamma_{\min}$, $\gamma_{\max}$, $\gamma_{\mathrm{br}}$, $p_{1}$, and $p_{2}$), and other three describe the properties of the emission region ($\delta$, $R_{\mathrm{size}}$ and $B$). The Doppler factor, the minimum and maximum energies of electrons are difficult to be constrained well by the current SED data. Therefore, in the model, we adopt the minimum Doppler factor of $\delta_{\min}=37$ and an appropriate minimum electron energy of $\gamma_{\min}=20$, which is constrained by the archival radio data and the observed X-ray data in some degree. In addition, based on the constraint of SED shape, the values of $\gamma_{\max}=5\times 10^{3}$ and $p_{2}=2.5$ are adopted in the model to further reduce the number of free parameters. The final model free parameters are reduced to five. We employ the Levenberg-Marquardt algorithm, which is provided by the Python lmfit package121212https://lmfit.github.io/lmfit-py/, to obtain the best-fitting parameters of the model for each observed SED. The best modeled SEDs are displayed in Figure 7, and the corresponding model parameters are listed in Table 4. The powers carried by the jet in the form of radiation ($P_{r}$), magnetic field ($P_{B}$), electrons ($P_{e}$) and cold protons ($P_{p}$, assuming one proton per emitting electron) are also given in Table 4. They are calculated as

[TABLE]

where $U_{i}$ is the energy density of the $i$-th component (see, e.g., Celotti & Ghisellini, 2008; Ding et al., 2017; Xue et al., 2019).

The results show that the radiative power of the jet ($P_{r}$) is higher than the kinetic power carried in relativistic electrons ($P_{e}$) and the power in Poynting flux ($P_{B}$), which is consistent with the results of typical high-power FSRQs (Celotti & Ghisellini, 2008). In the framework of the leptonic model, this result implies that the jet may be mainly loaded by hadrons, namely dynamically dominated by the bulk motion of cold protons, as both leptons and Poynting flux do not provide sufficient power to account for the observed emission (see Celotti & Ghisellini 2008; Ghisellini et al. 2010 for detailed discussions). In addition, the power in Poynting flux is roughly comparable to the kinetic power carried in relativistic electrons, and the magnetization parameter ($P_{B}/P_{e}$) ranges from 0.2 to 1.5, which indicates that the magnetic field energy and electronic kinetic energy of the jet are basically in equipartition.

The sizes of emission region obtained by the SED fittings range $(3-8)\times 10^{15}$ cm, which is slightly larger than that constrained by the doubling/halving time-scale. Considering the uncertainties in the time-scale measurement and the SED fittings, the sizes obtained by the two methods are basically consistent. The values of the magnetic field intensity $B$ and the electronic spectral index $p_{1}$ at different epochs are shown in Figure 8. $p_{1}$ shows a soft-hard-soft change. It should be pointed out that changing the values of $\delta$, $p_{2}$, $\gamma_{\min}$, and $\gamma_{\max}$ will not substantially affect the fitting results of $p_{1}$. In the Flare phase, the magnetic field intensity seems to decrease slightly, and then a harder electron spectrum with $p_{1}=1.07\pm 0.53$ appears in the Post-flare I epoch. Such a hard electronic energy spectrum usually can not be produced in the standard diffusive shock acceleration models, which generally forms an electronic energy spectrum with the electronic spectral index close to 2.0 (see, e.g., Drury, 1983; Bednarz & Ostrowski, 1998; Rieger et al., 2007). As thus, the fast $\gamma$-ray flare is probably not caused by shocks. Note that an electronic spectral index of 1.5 is generally considered as the lower limit still in agreement with the shock acceleration scenarios (e.g., Aharonian et al., 2006). Considering the uncertainty of $p_{1}$, the shock acceleration scenarios can not be eliminated completely in here. However, magnetic reconnection may be a more promising trigger for this fast gamma-ray flare. Three-dimensional numerical simulations show that the relativistic turbulence in AGN jet may trigger magnetic reconnection, which could further drive strong stochastic acceleration and form a hard electronic energy spectrum with an electronic spectral index of $\sim 1$ (see Guo et al., 2014, 2015).

5 Optical dominated outburst in 2017

As mentioned in Section 3.1, between May 2017 to March 2018, a significant optical dominant outburst occurred in PKS 1502+106. At the peak of the outburst, the long-term optical monitoring project at Steward Observatory monitored the highest optical polarization degree observed so far from PKS 1502+106, reaching $(47.4\pm 0.1)\%$, which is also one of the highest polarization levels ever observed for blazars (ATel #11047).

The distributions of optical polarization degree versus optical V-band fluxes and optical polarization angles are shown in Figure 9. Here we use 1-day time interval to cross-match the data of photometry and polarization measurements, and finally 153 data points are obtained. The data observed during the outburst are marked in red. There is a very tight positive correlation between the optical polarization degree and the optical fluxes during the outburst, with a spearman correlation coefficient of 0.82 (the chance probability $P=10^{-11}$). In addition, the optical polarization angles are within a narrow range during the outburst, with an average angle of $89.6^{\circ}$, which is close to the parsec-scale jet direction $\sim 81^{\circ}$ shown by the VLBI imaging (Karamanavis et al., 2016a). In non-outburst state (grey data points), the optical polarization degree and the polarization angles exhibit random distributions (see also Figure 10).

Kiehlmann et al. (2016) used a simple random walk polarization variability model to analyze the long-term polarization observation data of FSRQ 3C279. They found that the polarization variation is possibly dominated by a stochastic process during the low-brightness state, while during the outburst, the polarization variation is governed by a deterministic process (e.g., shock compression, non-axisymmetric magnetic field configuration, etc). This is consistent with what we see in PKS 1502+106. During the outburst, the polarization behaviors of PKS 1502+106 are in agreement with the expectations of the shock-in-jet model, where the magnetic field is compressed and aligned at the front of the shock so that the polarization degree and fluxes exhibit a close positive correlation (see, e.g., Hagen-Thorn et al., 2008). Interestingly, compared to outbursts with similar behaviors in other blazars discovered in previous studies (e.g., AO 0235+164 (Hagen-Thorn et al., 2008), 1ES 1959+650 (Sorcia et al., 2013)), the outburst in PKS 1502+106 lasted for a long time ($\sim$ one year) and has extremely high polarization degree.

5.1 Two-component decomposition

The polarization observations of Steward Observatory provide normalized (relative) Stokes parameters $q$ and $u$. The $q$–$u$ plane built from the data is shown in Figure 10, where the red data points represent the observed data during the outburst and the black star ($q_{c}=-0.106$ and $u_{c}=-0.068$) represents the average central point of $q$–$u$ calculated by iteratively discarding $>3\sigma$ outliers. The average central point are offset from the origin, and most of the red data points deviate from the average central point in the same direction. Therefore, we infer that the polarization of PKS 1502+106 could be composed of a constant or stable component associated with the jet configuration and a variable component that is related to the propagation of the shock (Holmes et al., 1984; Hagen-Thorn et al., 2008).

Based on the average central point, the polarization degree and polarization angle of the stable component are calculated as $P_{c}=\sqrt{q_{c}^{2}+u_{c}^{2}}=12.6\%$ and $\Theta_{c}=\frac{1}{2}\arctan(u_{c}/q_{c})=106.3^{\circ}$, respectively. For the variable component, following the method proposed by Hagen-Thorn et al. (2008), we use absolute Stokes parameter planes, i.e., $Q$–$I$ and $U$–$I$ planes, to determine its variability behaviors, where the absolute Stokes parameters are defined as $(Q,U)=(q,u)I$. If the variability of the variable composition is only caused by its flux variations, independent of its polarization degree and polarization angle, the relative Stokes parameters ($q_{v}$ and $u_{v}$) of the variable composition will remain unchanged. In this case, observed data points will lie on straight lines in the absolute Stokes parameter planes ($Q$–$I$ and $U$–$I$ planes), and the slopes ($r_{Q-I}$ and $r_{U-I}$) of these lines reflect the relative Stokes parameters of the variable component (i.e., $q_{var}=r_{Q-I}$ and $u_{var}=r_{U-I}$). The $Q$–$I$ and $U$–$I$ planes built from our observed data are shown in Figure 11, where $\alpha$ is a proportional term between $I$ and the V-band fluxes. In the $Q$–$I$ plane, the data points during the outburst (red points) show a very tight linear dependence, indicating that there is no change in the relative Stokes parameter $q_{v}$ during the outburst. In the $U$–$I$ plane, the data points during the outburst display three-stage linear relations. The blue, pink, and green data points correspond to the rising, peak and declining stages of the outburst, respectively. We perform linear fittings separately for the $U$–$I$ relations at the three different stages while performing only one linear fitting for $Q$–$I$ relation. The fitting results and the polarization parameters of the corresponding variable component are listed in Table 5. The variable component shows extremely high polarization degree, up to ($52.8\pm 1.7)\%$ in the peak stage. The polarization angle of the variable component is close to the parsec-scale jet direction at the rising and declining stages, and it has a change of $\sim 20^{\circ}$ at the peak stages.

5.2 The properties of the shock

Based on the variability behaviors seen above, below we briefly discuss the properties of the shock. As observed, this dramatic outburst seen in PKS 1502+106 occurred mainly in the optical band, suggesting that the emission in the shock zone is dominated by synchrotron radiation. In addition, the polarization angle of the variable component is close to the parsec-scale jet direction, which implies that the shock could be a transverse shock; the transverse shock orders the turbulent magnetic field along the front, which is perpendicular to the jet direction (Hagen-Thorn et al., 2008). In this scenario, the polarization degree depends on the viewing angle of the shock $\Psi$, the spectral index $\alpha_{\mathrm{o}}$ in optical band, and the ratio of densities of the shocked region to the unshocked region (compression ratio) $\eta=n_{\mathrm{shock}}/n_{\mathrm{unshock}}$ (Hughes & Miller, 1991):

[TABLE]

The viewing angle of the shock in the observer’s frame is subjected to relativistic aberration and determined by the bulk Lorentz factor $\Gamma$ and the viewing angle of the jet $\theta$:

[TABLE]

According to Eq. 5, the functions of the viewing angle of the shock and the compression ratio under the different polarization degrees are displayed in Figure 12. Meanwhile, according to Eq. 6 and the relation of $\delta=[\Gamma(1-(1-\Gamma^{-2})^{1/2}\cos\theta)]^{-1}$, the values of the Doppler factor corresponding to the viewing angles of the shock are marked in the figure. Here, we use the value of $\alpha_{\mathrm{o}}=1.8$ calculated from the average optical data during the outburst and the value of $\theta\approx 3^{\circ}$ constrained by the mm-VLBI observations (Karamanavis et al., 2016a). It can be seen that in a wide range of the Doppler factor, the polarization degree of the variable component wants to reach observed $\sim 50\%$, the compression ratio of the shock must be no less than $2.2$.

The time-scale of outburst relates to the thickness of the shock front, which is determined by the lifetime of relativistic electrons accelerated at the front. In the current situation dominated by synchrotron radiation, the lifetime of relativistic electrons in the observer frame is

[TABLE]

where $\nu_{\mathrm{GHz}}$ is observed synchrotron photon frequency in GHz (Hagen-Thorn et al., 2008). The functions of the magnetic field intensity and the cooling time under the different Doppler factors are displayed in Figure 13. As in Section 4.1, we use a double exponential function to fit the time profile of the outburst in the V band. The fitting result is shown in the inset of Figure 13. The rising and decaying time-scales of this outburst are $T_{r}=(73\pm 10)$ d and $T_{d}=(49\pm 5)$ d, and the corresponding fastest variability time (i.e., shortest doubling/halving time-scale) is $t_{\mathrm{var}}^{\mathrm{obs}}=(31.9\pm 2.7)$ d (marked as black line in Figure 13). In the case where the cooling time is equivalent to the fastest variability time, the magnetic field intensity is required to be $\sim 0.032$ G and weakly depends on the Doppler factor. Karamanavis et al. (2016b) performed a cross-correlation analysis for the radio light curves observed during the $\gamma$-ray outburst in 2008. Based on the observed time-delays between different bands, the structure of PKS 1502+106 in terms of synchrotron opacity was deduced, while the magnetic field intensity along the jet axis also was estimated (under the shock-in-jet scenario). They estimated the magnetic field intensity in the radio nucleus to be between 14 and 176 mG. Our estimated magnetic field intensity in the shock emission region is in harmony with their result. However, this magnetic field intensity is much smaller than that obtained in the $\gamma$-ray flare in 2015, which implies that the properties of the emission regions of the activities triggered by shock and magnetic reconnection are quite different.

6 Summary

After the $\gamma$-ray activities in 2008–2009, PKS 1502+106 entered a quiescent stage of up to six years. Until mid-2015, its prominent multi-band activities are re-detected. Using multi-band data from radio to $\gamma$-ray bands as well as optical polarization observations, we systematically explore the multiwavelength activities of high-redshift FSRQ PKS 1502+106 during 2014–2018. Two dramatic outbursts, a $\gamma$-ray dominated outburst in 2015 and an optical dominated outburst in 2017, are investigated in detail to explore the triggering mechanism of them and the physical properties of the emission regions. The main results are summarized as follows.

1. An hour-scale GeV $\gamma$-ray flare is discovered during the $\gamma$-ray dominated outburst in 2015. This fast flare shows obviously asymmetric time profile, and the $\gamma-$ray spectral index is $\Gamma_{\gamma}=1.82\pm 0.04$ at the peak of the flare, which is rarely seen in FSRQs. Based on the variability time-scale of the flare, the physical parameters of flaring region (e.g, minimum Doppler factor, emission region size, etc.) are constrained. See Section 4.1.

2. The $\gamma$-ray spectra at the different phases of the flare can be best fitted by the PLEC function. At the Flare phase, $\gamma$-ray spectrum emerges a curvature/break characteristic ($\sim 3.2\sigma$). The curvature/break characteristic is in line with the expectation of typical "cooling break", suggesting that it is most likely caused by the radiative losses of high-energy electrons. See section 4.2.

3. Based on a one-zone homogeneous leptonic model, multi-band SEDs at different epochs of the $\gamma$-ray flare are modeled. The results show that the multi-band radiation of PKS 1502+106 needs the involvement of EC process, and the soft photons in the EC process should mainly come from the BLR, which further confirms the result in Abdo et al. (2010a). In addition, SED modelings reveal the changes of the electron spectral index and magnetic field intensity in the emission region at different epochs of the flare. In the Flare phase, the magnetic field intensity seems to decrease slightly, and then a harder electron spectrum with the electronic spectral index of $p_{1}=1.07\pm 0.53$ appears. This result may imply that the fast $\gamma$-ray flare is generated by magnetic reconnection. See section 4.3.

The above conclusions are obtained in the framework of the leptonic scenario. In view of typical limitations in the lepto-hadronic models (see Section 4.3), the lepto-hadronic models are not further explored in our work. Nevertheless, it should be noted that the lepto-hadronic scenario may still a potential alternative. For example, some studies on the giant $\gamma$-ray flare of 3C 279 in June 2015 (similar to the time scale of the flare in PKS 1502+106) show that the lepto-hadronic scenario still has the potential to explain such fast $\gamma$-ray flare under certain conditions (see, e.g., Petropoulou et al., 2017; H. E. S. S. Collaboration et al., 2019). Whether the emission origin of blazars is leptonic scenario or lepto-hadronic scenario is still a controversial issue, which needs to be further explored in the future.

4. An optical dominated outburst occurred in 2017. During the outburst, the optical polarization degree and optical fluxes exhibit a very significant correlation. By analyzing Stokes parameters of polarization observations, our results show that the observed polarization could be composed of a stable component associated with the jet configuration and a variable component that is related to the propagation of shock. The polarization degree of the variable component is as high as $(52.8\pm 1.7)\%$ at the peak of the outburst. The outburst could be triggered by a transverse shock with a compression ratio of $\eta>2.2$, and the magnetic field intensity of the shock emission region is about $0.032$ G. See Section 5.

In PKS 1502+106, we see that both shock and magnetic reconnection may be the triggers of multi-band activities, and the multi-band activities triggered by the two may be significantly different. Perhaps the short-time high-energy flares are more likely to be triggered by magnetic reconnection, while the long-term outbursts in a low-energy band dominated by synchrotron radiation are more likely to be correlated with shocks. At present, there are few reports of flares triggered by magnetic reconnection. This may be attributed to the need for multi-band synergetic observation data to identify such flares. In addition, such flares perhaps have short time-scale, which makes it rare for such events to have multi-band simultaneous observations, thereby reducing the recognition rate of such flares. The coming era of multi-messenger time-domain astronomy will provide more high-quality multi-band synergetic observation data. By then, the studies of the multi-band activities for blazars will reveal a more complete picture of the energy dissipation mechanism in the jet (see, e.g., Rani et al., 2019; Burns et al., 2019).

We sincerely thank the anonymous referee for helpful suggestions. We acknowledge financial support from the National Key R&D Program of China grant 2017YFA0402703 (N.D., Q.S.Gu) and National Natural Science Foundation of China grant 11733002 (N.D., Q.S.Gu). Dingrong Xiong acknowledges financial support from the National Natural Science Foundation of China grant 11703078. This work has made use of Fermi data, obtained from the Fermi Science Support Center, provided by NASA’s Goddard Space Flight Center (GSFC). The data, software, and web tools obtained from NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC), a service of GSFC, were used. This work has made use of the Swift and XMM-Newton data. Also, data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. We also has made use of data from the OVRO 40 m monitoring program, which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.