IMPLICATIONS OF THE ANOMALOUS OUTBURST IN THE BLAZAR PKS 0208−512

The observing strategy of the Fermi Large Area Telescope (LAT) of scanning the sky every three hours—combined with supporting, multi-wavelength monitoring by numerous research groups—continues to provide flux and spectral variability information for a large sample of blazars in unprecedented detail (e.g., Nolan et al. 2012). In the second LAT active galactic nucleus (AGN) catalog, ≳95% of all sources are confirmed or candidate blazars (2LAC; Ackermann et al. 2012). Modeling of truly simultaneous time-variable spectral energy distributions (SEDs) from radio to γ-rays is the holy grail of blazar physics. Such analysis for a large sample of blazars is now possible with these data. In this Letter, we model the optical–near-infrared (OIR) to γ-ray SED of the blazar PKS 0208−512 during 2008–2011 to investigate the physical parameters related to its emission.

PKS 0208−512 is a flat spectrum radio quasar (FSRQ) at redshift z = 1.003 (Healey et al. 2008). It is detected at a 36σ level in the Fermi 2 yr catalog (2FGL; Nolan et al. 2012) and was regularly observed by the Yale/SMARTS optical/near-infrared monitoring program (Bonning et al. 2012; Chatterjee et al. 2012). Chatterjee et al. (2013, hereafter Paper I) showed that between 2008 August and 2011 September, PKS 0208−512 underwent three OIR outbursts of at least 1.3 mag and spanning one month or more. In contrast, at GeV energies the source shows similar flares only during intervals 1 and 3.

During the Fermi era, the GeV and OIR variability of blazars have been shown to be well correlated in most cases, which is consistent with the so-called "leptonic scenario." The OIR emission in blazar jets is believed to be due to synchrotron radiation from the relativistic electrons in the jet (Urry & Mushotzky 1982; Impey & Neugenbauer 1988; Marscher 1998). In the leptonic model, the γ-rays are generated by the same electron distribution by inverse-Compton upscattering of synchrotron photons from the jet itself ("synchrotron self-Compton" or SSC process; Maraschi et al. 1992; Chiang & Böttcher 2002) or external photons from the accretion disk, broad emission lines, or the torus ("external-Compton" or EC process; Sikora et al. 1994; Błażejowski et al. 2000; Dermer et al. 2009). In both SSC and EC scenarios, a strong correlation between the variations in OIR and γ-ray emission is predicted since they are produced by the same electrons. Hence, the OIR flare without any corresponding GeV variability observed for PKS 0208−512 appears to be anomalous and deserves a more detailed study. In this Letter, we model the OIR to γ-ray SED of the blazar PKS 0208−512 during the three OIR outbursts to identify the physical scenario causing the anomalous flare.

We analyzed data from Fermi/LAT using the standard Fermi Science Tools (v9r27p1) to derive the γ-ray spectra of PKS 0208−512. We included all sources within 15° of PKS 0208−512 extracted from the 2FGL. We kept their normalizations free and spectral indices frozen to their catalog values. We selected the good time intervals by using the logical filter "DATA_QUAL==1 & LAT_CONFIG==1 & ABS(ROCK_ANGLE) < 52." These data were analyzed with an unbinned likelihood analysis method using the standard analysis tool gtlike. We use the currently recommended set of instrument response functions (P7SOURCE_V6), Galactic diffuse background model (gal_2yearp7v6_v0.fits), and isotropic background model (iso_p7v6source.txt). To obtain spectra, we divide the photon energy range of 0.1–20 GeV into five spectral bins and model PKS 0208−512 with a simple power law in each bin, with the spectral index and normalization kept free. The spectral bins were decided by optimizing two quantities—namely, resolution and signal-to-noise—in each bin. The resultant spectral data are shown in Figure 1.

Zoom In Zoom Out Reset image size

The goal of the SED analysis in this Letter is to understand (1) what physical processes could explain the differences in the γ-ray/OIR luminosity ratio between the three flares, and (2) the relation between the OIR colors and the OIR fluxes. Hence, we calculate the spectra during these three intervals. The three intervals, as defined in Paper I, included some low-state data points. This is not suitable for the current analysis since our goal is to model the SEDs during the OIR outbursts. The intervals are now defined as (1) MJD 54710–54840, (2) MJD 55170–55230, and (3) MJD 55650–55850, such that the low-state data points are excluded. The length of the intervals are rather long which is useful in obtaining good signal-to-noise ratios in the spectral data. However, even with such long time intervals, the source is not significantly detected (with a detection criterion that the maximum-likelihood test statistic exceeds 25) in the highest energy spectral bin 8.1–20 GeV. In addition to these three intervals, we have calculated the OIR SED during a very low state MJD 55870–55915. Assuming that the contribution of the variable jet flux is extremely low during this time interval, the SED is dominated by thermal emission from the accretion disk and may thus yield a relatively precise estimate of the disk emission. That estimate is crucial in setting the scale of the broad-line region (BLR) and the dusty torus, which should be the main sources of external radiation for Comptonization into γ-rays.

