FERMI LARGE AREA TELESCOPE AND MULTI-WAVELENGTH OBSERVATIONS OF THE FLARING ACTIVITY OF PKS 1510-089 BETWEEN 2008 SEPTEMBER AND 2009 JUNE

We extracted light curves from the entire data set, to investigate the flaring activity. To take into account possible biases or systematics when the source flux is faint, we used two different time binnings of 1 day and 1 week. The light curves were extracted using gtlike, fitting the source spectrum by means of a power-law (PL) distribution $(dN/dE \propto E^{-\alpha _{\gamma }})$, where α_γ is the photon index, following the prescription given in the previous section. The flux was evaluated by integrating the fitted model above 0.2 GeV. The lower panel of Figure 1 clearly shows four major flaring episodes: between 2008 August 30 and 2008 September 26 (flare a), between 2009 January 4 and 2009 January 27 (flare b), between 2009 March 10 and 2009 April 9 (flare c), and between 2009 April 15 and 2009 May 12 (flare d). The source was almost quiescent between the end of 2008 September and the beginning of 2009 January (see Table 1 for a summary). The flux light curves with different temporal binning are compatible, with the daily integration binning showing better rapid flux variations that are smoothed in the weekly binning. A study of these variations based on the autocorrelation, Fourier analysis, and structure function is presented by Abdo et al. (2010d) together with other blazars.

Zoom In Zoom Out Reset image size

Figure 2 shows close-up light curves of the four flares and the green line represents the optical data in the R filter (see Sections 3 and 4). Since the statistics during the flares were high, it was possible to use also a 12 hr binning (blue points). Typically, the flares have a complex structure with peaks having durations from about 1 to 5 days and their moderately asymmetric profile can result from the overlapping of subsequent episodes. In a few cases, significant variations by a factor of two within 12 hr were detected. To have a better estimate of the rising and decaying timescales, we fit two rapid flares, with an almost regular shape, using an analytical law of the form A · exp^t/τ. In the case of the flare peaking at t≃ 54846 MJD (dashed black line in panel (b) of Figure 2), we find a rise faster than the decay: the rising timescale is τ ≃ 0.3 days and the decaying one is ≃1.4 days. The second event (flare peaking at t≃ 54962 MJD; dashed black line in panel (d) of Figure 2) followed the opposite behavior, having a best-fit rising timescale of ≃0.8 days, and a decay time of ≃0.25 days. Using a bin width of 6 hr, the shape of the flare is nearly symmetric, with rise and decay e-folding time of about ≃0.12 days. Such a fast variability can constrain the radiative region size R_rad by the well-known relation

$$\begin{equation} R_{\rm rad} \le \frac{c \Delta t \delta }{1+z}, \end{equation} \tag{ 1 }$$

where c is the speed of light, δ = 1/(Γ(1 − βcos θ)) is the beaming factor depending on the bulk Lorentz factor Γ and a viewing angle θ ≃ 1/Γ, and z is the cosmological redshift. Adopting the very fast superluminal velocity reported by Homan et al. (2002; β_app ≃ 20, in very good agreement with the results presented in Section 4.3) from δ ≃ Γ ≃ β_app we can estimate R_rad ≲ 9 × 10¹⁵ cm (in the case of τ ≃ 0.25 days). We will compare this result with other constraints derived in Section 5.1.

Zoom In Zoom Out Reset image size

We analyzed the γ-ray spectral shape of PKS 1510-089 during the whole period, the quiescent state, and the four flaring episodes using three spectral models: a PL, a log parabola (LP), $dN/dE \propto E/E_0^{-\alpha _{\gamma }-\beta \log (E/E_0)}$ (Landau et al. 1986; Massaro et al. 2004), and a broken power law (BPL). In the case of LP spectral law, the parameter β measures the curvature around the peak. The LP distribution has only three free parameters, and the choice of the reference energy E₀ does not affect the spectral shape; we fixed its value to 300 MeV. We performed the spectral analysis using an unbinned maximum-likelihood estimator (gtlike) and the same prescription given in Section 2. We used a likelihood ratio test ⁸⁴ (LRT; Mattox et al. 1996) to check the PL model (null hypothesis) against the LP model (alternative hypothesis). Since the PL is often rejected, we also test the LP model (null hypothesis) against the BPL model (alternative hypothesis). The results concerning the LRT are summarized in Table 2, and in Table 3 we report the details of the spectral analysis for each time range and spectral model. Due to the nonderivable character of the BPL law, we used also the loglikelihood profile method to determine the best-fit parameter for this model. The corresponding statistical uncertainty was estimated from the difference in the likelihood value with respect to its minimum such that −2ΔL = 1. The γ-ray spectrum of PKS 1510-089 is well described by an LP, with the only exception being the quiescent state. The value of the LRT reported in Table 2 shows that both for the flares (a, b, c, d) and for the full period, the LP model describes the spectrum better than PL with a probability higher than ≃99.6%. The BPL, in contrast, does not provide an improvement with respect to the LP model. The only exception is the flare b, but the false positive probability is about 27.33%, so there is no evidence to reject the null hypothesis. The LP model is then preferred because of the lower number of parameters. Moreover the curvature parameter β can be linked to physical processes such as the acceleration or the effects of the Klein–Nishina (KN) regime in the IC process, as we will discuss in Sections 5.1 and 5.2. For a better visualization of the SED shape and to show the departure from a PL trend, we produced an SED by performing an independent likelihood analysis starting form a grid of 20 energy bins logarithmically equispaced. The bins were then grouped in order to have at least 10 photons per bin, and the highest energy bin was chosen according to the maximum energy encircled within 95% of the point-spread function (PSF). The results are shown in Figures 3 and 4. In Figure 3, we show the full-period SED and the one extracted during the quiescent state. We plot by a dotted line the PL model, by a dashed line the LP model, and by a dot-dashed line the BPL model. With a red upward arrow, we indicate the highest energy event within 95% of the PSF for the whole period data set, corresponding to an energy of approximately 30 GeV. From the plot of the PL model residuals (lower panel), it is possible to clearly see the departure from a PL trend: the deviations, both at low and high energies, suggest for a spectral curvature confirming the LRT results. Moreover, it is possible to note that the BPL model does not deviate significantly from the LP trend, supporting again the LRT analysis. The SEDs of individual flares are plotted in Figure 4 and show that the spectrum was curved also during the single flaring episodes. We note that the flare-integrated spectral shape did not change significantly, despite the huge flux variations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Unbinned Likelihood Spectral Fit Results

