STOCHASTIC MODELING OF THE FERMI/LAT γ-RAY BLAZAR VARIABILITY

First, we modeled the Fermi/LAT light curves of the blazars in our sample with the OU process. We plot the light curves, the model realizations corresponding to the highest posterior likelihood, the standardized model residuals and their histograms compared to the expected standard normal distribution in Figures 2–6 (left). The sample posterior distribution functions of the characteristic timescale, τ₀ are presented in Figure 7. In Table 1, we list the median characteristic timescales and their uncertainties for all sources. The OU modeling gives the characteristic timescales in the τ₀ ≃ 0.7–43 day range (Table 1). In the case of 3C 454.3, the characteristic timescale is not constrained, and we derive the lower limit on its value, log [τ₀/day] > 2.0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is apparent that in general the OU process does not fit well the studied light curves. In particular, the model experiences systematic difficulties to fit the data during the quiet time periods when the assumed condition that TS ⩾ 25 forces the adaptive time bins Δt to be larger than a few days (see, e.g., Figures 3–5, left). To evaluate the quality of the fits, we compare the distributions of the standardized fit residuals with the standard normal distribution. In addition, for each source we calculate the auto-correlation function of the standardized residuals (Figure 8). It can be observed that while the distribution of the residuals approximates the standard normal distribution reasonably well, the auto-correlation functions deviate in several sources (e.g., B2 1633+38) from that expected for the white noise process. Thus, we proceed to testing the sup-OU model on the Fermi/LAT light curves of blazars.

Zoom In Zoom Out Reset image size

In the next step, we considered the superposition of the OU processes (sup-OU) and applied it to the Fermi/LAT blazar light curves. We observed a significant improvement to the quality of the fits in the low-flux periods, compared to the case of the OU modeling. The distributions of the standardized residuals assuming the sup-OU process can be compared with the OU model fit residuals and the standard normal distribution (Figures 3–6, right). The auto-correlation functions of the standardized residuals assuming the sup-OU process indicate that, indeed, in a majority of sources the sup-OU residuals approximate the white noise process better than the OU process residuals (Figure 8).

To formally evaluate what model is more successful in describing the Fermi/LAT blazar light curves, we computed the deviance information criterion (DIC; Spiegelhalter et al. 2002). DIC is used in Bayesian model selection problems where the posterior distributions of the models have been obtained by Markov chain Monte Carlo simulations. Models with smaller DIC are preferred to models with larger DIC. It is accepted that Δ(DIC) ≳ 10 is a difference substantial enough to prefer the model with smaller DIC. We use the following definition of the DIC

$$\begin{equation} {\rm DIC} = p_V + \bar{D}, \end{equation} \tag{ 3 }$$

where p_V = 0.5 × Var[D(θ)] describes the effective number of model parameters and hence the model complexity; while D(θ) is the deviance defined as D(θ) = −2log [p(y|θ)], where p(y|θ) is the model likelihood function, y is the data, and θ denotes the parameters in the model. The second term in Equation (3), $\bar{D}$, is the expected value of the deviance.

The comparison of DIC computed for the two models (Table 1) indicates that the sup-OU model is preferred over the OU model in all sources except 3C 454.3, 3C 66A, and PKS 2155-304, in which ΔDIC = DIC_supOU − DIC_OU < 10. Thus, for these three sources we accept the results of the OU modeling, while for the remaining sources we assume the results of the sup-OU modeling.

For each source, we plot the posterior probability distributions of the short and long characteristic timescales resulting from the sup-OU model in Figure 9. We notice that in many cases the posterior probability distributions of the characteristic timescales deviate from a Gaussian distribution. It can be inferred that while formally we can derive the median long and short characteristic timescales and their 90% confidence intervals in all sources, they may depend on our choice of τ_min and τ_max. The least deviated shapes are those representing the long characteristic timescale in PKS 1510−089; and the short characteristic timescales in 3C 454.3, and (to a lesser degree) in 3C 66A, PKS 0716+714, and BL Lac. Thus, we provide the median τ_L and/or τ_S and their uncertainties for these five cases, and the relevant limits for all sources (Table 1). The upper limits on τ_S can be compared with the smallest time spacing in the respective light curves, Δt (see Table 1), and the shortest sampled timescale, τ_min = Δt/10.

Zoom In Zoom Out Reset image size

