TIME SERIES ANALYSIS OF GAMMA-RAY BLAZARS AND IMPLICATIONS FOR THE CENTRAL BLACK-HOLE MASS

About 1% of galaxies have an active galactic nucleus (AGN), which emits 10⁹–10¹⁴ times the solar power over a wide range of the energy spectrum, from radio to gamma-ray energies. AGNs constitute one of the most violently variable and interesting classes of object in the universe. The activity of AGNs is believed to originate from the central supermassive black hole, with masses of 10⁶–10⁹ M_☉, and part of their energy is emitted in electromagnetic radiation from the surrounding region including the accretion disk formed around the black hole and relativistic jets ejected along rotation axes (e.g., Urry & Padovani 1995). A subclass of radio-loud AGNs are called blazars, in which the line of sight lies close to the jet axis, and the emission from relativistic jets is only visible in this class of AGNs due to the relativistic beaming effect, especially in the high-energy region. The electromagnetic spectra of blazars are dominated by non-thermal radiation produced in the jets. The popular scenario to explain these emission spectra assumes that the particles in the jets are accelerated to high energies by diffusive shocks in the jets and induce emission via interaction with surrounding matter/radiation (e.g., Fossati et al. 1998).¹

Observations of blazars at various wavelengths have revealed fast time variability which is most plausibly related to instabilities in the emission environment near the black hole (e.g., Ulrich et al. 1997). Past observations suggest the variability is larger at higher energies (the most extreme example is Mrk 501 Nowak et al. 2012), which may indicate that the higher energy emission comes from the region closer to the central black holes. The variability time scale reflects the size of the emission region, and thus the study of temporal behavior of gamma-ray flux is an excellent probe of the region close to the central engine, i.e., the supermassive black hole.

Blazars are known to show flaring activity which occurs randomly and continues for several days to months. In order to study their temporal variability precisely, blazars should be monitored continuously, or at least frequently. It is not an easy task for narrow field-of-view telescopes like optical, X-ray, and Cherenkov (TeV gamma-ray) instruments to monitor many blazars for long periods. Besides, their observations are limited as ground-based telescopes can only operate on dark, clear nights and X-ray satellites are used in pointing-mode observations.

Nevertheless, in the X-ray band, Hayashida et al. (1998) evaluated the central black hole masses in several AGNs based on their rather well-sampled X-ray light curves obtained with the Ginga satellite and suggested masses of one to two orders of magnitude smaller than previous estimates. Kataoka et al. (2001) studied time variability of three TeV-detected AGNs based on ASCA and/or RXTE observations and showed that (10⁷–10¹⁰) M_☉ black holes and internal shocks that start to develop at 100 times the Schwarzschild radii could explain the observed properties.

In the GeV gamma-ray band, the Fermi Gamma-ray Space Telescope has been monitoring the whole sky with the Large Area Telescope (LAT) since 2008. The LAT is a wide field-of-view gamma-ray imager that observes one-fifth of the sky at any instant and that scans the whole sky in a day (Atwood et al. 2009). The second Fermi-LAT catalog, which contains 1092 (28 identified and 1064 associated) AGNs among the 1873 detected sources in the 100 MeV to 100 GeV range (Nolan et al. 2012). Gamma-ray light curves of several tens of blazars are provided on a daily basis and this is a good database to study the time variability of blazars.

Abdo et al. (2010b) reported a detailed analysis of the variability of 106 objects in the Fermi-LAT Bright AGN Sample. They showed that the temporal behavior of gamma-ray fluxes of variable sources can be described by power-law power spectral density (PSD) in general, with a few blazars that showed strong activity exhibiting complex and structured temporal profiles. They examined whether it was possible to characterize blazar type with the PSD slope, but the results were not conclusive.

In this paper, we report the time series analysis of gamma-ray light curves of 15 blazars based on Fermi-LAT data and discuss the results in relation to the properties of the central engine. Our analysis is the first systematic study of long-term variability of blazar emission in the gamma-ray energy band, although the analysis method itself has been applied and reported previously (e.g., in the X-ray band; Lawrence et al. 1987; McHardy & Czerny 1987; Miyamoto et al. 1994; Hayashida et al. 1998; Kataoka et al. 2001).

