Fast γ-Ray Variability in Blazars beyond Redshift 3

Blazars, including flat spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs), are radio-loud active galactic nuclei (AGNs) whose relativistic jets are closely pointing to our line of sight (Blandford & Rees 1978; Urry & Padovani 1995). The jet-viewing angle θ_v is smaller than or comparable to the jet-beaming angle (1/Γ, where Γ is the jet bulk Lorentz factor), hence the jet emission is strongly boosted because of relativistic effects, making blazars so luminous that they are visible even at very high redshifts (e.g., Romani et al. 2004; Yi et al. 2014). High-redshift blazars are extremely powerful phenomena harboring supermassive black holes (SMBHs) heavier than one billion solar masses (e.g., Ghisellini et al. 2010). These sources are of major importance because they shed light on the formation and growth of the first generation of SMBHs and the role that relativistic jets played at that time (e.g., Volonteri 2010; Afonso et al. 2015). Studies of high-redshift blazars also provide information on whether and how the jet properties change with cosmic time, along with their potential impact on the evolution of AGNs, as well as their host galaxies (e.g., Volonteri et al. 2011; Fabian 2012; Morganti et al. 2013).

Since the relativistic jets are responsible for bright radiation of blazars, they are characterized by highly variable and polarized non-thermal continuum emission, and are generally detected in all observable bands, from radio to γ-ray (both GeV and TeV) regimes (Ulrich et al. 1997). Their spectral energy distribution (SED) reveals a two-bump shape, where one is likely due to synchrotron emission of relativistic electrons in magnetic fields and the other extending to γ-rays is usually explained as inverse Compton scattering of soft photons by the same population of relativistic electrons (e.g., Maraschi et al. 1992; Dermer & Schlickeiser 1993; Sikora et al. 1994; Błażejowski et al. 2000). One of the more extraordinary phenomena in blazars is their fast γ-ray variability. Imaging atmospheric Cerenkov telescopes (IACTs) like HESS and MAGIC have detected flux changes on a timescale of a few minutes in high-frequency peaked BL Lacs PKS 2155−304 (Aharonian et al. 2007) and Mkn 501 (Albert et al. 2007). Similar extreme behaviors have been also detected in other subtypes of jetted AGNs, including the low-frequency peaked BL Lac BL Lacertae (Arlen et al. 2013) and FSRQ PKS 1222+21 (Aleksić et al. 2011), strongly suggesting it is a common feature independent of the source type. In the GeV domain, flux variations on a timescale of several hours have been found in a few sources since the CGRO era (e.g., Mattox et al. 1997; Wehrle et al. 1998). The number of such detections has significantly increased thanks to Fermi, though the minimum variability timescale is generally limited to a few hours due to its routine survey observation mode (e.g., Foschini et al. 2011; Vovk & Neronov 2013; Liao & Bai 2015; Liao et al. 2016, but also see Foschini et al. 2013). Recently, taking advantage of a target of opportunity (ToO) repointing of the Fermi Large Area Telescope (LAT, Atwood et al. 2009), extremely fast GeV γ-ray flux variability with a doubling time of less than 5 minutes has been detected in the giant outburst 3C 279 (Ackermann et al. 2016). Fast γ-ray variability of blazars is important for investigating the speed, composition, and energetics of relativistic jets (e.g., Begelman et al. 2008).

The detected number of high-redshift blazar candidates (i.e., z  4) has significantly increased over the last decade, thanks to unprecedented powerful X-ray observatories (e.g., Swift; Gehrels et al. 2004), as well as wide-field surveys in the optical/UV (e.g., the Sloan Digital Sky Survey, SDSS; York et al. 2000), infrared (WISE; Wright et al. 2010), and radio bands (e.g., the NRAO VLA Sky Survey, NVSS; Condon et al. 1998). So far the farthest known blazar candidate is Q0906+6930 (z = 5.48; Romani et al. 2004), along with several others beyond redshift 5, including SDSS J102623.61+254259.5 (z = 5.2; Sbarrato et al. 2012), SDSS J013127.34−032100.1 (z = 5.18; Yi et al. 2014), and SDSS J114657.79+403708.6 (z = 5.0; Ghisellini et al. 2014). These sources are identified as blazars because of their high radio loudness and hard X-ray spectra (e.g., Romani 2006). From a γ-ray perspective, all γ-ray AGNs beyond redshift 3 are FSRQs and PKS 0537−286 (z = 3.1; Wright et al. 1978) stood as the farthest γ-ray blazar for a long time (Abdo et al. 2010a). Recently, 5 new γ-ray blazars with redshifts over 3.3 have been identified by Fermi-LAT Collaboration based on spatial overlap between the γ-ray sources and their optical counterparts, as well as their typical γ-ray FSRQ shape SEDs (i.e., high Compton dominance; Ackermann et al. 2017). However, variability information of these sources is still lacking.