We present the posterior probability distributions of the slopes, α, characterizing the intermediate part of the PSD between τ_L and τ_S in Figure 10. The distributions peak close to α ≃ 1 for the majority of the sources. The individual values and their uncertainties are listed in Table 1. The slopes do not depend on the blazar class (BL Lac versus FSRQ) or redshift (Figure 11, right and left). With the sup-OU model, we constrain the PSD slopes in all sources, except 3C 454.3, for which we find α < 0.43.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The main result of our study is the consistency of the four-year Fermi/LAT blazar activity with stochastic processes. The mixed OU process was preferred in 10 blazars from our sample, while in the 3 remaining blazars (2 BL Lacs, 3C 66A and PKS 2155-304; and 1 FSRQ, 3C 454.3) the single OU process and the mixed OU process resulted in fits of equally good quality.

The posterior probability distributions of the characteristic timescales from the OU process allowed us to constrain the Fermi/LAT characteristic timescales in two BL Lac sources ($\log [\tau _{\rm 0} / {\rm day}] = 1.39^{+0.43}_{-0.29}$ in 3C 66A, and $\log [\tau _{\rm 0} / {\rm day}] = 1.63^{+2.79}_{-0.35}$ in PKS 2155-304). This is the first time that the characteristic timescales have been observed in the Fermi/LAT light curves of blazars other than 3C 454.3. In 3C 454.3, the OU modeling resulted only in a lower limit, log [τ₀/day] > 2.0, which was inconsistent with the NM13 result. NM13 found that the PSD of 3C 454.3 can be characterized with α ≃ 1 and α ≃ 3 below and above the PSD break, respectively. The high frequency slope significantly steeper than ∼2 assumed in our modeling might have led to poor constraints on the characteristic timescales in this source. Characteristic timescales found through the OU modeling in 3C 66A and PKS 2155−304 fall outside the frequency range studied by NM13. We conclude that perhaps the method of the PSD construction used by NM13 was not sensitive enough to detect these BL Lac characteristic timescales in the Fermi/LAT light curves.

With regard to the sup-OU modeling, we reported the relevant limits on the characteristic timescales. In five cases with the posterior probability distributions resembling most closely a Gaussian distribution (τ_S in 3C 454.3, 3C 66A, PKS 0716+714, and BL Lac; τ_L in PKS 1510−089), we also provided the median values and their uncertainties. However, we cautioned that they might be affected by our choice of the lowest and highest sampled timescales, τ_min and τ_max. The upper limits on the short characteristic timescales obtained in our studies hint toward a sub-hour γ-ray time variability in at least four blazars: B2 1633+38, B2 1520+31, PKS 0454−234, and 3C 273, all belonging to the FSRQ class (Figure 11, left). Sub-hour variability timescales have been reported previously in another FSRQ, PKS 1510-089, in the context of the γ-ray flux doubling/halving timescales, rather than the characteristic timescales, in the study of the γ-ray flares in this source (Saito et al. 2013). Our adaptive light curves do not sample the sub-hour timescales due to a conservative TS constraint in each time bin. However, even with less constraining conditions on TS, it would be difficult to obtain sub-hour time resolution for a prolonged time required by our study, given the variable nature of blazars, as well as the characteristics of the Fermi/LAT instrument. Obtaining sub-hour time resolution in the Fermi/LAT light curves is feasible only for the highest γ-ray flux periods observed in the γ-rays. However, the results might be affected by the poorly understood systematic errors related to the Fermi/LAT data analysis on sub-orbital timescales. The lower limits on the long characteristic timescales indicate that observations conducted for a period of time longer than the four years considered here are needed in order to constrain τ_L.

Ackermann et al. (2010) computed a PSD and structure function of a 120 day long section of Fermi/LAT light curve of 3C 454.3, during which a powerful flare was recorded (MJD 55,120–55,260). They reported a 6.5 day timescale for this time period. We tested our models on the light curve of 3C 454.3 spanning the same range in MJD. Both the OU and sup-OU models resulted in fits of comparable quality (ΔDIC ≃ 3). However, the τ_S and τ_L distributions obtained through the sup-OU modeling were poorly separated, indicating at most a presence of one characteristic timescale in the light curve. The OU modeling was not successful in constraining the timescale and we only obtained a lower limit, log [τ₀/day] > 1.0. The OU limit on τ₀ derived for the 120 day light curve of 3C 454.3 and the 6.5 day timescale reported in Ackermann et al. (2010) are in disagreement, most probably due to the difference between the PSD assumed in the OU model and that found in Ackermann et al. (2010; α ≃ 1.40 and α ≃ 1.56 below and above the PSD break, respectively).

