CHARACTERIZING THE OPTICAL VARIABILITY OF BRIGHT BLAZARS: VARIABILITY-BASED SELECTION OF FERMI ACTIVE GALACTIC NUCLEI

Blazars are a relatively rare sub-class of active galactic nuclei (AGNs) in which a jet is aligned along the observer's line of sight, leading to the effects of relativistic beaming and unusual associated emission (Blandford & Rees 1978). Blazars are among the most variable extragalactic objects detected in time-domain optical surveys and have strong emission from radio to TeV energies (Ulrich et al. 1997). The central engine is believed to be accretion onto a supermassive black hole, driving relativistic outflows in a collimated jet with typical Lorentz factors on the order of Γ ∼ 10 that is pointed to within angle Γ⁻¹ of the observer. The term "blazars" usually encompasses both BL Lac objects and flat-spectrum radio quasars (FSRQs), which are believed to be jet-aligned Fanaroff–Riley type I and II AGNs, respectively (Urry & Padovani 1995). In this paper, we adopt this definition.

The canonical broadband spectral energy distribution (SED) of blazars typically includes several main components: (1) a synchrotron peak likely due to tangled magnetic fields in the jet that may extend from the radio to the soft-X-ray regime, (2) an inverse-Compton peak in the X-ray to GeV regime likely due to scattering of synchrotron or external photons off of relativistic electrons in the jet, (3) a possible inverse-Compton component in the soft X-rays due to a hot corona, and (4) occasional hints of an underlying accretion disk continuum or host galaxy emission. The optical and γ-ray observations we use in this paper are expected to be dominated by the synchrotron and jet inverse-Compton emission, respectively.

The strong high-energy inverse-Compton peak in the SEDs of blazars causes them to account for the vast majority of bright extragalactic γ-ray-emitting sources. The Fermi Space Telescope has surveyed the γ-ray sky since launch in 2008, and its Large Area Telescope (LAT; Atwood et al. 2009) instrument provides by far the deepest survey to date in the 100 MeV to 100 GeV regime. The current 2 year LAT source catalog (2FGL; Nolan et al. 2012) includes 1873 total sources, ∼44% of which are reliably associated with AGNs, and an additional ∼14% are candidate AGN associations. Of the reliably associated AGNs, the overwhelming majority are blazars. Approximately half of these Fermi blazars are BL Lac objects, and half are FSRQs (Ackermann et al. 2011).

The 2FGL catalog is a significant improvement over the 11 month Fermi LAT source catalog (Abdo et al. 2010) in number of sources, source detection methods, and source associations to known objects. However, ∼32% of sources in the 2FGL catalog still lack reliable associations with any known object, many of which may be blazars, especially at high Galactic latitudes (Ackermann et al. 2012). Although the angular resolution of the LAT is a dramatic improvement over previous all-sky γ-ray survey instruments such as EGRET, typical 95% confidence error ellipses of sources in the 2FGL catalog are still on the order of several arcminutes in size and can contain numerous candidate counterparts at other wavelengths, making association and identification of γ-ray sources difficult.

Recently, huge efforts have gone into identifying Fermi γ-ray sources by positional coincidences of candidate counterparts from multi-wavelength surveys (e.g., Stephen et al. 2010; Maeda et al. 2011), observations of contemporaneous variability at different wavelengths (e.g., Kara et al. 2012), and statistical methods based on observed γ-ray source properties (e.g., Ackermann et al. 2012). Notably, Massaro et al. (2011) applied a variant of the well-known method of selecting AGNs by their unique colors in the mid-IR (e.g., Lacy et al. 2004; Stern et al. 2005) to data from the WISE survey (Wright et al. 2010) and used their selection method to separate blazars from stars. Massaro et al. (2012) then further utilized their technique to find WISE-detected γ-ray blazar candidates in the 2FGL, achieving excellent completeness, although the efficiency of their method is unclear (see discussion in Section 5.1).

Aside from their distinct mid-IR colors, blazars are also unique in their strong variability at nearly all wavelength regimes. This will make blazars stand out in the flood of time-domain data from current and future large-scale optical imaging surveys such as Pan-STARRS (Kaiser et al. 2010), the Catalina Real-Time Transient Survey (Drake et al. 2009), the Palomar Transient Factory (PTF; Rau et al. 2009), and the Large Synoptic Survey Telescope (LSST; Ivezić et al. 2008). It is thus an auspicious time to explore the possibility of selecting and identifying γ-ray-emitting blazars en masse in Fermi error ellipses by their optical photometric variability.

Blazars have observed variability on timescales from hours to years, and their distinctiveness in time-domain surveys results from two effects. First, the relativistic jet strongly beams non-thermal emission along the line of sight toward the observer, boosting the luminosity by several orders of magnitude. Second, relativistic Doppler-boosting shortens the observed timescale of variability in comparison to that in the rest frame of the outflow. However, the physical mechanisms causing the variability are far less clear. Since blazar emission is largely dominated by non-thermal emission from the jet, it is possible that internal shocks from overtaking collisions of fluid shells with different velocities along the jet produce time-variable synchrotron emission (e.g., Böttcher & Dermer 2010), although other emission components such as the accretion disk may also be time variable.