All measurements in the B, V, R, J, and K bands were obtained with the ANDICAM instrument on the SMARTS 1.3 m telescope located at CTIO, Chile (Bonning et al. 2012). We calculate the OIR SEDs by averaging the flux over each of the three intervals including only nights when data in all five bands are available.

We use the Swift/X-ray telescope (XRT) data products generator (Evans et al. 2009) to determine the X-ray spectra in the 0.3–10 keV energy range during the three relevant intervals.

The selected SEDs of PKS 0208−512 are modeled with the leptonic radiative code Blazar (Moderski et al. 2003), which includes synchrotron emission and Comptonization of synchrotron and external radiation fields. We begin by determining the luminosity of the thermal component produced by the accretion disk. We use the composite spectrum of radio-loud quasars from Elvis et al. (1994) and determine its normalization to best fit the observed low-state OIR SED. As can be seen in Figure 1, and as was discussed in detail in Paper I, PKS 0208−512 shows a systematic softening of the OIR spectrum with increasing OIR flux level, which is the redder-when-brighter behavior observed in other FSRQs (e.g., Villata et al. 2006; Bonning et al. 2012). We checked whether this spectral transition can be explained by the superposition of synchrotron and thermal SED components. We identified the normalization of the composite spectrum, for which the observed spectral transition is accurately reproduced by solely varying the normalization of the synchrotron component. By integrating the luminosity of the UV-peaking thermal component of the normalized composite spectrum, shown in Figure 1, we find an accretion disk luminosity of L_d ≃ 8 × 10⁴⁵ erg s⁻¹. Our value is lower by a factor of ≃ 2.2 than the value adopted by Ghisellini et al. (2011), which can be explained by variability of the accretion disk emission on yearly timescales.

We estimate the characteristic radius of the BLR r_BEL ≃ 0.091 pc and the inner radius of the dusty torus r_IR ≃ 2.2 pc (assuming a dust sublimation temperature of T_IR = 1200 K) using scaling relations from Sikora et al. (2009). The energy densities of external radiation components depend on the assumed covering factors ξ_BEL and ξ_IR, for which we try different values.

We scan the OIR light curves for pairs of initial (t_i) and final (t_f) epochs such that F_f/F_i > 2, where F_i and F_f are the corresponding fluxes, and calculate the observed variability timescale t_{var, obs} = (t_f − t_i) × ln 2/ln (F_f/F_i). We find the shortest t_{var, obs} ≃ 2 days during flares 1 and 3. We constrain the location of the emitting region r, and the jet Lorentz factor $\Gamma _{\rm j}=(1-\beta _{\rm j}^2)^{-1/2}$, where β_j = v_j/c is the dimensionless jet velocity, corresponding to different values of the collimation efficiency parameter Γ_jθ_j, where θ_j is the jet opening angle (K. Nalewajko et al., in preparation; see Appendix A):

$$\begin{equation} \Gamma _{\rm j}(r,\Gamma _{\rm j}\theta _{\rm j}) \simeq \left(\frac{\mathcal {D}}{\Gamma _{\rm j}}\right)^{-1/2}\left[\frac{(1+z)\Gamma _{\rm j}\theta _{\rm j} r}{ct_{\rm var,obs}}\right]^{1/2} \,, \end{equation} \tag{ 1 }$$

where $\mathcal {D}=[\Gamma _{\rm j}(1-\beta _{\rm j}\cos \theta _{\rm obs})]^{-1}$ is the Doppler factor, and θ_obs is the viewing angle. There are strong theoretical and observational indications for Γ_jθ_j ⩽ 1 in AGN jets (Pushkarev et al. 2009; Komissarov et al. 2009). This constraint is illustrated in Figure 2, and one can see that it favors locations outside the BLR, although solutions located at r_BLR cannot be entirely excluded.

Zoom In Zoom Out Reset image size

A further constraint can be placed using the relation between r and Γ_j corresponding to a fixed value of the SSC luminosity L_SSC (K. Nalewajko et al., in preparation; see Appendix B):