Here, we perform a detailed analysis of Fermi-LAT γ-ray data of all seven known γ-ray blazars beyond redshift 3. This work is organized as follows. In Section 2, strategies in the Fermi-LAT data analysis procedure are introduced. The results of the analysis are reported in Section 3, in which we mainly focus on the results in the temporal domain. Finally, in Section 4 we summarize our results with some discussions.

In this paper, we used publicly released Fermi-LAT Pass 8 data (P8R2_SOURCE_V6, FRONT+BACK) and updated Fermi Science-Tools package of version v10r0p5 to perform the analysis. For each target, we considered the data set within a 10° region of interest (ROI) from 2008 October 27 to 2017 June 12 (i.e., Mission Elapsed Time, MET, between 246,823,875 and 518,983,875 s), along with the energy range from 100 MeV to 100 GeV. In order to reduce contamination from the Earth's limb, we removed γ-ray events with a zenith angle larger than 90°. The recommended quality-filter cuts (i.e., DATA_QUAL==1 && LAT_CONFIG==1) have been followed to ensure that the spacecraft keeps in good condition and hence the data set is valid for scientific use.

First, we used the script make3FGLxml.py⁴ to generate an initial background model for each target, which includes all point sources within 15° in the third Fermi-LAT source catalog (3FGL; Acero et al. 2015), as well as the latest galactic diffuse γ-ray emission model gll_iem_v06.fits and isotropic emission template iso_P8R2_SOURCE_V6_v06.txt⁵ (Acero et al. 2016). γ-ray locations and spectral templates (i.e., power-law function, dN/dE ∝ E^−Γ, where Γ is the spectral index) of the targets were set as the same as those in the literature (Acero et al. 2015; Ackermann et al. 2017). Then, an unbinned likelihood algorithm implemented in the gtlike task was adopted to extract the flux and spectrum. During the extraction, the flux and spectral parameters of sources within the 10° ROI, together with the normalization factors of the two diffuse components, were set free. The significance of the target is quantified with the test statistic (TS), defined as TS = − 2ln(L₀/L) (Mattox et al. 1996), where L and L₀ are the maximum likelihood values for the model with and without target source, respectively. Since 3FGL is based on 4-year Fermi-LAT data that has a narrower data time range than ours, we have checked whether there are newly emerging γ-ray sources beyond 3FGL by generating a residual TS map for each ROI. All new γ-ray background sources corresponding to excesses with TS > 25 (i.e., detection significance of 4.2σ) in the TS maps were added into the background model. For these new sources, their spectral models were also set to be power-law functions, and γ-ray positions were obtained by the gtfindsrc task. After all these steps, the final results were obtained by analyzing the updated background model.

In the temporal analysis, for background sources neither bright (i.e., TS value twice that of the target) nor close (<4°) to the targets, their spectral indexes were frozen to values from the global fit. Meanwhile, considering high-Galactic-latitude background point sources that are dominated by blazars with highly variable γ-ray emissions, we removed the weak background sources (i.e., TS < 1) from the source model. When a special γ-ray flare event was focused, first, an analysis of the entire flare epoch (i.e., several tens of days) was performed. Then, short-term light curves were extracted to search for evidence of intraday γ-ray variability. During the extraction, normalizations of two diffuse emission components, and spectral parameters for all background sources, were fixed to values from the analysis covering the whole flaring period. Note that this strategy has also been adopted for several similar studies (e.g., Saito et al. 2013; Liao & Bai 2015; Ackermann et al. 2016; Liao et al. 2016). Meanwhile, we also took into account the dependence of the effect area of Fermi-LAT due to its square shape, and then set phibins = 5 in the gtltcube task. When the TS value of the target is smaller than 4, the pyLikelihood UpperLimits tool was adopted to calculate the 95% upper limit of the flux.

3.1. Global Analysis

There are in total seven known γ-ray blazars beyond redshift 3 so far. Their basic information is listed in Table 1. Two sources (i.e., NVSS J053954−283956 and NVSS J080518+614423) are known as relatively strong γ-ray sources and have been included in Fermi source catalogs (e.g., Abdo et al. 2010b; Acero et al. 2015). γ-ray emissions for the other five sources have been detected from a recent search of high-redshift γ-ray blazars performed by the Fermi collaboration (Ackermann et al. 2017). The redshift range of these seven sources is between 3.0 and 4.3, especially NVSS J151002+570243, which is located beyond redshift 4. Except for NVSS J053954−283956, whose listed SMBH mass was obtained from modeling its big blue bump from the accretion disk (Ghisellini et al. 2010), all other estimations are based on optical spectroscopic observations (Torrealba et al. 2012; Alam et al. 2015).

Table 1.  Basic Information and γ-ray Properties of the Known 7 γ-ray Blazars beyond Redshift 3