Other classes of AGNs, including the far more numerous "normal" quasars (i.e., non-FSRQ, type 1 AGNs), are also highly optically variable, although the predominance of particular physical variability mechanisms may differ among AGN subclasses. Sesar et al. (2007) have shown that ∼90% of quasars in Sloan Digital Sky Survey (SDSS; York et al. 2000) Stripe 82 are variable at the >0.03 intrinsic rms variability level (defined in Equation (1)), and previous studies using large samples of individual quasar light curves have shown that quasar optical variability is well described by a damped random walk (DRW) model (Kelly et al. 2009; MacLeod et al. 2010; Kozłowski et al. 2010; Butler & Bloom 2011; Zu et al. 2012). MacLeod et al. (2011, hereafter MA11) further utilized the DRW model as a tool to separate stars and quasars by optical photometric variability in SDSS Stripe 82 and showed that quasar selection using this method can achieve efficiency E ⩾ 80% and completeness C = 90% by variability alone. This is particularly useful for selection of quasars in certain redshift regimes, where color selection may fail due to overwhelming contamination from foreground stars. Furthermore, studies of the long-term variability characteristics of AGNs may also provide a unique perspective on accretion disk and jet physics.

The optical variability of small numbers of individual blazars has been studied extensively in the literature through intensive long-term monitoring campaigns (e.g., Webb et al. 1988; Carini et al. 1992). However, no consistent picture of variability has emerged, at least partially due to the heterogenous nature of these observations and the small sample sizes. Blazar light curves from a large-scale flux-limited time-domain survey should be more homogenous, and the depth of current surveys can yield orders of magnitude more blazar light curves than possible through targeted monitoring campaigns. However, the rarity of blazars places significant constraints on the specifications of any time-domain survey in which a large number of blazar light curves can be obtained. Despite the success of the SDSS Stripe 82 for time-domain astronomy, the depth of the survey (∼22.5 mag in r band) does not compensate for the small sky coverage (∼290 deg²) when it comes to studying blazars. We instead employ the recalibrated LINEAR survey, which essentially covers the ∼10,000 deg² SDSS photometric footprint down to ∼17 mag in r band (Sesar et al. 2011, hereafter SE11). Bauer et al. (2009a) provided a previous study of blazar optical variability in the Palomar-Quest Survey, and Bauer et al. (2009b) conducted a subsequent variability-based search for new blazars within that survey. However, the Palomar-Quest survey provided a median of ∼5 epochs of observation for each object per filter and could only be used to study the ensemble variability of blazars rather than the individual sources. By contrast, the ∼200 epochs of observation per object in the LINEAR survey allow us to directly model the light curves of individual blazars in this study.

The structure of this paper is as follows: In Section 2, we describe the LINEAR survey and the construction of the Variable LINEAR catalog, from which our blazar and comparison normal quasar samples are drawn. In Section 3, we describe the DRW model of variability and our results of modeling the LINEAR light curves of blazars, normal quasars, and stars. In Section 4, we test the efficiency and completeness of blazar selection using a DRW model within Fermi error ellipses. In Section 5, we discuss this method in the context of its future applications, as well as possible implications of our work for AGN physics. We summarize and conclude in Section 6.

2.1. Photometric Data

2.2. The Bright LINEAR Catalog

From the full LINEAR catalog, we select a bright subset catalog suitable for variability science by imposing a cut on the minimum number of good observations in each light curve of ⩾30 and a LINEAR recalibrated magnitude cut of 14 < m_LINEAR < 17 (where photometric errors are ≲0.11 mag). At magnitudes >17 and <14, single-epoch photometric errors rise rapidly due to photon noise and saturation, respectively. These cuts provide us with ∼4.5 million objects in the Bright LINEAR catalog.

To create a sample of known blazars and quasars in the Bright LINEAR catalog, we positionally match the LINEAR sources to catalogs of known BL Lac objects and FSRQs. To find the known BL Lac objects, we match to 1371 BL Lac objects in the Véron-Cetty & Véron (2010) catalog, 501 BL Lac objects in the Plotkin et al. (2008) catalog, and 637 radio-loud BL Lac objects in the Plotkin et al. (2010) catalog (none of the 86 radio-quiet objects in Plotkin et al. (2010) matched to a LINEAR object), all using a 30 matching radius, resulting in 101 matched distinct BL Lac objects. Known quasars in the Bright LINEAR catalog were identified by matching to the 105,782 quasars in the Schneider et al. (2010) SDSS Data Release 7 catalog, also with a 30 matching radius, resulting in 1020 matched distinct quasars. We also separate out the FSRQs by positionally matching these 1020 quasars to the CRATES radio survey (Healey et al. 2007) of 11,131 flat-spectrum (α > −0.5, where S_ν∝ν^α) radio sources using a 30 matching radius, resulting in 42 distinct FSRQs. The 978 non-flat radio spectrum quasars remaining are classified here as normal quasars (i.e., unlikely to be strongly dominated in the optical by jet emission).

We note that there is ultimately an overlap of three objects between the BL Lac object and FSRQ samples due to double counting objects that appear in the BL Lac object catalogs, the Schneider et al. (2010) quasar catalog, and the Healey et al. (2007) catalog, which may result from the intermittent appearance of broad emission lines in their optical spectra. Since our definition of "blazars" here includes both BL Lac objects and FSRQs, we count these three objects as blazars (but not normal quasars due to their flat radio spectra). In summary, we identified a total of 140 blazars and 978 normal quasars in the Bright LINEAR catalog.