$$\begin{eqnarray} \Gamma _{\rm j}(r,L_{\rm SSC}) &\simeq & 1.27\left(\frac{\mathcal {D}}{\Gamma _{\rm j}}\right)^{-1}\left[\left(\frac{L_{\rm syn}}{L_{\rm SSC}}\right)\left(\frac{L_\gamma }{\zeta L_{\rm d}}\right)\right]^{1/8} \nonumber \\ &&\times \left[\frac{(1+z)r}{ct_{\rm var,obs}}\right]^{1/4}\,, \end{eqnarray} \tag{ 2 }$$

where L_syn is the synchrotron luminosity, L_γ is the γ-ray luminosity, and the parameter ζ is a function of ξ_BEL and ξ_IR. The observed X-ray luminosity of PKS 0208−512 is L_X ≃ (3.5 ± 0.3) × 10⁴⁵ erg s⁻¹, but the observed X-ray spectrum is too hard to be dominated by the SSC component. It is thus required that L_SSC < L_X. This condition is illustrated in Figure 2, and it strongly constrains the jet Lorentz factor. For the subsequent modeling, we choose a particular solution with r = 2 pc and Γ_j = 20, which is allowed by the constraints obtained from both flares 1 and 2. In many blazar studies, the Doppler-to-Lorentz factors ratio is often assumed to be $\mathcal {D}/\Gamma _{\rm j} \simeq 1$. However, in our model located at r = 2 pc, adopting $\mathcal {D}/\Gamma _{\rm j} = 1$ would require Γ_j > 20. Higher values of the ratio are possible up to $\mathcal {D}/\Gamma _{\rm j} \lesssim 2$ for very narrow jets pointing exactly at the observer. Here, we adopt an intermediate value—$\mathcal {D}/\Gamma _{\rm j} = 1.4$—which corresponds to a θ_obs ≃ 19, a jet opening angle of θ_j ≃ 067, and a collimation factor of Γ_jθ_j ≃ 0.235.

Our next step is to model the SED observed during flare 3, which shows the highest OIR flux and a strong indication of a spectral peak within the Fermi/LAT band, at ≃ 800 MeV. We inject ultrarelativistic electrons with a broken power-law spectrum $N_\gamma \propto \gamma ^{-p_i}$, with spectral indices p₁ for γ < γ_br and p₂ for γ > γ_br, where γ_br is the break Lorentz factor. The observed OIR spectrum can be used to precisely determine p₂, and the γ-ray spectral peak can be used to determine γ_br. There is no strong constraint on p₁, so we take the maximum value allowed by the observed X-ray emission. The observed SED places strong constraints on the Compton dominance parameter q = L_EC/L_syn, and the peak frequencies ratio w = ν_EC/ν_syn. It can be shown that—when the Compton scattering takes place in the Thomson regime—every component of external radiation has a covering factor that scales as ξ∝q/w², independent of B' or Γ_j (Sikora et al. 2009; Nalewajko et al. 2012). If we fix ξ, L_EC, and ν_EC, we find that $L_{\rm syn}\propto \nu _{\rm syn}^2$, a relation which is roughly perpendicular to the observed OIR spectrum, and thus the synchrotron peak can also be robustly determined. We attempted to fit the observed SED for different values of ξ_IR. We find that ξ_IR ≲ 0.1, otherwise the value of w is too small as compared with observations. This value of the covering factor is typical for blazar emission models.

Finally, we model the SEDs of flares 1 and 2, assuming that they also are produced in the IR region. These flares have very similar average OIR spectra, but they differ by almost an order of magnitude in γ-ray luminosity. The γ-ray spectrum of flare 1 shows hints of a complex structure, which might require two spectral bumps, one peaking below 100 MeV, and another peaking at ∼1 GeV. We attempted to model this SED using a complex electron energy distribution, but we did not obtain a solution that would be acceptable both in the γ-ray and OIR bands. We note that a two-zone model, with differences in either location, magnetic field strength, and/or bulk Lorentz factor, may be required to explain the details of this SED. Instead, we present a model using a broken power-law electron distribution with a cutoff at γ_max, which reproduces the general shape of the γ-ray spectrum. The γ-ray spectrum of flare 2 is poorly constrained; therefore, we focus on the OIR spectrum. Our best-fit models of flares 1, 2, and 3 are shown in Figure 1 and the corresponding parameters are presented in Table 1.

Table 1. Parameters of the SED Models of Three Flaring States of PKS 0208−512