We use the "monitored source light curves" provided by the Fermi Science Support Center (FSSC)² for bright and transient gamma-ray sources. They are regularly updated throughout the mission. In this paper, we analyze the daily light curves of 15 AGNs (Table 1) in the energy range 100 MeV–300 GeV. The data period is between 2008 August 9 and 2012 April 26. These sources are selected because a large fraction of data points are detections, not upper limits, so that we can extract useful information on time variability. Note that the usage of daily light curves spanning 44 months naturally limits our time series analysis to the 10⁻⁸–10⁻⁵ Hz range. Faster variability observed in the case of 3C 454.3 (Abdo et al. 2011), for example, is out of the scope of the present analysis.

As we noted in Section 1, the variability time scale of AGN flares reflects the size of the emission region, and thus the study of temporal behavior of gamma-ray flux can be a good probe to explore the physical environment close to the central engine. However, characterizing the time scale is not a simple task since we know that the intensity of gamma-ray emission from AGN flares varies very irregularly. The fastest doubling time has been widely used (see Barr & Mushotzky 1986, for example) as a variability measure, but it depends on data quality and coverage. Here we adopt a spectral analysis, the normalized power spectrum density (NPSD), to evaluate the characteristic timescale of light curves which fluctuate chaotically, after Miyamoto et al. (1994).

2.1. Normalized Power Spectral Density

The PSD shows the degree of variation at every frequency (or cycle) by calculating the Fourier transform of time variable data (Lawrence et al. 1987; McHardy & Czerny 1987).

The NPSD, which is obtained by dividing the PSD by the average source intensity squared, has proven to be useful to compare variability at each frequency even if the brightness changes (Miyamoto et al. 1994; Kataoka et al. 2001). It is defined as

$$\begin{eqnarray} P(f) &=& \frac{[a^2(f)+b^2(f)-\sigma ^2_{\mathrm{stat}}/n] T}{F^2_{\mathrm{av}}} \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} a(f) &=& \frac{1}{n} \sum ^{n-1}_{j=0}F_j \cos (2 \pi f t_j) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} b(f) &=& \frac{1}{n} \sum ^{n-1}_{j=0}F_j \sin (2 \pi f t_j), \end{eqnarray} \tag{ 3 }$$

where F_j is the source count rate at time t_j (0 < t_j < n − 1), T is the total time length, F_av is the mean value of source count rates, and σ_stat is the error due to counting statistics. In our analysis, we calculated the power P(f) for some discrete frequencies given by f = k/T (k is an integer and 1 < k < n/2) and averaged. The Fermi-LAT light curves are given with flux errors (e_j) and their standard deviation ($\ssty \sqrt{\sum _j e_j^2/N}$) is substituted for σ_stat as a rough estimate of error of counting statistics (see Section 2.2 for more discussion). The error bars of the NPSD are standard deviations of powers in each frequency bin (see Hayashida et al. 1998; Kataoka et al. 2001, for more discussion of NPSD analyses). In our calculation of the NPSD, we did not use upper limits in the light curves. In addition, we did not interpolate any blank (i.e., missing) data contained in the light curves.

2.2. Poisson Noise

We calculated the NPSD for the Fermi-LAT daily light curves of 15 AGNs using nine frequency bins divided logarithmically from 10⁻⁷ Hz to 10^−5.2 Hz. Plots are shown in Figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One may note that data points at high frequencies are missing for most of the sources in Figure 1. One reason is that many observations in the Fermi-LAT light curves yield only upper limits. Another reason is that large flares which last for ten of days or more are rare: we cannot have points above ∼10⁻⁶ Hz without such flares.

We applied least-square fits to these points assuming a power-law

$$\begin{equation} f(\nu) \propto \nu ^{\gamma } \end{equation} \tag{ 4 }$$

and/or a broken power-law