In 10 of 13 Fermi/LAT blazar light curves in our sample, the sup-OU model was preferred over the OU model based on DIC. Thus, the second important result of our study is that for the range of timescales sampled with the sup-OU model (Figure 11, left), the underlying PSDs of these 10 sources resemble a power-law function proportional to 1/f^α with 0 < α < 2, with majority of the slopes clustering around unity (Figure 10). We were able to constrain α in all but one source, 3C 454.3, for which we found α < 0.43.

Abdo et al. (2010) reported on an average Fermi/LAT PSD slope for their brightest FSRQs, α_FSRQ = 1.4 ± 0.1, and BL Lacs, α_BLLac = 1.7 ± 0.3. In general, the slopes derived for the individual sources in our sample through the sup-OU modeling are shallower than those respective averages. In the FSRQ class, slopes of all our sources are inconsistent with the Abdo et al. average. In the BL Lac class, only PKS 2155−304 and (marginally) BL Lac and Mkn 421 have slopes consistent with the Abdo et al. average at the 2σ confidence level, but these are the three sources with the broadest posterior slope distributions in our sample (see Figure 10).

We also compared our derived PSD slopes to those found by NM13 (Table 1). We found that in PKS 0454−234, 3C 279, 3C 66A, and PKS 2155−304 the two methods are consistent with each other in terms of the derived PSD slopes within the provided error bars. In 3C 454.3, PKS 1510−089, 3C 273, BL Lac, and Mkn 421 our slopes were shallower than those derived by NM13. In 3C 454.3 we suspect that the reason for the discrepancy is in our assumption that the high frequency slope equals 2. If this slope is as steep as 3, as reported in NM13, then our modeling might result in a slope below the high frequency break that is flatter than the actual value in order to maintain a satisfactory fit.

The origin of blazar γ-ray variability seen in the Fermi/LAT light curves remains under discussion. In our studies, we followed a stochastic approach to describe this variability motivated by its success in explaining the optical and/or X-ray variability of nonblazar accreting black holes. We demonstrated that the OU and/or mixed OU processes approximate well the blazar γ-ray light curves. Our modeling implies the nonlinearity of these blazar light curves, similar to the case of the black hole binary and AGN X-ray light curves, as argued by Uttley et al. (2005).

The interpretation of our models with regard to the optical and/or X-ray black hole variability is along the line of the propagation class of models (e.g., Lyubarskii 1997, Kotov et al. 2001) where perturbations originating at the outer radii of an accretion disk propagate inward and couple with each other producing the observed broad PSD shapes. In particular, Kelly et al. (2011) showed that the mixed OU process is the solution to the linear diffusion equation perturbed by a spatially correlated noise field. Within this framework the low-frequency break in the PSD corresponds to the diffusion timescale in the outer region of the accretion flow, and the shape of the PSD above the low-frequency break depends on the viscosity of the flow. An additional high-frequency break in the PSD may exist if the noise field is spatially correlated. The high-frequency break then corresponds to the time it takes a perturbation traveling at the viscous speed to cross the characteristic spatial scale of the noise field.

The short characteristic timescale corresponding to the high-frequency break in the black hole binary and AGN X-ray PSD is an important variability parameter, due to its correlation with the black hole mass and mass accretion rate that holds over many orders of magnitude (McHardy et al. 2006, Kelly et al. 2011). Kelly et al. (2011) demonstrated that τ_S obtained through direct modeling of the 2–10 keV X-ray light curves of a sample of 10 nearby Seyfert galaxies (Sobolewska & Papadakis 2009) with the sup-OU process are in agreement with those derived using the PSD construction methods (Markowitz et al. 2003, McHardy et al. 2007, González-Martín & Vaughan 2012). However, blazar γ-ray emission is dominated by the processes related to the jet physics and beamed due to relativistic effects, unlike the X-ray coronal emission of the radio-quiet AGN and black hole binaries. Nevertheless, it could be envisaged that the accretion disk variability structure is carried outward in the jet, leading to the jet variability with properties inherited from the innermost radii of the accretion disk/X-ray corona system (e.g., Lohfink et al. 2013). It would then be expected that the X-ray and Fermi/LAT short characteristic timescales of blazars were consistent with those found in X-rays in radio-quiet and radio-loud AGNs. Thus, it is interesting to compare the blazar and nonblazar variability properties across different energy bands in order to search for evidence in favor of or against a common variability mechanism (Figure 11, middle; note differing frequency bands sampled by the respective studies).