Flare	r	Γj	ξIR	B'	p1	p2	γbr	γmax	$u_{\rm e}^{\prime }$	Ne(1054)
	(pc)			(G)					(erg cm−3)
1	2	20	0.1	0.17	1.5	2.55	850	12000	0.0145	1.03
2	2	20	0.1	0.34	2.05	3.7	3000	...	0.0103	3.47
3	2	20	0.1	0.275	1.6	4.1	4400	...	0.0063	0.52

Download table as:  ASCIITypeset image

In Table 1, we also report the comoving energy density of electrons $u_{\rm e}^{\prime }$ and the total number of electrons N_e. The numbers suggest that flare 1 requires the highest $u_{\rm e}^{\prime }$, and flare 2 requires the highest N_e. These numbers are, however, strongly dependent on the low-energy slope of the electron distribution p₁, which is the most uncertain parameter of SED models, especially for flare 2.

The main problem that we address in this work is the difference between the Compton dominance q among the three flares of PKS 0208−512. In the framework of EC models, we expect that $q = L_\gamma /L_{\rm IR} \simeq L_{\rm EC} / L_{\rm syn} \simeq u_{\rm ext}^{\prime }/u_{B}^{\prime }$. Our modeling results correspond to q₁ ≃ 20.9, q₂ ≃ 5.2, and q₃ ≃ 8.0 for flares 1, 2, and 3, respectively. Thus, the difference between q values for flares 1 and 2 is a factor of ≃ 4. In order to explain the differences between the three flares of PKS 0208−512, the magnetic field strength should vary by a factor of ≃ 2, assuming $u_{\rm ext}^{\prime }$ to be constant. Alternatively, we can assume that the magnetic field profile is fixed, with B'∝r⁻¹, which means that $u_{B}^{\prime } \propto r^{-2}$. Since the comoving energy density of external radiation fields is expected, in general, to have a different distance dependence, $u_{\rm ext}^{\prime } \propto r^{-a^{\prime }}$, we can explain these differences by varying the location of the emitting region. For example, in the case that a' = 0 or a' = 4, distance variations by a factor of ≃ 2 are required.

Either requirement for the differences in q between the three flares of PKS 0208−512 is relatively easy to satisfy. So, it is interesting to consider how common such variations in q might be between different flares produced by the same blazar, or between different blazars. Although no systematic study of the relative amplitudes of correlated γ-ray and optical flares in FSRQs has been conducted to date, our experience in working with the SMARTS and Fermi/LAT data suggests that flare 2 is indeed anomalous. This implies that either the distance at which such flares are produced is finely tuned or most flares are produced in jet regions where a' ≃ 2. In light of the generally irregular behavior of blazars, we are inclined to favor the latter scenario.

The spatial distribution of the external radiation fields that are Comptonized in blazars remains poorly understood. Typically, it is assumed that the reprocessing medium is concentrated at a characteristic radius r_ext (standing for r_BEL for the BLR and r_IR for the dusty torus). Then, the distribution of external radiation in the external frame can be approximated by a broken power-law distribution, with a ≃ 0 for r < r_ext and a ≃ 2 for r > r_ext. In fact, in the external frame where external radiation sources are non-relativistic and external radiation fields are not subject to extinction, we have a ⩽ 2, and this limit is realized for point sources. In the comoving frame, the distribution of external radiation can be steeper due to the Lorentz transformation $u_{\rm ext}^{\prime } \simeq \mathcal {D}^{\prime 2}u_{\rm ext}$, where $\mathcal {D}^{\prime } = \Gamma _{\rm j}(1-\beta _{\rm j}\cos \psi)$ is a Doppler factor defined in the comoving frame, and ψ ≃ r_ext/r for r ≫ r_ext is the angle between the velocities of the external radiation photons and the emitting region in the external frame. For large Γ_j and small ψ, one can approximate $\mathcal {D}^{\prime } \simeq (\Gamma _{\rm j}/2)(1/\Gamma _{\rm j}^2 + r_{\rm ext}^2/r^2)$, which highlights that there are, in fact, two regimes. For r ≪ Γ_jr_ext, one has $\mathcal {D}^{\prime } \simeq \Gamma _{\rm j}r_{\rm ext}^2/(2r^2)$, and thus a' ≃ 6. For r ≫ Γ_jr_ext, one has $\mathcal {D}^{\prime } \simeq 1/(2\Gamma _{\rm j})$, and thus a' ≃ 2. The latter regime is not very practical, because at that point the external radiation is strongly deboosted. In the former regime, one can obtain a' < 6 by assuming a significant spatial distribution (stratification) of the external medium emissivity j_ext(r). A distribution of external radiation with a' ≃ 2 can be obtained by requiring (approximately) that j_ext∝r⁻².