NVSS Name	b	z	F1.4 GHz	Rmag	MBH, ⊙	Fγ	Γγ	TS	σvar
J053954−283956	−273	3.1	0.86	19.0	9.3	4.92 ± 0.23	2.80 ± 0.04	1737.8	15.7
J064632+445116	175	3.4	0.45	18.5	9.1	1.35 ± 0.19	2.95 ± 0.13	103.3	1.2
J080518+614423	324	3.0	0.83	19.6	9.07	2.21 ± 0.13	2.83 ± 0.05	467.2	10.8
J135406−020603	501	3.7	0.73	19.2	8.9	1.0 ± 0.16	2.83 ± 0.13	60.1	5.0
J151002+570243	503	4.3	0.20	19.9	8.5	0.42 ± 0.13	2.72 ± 0.20	32.7	3.2
J163547+362930	421	3.6	0.15	20.6	8.7	1.93 ± 0.29	3.21 ± 0.12	129.9	6.2
J212912−153841	−419	3.3	0.59	16.5	9.8	1.62 ± 0.19	3.04 ± 0.12	108.9	2.4

Note. (1) NVSS name of the object; (2) Galactic latitude; (3) redshift; (4) NVSS radio flux density at 1.4 GHz in Jy; (5) apparent R band magnitude; (6) logarithm of black hole mass; (7) average γ-ray flux of a 105-month LAT data analysis on a scale of 10⁻⁸ ph cm⁻² s⁻¹; (8) γ-ray spectral index corresponding to column (7); (9) TS value corresponding to column (7); (10) significance of the variability estimated from the monthly γ-ray light curve. The Galactic latitudes are derived from NED. The redshifts and the R magnitudes are obtained from the Half Million Quasars catalog (Flesch 2015). The NVSS 1.4 GHz flux densities are from Condon et al. (1998). Note that the ${\sigma }_{\mathrm{var}}$ here for NVSS J163547+362930 could be influenced by the neighbor 3FGL J1635.2+3809. For more details see Section 3.2.1.

Download table as:  ASCIITypeset image

Before performing the temporal analysis, a fit of the entire ~105-month LAT data for each source was performed. All seven sources are found to be significant γ-ray emitters with soft spectra (i.e., Γ_γ > 2.7), which is consistent with the literature (Acero et al. 2015; Ackermann et al. 2017). Averaged γ-ray fluxes and corresponding γ-ray spectral indexes and TS values are also provided in Table 1. NVSS J053954−283956 is the brightest one among these sources, and it has the highest TS value. Note that for sources selected from Ackermann et al. (2017), our averaged fluxes are generally lower than their values, because we chose a different energy range of LAT data (i.e., 0.1–100 GeV) from theirs (i.e., 0.06–300 GeV). For sources with TS > 100, a sophisticated spectral model (i.e., Log-parabolic function) is adopted to fit all the LAT data, but no significant improvement is found compared to the initial usage of a power-law function.

3.2. Temporal Behaviors

3.2.1. Monthly γ-ray Light Curves

Since the majority of our sources are not bright in γ-rays (i.e., TS_{105 month}  100), firstly we evenly divide the total LAT data into 21 time bins (each with a bin of about 5 months) to extract a γ-ray light curve for each source; see Figure 1. For NVSS J053954−283956 and NVSS J080518+614423, significant γ-ray variability is apparent, confirming the results in Fermi source catalogs (e.g., Nolan et al. 2012). For the other five sources, their error bars are relatively large and there are many upper limits in the light curves. So we use the "variability index" test (Nolan et al. 2012) to quantify the significance of the variability of the light curves, from which information on the upper limits can be properly considered. The null hypothesis of this test is that source flux is constant across the full time period. The variability index is derived using the following expression in Nolan et al. (2012):

$$\begin{eqnarray}&&{\mathrm{TS}}_{\mathrm{var}}=-2\,\displaystyle \sum _{i}\displaystyle \frac{{\rm{\Delta }}{F}_{i}^{2}}{{\rm{\Delta }}{F}_{i}^{2}+{f}^{2}{F}_{\mathrm{const}}^{2}}\mathrm{ln}\left(\displaystyle \frac{{{ \mathcal L }}_{i}({F}_{\mathrm{const}})}{{{ \mathcal L }}_{i}({F}_{i})}\right),\end{eqnarray} \tag{ 1 }$$

where for the ith time bin, ${{ \mathcal L }}_{i}$ is the likelihood value, F_i is the observed photon flux, ΔF_i is the statistical uncertainty of F_i, F_const is the assumed constant flux, and f is the systematic correction factor, from which we take a value of 2% following Nolan et al. (2012). In our analysis, the optimized constant flux for each source is close to the average flux from the global analysis (within 1σ statistical uncertainty). If the null hypothesis is correct, the derived ${\mathrm{TS}}_{\mathrm{var}}$ in case of a light curve with 21 time bins should follow a χ² distribution with 20 degrees of freedom and hence the variability significance in σ units is obtained (also listed in Table 1). It is not surprising that the σ_var values of NVSS J053954−283956 and NVSS J080518+614423 are high (i.e. >10). Interestingly, we find that the γ-ray emissions of NVSS J135406−020603, NVSS J151002+570243, and NVSS J163547+362930 are significantly variable (i.e., σ_var > 3). Due to the limited angular resolution of Fermi-LAT, strong variability in nearby background sources can cause artificial variability for the target. Therefore, we have checked whether there are any such neighbors around these five sources. We find that this situation only happens to NVSS J163547+362930, for which there is one bright and highly variable background source: 3FGL J1635.2+3809, which is about 16 away. Since the 68% C.L. contamination angle of LAT for 500 MeV photons is about 15, to avoid significant impact from the neighbor, individual light curves with a lower-energy cut of LAT data from 100 to 500 MeV, both for the target and the neighbor, are extracted (see Figure 2). There are two major γ-ray flares in the >100 MeV light curve of NVSS J163547+362930, with peaking times around MJD 56192 and MJD 56792, respectively. The former flare coincides with a high flux state of the neighbor 3FGL J1635.2+3809 and disappears in the >500 MeV light curve, indicating that it is probably artificial. However, the other flare corresponds to a low flux state of the neighbor and remains to significant in the >500 MeV light curve, which suggests an intrinsic link between this flare and the target. The variability index of the >500 MeV light curve of NVSS J163547+362930 is also calculated, given as σ_var = 4.9. Although it is smaller than that in >100 MeV case, the γ-ray emission of NVSS J163547+362930 is still proven to be significantly variable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. Detecting Fast γ-ray Variability