Time Range	PL			LP				BPLa
	αγ	F100b	loglike	αγ	β	F100b	loglike	αγ	$\alpha _{\gamma _1}$	Ebc	F100b	loglike
Full	2.44 ± 0.01	1.32 ± 0.03	323082.6	2.23 ± 0.02	0.09 ± 0.01	1.12 ± 0.03	323056.4	2.30 ± 0.02	2.63 ± 0.04	980 ± 130	1.37 ± 01	83057.6
								...	...	...	...	...
Quiescent	2.43 ± 0.07	0.27 ± 0.03	83057.7	2.3 ± 0.1	0.03 ± 0.05	0.25 ± 0.04	83057.5	2.5 ± 1.5	2.4 ± 1.5	700 ± 8000	0.3 ± 02	28903.5
								...	...	...	...	...
Flare a	2.55 ± 0.01	1.55 ± 0.02	28908.9	2.2 ± 0.1	0.19 ± 0.06	1.16 ± 0.12	28902.4	2.39 ± 0.01	3.22 ± 0.04	1500 ± 40	1.37 ± 0.01	28903.5
								...	...	...	...	...
Flare b	2.35 ± 0.04	2.1 ± 0.1	23932.7	2.0 ± 0.1	0.10 ± 0.03	1.7 ± 0.1	23928.5	2.26 ± 0.05	3. ± 0.03	3400 ± 800	1.9 ± 0.1	23927.9
								...	...	...	...	...
Flare c	2.37 ± 0.3	3.9 ± 0.1	38328.9	2.13 ± 0.06	0.1 ± 0.02	3.3 ± 0.2	38318.5	2.28 ± 0.04	2.9 ± 0.2	1918 ± 593	3.7 ± 0.1	38319.9
								...	...	...	...	...
Flare d	2.44 ± 0.04	2.9 ± 0.1	31326.1	2.24 ± 0.08	0.09 ± 0.03	2.5 ± 0.1	31322.0	2.31 ± 0.08	2.6 ± 0.1	1000 ± 300	2.6 ± 0.1	31323.3

Notes. ^aBLP fit by means of loglikelihood profile. ^b10⁻⁶ photons cm⁻² s⁻¹. ^c MeV.

Download table as:  ASCIITypeset image

2.2.1. Spectral Evolution

To complete the analysis of the spectral behavior of PKS 1510-089, we investigate whether spectral changes are seen between the quiescent and the flaring state. Since we are mainly interested in the search of possible trends rather than in the best description of the spectral distribution, we simply evaluated the PL spectral indices in the various brightness states, which can be considered representative of the mean slope. In the upper panel of Figure 5, we plot the photon index against the flux above 0.2 GeV, resulting from the same spectral analysis showed in Figure 1. In the case of daily integration (green circles), and more marginally for the weekly integration (red circles), this plot seems to show a softer when brighter trend, up to a flux level of F(E > 0.2 GeV)≃ 2.4 ×10⁻⁷ photons cm⁻² s⁻¹. Above this value of the flux, the source has a harder when brighter trend. We analyze the correlation of the harder when brighter trend, for the weekly binning, using a Monte Carlo method that takes into account the dispersion of flux and index measurements. In detail, we re-sample the flux and index values for each observed pair, extracting the data from a normal distribution centered on the observed value and with a standard deviation equal to the 1σ error estimate. We find a correlation coefficient of r = 0.43 with a 95% confidence limit 0.24 ⩽r⩽ 0.58. The trend for F(E > 0.2 GeV)≳ 2.4 ×10⁻⁷ photons cm⁻² s⁻¹ is reported in the inset of the upper panel of Figure 5.

Zoom In Zoom Out Reset image size

Although for some EGRET blazars, Nandikotkur et al. (2007) observed a similar flux-hardness anticorrelation at low fluxes, we need to take into account possible effects coming from the poor statistics when the source flux is low. As first, we note that moving to the weekly binning the trend is less evident, although the flux range is the same as that of the daily binning. As a further check, in the lower panel of Figure 5 we plot the photon index against the number of photons predicted by the best-fit model. It is clear that the dispersion of the photon index is related closely to the number of predicted photons, above N ≃ 20 the trend is the same for both the two integration timescales, and the photon index clusters around 2.5, without showing the very soft and very hard index values present at low number of events.

In conclusion, we cannot exclude the presence of a softer when brighter trend, for low flux levels and short timescales, but the statistical effects we present do not allow to obtain a purely physical interpretation. Similar results have been found for other Fermi Blazars (Abdo et al. 2010c), independently of their redshift or class (BL Lacs/FSRQs).

Even if we did not find a strong evidence for a index–flux correlation, the variation of the photon index returned by the gtlike fit as a function of time (see upper panel of Figure 1), shows that the dispersion on the photon index is larger until mid-March roughly, and gets narrower after. This feature is emphasized in Figure 6, where in the upper panel we plot the histogram of the photon index for the weekly integration, before MJD 54905 (corresponding to 2009 March 15, blue shaded histogram), and after MJD 54905 (red empty histogram). In the lower panel, we plot the same analysis for the case of daily integration. The distributions before and after MJD 54905 have the same mean, in both daily and weekly integrations (≃2.5), but very different standard deviations. In the case of daily integration, we have 0.45 and 0.21 before and after MJD 54905, respectively. In the case of the weekly integration, we have a standard deviation of 0.28 before MJD 54905 and 0.07 after. To test more quantitatively whether or not the distributions are different, we use a Kolmogorov–Smirnov (KS) test. We applied the test to the distributions of the spectral indices before and after MJD 54905, for the daily and weekly integrations. The test returns a p-value of ≃0.13 and ≃0.35, for the case of weekly and daily binning, respectively. The KS test gives only a marginal indication that the data sets, before and after MJD 54905, may not be drawn from the same distribution.

Zoom In Zoom Out Reset image size

In conclusion, the typical γ-ray photon index of PKS 1510-089 is around 2.5 and values largely different from this were never observed in high states. LP best fits indicate a significant but quite mild spectral curvature, that in the brightest flares (b, c, and d) was always very close to β = 0.1.