The potential of the LINEAR survey for blazar time-domain studies can be seen in Figure 1, which shows normalized histograms of the number of LINEAR epochs of observation for each of the 140 objects in the blazar sample, the 987 objects in the normal quasar sample, and all other objects (mainly stars) in the Bright LINEAR catalog. The histograms are all generally consistent with each other, peak at ∼200 epochs, and lack objects with <30 good observations due to the previously imposed cut.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before attempting to model the light curve of every object, we would first like to compare the general level of variability of blazars, normal quasars, and other objects in the Bright LINEAR catalog. We follow Sesar et al. (2007) and estimate the intrinsic rms variability σ for each light curve, defined as

$$\begin{equation} \sigma = [\Sigma ^2 - \psi (m_{{\rm LINEAR}})^2]^{1/2} \end{equation} \tag{ 1 }$$

for all objects with Σ > ψ(m_LINEAR) and σ = 0 otherwise, where Σ is the rms scatter of m_LINEAR for all observations of each object, and ψ(m_LINEAR) is the median photometric error of LINEAR objects as a function of magnitude found by SE11 (see their Equation (7)). Since the photometric errors in LINEAR are relatively large, Σ mainly reflects photometric noise rather than variability for the vast majority of objects. Equation (1) attempts to remove the effects of the median photometric error, and thus σ is a better measure of intrinsic variability. Figure 2 shows the median σ and ψ(m_LINEAR) of all objects in the Bright LINEAR catalog as a function of m_LINEAR magnitude.

Zoom In Zoom Out Reset image size

We compare the variability of blazars to normal quasars and other LINEAR objects, as it is expected that blazars should be systematically more variable. Figure 3 shows the integrated distribution of blazars, normal quasars, and other LINEAR bright objects in the Bright LINEAR catalog as a function of the estimated intrinsic rms variability σ. We estimate that ∼36% of blazars, as compared to only ∼11% of quasars, are variable with σ > 0.15 mag in this passband. At this level of variability, the estimate of σ is not significantly affected by uncertainties in ψ(m_LINEAR). Only ∼4% of the other bright LINEAR objects (non-blazars and non-quasars) are variable above this σ > 0.15 mag level, and they are likely a combination of variable stars (e.g., RR Lyrae), eclipsing binaries, underestimates of the photometric errors ψ, source blending, or other photometric issues.

2.3. The Variable LINEAR Catalog

MA11 designed and tested optical variability-based selection of quasars in SDSS Stripe 82 by modeling the individual light curves of known quasars as a DRW. We perform a similar analysis of the LINEAR light curves of known blazars and normal quasars and show that blazar light curves lie in a distinct region of DRW variability parameter space. The DRW model statistically parameterizes variability using three parameters: a mean light-curve magnitude, a characteristic damping timescale of variability τ, and a driving amplitude of short-term stochastic variability $\hat{\sigma }$ (see Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010; Zu et al. 2012, for further discussion). The fitted parameters can be expressed in terms of the structure function (SF), defined as the rms magnitude difference between all pairs of observations in each individual light curve as a function of the time lag. While this model is phenomenological, it is likely that the fitted parameters reflect the physical processes that cause the variability.

Following Kozłowski et al. (2010) and MA11, we model each individual light curve in our sample as a stochastic process described by the exponential covariance matrix

$$\begin{equation} S_{ij} = \sigma ^{2}\mathrm{exp}(-|t_{i}-t_{j}|/\tau) \end{equation} \tag{ 2 }$$

for each pair of observations at time t_i and t_j in the light curve. This describes a DRW process with a characteristic damping timescale τ beyond which the structure function will asymptote to a constant value of SF$_\infty = \hat{\sigma }\sqrt{\tau }$ and a long-term standard deviation of variability $\sigma = \sqrt{2}$SF_∞. The short-term driving amplitude of variability $\hat{\sigma } = \sigma \sqrt{2/\tau }$ determines the rise of SF($\Delta \mathit {t}$) for Δt ≪ τ. We model each individual light curve and calculate $\hat{\sigma }$ and τ, along with their likelihood distributions, using the method of Press et al. (1992), its generalization in Rybicki & Press (1992), and the fast computational implementation in Rybicki & Press (1995). The corresponding structure function for the model for time lag Δt is

$$\begin{equation} {\rm SF} (\Delta \mathit {t})= \rm {SF}_{\infty } (1-{\rm e}^{-|\Delta \mathit {t}|/\tau })^{1/2}. \end{equation} \tag{ 3 }$$

In fitting each individual light curve with the DRW model, we calculate the likelihood of a DRW solution L_DRW, as well as the likelihood of a pure white noise solution L_noise. Light curves that are better described by the DRW model than white noise at a 5σ level will have ΔL_noise = ln(L_DRW/L_noise) > 5. All our FSRQ and BL Lac object light curves in the Variable LINEAR catalog fit the DRW model at above this 5σ level. To remove light curves for which the survey length is shorter than τ (thus leaving τ unconstrained), we also calculate the likelihood of a runaway timescale, L_∞. Objects for which τ is constrained in the DRW model will have ΔL_∞ = ln(L_DRW/L_∞) > 0.05 (MacLeod et al. 2010). MacLeod et al. (2010) showed that the DRW is an excellent fit to quasar light curves and used cuts on minimum ΔL_noise and minimum ΔL_∞ to select quasars with high-quality light curves. In the higher stellar density regions of Stripe 82, MC11 used cuts on maximum $\hat{\sigma }$, minimum τ, and minimum ΔL_noise to achieve impressive completeness C = 93% (defined as percentage of the total confirmed quasars in the sample that are selected) and efficiency E > 78% (percentage of quasar candidates selected that are confirmed quasars) in quasar selection. Furthermore, Kozłowski et al. (2010) showed that the method works even in the high stellar density regions of the Magellanic Cloud fields.