According to the monthly γ-ray light curves, fluxes of several time bins are significantly higher than the averaged flux. Together with relatively large TS values (i.e., >50), this allows us to further search for any possible fast γ-ray variability. Except for NVSS J135406-020603, for which no significant short-term variability is found, detailed temporal analyses of NVSS J053954−283956, NVSS J080518+614423, and NVSS J163547+362930 are described below.

3.2.2.1. NVSS J053954−283956

NVSS J053954−283956, also known as PKS 0537−286, is one of the most luminous high-redshift quasars (z = 3.104, Wright et al. 1978). Its first detection was at radio frequencies (e.g., Bolton et al. 1975). It is also a bright and well-studied source in X-rays (e.g., Zamorani et al. 1981; Sambruna et al. 2007; Bottacini et al. 2010), showing an extremely high X-ray luminosity (~10⁴⁷ erg s⁻¹ in the 0.1–2 keV range) and a particularly hard spectrum (Γ_X ~ 1.2), indicative of a significant contribution of non-thermal jet emission. From a temporal perspective, modest optical-NIR and X-ray continuum variations have been observed (Bottacini et al. 2010). Here we present its γ-ray temporal characteristics. As shown in the monthly γ-ray light curve of NVSS J053954−283956 (Figure 1), there are three γ-ray flares. Enlarged 15-day and 3-day time bin γ-ray light curves corresponding to these flares are presented in Figure 3. For the first two flares (i.e., flare-A and flare-B), no significant intraday γ-ray variability is found, confirming a previous study that suggests a minimum γ-ray variability timescale of ~18 days (Vovk & Neronov 2013). However, the third flare (i.e., flare-C) detected in 2017 May,⁶ exhibits totally different behavior. In the 15-day time bin light curve of the flare-C, the flux of the eighth bin is significantly higher than the averaged flux, whereas the target maintains the low flux state for the remaining time. Moreover, this behavior is confirmed by the further 3-day time bin γ-ray light curve, where an intense γ-ray outburst suddenly appears. The flux quickly rises to the maximum value within about 6 days (1.46 days at the source frame) and the descent time is as short as the ascent time. In addition to the short variability timescale, the variation amplitude of this outburst is one order of magnitude larger than the average flux. The 3-day peaking flux reaches (1.2 ± 0.1) × 10⁻⁶ ph cm⁻² s⁻¹, with a very high TS value (426), while the averaged fluxes in the epoch of flare-C and all 105 months are (1.0 ± 0.1) × 10⁻⁷ ph cm⁻² s⁻¹ and (4.9 ± 0.2) × 10⁻⁸ ph cm⁻² s⁻¹, respectively. These results put a tight constraint on the doubling timescale at the source frame, ${\tau }_{\mathrm{doub},\mathrm{source}}={\rm{\Delta }}t\times \mathrm{ln}2/\mathrm{ln}({F}_{1}/{F}_{2})/(1+z)\lt 0.5\,\mathrm{day}.$ A power-law function provides an acceptable description of the burst SED, while the log-parabolic function does not bring any significant improvements. The spectral index of the burst SED (Γ = 2.53 ± 0.09) is slightly harder than the averaged SED (Γ = 2.78 ± 0.04), consistent with Cheung (2017). Meanwhile, a subsequent γ-ray localization analysis finds that the optimized location at this time is R.A. 8504 and decl. − 2864. The corresponding 95% C.L. error radius is 010, which is consistent with the value (008) listed in 3FGL (Acero et al. 2015). Since the angular separation between the γ-ray position and radio location of NVSS J053954−283956 is 006, it still locates within the 95% C.L. error radius. Note that in the normal observation mode LAT performs a complete and uniform coverage of the sky in 3 hr, thus we limit the minimum time bin in our analysis to 3 hr. Therefore, 12, 6, and 3 hr time bin γ-ray light curves are extracted to perform further investigations; see Figure 4. Since in the epoch of flare-C the target is at a high flux state in which most of the time bins are not upper limits, we adopt a simple χ² test to check whether the source is significantly variable. By optimizing the assumed constant flux, we find evidence of significant variability on intraday γ-ray light curves: (p, χ²/dof) = (7.0 × 10⁻⁸, 61.4/14) for a 12 hr light curve and (1.6 × 10⁻⁵, 65.7/25) for a 6 hr light curve, respectively. However, variability for the 3 hr light curve is not statistically significant due to large uncertainties, (p, χ²/dof) = (0.11, 26.7/19). Variability timescales are also estimated (listed in Table 2) by fitting data in the ascent phase with the exponential function:

$$\begin{eqnarray}&&F(t)=F({t}_{0})\cdot {2}^{-(t-{t}_{0})/\tau },\end{eqnarray} \tag{ 2 }$$

where F(t) and F(t₀) are the fluxes at times t and t₀, respectively, and τ is the characteristic timescale. In the 12 hr light curve, a quick rise begins in the sixth time bin (i.e., at MJD 57877.5, with a flux of (4.3 ± 1.6) × 10⁻⁷ ph cm⁻² s⁻¹), one day later the flux reaches the peak (i.e., at MJD 57878.5, (1.6 ± 0.3) × 10⁻⁶ ph cm⁻² s⁻¹). Since the eighth time bin, the target maintains a high flux state (10⁻⁶ ph cm⁻² s⁻¹) but with a relatively modest descent that costs about 2.5 days. Then the target is back to a quiescent flux state. This variability trend is confirmed by the 6 hr light curve, from which the ascent time is constrained as 18 hr (as short as 4.4 hr at the source frame). Together with the variation amplitude (i.e., F_bin15/F_bin12  8.4), τ_doub,source can be estimated as short as 1.4 hr. A similar τ_doub,source (i.e., 1.3 hr) can be also derived from the ascent phase in the 3 hr light curve despite the large error bars. Meanwhile, the 6 hr light curve reveals that the entire outburst may constitute of several sub-flares. Interestingly, violent variability may appear in these sub-flares. For example, in the 3 hr light curve the flux rises from (3.8 ± 2.7) × 10⁻⁷ ph cm⁻² s⁻¹ at MJD 57879.70 to (2.5 ± 0.9) × 10⁻⁶ ph cm⁻² s⁻¹ at MJD 57879.82, leading to a very short variability timescale of τ_doub,source ~ 16 minutes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Flux Doubling Timescales

NVSS Name	Epoch	Time Bin	τdoub,source (hr)	p-value
J053954−283956	flare-C	3 day	8.3 ± 0.2	2.9 × 10−50
J053954−283956	flare-C	12 hr	2.4 ± 0.6	7.0 × 10−8
J053954−283956	flare-C	6 hr	1.3 ± 0.3	1.6 × 10−5
J053954−283956	flare-C	3 hr	1.5 ± 0.6	0.11
J080518+614423	flare-A	3 day	19.6 ± 2.5	1.0 × 10−27

Note. (1) NVSS name of the object; (2) data time epoch given in Figure 1; (3) estimated doubling timescales at the source frame, along with 1σ errors; (4) probability that the null hypothesis stands.

Download table as:  ASCIITypeset image

3.2.2.2. NVSS J080518+614423