We analyzed XRT (Burrows et al. 2005; Gehrels et al. 2004) data using the xrtpipeline tool provided by the HEADAS v6.7 software package, for data observed in photon counting mode. Events in the 0.3–10 keV energy band were extracted, selecting grades in the range 0–12, and default screening parameters to produce level 2 cleaned event files were applied. Due to the low count rate (<2 counts s⁻¹), we did not find any signature of the pile-up effect.

We used data from the Burst Alert Telescope (BAT) on board the Swift mission to derive the spectrum of PKS 1510-089 in the 14–195 keV band. The spectrum is constructed by averaging the spectra of the source extracted over short exposures (e.g., 300 s) and it is representative of the sources emission over the five year time range spanned by the observations. These spectra are accurate to the mCrab level and the reader is referred to Ajello et al. (2008, 2009a, 2009b) for more details. PKS 1510-089 is bright in BAT and the approximate significance of the BAT spectrum used for this analysis is ≃13σ.

We followed the steps outlined in the UVOT User's Guide, to perform UVOT data reduction and analyses in all the six available filters (V, B, U, UVW1, UVM2, and UVW2). We started from the raw data stored in the HEASARC archive and we made sure that the sky coordinates were updated, the modulo-8 correction was applied, duplicated FITS extensions have been removed, and the aspect correction was calculated. Based on the active galactic nucleus (AGN) intensity, the optimal source extraction region is a 5'' circle. The background region is an annulus with inner–outer radii of 15''–27'', 27''–35'', depending on the filter used. In order to improve the astrometry, the NASA/IPAC Extragalactic Database (NED) position has been adjusted using the uvotcentroid task, and field of view sources have been excluded from the background region. The standard output of the uvotsource task has been used to extract the photometric light curves. We corrected the magnitudes for Galactic extinction assuming E(B − V)_Gal = 0.097 mag. This value was calculated from Schlegel et al. (1998) tables using tools provided by the NASA/IPAC archive.⁸⁵ The absorption for the other filters was calculated according to the extinction laws of Cardelli et al. (1989). The de-reddened magnitudes were converted into fluxes in physical units taking into account the zero points by Poole et al. (2008).