$$\begin{equation} f(\nu) \propto \lbrace \begin{array}{@{} l}\nu ^{\gamma _1} \ (\nu < \nu _{b}) \\ \nu ^{\gamma _2} \ (\nu \ge \nu _{b}), \end{array} \end{equation} \tag{ 5 }$$

where ν is the frequency and ν_b is the "turnover" frequency. Fit parameters are summarized in Table 1. The fit lines are overplotted in Figure 1, where broken power-law lines are plotted only when the reduced χ² values are smaller than single power-law values. We see NPSDs for four sources, PKS 0537−441, 3C 279, 3C 454.3 and PKS 2326−502, are better fitted by broken power-laws than by single power-laws.

The NPSD plots for PKS 0537−502, 3C 279 and PKS 2326−502 show upward turnovers above 10^−6.18 Hz, 10^−5.88 Hz and 10^−6.10 Hz, respectively, but the slopes above these frequencies (γ₂) have large uncertainties and are consistent with zero: thus they may have reached a constant Poisson-like noise level which is not removed by our rough estimate of counting statistical error (see Section 2.1). We checked the difference of NPSD values before and after removing the Poisson noise, assuming the value of Poisson noise is the normalized square-root of sum of squares of several flux's error: $\ssty \sqrt{\sum e^2_i/N}$. With this procedure, the NPSD values at high frequencies showed smaller values than those before removal, but the slope above the turnover of the NPSD plot remained flat (consistent with zero slope). Thus, even if we could not remove the effect of Poisson noise completely, it seems this flattening behavior does not have a physical origin.

On the other hand, the NPSD plot for 3C 454.3 shows a turnover at 10^−5.83 Hz and the slope above it, −3.08 ± 0.83, is well determined. The reduced χ² value decreases from 1.29 for single power-law fit, which is not at acceptable level, to 0.23 for broken power-law fit, which is acceptable. Thus, only the plot for 3C 454.3, which exhibited an extraordinary large flare in 2010 November (Abdo et al. 2011), seems to show a physically meaningful turnover, at 10^−5.83 Hz, which we discuss further in the next section.

The internal shock model is a popular scenario of blazar emission as it can explain spectral energy distributions and time-lag features (Böttcher & Dermer 2010 and references therein). Wehrle et al. (2012) studied multiwavelength variations of the 3C 454.3 outburst from 2010 November to 2011 January with observations by Herschel, Swift, Fermi-LAT, optical telescopes, and submillimeter arrays. They proposed a model in which turbulent plasma crosses a conical standing shock in the parsec-scale region of the jet, based on time-resolved spectral energy distributions for this outburst. Thus, here we assume the internal shock model as the emission mechanism for gamma-ray flares and we interpret the characteristic time scale that we found in terms of this model.

Kataoka et al. (2001) studied the X-ray variability of three TeV blazars, Mrk 421, Mrk 501, and PKS 2155−304, using ASCA and RXTE data. In order to interpret the observed characteristic time scale which they found in their NPSD plots, they assumed a simple model based on the internal shock model. They considered two relativistic blobs with bulk Lorentz factors Γ and a₀Γ (a₀ > 1) ejected at the times t = 0 and t = τ₀ (> 0), respectively, and when the second, faster blob catches up and collide with the first, slower blob, the resulting shock generates a high-energy flare. In this model, the mass of the central black hole, M_CBH, is derived from the variation time-scale, t_var, as

$$\begin{equation} M_{\mathrm{CBH}} \simeq 9 \times 10^8\, M_{\odot } \frac{t_{\mathrm{var}}}{\mathrm{day}}\frac{10}{k}\frac{a^2_0-1}{2a^2_0}, \end{equation} \tag{ 6 }$$

where M_☉ is the solar mass and k = cτ₀/R_g ⩾ 3 with the Schwarzschild radius R_g. They derived (10⁷–10¹⁰) M_☉ as the masses of the central black holes of these blazars.