Blazar PSDs have only been studied in several sources in energy bands other than the γ-rays. Soft X-ray blazar PSDs (here below 20 keV) were computed for three BL Lac sources (Mkn 421 and PKS 2155−304 also included in our Fermi/LAT sample, and Mkn 501) by Kataoka et al. (2001), and for FSRQ sources by McHardy (2008; 3C 273), Chatterjee et al. (2008; 3C 279). Characteristic timescales corresponding to the high-frequency PSD break were detected in all of these sources with the exception of 3C 279. The sampled frequency bands and detected PSD breaks are indicated in Figure 11 (middle). In general, these studies demonstrated that in soft X-rays blazars show high-frequency breaks and PSD shapes below these breaks that are largely consistent with those of Seyfert galaxies and radio-loud AGNs (e.g., 3C 120, Marshall et al. 2009; 3C 390.3, Gliozzi et al. 2009), despite a possible contribution from a jet emission in the latter class. However, at the frequencies above the high-frequency break, the soft X-ray PSDs of BL Lacs seem to drop more rapidly than it is observed in Seyfert galaxies. Hard X-ray (14–150 keV; Swift/BAT) PSDs have recently been constructed by Shimizu & Mushotzky (2013) for a sample of 30 AGN, which included two FSRQs (3C 273 and 3C 454.3) and one BL Lac (Mkn 421). This study covered longer timescales than the soft X-ray study of Kataoka et al. (2001). A characteristic timescale was only detected in 3C 273, preventing a detailed cross-band comparison. The PSDs of all the other AGNs in their sample were described with unbroken power-law functions with a typical slopes, 0 < α < 2, when constrained.

The γ-ray blazar PSD slopes derived in our study are, in general, consistent with those observed in X-rays in nonblazar AGNs. However, we were able to derive only the upper limits on the blazar short characteristic timescales. Given that they are consistent with the low end of τ_S detected in the radio-quiet AGNs (Figure 11), we cannot discard the possibility that the γ-ray τ_S in blazars and the X-ray τ_S in nonblazar AGNs are comparable.

However, the OU γ-ray characteristic timescales detected in PKS 2155−304 and 3C 66A are much longer than those of the BL Lac sources (and also FSRQs) in the soft X-rays. In addition, closer inspection of a few available blazar PSDs across different energy bands suggests that they are rather different. In particular, in PKS 2155−304 the value on τ_S derived in our work is in disagreement with the X-ray τ_S found in Kataoka et al. (2001). The PSD slope derived here for 3C 279, α ≃ 1.1, is in disagreement with that reported in the soft X-rays by Chatterjee et al. (2008, α ≃ 2.3). We do not detect a γ-ray characteristic timescale in 3C 273, while the soft and hard X-ray detections were reported at the values well within our range of sampled timescales (McHardy 2008; Shimizu & Mushotzky 2013). These discrepancies imply that either the blazar variability process is energy dependent, or the origins of the X-ray and γ-ray blazar variability are different.

The processes intrinsic to the jet physics such as the internal shock model (e.g., Ghisellini 1999) or the minijets-in-a-jet model (e.g., Begelman et al. 2008, Giannios et al. 2009, Narayan & Piran 2012) may well dominate the blazar variability structure leading to the γ-ray variability properties that are largely different from those implied by the accretion disk physics and propagation models. Alternatively, the γ-ray variability pattern may be a combination of the baseline stochastic disk/corona variability carried out in a jet, and variability intrinsic to the jet. If the propagation model of variability envisaging coupling of perturbations at different radii continues toward the jet, and the disk X-ray perturbations couple with the X/γ-ray perturbations intrinsic to the jet, then the resulting Fermi/LAT variability structure may differ from that known from the X-ray band. Additional complications may be introduced by altering the variability timescales due to relativistic effects.