In Paper I, it was speculated that the OIR outburst during interval 2 was caused by a change in the magnetic field without any change in the Doppler factor of the emitting region or total number of emitting electrons. This will cause a synchrotron flare without any variability in the EC emission. Similarly, if the magnetic field is stronger at the location of flare 2, this will result in an OIR outburst much more significant than the GeV flare and the latter may not be detected. Our analysis in this Letter has demonstrated that such a scenario is a possible explanation for the anomalous flare.

An alternative scenario was discussed in Paper I. A smaller Γ_j at the location of flare 2 will result in a smaller Compton dominance causing the OIR to be more dominant than the GeV emission. The constraints on r and Γ_j, shown in Figure 2, are inconsistent with such a scenario in which flare 2 is produced in the inner region of the jet, where the jet is still undergoing acceleration and Γ_j is lower than that for the other two flares. The emitting region cannot be located at r < 0.1 pc for any value of Γ_j, unless the variability timescale is much shorter than observed (i.e. hours). However, this scenario may still be valid if the jet continues to accelerate at or beyond ∼ pc scale such that Γ_j at the location of flare 2 can be significantly smaller than those corresponding to flares 1 and 3.

This work made use of data supplied by the UK Swift Science Data Center at the University of Leicester, and of SMARTS optical/near-infrared light curves which are available at www.astro.yale.edu/smarts/glast/home.php. R.C. and A.D.M. were partially supported by NASA through ADAP award NNX12AE38G and EPSCoR award NNX11AM18A. K.N. acknowledges support by the NSF grant AST-0907872, the NASA ATP grant NNX09AG02G, the NASA Fermi GI program, and the Polish NCN grant DEC-2011/01/B/ST9/04845.

We assume that the emitting region has characteristic size R, which is related to the observed variability timescale via $R \simeq c \mathcal {D}t_{\rm var,obs}/(1+z)$. We can also relate R to the location of the emitting region via R ≃ θr, where θ is the opening angle of the emitting region. We distinguish this parameter from the jet opening angle θ_j, demanding that θ ⩽ θ_j. It is convenient to combine θ with the Lorentz factor, and use the collimation factor Γθ. We can now derive direct constraints on the source Lorentz factor as a function of Γθ:

$$\begin{equation} \Gamma (r,\Gamma \theta) \simeq \left(\frac{\mathcal {D}}{\Gamma }\right)^{-1/2}\left[\frac{(1+z)\Gamma \theta r}{ct_{\rm var,obs}}\right]^{1/2} \,. \end{equation} \tag{ A1 }$$

We calculate the luminosity of the SSC component using the relations $L_{\rm SSC}/L_{\rm syn} \simeq g_{\rm SSC}(u_{\rm syn}^{\prime }/u_{B}^{\prime })$, where g_SSC is a correction factor (mainly due to spectral shape and source geometry) and $u_{\rm syn}^{\prime } \simeq L_{\rm syn}^{\prime }/(4\pi c R^2)$, where $L_{\rm syn}^{\prime } = L_{\rm syn}/\mathcal {D}^4$. The Compton dominance parameter $q = L_\gamma /L_{\rm syn} = g_{\rm EC}(\mathcal {D}/\Gamma _{\rm j})^2(u_{\rm ext}^{\prime }/u_{B}^{\prime })$, where g_EC is a correction factor (mainly due to Klein–Nishina effects). Energy density of external radiation in the jet comoving frame is $u_{\rm ext}^{\prime } \simeq \zeta \Gamma _{\rm j}^2L_{\rm d}/(3\pi cr^2)$, where ζ is related to the covering factors of the sources of external radiation (Sikora et al. 2009), and L_d is the accretion disk luminosity. Putting these relations together, we obtain

$$\begin{eqnarray} \Gamma _{\rm j}(r,L_{\rm SSC}) &\simeq& \left(\frac{\mathcal {D}}{\Gamma _{\rm j}}\right)^{-1}\left[3\left(\frac{g_{\rm SSC}}{g_{\rm EC}}\right)\left(\frac{L_{\rm syn}}{L_{\rm SSC}}\right)\left(\frac{L_\gamma }{\zeta L_{\rm d}}\right)\right]^{1/8}\nonumber\\ &&\times\left[\frac{(1+z)r}{2ct_{\rm var,obs}}\right]^{1/4}\,. \end{eqnarray} \tag{ B1 }$$