GASP is performing a long-term monitoring of 28 γ-ray loud blazars in the optical, near-infrared, millimeter, and radio bands (Villata et al. 2008, 2009b). The GASP has been following PKS 1510-089 since 2007 January, and contributed to MW studies involving γ-ray data from AGILE (Pucella et al. 2008; D'Ammando et al. 2009a). The optical and near-infrared GASP data for the present paper were acquired at the following observatories: Abastumani, Armenzano, Calar Alto, Campo Imperatore, Castelgrande, Crimean, Kitt Peak (MDM), L'Ampolla, Lowell (Perkins), Lulin, Roque de los Muchachos (KVA and Liverpool), Sabadell, San Pedro Martir, St. Petersburg, Talmassons, and Valle d'Aosta. Magnitude calibration was performed with respect to a common choice of reference stars in the field of the source from the photometric sequence by Raiteri et al. (1998). Conversion of magnitudes into de-reddened flux densities was obtained by adopting the Galactic absorption value A_B = 0.416 from Schlegel et al. (1998), consistent with the E(B − V) color excess, the extinction laws by Cardelli et al. (1989), and the mag-flux calibrations by Bessell et al. (1998).

The GASP millimeter-radio data were taken at Medicina (5, 8, and 22 GHz), Metsähovi (37 GHz), Noto (43 GHz), SMA (230 GHz), and UMRAO (4.8, 8.0, and 14.5 GHz).

The 2 cm VLBA/MOJAVE program (Lister et al. 2009b and references therein) has been monitoring PKS 1510-089 at 15 GHz with the Very Long Baseline Array (VLBA) since 1995. Method of observations, data processing, and imaging is discussed by Lister et al. (2009b). Typical resolution of these images is about or better than 3 pc.

In addition to the 15 GHz MOJAVE VLBA monitoring, single-epoch simultaneous multifrequency 5–43 GHz VLBA measurements were done on 2009 April 9, in support of the first year Fermi observations (Sokolovsky et al. 2010). Accuracy of flux density measurements is dominated by calibration uncertainties: about or less than 5% at 5, 8, and 15 GHz, about our less than 10% at 24 and 43 GHz.

We performed the spectral analysis of Swift-XRT data after grouping the photons to have a minimum number of 10 photons per bin, and we fitted the spectra by means of a photon PL distribution $F(E)=K E^{-\alpha _X}$, plus a Galactic absorption with an equivalent column density N_H = 7.88 × 10²⁰ cm⁻² (Lockman & Savage 1995). We extracted the X-ray spectrum for each pointing. The corresponding spectral analysis results are reported in Table 4.

The X-ray light curve obtained from the fluxes reported in Table 4 shows a modest variability if compared to optical and γ-ray flares (see Figure 7). We prefer to present light curves in terms of νF(ν) to make easier the comparison between the various bands and the SED changes. The average flux integrated in the 0.3–10.0 keV range is around 10⁻¹¹ erg cm⁻² s⁻¹, the lowest flux, recorded on 2009 January 16, was (6.5 ± 0.8) ×10⁻¹² erg cm⁻² s⁻¹. The highest flux, recorded on 2009 April 28, was (15.0 ± 1.5)×10⁻¹² erg cm⁻² s⁻¹. In this case, the flux increased by a factor of 2 within a day, and the spectrum reached the hardest state (α_X = 1.13 ± 0.13). This is the most relevant X-ray flaring episode in our data set and looking at the MW light curve in Figure 7 it seems to have no counterparts in other wavelengths. Since the statistics are low, we performed the spectral analysis using the Cash method (C-stat; Cash 1979) based on the use of a likelihood function. This method returns flux and photon index values that are compatible with those coming from the χ² method. Even if the two methods results are compatible, the significance of this flare is low (≃2σ), so we do not investigate possible physical implications.

Zoom In Zoom Out Reset image size

During our observations, the source spectrum was always hard, with a photon index ranging between about 1.3 and 1.6. The plot of the flux in 0.3–10.0 keV range versus the photon index (see Figure 8) is compatible with a harder when brighter trend. Using the Monte Carlo method described in Section 2.2.1, we get a correlation coefficient r = −0.31 with a 95% confidence limit of −0.55 ⩽ r ⩽ −0.05. This spectral trend is consistent with the same analysis performed by Kataoka et al. (2008). We note also that Kataoka et al. (2008) found a soft-X-ray excess in the Suzaku data, but the statistics of the individual pointings in our data set are not sufficient to detect such a feature.

Zoom In Zoom Out Reset image size

In order to increase the statistics and to look for differences between the different γ-ray flares, we produced X-ray SEDs averaged during the b, c, and d γ-ray flaring intervals, and during the post-d flaring period, reported in Figure 9. These SEDs show that the average state of the X-ray emission was almost steady, without drastic differences between the flares and the post-flare integration period. We also note that a possible soft-X-ray excess is visible in the post-b flare-averaged SED. The scatter plot in Figure 10 shows no correlation between the XRT flux and the LAT flux. This absence of correlation is relevant to the understanding of the emission scenario that we will discuss in Section 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The five year integrated BAT SED is plotted in Figure 9. The photon index, in the 14–150 keV band, is 1.37^+0.08_−0.19. Despite the long integration time of the BAT data, the photon index value is almost compatible with the range of values observed in the XRT data in our data set, and in other historical observations. This suggests that the X-ray and hardX-ray flux and spectral shape of this source are quite stable, or at least that our X-ray sampling is representative of the X-ray and hard X-ray shape on timescales of years.

Simultaneous Swift-UVOT and GASP observations provide a valuable data set to study the low-energy bump of the SED. In Figure 11, we plot the SEDs obtained from UVOT and GASP data, simultaneous within a daily timescale. These data show a minimum around the frequency of 5 × 10¹⁴ Hz, which does not seem to vary when the source is flaring. The frequency of the high-energy peak can be estimated well and is close to 10¹⁵ Hz, while that of the peak at low frequency cannot be well established and should be roughly estimated around 10¹³ Hz. The UV peak, as in many quasars, is likely due to the BBB that usually is understood as thermal emission from the accretion disk surrounding the BH. Assuming that the disk luminosity is steady or slowly variable compared to the synchrotron emission, we expect that the UV excess gets less and less evident, as the optical flux increases. To test this scenario, in Figure 12 we plot the ratio of the flux in the highest energy UVOT filter (UVW2) to the flux in the B UVOT filter, as a function of the optical R flux. This plot shows that the UV spectrum gets harder when the optical R flux is lower. This trend is consistent with the UV excess decreasing as the R flux is increasing, as expected in a scenario in which the UV bump originates in the thermal emission of the accretion disk. Despite this trend, we note that even in the highest state plotted in Figure 11, the corresponding UV bump is still prominent. For other FSRQs, such as 3C 454.3, the BBB is not visible during the high optical states (Giommi et al. 2006; Villata et al. 2006), since it is dominated by the synchrotron flux, as we would expect here in the case of PKS 1510-089. In panels (a) and (b) of Figure 13 UV light curves for all the six filters are shown, and in panel (c) we plot the ratio of the UVOT V filter flux to the UVOT W2 filter. The trend of the hardness ratio is clearly anticorrelated with that of the fluxes. Again, this supports the BBB scenario, with the UV spectrum harder when the BBB is more evident, namely, when the synchrotron flux is lower.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 14, we report the scatter plot of the Fermi-LAT flux F(E >200 MeV) versus the UVOT νF(ν) in the UVW2 filter, to check for a possible significant correlation between the γ-ray and UV fluxes. The correlation coefficient of the logarithms of the UV and γ-ray fluxes, obtained through the Monte Carlo method described in Section 2.2.1, is r = 0.2 with a 95% confidence interval of 0.05 ⩽ r ⩽ 0.34. If we exclude the point with the lowest γ-ray flux, the correlation coefficient is r = 0.05, suggesting that the overall correlation is not significant.

Zoom In Zoom Out Reset image size

This lack of correlation hints that the spectral evolution of the UV spectrum results from the contamination of the high-energy branch of the synchrotron emission, and that a change in the BBB luminosity is not the driver of the γ-ray luminosity variations. We will discuss this in Section 5.

The optical (R band) and near-IR (J, H, and K bands) observations span the period from 2008 May until the end of 2009 June. In Figure 15, we report the light curves of the optical data. The best sampled band is the R one, with more than 600 observations. The optical flaring activity increased dramatically after mid-2009 March, corresponding to the γ-ray flare c. On 2009 May 8, the source was very bright with an R magnitude of 13.60 ±  0.02 (Larionov et al. 2009a), and in 2 days it reached its historical peak at R = 13.07 ±  0.02, about 1.3 mag brighter than the previous record level, observed on March 27 (Larionov et al. 2009b) of the same year. This dense monitoring is fundamental to understand the correlation between the optical and the γ-ray activity, both for the long-term trends and the single flares. In Figure 16, we show the scatter plot of the Fermi-LAT flux F(E >200 MeV) versus the optical νF(ν) in the R filter. The correlation coefficient of the logarithm of the optical and γ-ray fluxes, evaluated through the Monte Carlo method described in Section 2.2.1, is r = 0.42 with a 95% confidence interval 0.36 ⩽ r ⩽ 0.46, higher than that found for the UV band. This finding hints that the driver of the flaring activity can be a change in the high-energy branch of the electron distribution. We will investigate this scenario more accurately in Section 5. Despite the statistical significance of this result, there is an evident dispersion in the scatter plot that could be related to the inter-band time lags. Indeed, temporal lags could be related to the internal source photon absorption, to the cooling time of the radiating particles, or to inhomogeneities in the emitting region. To search for possible time lag between the optical (R band) and the γ-ray band, we used the discrete cross-correlation function (DCCF) method. Taking into account the whole data set, the DCCF analysis returns a lag of 13 ± 1 days, with the γ-ray band leading the optical one (see Figure 17, top panels). Since the R flux increased dramatically after the end of 2009 February, we divide the data set in two sections, before and after MJD 54890 (corresponding to 2009 Feburary 28). The DCCF results, plotted in the middle panels of Figure 17, show that the 13 day delay is still apparent. The same result holds if we analyze individually the flaring sequences d and b (Figure 17, bottom panels). Looking at the flare light curves, it is clear that the 13 day lag is caused by a change in the relative flux of the two bands such that the γ-ray flux is stronger in the first part of the flare and the R-band flux gets stronger in the second half. We plot in the lower panel of Figure 18 the LAT light curve above 200 MeV, and in the upper panel the optical light curve, backward shifted by 13 days: the four γ-ray flares a, b, c, d, seem to have an optical counterpart. In particular, the peaks of the bright events in May would be very close in time to those in the LAT daily light curve at MJD ≃ 54948 (see Figure 2). To check the effect of this lag, we evaluate the correlation between the logarithms of the optical and γ-ray flux after applying the above time shift (see Figure 19) and find r = 0.62 with a 95% confidence interval 0.57 ⩽ r ⩽ 0.66. The values of the correlation coefficient and of its 95% confidence interval, after applying the 13 days time shift, are significantly higher than those obtained without the time shift. Moreover, Figures 19 and 16 show that time shift reduces the dispersion in the scatter plot.

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

This 13 day correlation between the optical and γ-ray emission appears in the different outbursts during this flaring activity but we do not have any indication that it is a characteristic behavior related to the temporal and energetic evolution of the flares. Only through long-term simultaneous γ-ray and optical observations we may better understand the actual level of randomness of this correlation and its possible physical meaning.

Radio data have a much lower sampling with respect to optical and γ-ray, but at 14.5 GHz and 37 GHz it is possible to follow the overall trend of the radio flaring activity. Even if radio fluxes do not show a correlation with the optical and the Fermi-LAT flaring trend, it is possible to identify a large radio flare, starting at about MJD 54900 (see Figure 7). This radio flare, visible at 14.5, 37, and 230 GHz, seems to start quasi-simultaneously with the γ-ray flare c at 37 GHz, and keeps increasing until the end of flare d. The 230 GHz light curve shows a structure similar to a plateau, starting when the γ-ray flare c has ended, and a possible plateau is present also at 37 GHz, starting at the end of flare d. The 14.5 GHz fluxes seem to lag behind the 37 and 230 GHz. Analysis in the radio band is much more complex than optical and γ-ray data, because of the synchrotron self-absorption and the longer cooling times, and the possibility of the overlapping of different flares.

In 2009 June and July, the integrated parsec-scale flux density had reached its historical maximum since 1995 (≃ 4 Jy), as the 15 GHz MOJAVE VLBA measurements show. The VLBI core flux density has also shown the highest value (see Figure 20, upper panel). If this major radio flare, observed in the core of the parsec-scale jet, is connected to the huge γ-ray flares which have happened in the first half of 2009, this determines the source of the high-energy emission to be located around the base of the parsec-scale jet following causality arguments, as suggested by Kovalev et al. (2009). The delay between the peaks of light curves in the γ-ray and radio bands can be, at least partly, explained by the synchrotron self-absorption of the emission at radio frequencies.

Zoom In Zoom Out Reset image size

In Figure 21, we present the three highest frequencies Stokes I parsec-scale images, from the 5–43 GHz VLBA measurements performed on 2009 April 9. The size of the bright parsec-scale core at 24 and 43 GHz is estimated to be about 60–70 μas or 0.3–0.4 pc. The core shows a flat radio spectrum (see Figure 22) indicative of a synchrotron self-absorbed region, while the first well-resolved jet feature is already optically thin radio spectral index s_r = −0.9 $(F(\nu)\propto \nu ^{s_r})$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The highest apparent speed of a component motion in the jet of PKS 1510-089 observed at 15 GHz by Lister et al. (2009c) is v_app = 24c which makes its jet highly relativistic and Doppler boosted, a typical case for γ-bright blazars (Lister et al. 2009a). Marscher et al. (2010a) report on two new knots, observed in VLBA images at 43 GHz. The first with an apparent speed of 24 ± 2c, passed the core on MJD 54674.5 ± 20 (2008 July 27). The second knot passed the core on MJD 54958.5 ± 4 (2009 May 7) which is consistent, within the estimated uncertainties, with the huge optical flare observed on 2009 May 8. A feature, coincident to the second knot, is also seen to emerge in the MOJAVE 15 GHz VLBA images from 2009 June through 2009 December. The fitted speed at 15 GHz (675 μas yr⁻¹) is somewhat slower than the Marscher et al. (2010a) speed (970 ± 60 μas yr⁻¹), but this could be due to blending with a third feature seen in the 43 GHz data to emerge sometime in 2009 June (Marscher et al. 2010b). Despite the high resolution of the 43 GHz images, a more accurate ejection date for the latter feature could not be obtained due its rapid angular evolution and close proximity (0.1 mas) to the feature ejected in 2009 May.

PKS 1510-089 was one of the brightest and most active blazars observed by Fermi-LAT during the first year of survey. Flaring episodes with timescales from weeks to months in addition to rapid and intense outbursts, were observed both at optical and γ-ray energies. The estimated isotropic luminosity above 100 MeV during the flare c was about 8 × 10⁴⁷ erg s⁻¹. This value represents the brightest flare-averaged state of the source. The brightest daily time-resolved luminosity, recorded during the flare c on 2009 March 26, was of ≃2 × 10⁴⁸ erg s⁻¹. During the flare d, variability timescales down to a fraction of a day were observed both at optical and LAT energies.

The rapid variability and the powerful γ-ray luminosity raise the problem of the pair production opacity. Indeed, without beaming effects, the source size estimated from the observed variability timescale (R_rad = cΔt/(1 + z)) makes the source opaque to the photon–photon pair production process, provided that γ-ray and X-ray photons are produced cospatially.

If the component emerged on 2009 May 7, observed also in the MOJAVE 15 GHz VLBA, is related to the flares a and b, then we can use the value of the beaming factor derived from the motion of radio knots, namely, δ ≃ Γ ≃ β_app ≃ 21. This estimate is compatible with other VLBI estimates (Homan et al. 2002), and is slightly larger than an alternative estimate based on the variability observed at 22 and 37 GHz, δ_var = 16.7 (Hovatta et al. 2009). Taking into account the most rapid timescale estimated in Section 2.1, Δt ≃ 0.25 days, the actual emitting region size results of the order of R_rad ⩽ δΔtc/(1 + z) ≃ 1 × 10¹⁶ cm. Using the above fastest timescale and the quasi-simultaneous observed X-ray flux F_X ≃ 8 × 10⁻¹² erg cm⁻² s⁻¹ (observed at a typical frequency of 10¹⁸ Hz), one can impose a limit on the minimum value of the beaming factor resulting in a source transparent to the photon–photon annihilation process (Maraschi et al. 1992; Mattox et al. 1993; Madejski et al. 1996). Combining the source size from Equation (1), the X-ray photon energy in the source frame, and the intrinsic X-ray luminosity with the optical depth expression, we get a model independent estimate of δ ≳ 8. This lower estimate indicates that the values derived from the VLBI images are well above the pair production transparency limit.

X-ray fluxes, spectral indices, and trends in our data set are compatible with those reported in previous analysis with several X-ray telescopes (BeppoSAX, ASCA, Suzaku; Kataoka et al. 2008). The harder when brighter trend and the lack of X-ray/γ-ray correlation are useful to constrain both the low-energy tail of the electron distribution and the emission scenario. The typical value of the soft-to-hard-X-ray photon index α_X is close to 1.4, with a quite narrow dispersion (≃0.1), and is similar to that observed in FSRQs with z > 2 (Page et al. 2005). Since in FSRQ objects the X-ray band samples the low-energy tail of the ERC component, the X-ray energy spectral index (s_X) constrains the slope p of the low-energy tail of the electron distribution in the range 1.6–2.0 (p = 2s_X + 1; e.g., Rybicki & Lightman 1979).

The high-energy spectral index of the electron distribution can be estimated from the optical/UV spectral shape. In Section 4.2, we showed that the optical/UV spectral index depends on the relative contamination between the synchrotron component and the BBB emission. A reliable estimate of the synchrotron spectral index could be achieved from the optical data, but, unfortunately, our spectral coverage with two or more simultaneous frequencies is only in a limited time window. In Figure 23, we plot the spectral index in the R to H bandpass versus the γ-ray flux integrated above 200 MeV. If we take into account the blue points which refer to the subflare peaking at about MJD 54909 (see Figure 2, panel (d)), we note that the γ-ray flux increased by about a factor of 2.5 with the optical spectral index almost constant around ≃1. This value of the optical spectral index corresponds to an electron energy distribution index of about 3. We note that the steadiness of the optical R − H spectral index, during this γ-ray subflare, is consistent with the spectral evolution of the γ-ray emission. Indeed, extracting the γ-ray photon indices simultaneous within 1 day with the optical ones, we note that also these values are almost stable around the value of 2.4. The decay of the flare occurred with a significant spectral softening, consistent with a cooling dominated regime, with a corresponding electron spectral index in the range between 3.5 and 4.0.

Zoom In Zoom Out Reset image size

The observed MW SEDs of PKS 1510-089, reported in Figure 24, show that during the flaring state, the IC peak dominated over the synchrotron one by more than 1 order of magnitude. Indeed, even if both the synchrotron and IC peak frequencies are not sampled in our data set, we can estimate their typical peak flux by extrapolating the LP best-fit model for the γ-ray data, and by fitting the radio, optical, and UV SED points by means of cubic function. In the case of the synchrotron component, simultaneous millimeter and optical/UV data indicate that the peak flux should not be higher than a few times 10⁻¹¹ erg cm⁻² s⁻¹. In the case of the IC component, since the spectrum is mildly curved with the peak energy below the threshold of our analysis (200 MeV), the peak flux must be a few times 10⁻¹⁰ erg cm⁻² s⁻¹. The presence of a pronounced BBB and the Compton dominance of about 10, during the flaring states, suggest that the contribution from photons originating outside the jet is appropriate to model the broadband SED.

Zoom In Zoom Out Reset image size

The absence of X-ray/γ-ray flux correlation hints that the ERC flux variations depend on a change in the high-energy spectral index of the electron distribution, instead of a change in the external radiation field. In Section 4.2, we noted that the lack of γ-ray/UV correlation suggests that BBB variations are excluded as a main driver of the γ-ray variability. We tested the change of electron distribution by studying the correlation between the logarithms of the γ-ray and of the optical R fluxes (see Figure 16). The relation between the γ-ray and optical flux, fitted by means of PL, returns an exponent of ∼0.5 (or ∼0.6 if we use the 13 day shifted optical light curve discussed in Section 4.2). Despite the large scatter, this value is very interesting, because synchrotron and IC fluxes correlate very differently in the case of SSC, ERC/BLR, and ERC/DT. Moreover, the correlation depends on the low- and high-energy range chosen, respectively, for the synchrotron and IC component (Katarzyński et al. 2005). To have an estimate of the PL exponent for this correlation, we reproduced numerically a set of SEDs comparable to that observed for PKS 1510-089, using a well-tested code (Tramacere 2007; Tramacere et al. 2009; Tramacere & Tosti 2003). We calculated an ERC-dominated model, where the γ-ray emission is dominated by the BLR contribution and the electron distribution is a BPL. We fixed all the parameters and we increased the high-energy index of the electron distribution from −3.5 to −2.5 to evaluate the correlation between the optical energy emission (νF(ν) at 10¹⁴ Hz) and the integrated F(E > 200 MeV) γ-ray flux, taking separately into account the three IC components SSC, ERC/BLR, and ERC/DT, as shown in Figure 25. The resulting flux relations were fitted by simple PLs and the resulting exponents were 3.1 for the SSC, 1.8 for ERC/DT, and 0.4 for ERC/BLR. The last one was the only one found to be consistent with that observed for PKS 1510-089 (≃ 0.5–0.6) and therefore, we can disfavor both SSC and ERC/DT.

Zoom In Zoom Out Reset image size

The SEDs of the ERC/BLR model in the upper panel of Figure 25 have a curvature more pronounced than the ERC/DT ones. This is due to the KN cross section, because UV photons emitted from the accretion disk and reflected toward the jet by the BLR are blueshifted roughly by a factor of Γ. For Γ≃ 10, the typical energy of these photons in the emitting region rest frame is then ≃ 10¹⁶ Hz, hence the IC scattering with the electrons with γ≃ 1000 occurs under a mild KN regime. The resulting smooth curvature has a value of ≃0.1 that is compatible with the observed one.

We attempt a leptonic ERC/BLR-oriented SED modeling. We aim to reproduce the SEDs for the three flares with simultaneous data from radio to γ-ray energies, namely, flares b, c, and d. Moreover, we try also to fit the quiescent state. Although variability timescales reached values of a fraction of day, LAT data required longer integration times to produce an SED enough good for spectral modeling. The intermediate fluxes of the synchrotron and BBB components, observed during the γ-ray integration period (red filled circles in Figure 24), are used in the fitting procedure as representative of the flare-averaged state of the low-energy bump.

Since the correlation between the γ-ray and the optical fluxes (see Section 5.1 and Figure 25) favors an ERC/BLR scenario, we assume that the dissipation zone is in the subparsec scale. This is also consistent with the mild curvature observed in the γ-ray spectra. Indeed, as shown in the previous section, the ERC/BLR process occurs under the KN regime leading to the curved MeV/GeV spectral shape that matches the one observed in the Fermi spectra.

We assume a jet viewing angle of θ = 25 for both the flaring and the quiescent states. During the flaring states, we choose a bulk Lorentz factor in the range [14–16], resulting in a beaming factor range of [20–21.5] that is compatible with the VLBI observations. During the quiescent state, we use a bulk Lorentz factor of 12.0, corresponding to a beaming factor of about 18.8.

As a further step, we estimate the accretion disk physical characteristics using the UV data. We use the UV observations during the lowest synchrotron state of our data set as an estimate for the upper limit of the accretion disk luminosity (L_d). According to the observed UV flux and to the luminosity distance, we set a reference value of L_d ≃ 3 × 10⁴⁵ erg s⁻¹. We model the accretion disk, following the prescription in King (2008) and rearranging the expression as in Ghisellini et al. (2009), using a multi-temperature blackbody with a temperature profile given by

$$\begin{equation} T_{\rm disk}^4(R)=\frac{3 R_S L_d }{16 \pi \epsilon \sigma _{\rm SB} R^3} \Big [1-\Big (\frac{3 R_S}{R}\Big)^{1/2} \Big ], \end{equation} \tag{ 2 }$$

where σ_SB is the Stefan–Boltzmann constant, R_S = 2 GM_BH/c² is the Schwarzschild radius for a BH mass M_BH, 3R_S is the last stable orbit in the case of Schwarzschild BH, and is the accretion efficiency that is linked to the bolometric luminosity and to the accretion mass rate $\dot{M}$, by $L_d=\epsilon \dot{M}c^2$. We assume that the radiative region of the disk extends from ≃3R_S to ≃500R_S. Since T_disk(R) peaks at R ≃ 4R_S, we can use the UV data to constrain the peak of T_diks(R) and we can use Equation (2) to estimate R_S. From our data, we get a value of T^peak_disk ≃ 4 × 10⁴ K. Assuming as a reference value for the accretion efficiency ≃ 0.1, we get R_S ≃ 1.6 × 10¹⁴ cm. The BH mass M_BH is then about 5.4 × 10⁸ M_☉ and the accretion rate of about 0.5 M_☉ yr⁻¹ that corresponds to ≃0.04 times the Eddington accretion rate $(\dot{M}_{\rm Edd})$. The value of the BH mass estimated from our UV data is compatible with that obtained by Oshlack et al. (2002) using the virial assumption with measurements of the Hβ FWHM and luminosity (M_BH≃ 3.9 ×10⁸ M_☉), and by Xie et al. (2005; M_BH≃ 2 × 10⁸ M_☉).

The radius of the BLR can be estimated using the empirical relation by Kaspi et al. (2005) and Bentz et al. (2006):

$$\begin{equation} R_{\rm BLR}\simeq 20\Big [ \frac{\nu ^* L(\nu ^*)}{10^{44}\, \hbox{erg}\, \hbox{s}^{-1}} \Big ]^{0.5} \hbox{lt-days}, \end{equation} \tag{ 3 }$$

where ν* ≃ 5.878 × 10¹⁴ Hz. According to the UVOT data, we estimate ν*L(ν*) < 2 × 10⁴⁵ erg s⁻¹, and we get a value of R_BLR < 2.3 ×  10¹⁷ cm. We use BLR radius of about 1.6 ×  10¹⁷ cm with reflectivity value of 0.1 (Ghisellini et al. 2009). The DT is modeled as a BB with a temperature T_DT of about 100 K, and a reflectivity of 0.3 (Ghisellini et al. 2009). The distance is set to a typical value of ≃10¹⁸ cm, and it is fine tuned in the fit in order that the BLR/DT+BLR/ERC correctly match the X-ray and hard X-ray data.

We use an emitting region size (R_rad) in the range [2–5] × 10¹⁶ cm that is almost twice the value estimated from the fast variability. In this regard, we note that since we are describing flare-averaged states, with integration times of the order of a few weeks, the discrepancy with the fast variability estimate is not problematic. As electron energy distribution N(γ) we use a BPL:

$$\begin{equation} N(\gamma)\propto \left\lbrace \begin{array}{@{}ll} \gamma ^{-p} \quad \mbox{for}\quad \gamma _{\rm min}\le \gamma <\gamma _{br}\\ \gamma ^{-p1} \quad\mbox{for\quad} \gamma _{br}\le \gamma \le \gamma _{\rm max} \end{array}\right. \end{equation} \tag{ 4 }$$

with the number density of the electrons

$$\begin{equation} N_e=\int _{\gamma _{\rm min}}^{\gamma _{\rm max}} N(\gamma)d\gamma. \end{equation} \tag{ 5 }$$

The low-energy spectral index (p) is chosen ≃1.9, as hinted by the typical X-ray photon index (see Section 5.1); the value of the break energy (γ_br) is chosen to match the position of the peak energy of the IC bump. It is of the order of [250–300] and it is tuned in the three flares to fit the data. The high-energy spectral index (p₁) is chosen according to the average spectral index of the LAT data, and the reference value is ≃3. Magnetic field intensity B is chosen to be 1.0 G for flares b and c, while for the flare d we need to use a larger value of B and a more compact region size.

All the flares are assumed to occur at dissipation distance from the disk (R_diss) that is within the BLR and the DT (R_diss < R_BLR, R_diss < R_DT). The resulting best fit of the MW SEDs is reported in Figure 24, and the corresponding values of the best-fit parameters for the three flares are reported in Table 5.

To have an indication of the change in the energetic contents of the jet as a function of the flaring activity, we try to fit also the quiescent γ-ray SED. Considering the lack of MW data during this time interval, we assume that the X-ray SED is close to the state observed during flare b, as supported by the low X-ray variability. From the MW light curves reported in Figure 7, we infer that the optical/UV flux is at a level of few times 10⁻¹² erg cm⁻² s⁻¹. We can have an acceptable description of the MW SED during the quiescent γ-ray period using parameters comparable to those chosen for the flare b, but we need to decrease significantly γ_br to ≃65. The best-fit model is represented by the black dash-dotted line in the top panel of Figure 24.

It is interesting to study the evolution of the energetic content of the jet resulting from the SED modeling discussed above. We evaluate the total kinetic power of the jet as

$$\begin{equation} L_{\rm jet}=\pi R_{\rm rad}^2 c \Gamma ^2(U^{\prime }_e + U^{\prime }_B + U^{\prime }_p), \end{equation} \tag{ 6 }$$

where the rest-frame magnetic energy density is given by U'_B = B²/8π, the electron energy density is U'_e,⁸⁶ and U'_p = 0.1N_em_pc² is the cold-proton energy density, assuming that we have 1 cold proton per 10 electrons (Sikora et al. 2009). The power carried by the jet in terms of radiation is given by L_rad ≃ L'Γ²/4. We evaluate L' summing up the numerical integration of each radiative component (synchrotron, SSC, ERC/BLR, ERC/DT) as observed in the jet rest frame. In Table 6, we report the results for flares b, c, and d and the quiescent state, corresponding to the physical parameters reported in Table 5.

The total kinetic power of the jets is almost steady, except during the flare d when it increased by a factor of two. Indeed, for the case of flares b and c, we find values of L_jet that are comparable with that of the quiescent state. In contrast, we note that the value of the electron energy density during all the flaring states is much larger than that estimated during the quiescent state, and, during the flare d, it increased with respect to the quiescent state by about a factor of 2.7. The low values of U'_e/U'_B reflect the low SSC contribution. Indeed, for the choice of our model parameters, the SSC is always negligible compared to the other radiative components. As final remark, we note that the radiative efficiency of the jet (L_rad/L_jet), is ≃0.03 during the quiescent state, and increases up to ≃0.2 during flare c.

We now compare the jet total kinetic luminosity with the accretion disk luminosity. According to the analysis reported in Ghisellini et al. (2010), if the jet power comes from the accretion process, then the accreting mass has to account both for the disk luminosity and for jet power. We can write

$$\begin{equation} \epsilon _{\rm tot}\dot{M}c^2=\epsilon _{D}\dot{M}c^2+\epsilon _{\rm jet}\dot{M}c^2. \end{equation} \tag{ 7 }$$

The ratio of L_jet to L_D, reported in Table 6 is equal to _jet/_D. Using as _D the typical value of 0.1, we need a total efficiency _tot ≳ 0.26 that is still compatible with the maximum efficiency in the case of a Kerr BH.

It is interesting to compare the flaring activity of PKS 1510-089 with other FSRQs observed by Fermi. 3C 454.3 (z = 0.859), reached an isotropic luminosity above 100 MeV L(E > 100) ≃ 8 × 10⁴⁸ erg s⁻¹ during the period of 2008 August/September (Abdo et al. 2009a). It became the highest ever recorded object in the gamma rays during the flare in 2009 December (Escande & Tanaka 2009), with a luminosity L(E > 100) ≃ 1 × 10⁴⁹ erg s⁻¹, about five times larger than the maximum flare resolved luminosity of PKS 1510-089. The source during the flares exhibited rapid variability down to sub-daily timescales that is comparable to that observed in PKS 1510-089. The γ-ray spectra of 3C 454.3 object show a spectral break around 2 GeV (Abdo et al. 2009a). In this case, we see a relevant difference with respect to PKS 1510-089 that exhibited a mild curvature spectrum. A possible reason, as explained in the previous section, is that in the case of 3C 454.3 the dominant external photon field originates in the DT and not in the BBB. Thus, the IC process occurs under the TH regime, and the break reflects the break in the electron distribution. The BH mass is estimated to be M_BH ≃ 4 × 10⁹ M_☉, and the disk luminosity of the order of L_D ≃ 2 × 10⁴⁶ erg s⁻¹, both these values are 1 order of magnitude larger than PKS 1510-089. Another interesting difference is that the 3C 454.3 exhibits no BBB features during the high optical and γ-ray outburst.

The FSRQ PKS 1502+106 (z = 1.839), during the outburst in 2008 August (Abdo et al. 2010b), reached a luminosity L ≃ 1 × 10⁴⁹ erg s⁻¹, showing a very large Compton dominance up to 100. As in the case of 3C 454.3, this object exhibited no BBB feature at UV energies. The flare of this object was an isolated episode, rather than a flaring sequence. The γ-ray spectrum showed a curved shape, similar to PKS 1510-089, during the full integration period, and during the post-flare state. The curvature values, 0.1–0.2, are comparable to those observed in PKS 1510-089, hence also in this case an origin in the KN effect is possible, as well as a curved electron distribution with the IC process in TH regime. The MW variability of this object is very different from that observed in PKS 1510-089. Indeed, the X-ray flux followed both the γ-ray and the optical flux. Moreover, there was a strong correlation between the UV and the γ-ray. Probably this can be explained with a sharp change in the density of the radiating electrons, or in the beaming factor.

PKS 1454-354 (z = 1.424), another distant FSRQs discovered in the γ-ray by Fermi during the 2008 August/September flaring activity (Abdo et al. 2009b), reached a luminosity L(E > 100) ≃ 7 × 10⁴⁸ erg s⁻¹, and exhibited a rapid variability with hours/days timescales.

We conclude that the most rapid timescale for the γ-ray variability (down to 6/12 hr) is common to this class of objects, despite the bolometric γ-ray luminosity. This timescale is limited by the minimum LAT integration time required to extract a flux for object with F(E > 200 MeV) of the order of 10⁻⁵ photons cm⁻² s⁻¹, hence a faster variability could be possible. Compared to those listed above, PKS 1510-089 seems to be a less powerful object in terms of both γ-ray luminosity and thermal luminosity. The curvature in the γ-ray spectrum, when compared to the spectral break observed in other sources such as 3C 454.3, may be understood as a signature of the KN regime, hence of a dissipation region within the BLR.