Though this model was devised to explain the X-ray time variability of blazars, here we assume the same mechanism works for the gamma-ray time variability³, and we applied the above equation to estimate the central black hole mass of 3C 454.3. We take the variation time-scale as the inverse of the turnover frequency which we observed in the NPSD plot of 3C 454.3, t_var = 1/(10^−5.834 Hz) = 6.82 × 10⁵ s =7.89 days. Figure 2 shows the result for several values of model parameters: a₀ = 1–100 and k = 5, 20, 100 following Kataoka et al. (2001). We can infer that the central black hole mass is in the range (10⁸–10¹⁰) M_☉ from this plot in most of the parameter space (a₀ ≳ 2, which means the Lorentz factors of colliding blobs differ significantly).

Zoom In Zoom Out Reset image size

Alternatively, we can infer the range of the unknown parameter k, or the light crossing time in units of the Schwarzschild radius, by assuming the central black hole mass estimated by other methods.

3C 454.3 is one of the most well-known and well-studied gamma-ray sources. Recently it showed two large flares, in 2009 November–December (Ackermann et al. 2011; Striani et al. 2010) and 2010 December (Abdo et al. 2011). Bonnoli et al. (2011) analyzed the 2009 flare and estimated the central black hole mass by refining the discussion of Gu et al. (2001), who used the broad line width and the distance of the broad line region (BLR) from the center. Assuming the broad emission lines being produced in clouds which are gravitationally bound and orbiting with Keplerian velocities (Dibai 1981), the central black hole mass can be given by M_CBHR_BLRV²G⁻¹, where R_BLR is the radius of the BLR and V is the velocity of the clouds in the BLR. Gu et al. (2001) and Bonnoli et al. (2011) estimated the central black hole mass of 3C 454.3 is 4 × 10⁹ M_☉ and 5 × 10⁸ M_☉, respectively.

Figure 3 shows the variation of k = cτ₀/R_g as a function of the ratio of blob Lorentz factor, a₀, for two estimated values of the central black hole mass. We see k is in the range of 7 to 70 for a₀ ≳ 2. This parameter range is more restrictive than in the case of general consideration for blazar X-ray flares by Kataoka et al. (2001), where k in the range of 5–100 for t_var = 1–10 days and M_CBH = (10⁷–10¹⁰) M_☉.

Zoom In Zoom Out Reset image size

The GeV gamma-ray flux during the strong outburst of 3C 454.3 in 2010 November exhibited a very fast variability with the rise time of 4.5 ± 1 hr and the fall time of 14 ± 2 hr (Abdo et al. 2011). We tried to apply the Equation (6) assuming the variation time scale as this rise time, and the results are shown in Figures 2 and 3. The central black hole mass inferred from this time scale (<2 × 10⁸ M_☉; see Figure 2) is much smaller than that inferred from the characteristic time scale in NPSD, and the light crossing time (k < 2; see Figure 3) is less than the time to cross the last stable orbit of the black hole. Thus it seems unreasonable to assign the short time scale of 4.5 hr to the internal shock model under consideration in this paper: it should be interpreted as, e.g., shocks forming with strong anisotropic geometries, albeit a low duty cycle (e.g., Salvati et al. 1998), or existence of small active regions, inside a larger jet, moving faster than the rest of the plasma, occasionally pointing toward us (Ghisellini & Tavecchio 2008).

Using the gamma-ray daily light curves observed by Fermi-LAT over 3.6 yr, we studied the temporal behavior of 15 AGNs by calculating the NPSDs for each sources. One source, 3C 454.3, showed a clear turnover in the NPSD curve which corresponds to a characteristic time scale of 6.82 × 10⁵ s.

This time scale can be interpreted as a result of an collision of blobs in the internal shock of the blazar jet, as discussed by Kataoka et al. (2001) for the NPSD of X-ray data on blazars. The time variation scale we found indicates the central black hole mass of 3C 454.3 is in the range of (10⁸–10¹⁰) M_☉.

Alternatively, if we assume the central black hole mass of 3C 454.3 is (0.5–4) × 10⁹ M_☉ after Gu et al. (2001) and Bonnoli et al. (2011), we can infer the time interval emitted the internal shock wave is 7–70 times the light crossing time of the Schwarzschild radius of the central black hole.

We thank the anonymous referee for useful comments and suggestions. This work is supported in part by the Grant-in-Aid from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, No. 22540315.