We calculate the best-fit DRW variability parameters of the individual light curves of the 60 blazars, 155 normal quasars, and a random sample of 6000 other objects (mainly foreground Galactic stars, representative of the other >180,000 objects) from the Variable LINEAR catalog. To first understand the underlying distribution of DRW variability parameters for these different populations, we redshift-correct our parameters as outlined in Kelly et al. (2009) using spectroscopic redshifts from the Schneider et al. (2010) catalog for all FSRQs and normal quasars. For BL Lac objects, redshifts are much less accurate (and sometimes not possible) due to the lack of strong spectral features. Since our sample of BL Lac objects is drawn from a variety of catalogs, there are sometimes discordant redshifts reported; in such cases, we preferentially adopt the spectroscopically derived redshifts (or the lower limits) from the Plotkin et al. (2008) and Plotkin et al. (2010) catalogs. For three BL Lac objects, there is no redshift reported in any of the catalogs, and so for these three objects we adopt the mean redshift found for all other LINEAR BL Lac objects, z = 0.28. While this assumption is not ideal, very few extragalactic objects in the LINEAR survey will be at high redshifts due to the shallow optical flux limit, so uncertainties on these three redshifts will not strongly impact our results.

Imposing the restrictions on ΔL_∞ > 0.05 and ΔL_noise > 5 leaves us with 51 blazars and 121 normal quasars in our sample. Figure 4 shows the best-fit DRW parameters τ against SF_∞ of blazars, normal quasars, and other objects (mostly Galactic stars), along with normalized histograms of their distributions. Although we use separate symbols for BL Lac objects and FSRQs in Figure 4, we do not separate these two blazar sub-populations in our subsequent analysis due to the small sample size and unknown biases in our blazar sample stemming from possible correlations between variability and physical properties of the blazars. Future surveys yielding larger samples may be able to robustly constrain differences between these two sub-populations and provide further insight into their jet properties. Figure 4 is directly comparable to the analysis of Kozłowski et al. (2010) and MA11 and confirms that the structure functions of normal quasar light curves tend to have longer characteristic timescales of variability τ and slightly larger SF_∞ than stars. More importantly, our analysis also shows that blazars have characteristic timescales τ in between those of normal quasars and other objects (mainly variable stars), as well as larger values of SF_∞. The driving amplitude on short timescales $\hat{\sigma }$ is related to these two parameters by the relation $\hat{\sigma }$ = SF_∞/$\sqrt{\tau }$, as shown in Figure 4. Figure 5 shows the distribution for blazars and normal quasars in $\hat{\sigma }$, where normal quasars tend to peak at log ($\hat{\sigma }$/mag yr^−1/2) ∼−0.75, while blazars tend to peak at log $\hat{\sigma }$ ∼0.0. A two-dimensional Kolmogorov–Smirnoff test between the distributions of $\hat{\sigma }$ for blazars and normal quasars in Figure 5 gives a p-value of 2.04 × 10⁻⁵. For the distributions of τ in Figure 4, the p-value is 0.032.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While MA11 performed quasar selection using this approach by placing cuts on maximum $\hat{\sigma }$ and minimum τ, we can also see from Figure 4 that a minimum τ cut would separate blazars from other objects (mostly Galactic stars), while a minimum $\hat{\sigma }$ cut can help separate blazars from quasars. However, it is clear from Figure 4 that highly efficient and complete quasar and blazar selection purely by optical variability in the LINEAR survey is difficult, as there will be either heavy contamination or a low recovery rate. The main cause of this is the bright magnitude range (14 < m_LINEAR < 17) of LINEAR, which overwhelmingly probes Galactic stars and very few actual AGNs (e.g., there are >180,000 sources in the Variable LINEAR catalog, only 214 of which are known normal quasars and blazars), and so deeper survey data containing a substantial fraction of extragalactic sources are necessary. However, currently available data from deeper time-domain optical surveys such as SDSS Stripe 82 have comparatively small sky coverage and do not yield a large enough sample of the much rarer blazars.

Despite the overwhelming stellar contamination in the bright magnitude range of the LINEAR survey, our proposed blazar selection method is still useful with LINEAR if we restrict our search to the small regions of sky associated with the positions of Fermi sources. This greatly reduces the stellar contamination, while simultaneously providing a test of a novel method (based on long-term photometric variability characteristics) to associate γ-ray sources with counterparts at optical wavelengths. We will show the DRW parameters of known Fermi AGNs (overwhelmingly blazars), quantify the completeness and efficiency of blazar selection using this method, and identify new variability-selected candidate Fermi blazar counterparts.

We first positionally match the Variable LINEAR catalog described in Section 2.3 to the 1873 sources in the 2FGL catalog, approximating each error ellipse as a circle with radius equal to the semimajor axis of the 95% confidence ellipse. This initial match yields 480 variable LINEAR objects in 173 Fermi error circles. These 480 variable LINEAR optical objects in Fermi error circles are then positionally matched using a 30 matching radius to the 929 associated AGNs in the 2nd Fermi AGN catalog (Ackermann et al. 2011), restricting consideration to those that have Bayesian association probabilities >0.8 (this association method is described in the Appendix of Abdo et al. 2010). This yields a sample of 47 known (i.e., confidently associated) Fermi γ-ray-emitting AGNs. We calculate DRW parameters for the LINEAR light curves of all 480 variable objects and make nominals cut on ΔL_noise > 5 and ΔL_∞ > 0.05 as in Section 3.