NVSS J080518+614423 (z = 3.033, Sowards-Emmerd et al. 2005) was also first known as a radio emitter (Becker et al. 1991) until its optical counterpart was identified (Snellen et al. 2002). Similar to NVSS J053954−283956, it was detected by WISE and Swift-BAT (D'Abrusco et al. 2012; Baumgartner et al. 2013). In its monthly γ-ray light curve, after two flares (i.e., flare-A and flare-B) in the first and second year of Fermi observation, it maintained low flux state for several years. These two flares were confirmed by enlarged 15-day time bin γ-ray light curves; see Figure 5. More importantly, evidence of intraday γ-ray variability as found from the 3-day time bin γ-ray light curve of flare-A, of which the variability index was given as σ_var = 10.8σ. But no similar behavior can be found in flare-B due to large uncertainties; also see Figure 5. In consideration of the ascent time of 6 days, together with the variability amplitude of 3.5, the intrinsic doubling timescale in flare-A can be estimated as τ_doub,source  19.6 hr. Note that the 3-day peaking flux ((3.9 ± 0.7) × 10⁻⁷ ph cm⁻² s⁻¹) is roughly 15 times the 105-month averaged flux. Intraday γ-ray light curves (i.e., 12 and 5 hr time bins) corresponding to this epoch have also been extracted; see Figure 6. No further constraints of the variability timescale are obtained. Similar to NVSS J053954−283956, the spectral index in flare-A of NVSS J080518+614423 (Γ = 2.40 ± 0.12) is harder than the averaged SED (Γ = 2.82 ± 0.05). Meanwhile, the optimized γ-ray positions at this time are R.A. 12125 and decl. 6173, with a 95% C.L. error radius of 012. The angular separation between the γ-ray position and the radio position of NVSS J080518+614423 is only 004, supporting the association.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2.3. NVSS J163547+362930

NVSS J163547+362930 was discovered by the MIT-Green Bank 5 GHz survey (Griffith et al. 1990). It was then included in the DR10 SDSS quasar catalog, with a redshift estimation of 3.647 (Pâris et al. 2014). As shown in Figure 2, there is a flare in the monthly >500 MeV γ-ray light curve of NVSS J163547+362930 that is probably not from the neighbor. Therefore, an enlarged 15-day time bin >500 MeV light curve has been extracted, together with another one from the neighbor 3FGL J1635.2+3809; see Figure 7. Interestingly, though the variability at the overall period is not significant (i.e., σ_var = 1.5σ) due to large uncertainties, the flux of the tenth bin (TS = 39) is roughly three times the averaged flux during this epoch. The neighbor is not detectable (TS < 4) for Fermi-LAT at the same time, confirming what we find in the monthly light curve. The optimized γ-ray location for this time bin is R.A. = 24889 and decl. = 3646. Since the angular separation is 006 while the 95% C.L. γ-ray error radius is 010, the association between NVSS J163547+362930 and the γ-ray source is confirmed. Similar to two former sources, a sign of bluer-when-brighter spectral variability behavior has been found for NVSS J163547+362930 (i.e., Γ_{flare,>500 MeV} = 2.58 ± 0.25 while Γ_{105 month,>500 MeV} = 2.73 ± 0.23). From the 15-day light curve, intrinsic doubling timescales can be constrained as ~2.2 and 1.4 days for the ascent and descent phases, respectively. Unfortunately, no evidence of intraday γ-ray variability was found in further temporal analyses.

Zoom In Zoom Out Reset image size

Benefiting from a large effective area and a wide field of view of Fermi-LAT, together with its routine survey mode covering the entire sky every 3 hr, our understanding of the γ-ray variability properties of blazars has been profoundly improved. Significant γ-ray variability has been accepted as a common feature of γ-ray blazars, with a variability amplitude of up to more than two orders of magnitudes and a variability timescale ranging from minutes to years (e.g., Nolan et al. 2012; Liao & Bai 2015; Ackermann et al. 2016; Liao et al. 2016). Over a long time period, evidence of quasi-periodic modulation of the γ-ray emission of several blazars has been reported (e.g., Ackermann et al. 2015; Zhang et al. 2017). In addition to flux variation alone, changes to the γ-ray spectrum are often observed during different flux states (Abdo et al. 2010c). Moreover, multi-wavelength campaigns that include γ-ray observation as well as complementary observations from radio to X-rays, have become a regular approach to investigating the physical processes of AGN jets, and γ-ray emission of blazars is observed to be tightly connected with emissions in other windows of electromagnetic radiation (e.g., Marscher et al. 2008; Abdo et al. 2010d; Liao et al. 2014, 2016). In high-redshift regimes (i.e., z > 2), γ-ray variability of several blazars has been studied in detail, including 0836 + 710 (z = 2.22; Akyuz et al. 2013, TXS 0536 + 135 (z = 2.69; Orienti et al. 2014), PKS 1830−211 (z = 2.51; Abdo et al. 2015), CGRaBS J0225+1846 (z = 2.69; Paliya et al. 2016), and PKS 2149−306 (z = 2.35; D'Ammando & Orienti 2016). Among these sources, PKS 2149−306 and PKS 1830−211 are the most luminous ones, with peaking γ-ray luminosities of 1.5 × 10⁵⁰ and 2.9 × 10⁵⁰ erg s⁻¹ (Abdo et al. 2015; D'Ammando & Orienti 2016), respectively. Meanwhile, evidence of fast γ-ray variability has been claimed for 0836+710 and PKS 1830−211 (τ_doub,source ~ 2 hr; Akyuz et al. 2013; Abdo et al. 2015). By comparison, blazars beyond redshift 3 are the focus of this study. On one hand, violent γ-ray variability, with large amplitudes (i.e., over one order of magnitude) for NVSS J053954−283956 and NVSS J080518+614423, are reported. In particular, the former has a peaking γ-ray luminosity of 1.1 × 10⁵⁰ erg s⁻¹ (in this work we take a ΛCDM cosmology with H₀ = 67 km s⁻¹ Mpc⁻¹, Ω_m = 0.32, and Ω_Λ = 0.68; Planck Collaboration et al. 2014), becoming the third most luminous source among the high-redshift blazars and the most distant γ-ray flaring source so far. Meanwhile, the intrinsic doubling variability timescale of NVSS J053954−283956 is constrained as short as 1.4 hr, which also makes it the most distant known source with intraday γ-ray variability . The energy dissipation mechanism corresponding to this extreme phenomenon could be the magnetic reconnection process and "minijets" scenario (e.g., Giannios et al. 2009; Cerutti et al. 2012; Blandford et al. 2017). On the other hand, our study reveals significant long-term γ-ray variability in three of these sources. The significant γ-ray variability, together with the results of γ-ray localization analyses in the flaring epoch, strongly indicate a blazar nature for these γ-ray sources. Meanwhile, the observed bluer-when-brighter spectral variability behaviors suggest that the identification of flaring epochs of high-redshift blazars is helpful for searching the distant high-energy γ-ray photons and may be useful for constrianing extragalactic background light (EBL) models.

Since variation of optical emissions of FSRQs has always been observed simultaneously with variation of their γ-ray emission, we also look into archival Palomar Transient Factory (PTF)/intermediate PTF (iPTF) data to search for additional evidence of fast variability for these sources. A detailed description of the PTF/iPTF project can be found in Rahmer et al. (2008), and data reduction pipelines and photometric calibration procedures can be found in the literature (e.g., Law et al. 2009; Ofek et al. 2012a, 2012b; Laher et al. 2014; Surace et al. 2015). Catalog Mould R (i.e., R_PTF) band SExtractor (Bertin & Arnouts 1996) data from the IPAC pipeline for all seven high-redshift blazars have been downloaded from the PTF/IPAC data archive hosted at the NASA/IPAC Infrared Science Archive (IRSA).⁷ The R_PTF mag for each source is extracted by matching the detected sources in the catalogs to the input high-redshift blazar with a match radius of two arcseconds. Unfortunately, except for NVSS J151002+570243 and NVSS J163547+362930, the PTF/iPTF data of these sources are rather sparse (i.e., ≤20 nights for 105 months). Meanwhile, the optical emission of NVSS J151002+570243 is faint (i.e., r_sdss  20.3 mag), close to the detection limit of a routine 60 s R_PTF band exposure (~21 mag; Ofek et al. 2012a). Thus, only the light curve of NVSS J163547+362930 is presented. The light curve is extracted by a python script⁸ with two steps. First, several comparison stars (i.e., CLASS_STAR > 0.95 and FLAGS = 0) are picked under photometric conditions (i.e., PHTCALFL = 1 and PCALRMSE < 0.04, 19 nights) based on their stable fluxes (i.e., mag_std < 0.04) throughout the entire time range. Then, differential photometry based on these comparison stars is adopted to correct outliers affected by bad weather in the "raw" light curve. As shown in Figure 8, significant variability can be directly seen in the daily averaged PTF/iPTF light curve of NVSS J163547+362930 that includes 92 nights with a time range between MJD 55635.5 and MJD 56847.2. There are four main optical flares in the light curve. Interestingly, in one flare that has good data coverage, the optical flux quickly rises from MJD 56060.4 with R_PTF = 20.6 to MJD 56065.3 with R_PTF = 19.1, indicating that the doubling time at the source frame then is ~12 hr; see the zoomed-in panel of Figure 8. This intraday optical variability behavior indicates that the central engine is highly active, consistent with the signs of fast γ-ray variability based on the 15-day γ-ray light curve, though no simultaneous iPTF observation is accessible when the γ-ray flare appears.

Zoom In Zoom Out Reset image size

Based on the observed fast γ-rays variability, values of the Doppler factor of the emitting jet blob should be high to avoid heavy absorption on γ-rays from soft photons via the γγ process. The optical depth of γγ absorption can be calculated as (Dondi & Ghisellini 1995):

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }({x}^{{\prime} })=\displaystyle \frac{{\sigma }_{{\rm{T}}}}{5}{n}^{{\prime} }({x}_{{\rm{t}}}^{{\prime} }){x}_{{\rm{t}}}^{{\prime} }{R}^{{\prime} },\end{eqnarray} \tag{ 3 }$$

where σ_T is the scattering Thomson cross section, n'(x') is the number density of the target photon, ${x}_{{\rm{t}}}^{{\prime} }$ is the energy of the target photon in dimensionless units, and R' is the absorption length. The soft photons can be from the jet radiation itself and external emission is from the accretion system (e.g., from an accretion disk or broad emissions lines). Since γ-rays with energy 3 GeV are detected during the flaring epoch of NVSS J053954−283956 and NVSS J080518+614423, energies of absorption soft photons could be at several keV and several tens of eV for the internal and external absorptions, respectively. Only the former is considered here because little information is available for emission at extreme ultraviolet wavelengths. Adopting variability timescales of 1.4 and 19.6 hr for these two sources and setting ${L}_{\mathrm{soft}}$ as 10⁴⁷ erg s⁻¹ (Ghisellini et al. 2010), we have δ  11 and δ  7. A similar calculation can be applied to NVSS J163547+362930 (δ  7 while E_γ  2 GeV) by assuming the optical and γ-ray emissions share the same radiation region. Therefore, the radii of the radiation region (i.e., ${R}_{\gamma }\lt \delta c{\tau }_{\mathrm{doub},\mathrm{source}}$) can be constrained to be smaller than 1.7 × 10¹⁵ cm, 1.5 × 10¹⁶ cm, and 9.1 × 10¹⁵ cm for NVSS J053954−283956, NVSS J080518+614423, and NVSS J163547+362930, respectively. The corresponding characteristic distance scale of the radiation region along the jet for a conical geometry is r_γ  R_γ/θ  R_γΓ  R_γδ, where θ is the jet-opening angle. Compared with the typical size of the broad line region (BLR; i.e., ~0.1 pc; Tavecchio et al. 2010), the locations of the γ-ray emission regions of these three sources (i.e., <0.03 pc) could be within the BLR. Meanwhile, the total jet power required to produce such high apparent γ-ray luminosities can be estimated as ${L}_{\mathrm{jet}}\simeq {L}_{\gamma }/({\eta }_{\mathrm{jet}}{{\rm{\Gamma }}}^{2})$, where η_jet is the jet radiative efficiency, typically ~0.1 (Nemmen et al. 2012), and L_γ are ~10⁵⁰ erg s⁻¹ and 10⁴⁹ erg s⁻¹, corresponding to NVSS J053954−283956and NVSS J080518+614423, respectively. Since their Eddington luminosities are given as ~3 × 10⁴⁷ erg s⁻¹ and 2 × 10⁴⁷ erg s⁻¹ (Ghisellini et al. 2010), their jet powers will exceed the Eddington luminosities if the Γ of the former source is smaller than 60, while the upper limit is 25 for the other source.

It is interesting to compare our results with those from SED modelings (Ghisellini et al. 2010; Ackermann et al. 2017), though the majority of the data used there are non-simultaneous. The Doppler factor values of these high-redshift blazars derived from the SED modelings are ~11–15, consistent with our results. Meanwhile, similar to our results, the locations of γ-ray emission regions are also found within the BLR from SED modelings. In addition to these approaches, direct measurements of the ejection speed of jet blobs and the observed brightness temperature using the shortest radio variability timescale (e.g., Hovatta et al. 2009), and comparing these approaches with the theoretically expected brightness temperature assuming equipartition (i.e., T_B = 5 × 10¹⁰ K; Readhead 1994), can also shed light on the Doppler factor. However, these high-redshift blazars are not included in current VLBA and radio flux monitoring projects (e.g., Lister et al. 2009; Richards et al. 2011; Jorstad et al. 2017). Nevertheless, let us make a comparison between our high-redshift blazars and those at low redshifts. MOJAVE parsec-scale, kinematic VLBA observations give an averaged apparent jet speed of β_app  9 of their sample (Lister et al. 2013). Particularly, observed β_app bright sources with detections of fast γ-ray variability are larger than 10 (Jorstad et al. 2005). Based on multi-wavelength radio light curves, a mean value of ~12 for the Doppler factor for F-GAMMA FSRQs is suggested (Liodakis et al. 2017). These indications from radio observations are confirmed by SED modeling studies (e.g., Ghisellini & Tavecchio 2015; Zhang et al. 2015). Meanwhile, studies of the luminosity functions of blazars and their parent populations allow for a constraint of the Doppler factor, which is in agreement with the kinematic radio observations (e.g., Ajello et al. 2012). In conclusion, no significant difference of the Doppler factor between these seven high-redshift blazars (i.e., z > 3) and nearby ones is found, and future simultaneous multi-wavelength campaigns will be helpful for understanding these highly active phenomena.