Due to the bright magnitude range of the LINEAR survey, it is likely that the flux limit of the LINEAR survey is too bright to actually detect the optical counterparts of the faintest Fermi γ-ray sources, leaving us a population of orphan error circles that contain only LINEAR contaminants. Since these fainter sources will also preferentially have larger error circles, we can remove them by considering further cuts on the maximum radii of the error circles. In Figure 6, we show the distribution of the total number of LINEAR variable objects in Fermi error circles as a function of error circle radius r, for all 173 error circles that contain at least one LINEAR variable object, as well as the distribution for the 47 error circles containing a LINEAR variable object that matched to an AGN in the 2nd Fermi AGN catalog. Figure 6 shows that the number of variable LINEAR objects in error circles with known γ-ray AGN counterparts peaks at 2' and drops rapidly; beyond 10', there is only one Fermi error circle that contains a known γ-ray AGN counterpart. However, the distribution of the total number of variable LINEAR objects in all error circles as a function of radius has a long tail to large r, as there are a small number of large error circles (r > 15') with large numbers of variable LINEAR objects, as expected. This suggests that a cut on the error circle radius of ∼10' is appropriate, but to gauge the impact of such a cut, we will consider in detail a range of 8', 10', and 12' cuts on the maximum error circle size. A cut at r < 10' removes from consideration only one Fermi error circle containing a known γ-ray AGN, while a cut at r < 8' removes four. Scatter in the negative relation between the optical flux and γ-ray error circle size of Fermi blazars may introduce some biases in our resulting sample toward blazars with higher than average optical-to-γ-ray flux ratios when we place a cut on the error circle size. However, this bias is expected to be small, as only ∼7% of reliably associated AGNs in the 2nd Fermi AGN catalog with optical magnitudes <17 reside in Fermi error circles >10' in radius. We note that although the Fermi AGN association success rate of Ackermann et al. (2011) as a function of r will also affect Figure 6, this effect is secondary to the LINEAR flux limit for the high association probabilities considered here.

Zoom In Zoom Out Reset image size

Figure 7 (left panel) shows the distribution in τ and SF_∞ of LINEAR optical counterparts to Fermi γ-ray AGNs, as well as all other variable objects in Fermi 2FGL error circles of r < 10'. We do not apply cosmological redshift corrections here for blazar selection. Similar to the conclusions drawn from Figure 4, there is clear separation between γ-ray-emitting AGNs and other LINEAR variables in DRW parameter space. The significant normal quasar population seen in Figure 4 is minimized when looking only in Fermi error circles. As expected, almost all γ-ray AGNs recovered in our sample have blazar-like variability. Figure 7 (right panel) shows that variable objects in error circles with radius r > 10' largely have SF_∞ and τ similar to that of stars and are thus mostly contaminants, rather than new blazar candidates that were removed by the cut on r < 10'.

Zoom In Zoom Out Reset image size

We can select γ-ray AGNs by placing cuts on minimum τ and $\hat{\sigma }$. For different log τ and log $\hat{\sigma }$, we estimate efficiency and completeness of γ-ray AGN selection. The results for a 10' cut on the error circles are shown in Figure 8. Selection cuts at large minimum τ and minimum $\hat{\sigma }$ yield efficiencies of 0 as there are no such objects. Figure 9 shows the maximum completeness achievable by tweaking the selection cuts on minimum τ and $\hat{\sigma }$, as a function of efficiency. This shows how the completeness changes as a function of efficiency in γ-ray AGN selection, for error circles with a maximum radius of 8', 10', and 12'. As expected, the efficiency decreases as the maximum size of error circles considered increases, at a set completeness.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By jointly maximizing efficiency and completeness via the quantity $\sqrt{E^2+C^2}$ from the curves in Figure 9, we can achieve E ⩾ 88% and C = 88% for Fermi error circles of r < 8' using cuts on log τ > 1.35 and log $\hat{\sigma } > -1.20$. For error circles of r < 10', E ⩾ 76% and C = 86% using cuts on log τ > 1.53 and log $\hat{\sigma } > -0.90$, and for error circles of r < 12', E ⩾ 70% and C = 86% using cuts on log τ > 1.53 and log $\hat{\sigma } > -0.90$. This demonstration verifies that we are able to select (in this case, recover in double-blind fashion) γ-ray-emitting AGNs in Fermi error circles with high completeness and efficiency by modeling their optical light curves as a DRW and imposing selection cuts on the variability parameters.

Our definition for completeness is similar to the convention of many other recent time-domain AGN studies (e.g., MA11). The completeness we calculate here is specifically for the selection of variable LINEAR objects matched to Fermi AGNs in Ackermann et al. (2011) with Bayesian association probabilities >0.8, and the efficiency is calculated assuming that these associations are all correct. However, this efficiency is a lower bound, as it will increase if some of the contaminants we encounter in selection are actually blazars. This may occur in our analysis if some LINEAR variables that pass all DRW blazar selection cuts and did not match to an AGN in the 2nd Fermi AGN catalog may actually be additional γ-ray AGN counterparts, not already recognized as such. Indeed, this is likely the case for at least one object, lying in the error circle of Fermi source 2FGL J1649.6+5238.

In the 2nd Fermi AGN catalog, Fermi source 2FGL J1649.6+5238 is listed as associated with the radio source 87GB 164812.2+524023 at a 0.0 Bayesian probability, calculated based on the local density of sources from catalogs of likely counterparts. This radio source is positionally coincident with a LINEAR object at R.A. = 16^h49^m24^s.99, decl. = 52°35'1505, which has DRW parameters calculated from its LINEAR light curve of log τ = 1.936, log $\hat{\sigma } = -0.026$, and log SF_∞ = 0.942, very similar to that of confirmed γ-ray AGNs in Figure 7. This LINEAR source was not counted as a γ-ray AGN in our calculations of efficiency and completeness because its Bayesian probability of association was below 0.8 in Ackermann et al. (2011). Furthermore, this LINEAR counterpart also did not match to any known blazar in the BL Lac object and quasar catalogs used in Section 2.2.

Based on these blazar-like DRW parameters calculated from the light curve, we suggest this LINEAR variable object as the plausible optical counterpart to Fermi source 2FGL J1649.6+523, as well as increased confidence in associating 87GB 164812.2+524023 as the radio counterpart. We note that this Fermi source is also associated with this radio source with a 0.82 probability using the log N–log S method in Ackermann et al. (2011), based on observed properties of candidate radio counterparts. The origin of this disparity lies in the differences in the approaches used to calculate the Bayesian and log N–log S (radio) association probabilities. We opted to use the Bayesian probability >0.8 criterion for known γ-ray AGNs, but the use of the log N–log S probability also recovers this radio source as a likely γ-ray AGN. It is not our intention to compare the two statistical approaches, but rather to point out that the discrepancies in the results of these two methods may lead to different results in analysis. In cases where these two methods lead to highly discrepant results, our completely independent optical variability-based method can provide valuable additional information.

Finally, we note that our variability-based approach may be further applicable to mass identification of unassociated γ-ray sources using deeper time-domain optical surveys. Data from current surveys such as the PTF, with single-epoch depth of r ∼ 20.6 mag, 30 epochs of observation over 8200 deg² of sky, and ∼0.01 mag repeatability, can be used to identify unassociated Fermi blazars using the method we have presented here. Future surveys such as LSST, with single-epoch depth of r ∼ 24.5 mag over 20,000 deg² of sky and ∼10³ epochs per source, can be used to vastly increase the existing sample size.

5.1. Efficiency and Completeness

In Section 3, we noted that blazar selection using the Variable LINEAR catalog across the full field of the survey (as contrasted to blazar selection only within Fermi error circles) was not efficient and complete. At first, this may seem to be in stark contrast to the success of quasar variability selection using a similar DRW modeling approach (e.g., MA11; Butler & Bloom 2011; Kozłowski et al. 2010). This is primarily a consequence of the shallow depth of the LINEAR survey rather than a reflection of the ultimate performance expected from our approach (especially for fainter blazars). At the depth of SDSS Stripe 82 (flux limit of g ∼ 20.5), objects with low-redshift quasar-like colors account for ∼63% of all variable objects, with Galactic stars making up the vast majority of the remaining variable objects (Sesar et al. 2007). In contrast, at the shallow depth of the LINEAR survey (flux limit of r ∼ 17), only ∼0.08% of objects in the Variable LINEAR catalog are quasars (see Section 2.3). As such, quasar selection by imposing appropriate maximum $\hat{\sigma }$ and minimum τ cuts (as done in MA11) on the LINEAR quasar sample in Figure 4 will not yield E beyond a few percent for any value of C > 50% due to scatter in the parameters of the enormously larger population of variable stars. Blazar selection in LINEAR will be similarly inefficient. These issues are exacerbated by the large errors on the estimated DRW parameters as compared to SDSS, caused by the larger photometric uncertainties in LINEAR. In any case, current and future time-domain surveys such as PTF and LSST will alleviate both these issues by going many magnitudes deeper while providing ∼1% photometry.

As noted in Section 4, our calculated efficiencies of blazar selection are lower limits, since calculating the true efficiency requires correct identification of every object as either a blazar or contaminant. There are many contaminating variable LINEAR objects that fall in the Fermi error circles we considered in our analysis in Section 4 that have blazar-like variability (i.e., have DRW parameters similar to blazars and are counted as contaminants when jointly optimizing completeness and efficiency for known γ-ray blazar selection). Some of these contaminants are certainly blazars and may actually be the γ-ray source. Securing identification of these blazar-like variable objects will require spectroscopic follow-up. We list in Table 1 the positions, parent 2FGL error circle, and DRW parameters of 12 such new candidate LINEAR variable objects with blazar-like optical variability that fall in Fermi 2FGL error circles with r < 10'. We also list any previously associated counterparts of each parent 2FGL source from Ackermann et al. (2011) with Bayesian association probabilities >0.8. The LINEAR blazar candidates are not positionally coincident with these previously associated counterparts; rather, they are additional γ-ray blazar candidates in their respective error circles. Since it is possible that some Fermi error circles may contain more than one γ-ray-emitting object, these LINEAR blazar candidates are still worthy of spectroscopic follow-up.

Table 1. J2000 Positions of New Candidate LINEAR Blazars in Fermi Error Circles with r < 10' Selected by Variability

R.A.	Decl.	Fermi 2FGL Source	log	log	2nd Fermi AGN Catalog Association
(deg)	(deg)		(τ/days)	($\hat{\sigma }$/(mag yr−1/2))
128.5027	42.2570	2FGL J0834.3+4221	1.89	−0.20	OJ 45
153.3687	34.4883	2FGL J1013.6+3434	1.58	−0.28	Unassociated
172.2141	−5.6568	2FGL J1129.0−0532	2.20	0.02	Unassociated
189.5802	−20.0237	2FGL J1238.1−1953	2.07	−0.11	PMN J1238−1959
219.1272	23.2824	2FGL J1436.9+2319	1.66	−0.29	PKS B1434+235
220.4843	43.9407	2FGL J1442.0+4352	2.00	−0.26	BZB J1442+4348
232.8992	57.4483	2FGL J1531.0+5725	2.50	−0.04	Unassociated
232.9009	57.3541	2FGL J1531.0+5725	2.36	−0.58	Unassociated
246.9883	32.4363	2FGL J1627.8+3219	2.30	−0.34	Unassociated
250.2259	11.7870	2FGL J1641.0+1141	1.89	−0.20	MG1 J164058+1144
252.3542	52.5875	2FGL J1649.6+5238	1.93	−0.02	Unassociated
331.7412	−0.4681	2FGL J2206.6−0029	1.80	0.27	PMN J2206−0031

Download table as:  ASCIITypeset image

We further note that the efficiencies we have calculated in γ-ray blazar selection in Section 4 are conservative, as we have considered all variable LINEAR objects lying in all Fermi error circles below a certain size. However, ∼10% of sources in the 2FGL are associated with non-AGN objects (pulsars, supernova remnants, etc.), while ∼30% of sources are unassociated. We chose to include all these error circles in our analysis (and not just those containing known Fermi blazars) in a double-blind test, as it is a more faithful demonstration of the efficacy of this approach.

We would like to compare the efficiency and completeness of our variability-based Fermi blazar selection method to other methods. However, a direct comparison is difficult; few studies have attempted to use characteristic properties of blazars to systematically select Fermi blazars in double-blind fashion as we have done. Furthermore, the completeness and efficiency of those methods are almost never jointly optimized (as done in Section 4). For example, Massaro et al. (2011) used mid-IR color selection to recover known blazars in the WISE survey and extended their approach to find new Fermi blazar candidates (Massaro et al. 2012). By parameterizing the similarity of the mid-IR colors of WISE sources to WISE-detected Fermi blazars, Massaro et al. (2012) imposed a mid-IR color-based selection criterion that is able to select known Fermi blazars with 87% completeness over the full region of sky surveyed by WISE in its first year. However, Massaro et al. (2012) did not jointly calculate the completeness and efficiency of their approach. Their mid-IR color selection method has the advantage of the all-sky coverage of the WISE survey. However, mid-IR color selection is dependent on the degree to which non-thermal jet emission dominates the mid-IR emission of individual blazars. Indeed, Massaro et al. (2011) find that FSRQs generally have WISE colors closer to normal quasars than to BL Lac objects. This may be due to significant thermal mid-IR emission in the SEDs of FSRQs, similar to that in normal quasars (e.g., from dust emission; see Plotkin et al. 2012), thus making blazar selection more difficult.

Massaro et al. (2012) apply their selection method to WISE sources in unidentified Fermi error circles to find 297 γ-ray blazar candidates in 156 out of 313 unidentified Fermi error circles. However, their estimates of the efficiency of their method by systematically offsetting the position of each Fermi error circle to assess the number of random associations gave 262 false blazar candidates (using the same mid-IR color selection criteria), suggesting that the efficiency of this approach may be low. Nevertheless, mid-IR color- and variability-based selections of blazars are independent and complementary methods, with different selection biases. A combination of both approaches, using the full WISE survey as well as data from current time-domain optical imaging surveys, may be highly efficient and complete.

5.2. Implications for Blazar Variability

Our finding that optical blazar variability is well described as a DRW process with distinct variability parameters may have interesting potential implications for blazar jet physics. The SF_∞ and τ DRW parameters should both be affected by relativistic effects. For example, the observed blazar characteristic damping timescale τ_{blz, obs} (after correcting for cosmological redshift) should be shortened in comparison to the interval τ_{blz, rest} in the rest frame of the emitting material due to the relativistic motion, with τ_{blz, rest} = δτ_{blz, obs}, where δ is the kinematic Doppler factor. While τ_{blz, rest} is uncertain and difficult to measure, we can place constraints on it through comparisons to τ_qso, the characteristic DRW timescale of normal quasars (which are not affected by relativistic effects), and independent measurements of δ.

In Figure 4, the distribution of τ peaks at τ_{blz, obs} ∼ 80 days for blazars and τ_qso ∼ 320 days for normal quasars. While these distributions have large tails, their relative peaks provide an adequate basis for comparison. If the underlying variability in quasars ultimately also drives blazar variability, then δ ∼ τ_qso/τ_{blz, obs} ∼ 4. Figure 10 shows the distribution of this estimated δ for blazars in the Variable LINEAR catalog, assuming τ_qso = 320 days. This is smaller than estimates of the typical blazar Doppler factor from radio observations of Fermi-detected blazars (δ ∼ 10–30; Savolainen et al. 2010) and may imply that the variability of the non-thermal jet emission in the jet rest frame (e.g., due to shocks in the jet) occurs on longer timescales than the thermal disk variability of normal quasars.

Zoom In Zoom Out Reset image size

We can instead adopt the observed typical Doppler factors for blazars and calculate τ_{blz, rest} for each individual blazar. Ideally, we would use a measured value of δ for each individual blazar in this calculation, but direct measurements of δ from the literature are sparse and inhomogeneous. In Figure 10 (bottom panel), we show the distribution of τ_{blz, rest} assuming a reasonable δ = 10 for all blazars. The peak in Figure 10 occurs at τ_{blz, rest} ∼ 3 years, which is longer than the τ_qso ∼ 320 days we measure for LINEAR normal quasars. The discrepancy between δ for Fermi blazars implied by radio observations and δ calculated assuming intrinsic variability similar to normal quasars suggests different stochastic mechanisms driving the variability in the disk and the jet. While a detailed physical interpretation of this is currently unclear, our long-term optical variability results may provide additional constraints on models of blazar jets.

5.3. Short-timescale Variability of Quasars and Blazars

We have used light curves from the LINEAR optical imaging survey to study the variability of blazars and to compare their variability to that of normal quasars. In contrast to many earlier studies of blazar variability, this study is the first to use a reasonably large, homogeneous sample of individual light curves from a wide-field survey, where the distribution of various variability properties for blazars and quasars can be calculated and compared.

In a time-domain survey, blazars are among the most violently variable objects detected and are generally more variable than normal quasars. We estimate that ∼36% of blazars but only ∼11% of quasars are variable at the extreme >0.15 mag intrinsic rms variability level. We show that, like quasars, blazar light curves are well fit by the DRW model for variability, but the blazars lie in distinct variability parameter space with higher τ than stars and higher $\hat{\sigma }$ than normal quasars. This suggests that blazars can be selected with high efficiency E and completeness C by imposing minimum selection cuts on τ and $\hat{\sigma }$.

Due to the overwhelming numbers of variable stars in the bright magnitude regime probed by LINEAR, it is difficult to select blazars (and quasars) at high efficiency and completeness in the full LINEAR sample using optical variability methods alone. We instead examine a more focused test to select γ-ray-emitting blazars in Fermi error ellipses from the 2FGL catalog. We calculate DRW parameters for light curves of all LINEAR variable objects within Fermi error ellipses and place minimum selection cuts on τ and $\hat{\sigma }$. Using this approach, we are able to recover the corresponding γ-ray-emitting AGN counterparts in the 2nd Fermi AGN catalog with E ⩾ 88% and C = 88% for 95% confidence error ellipses with semimajor axis r < 8' and E ⩾ 70% and C = 86% for r < 12'.

Our γ-ray blazar selection has uncovered a variable LINEAR object, coincident with radio source 87GB 164812.2+524023, in the error ellipse of Fermi source 2FGL J1649.6+5238. This object was not associated with the Fermi source in the 2nd Fermi AGN catalog using Bayesian probabilities. Our analysis shows that this object has optical variability properties consistent with γ-ray-emitting blazars and is likely to be the γ-ray source. We find a total of 12 objects with LINEAR variability parameters similar to blazars lying in Fermi error circles with r < 10' (see Table 1). Confirmation of these variability-selected blazar candidates will require spectroscopic follow-up.

Our results suggest that the variability of the non-thermal jet emission in blazars is stochastic in nature, with higher amplitudes and shorter observed timescales in comparison to the variability of normal quasars, likely due to the effects of relativistic beaming. Assuming a reasonable Doppler factor of 10, we estimate that blazars are characteristically variable on timescales of ∼3 years in the rest frame of the emitting material, longer than the ∼320 day characteristic disk flux timescale for quasars. This suggests that different physical mechanisms dominate the observed variability from blazars and quasars, likely connected to the jet and to the disk, respectively. The fitted parameters imply that blazars have a typical intraday variability amplitude of ∼0.05 mag, compared to ∼0.01 mag for normal quasars. Furthermore, there is a significant fraction (∼15%) of normal quasars that exhibit large intraday variability of >0.1 mag, detectable above the photometric uncertainty. The source of the extreme variability of these quasars is unclear.

We argue that variability-based blazar selection is likely to be highly efficient and complete in deeper optical time-domain imaging surveys, and that the variability-based blazar selection method presented in this paper is capable of greatly increasing the number of known blazars. Our approach, based on long-term photometric variability characteristics of a large sample of individual blazar light curves, may also bring a new perspective on accretion disk and jet physics. Finally, this work typifies the fruitful ancillary science made possible by combining data from surveys as diverse in all respects as LINEAR, SDSS, and Fermi.

The authors thank Eric Agol (University of Washington) and Ying Zu (Ohio State University) for helpful discussions regarding AGN variability. Ž.I. and C.L.M. acknowledge support by NSF grant AST-0807500 to the University of Washington, NSF grant AST-0551161 to LSST for design and development activity, and Croatian National Science Foundation grant O-1548-2009. R.M.P. acknowledges support from a Netherlands Organization for Scientific Research (NWO) Vidi Fellowship. C.S.K. is supported by NSF grant AST-1009756. B.S. acknowledges NSF grant AST-0908139 awarded to Judy Cohen for partial support.

The LINEAR program is sponsored by the National Aeronautics and Space Administration (NRA No. NNH09ZDA001N, 09-NEOO09-0010) and the United States Air Force under Air Force Contract FA8721-05-C-0002. Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the United States Government.