According to the Swift-BAT detected high-redshift blazars (Ajello et al. 2009), the number of SDSS-FIRST detected jetted AGNs with z > 3 is fairly less than predictions of the orientation-based AGN unified scheme (Volonteri et al. 2011). One possible explanation is that the averaged Lorentz factor of these Swift-BAT high-redshift blazars (i.e., Γ ~ 5) is generally lower than routine values (i.e., Γ ~ 15). Violent and rapid variability in blazars beyond redshift 3 are found in this study. Note that in addition to the two detections of FSRQ γ-ray variability with timescales of several minutes, based on either IACT observations or ToO repointing observations from Fermi-LAT (Aleksić et al. 2011; Ackermann et al. 2016), there are only a few FSRQs whose minimum variabilities are ~1–2 hr (e.g., Foschini et al. 2013; Saito et al. 2013; Hayashida et al. 2015; Liao & Bai 2015; Liao et al. 2016). NVSS J053954−283956 is one of the most violently active γ-ray FSRQs. Meanwhile, the light curves extracted here are under the survey mode operation of Fermi-LAT, due to the limited exposure time, and the reported doubling timescales should be treated as upper limits only. Moreover, in addition to NVSS J053954−283956, intraday variations are detected in two other high-redshift blazars (i.e., three in all seven), which suggests that such phenomena should not be rare for these sources. Therefore, this deficit could be due to selection effects rather than the gradient descent of the averaged Lorentz factor of blazars at higher redshifts.

We appreciate the helpful suggestions from the anonymous referee. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA's Goddard Space Flight Center. This research has also made use of the NASA/IPAC Extragalactic Database and the NASA/IPAC Infrared Science Archive, which are operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research makes use of the SIMBAD database, operated at CDS, Strasbourg, France.

This work was supported in part by the National Basic Research Program of China (No. 2013CB837000), and the National Natural Science Foundation of China under grants 11525313 (i.e., Funds for Distinguished Young Scholars), 11703093, and U1738136.