INFRARED TIME LAGS FOR THE PERIODIC QUASAR PG 1302-102

Optical reverberation mapping, measuring the time delay between variations in the accretion disk optical continuum luminosity of an active galactic nucleus (AGN) and the delayed response of the broad emission line luminosity, has become a standard technique for studying the physics and structure of the inner regions of AGNs (e.g., Blandford & Payne 1982; Peterson 1993; Bentz et al. 2006). Such studies enable measurement of the physical size of the broad-line region, opening a way to estimate the virial black hole masses of distant AGNs (e.g., McLure & Jarvis 2002; Vestergaard & Peterson 2006; Jun et al. 2015) through the relation between broad-line region size and optical luminosity (${R}_{{\rm{BLR}}}\sim {L}^{0.5};$ e.g., Kaspi et al. 2000; Bentz et al. 2006).

Analogous to optical broad-line reverberation analyses, light-travel time lags between the ultraviolet/optical and near-infrared continuum provide a measurement of the physical size of the infrared-emitting dusty tori that are believed to surround most AGNs (e.g., Penston et al. 1974; Clavel et al. 1989; Suganuma et al. 2006). The inferred time-lag-based sizes from this dust reverberation, supported by direct interferometric K-band measurements for a few nearby AGNs (e.g., Swain et al. 2003; Kishimoto et al. 2011a), shows that the dust-emission region scales with the ultraviolet/optical luminosity (${R}_{{\rm{dust}}}-L$ relation; e.g., Oknyanskij & Horne 2001; Koshida et al. 2014) with a similar square root power as the ${R}_{{\rm{BLR}}}-L$ relation, but with a larger normalization. These radius–luminosity relations indicate that the geometric components of an AGN are spatially structured, such that illumination of the accretion disk propagates first through the broad-line region and later through the cooler, dusty region.

Recently, Graham et al. (2015a) reported an unusual example of an optically variable quasar, PG 1302-102 at z = 0.278. The nine-year optical light curve from the Catalina Real-time Transient Survey (CRTS; Drake et al. 2009), combined with historical data going back ∼20 years, shows periodic variability with an observed period of 1884 ± 88 days. Additional examples of similarly periodic quasars have recently been reported by Graham et al. (2015b) and Liu et al. (2015). This unusual phenomenology is distinct from the typical variability of most quasars, which is akin to a damped random walk (e.g., MacLeod et al. 2012; Graham et al. 2014). The most likely interpretations of the periodic quasars involve close, subparsec binary supermassive black hole systems (e.g., D'Orazio et al. 2015a, 2015b; Graham et al. 2015b).

Taking advantage of the distinct periodic signals in the optical, we investigate the infrared light curves of PG 1302-102 using data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and Akari (Murakami et al. 2007) missions. We also include more recent data from the Near-Earth Object WISE Reactivation mission (NEOWISE-R; Mainzer et al. 2014). We measure the time lags between the optical and infrared emission, thereby confirming that the optical periodicity is reproduced in the dust emission and provides a size measurement of the dust-emission region. Furthermore, using multiband 3.4 μm (W1) and 4.6 μm (W2) data, we investigate wavelength-dependent trends in the light curve time delays, which provide further information about the size and structure of the dusty torus. In a separate, related paper, D'Orazio et al. (2015b) report on the UV light curve of PG 1302-102, showing that it is consistent with relativistic boosting from a supermassive black hole binary causing the observed UV and optical periodicity.

In this Letter, we describe the data collection and processing in Section 2, the measurement of time lags and amplitudes in Section 3, and their implication for the binary supermassive black hole interpretation in Section 4. Throughout, photometry is reported in the Vega system, and we adopt a flat ΛCDM cosmology with ${H}_{0}=70\;\;\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{Mpc}}^{-1},$ Ω_m = 0.3, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7.$

PG 1302-102 is bright enough to have been detected in single exposures throughout the WISE mission. Collecting both cryogenic and post-cryogenic multiepoch photometry from the AllWISE and NEOWISE-R data releases, as well as the most recent (pre-release) NEOWISE-R observations of PG 1302-102, we measure W1 and W2 light curves spanning 5.5 years (MJD 55208-57207). The data consist of seven groups spaced at intervals of roughly six months, albeit with a three-year gap when WISE was in hibernation (MJD 55572-56668). Each group has 10–15 observations, with the latest data taken in 2015 July. NEOWISE-R is expected to obtain ∼3 additional groups of observations over the next ∼1.5 years. We do not include longer-wavelength data at 12 μm (W3) or 22 μm (W4) which are restricted to the cryogenic period of the WISE mission as the temporal coverage of those observations consist of only two groups. We further removed the data suffering from poor quality frames (qual_frame = 0), charged particle hits (saa_sep < 0), or scattered moonlight (moon_masked = 1).

Supplementing the WISE observations of PG 1302-102, we also analyzed archival 2.5–5.0 μm grism spectra obtained by Akari six months prior to the WISE observations (Ohyama et al. 2007). Out of three post-cryogenic Akari spectra of PG 1302-102 taken on MJD 55022-55023, we removed the second spectrum as it is significantly affected by hot pixels. We convolved the processed and extracted spectra from the other two observations with the WISE filter responses (Wright et al. 2010) to obtain W1 and W2 equivalent magnitudes. Because the W2 filter response extends longward of 5 μm, the wavelength cutoff of post-cryogenic Akari, we extrapolated the spectrum redward based on a power-law fit, taking into account the 1σ uncertainties in the fitted parameters. Finally, we averaged the two Akari magnitudes, separated by one day, to improve the sensitivity. For ease of description, in the following, we discuss the W1 and W2 light curves, though those light curves also include these earlier epoch Akari observations. Figure 1 presents the CRTS, WISE, and Akari light curves with the periodic functional fits (Section 3) overplotted.

Zoom In Zoom Out Reset image size

Assuming that the optical variability of PG 1302-102 arises from the accretion disk(s) of the (binary) AGNs, we expect the fluctuations in the optical luminosity to process through the inner dusty torus and produce a distribution of time lags in the hot dust emission corresponding to the light-travel time to the inner portions of the torus. We aim to measure the cross-correlation signal between the optical and infrared magnitudes and test if the infrared light curves trace the periodic variability appearing in the optical time series data.

Because the optical light curves in Graham et al. (2015a) show a well-defined periodic behavior, we first cross correlate their sinusoidal model fit with the W1 and W2 band data. In Figure 2(a), we plot the cross-correlation amplitudes for a range of time lags. The cross-correlation amplitudes peak at 0.69 for W1 and 0.63 for W2, indicating that the WISE data are moderately correlated to the periodic optical light curve model.

Zoom In Zoom Out Reset image size

To give proper treatment of the uncertainties into the cross-correlation analysis (e.g., Peterson et al. 1998), we performed a set of Monte Carlo simulations, including the observational uncertainty from the CRTS data (0.06 mag) and the uncertainty in the optical period (88 days). We performed 100,000 realizations adding random Gaussians to the model magnitudes with σ = 0.06 mag and modifying the period by adding random Gaussians with σ = 88 days. We measured the centroid of each simulated cross-correlation function where the amplitudes were above 0.8 times the peak value, while limiting the range of time lags to [−500, 1500] days to do away with the other periodic solutions. In Figure 2(b), we plot the resultant distribution of the time lags from the Monte Carlo simulations. Adopting the fiducial time lag as the mean of the distribution of simulated time lags and the uncertainty as the difference between the −34.13% and 34.13% percentile in the time lag distribution deviated from the mean (Peterson et al. 1998), the W1 and W2 time lags to the optical model (τ, observed-frame days) become

$$\begin{eqnarray}{\tau }_{W1} & = & (335\pm 153)+1884n\\ {\tau }_{W2} & = & (524\pm 148)+1884n\end{eqnarray} \tag{ 1 }$$

with any integer (n = 0, ±1, ±2, ...) multiples of 1884 days providing equally probable solutions (at least mathematically, if not physically).

The sizable cross-correlation amplitudes to the periodic model imply that the WISE light curves are also periodically variable. We checked this by fitting the infrared light curves of PG 1302-102 with periodic functions and calculating the signal-to-noise ratio (S/N) of the periodic modulation above white noise. Adopting the sinusoidal period from Graham et al. (2015a), we fit the magnitude m(t) at the observed-frame time t as $m(t)=A\mathrm{cos}(2\pi t/T)+B,$ solving for amplitude A and magnitude zero-point B and fixing the period T = 1884 days. Applying the best-fit values into the periodic S/N indicator of Horne & Baliunas (1986), we get values of 5.7, 4.1, and 3.4 in the optical, W1, and W2, respectively. The false-alarm probability F that a periodicity in the data arises from pure noise, Equations (13) and (23) from Horne & Baliunas (1986) yield F ∼ 3 × 10⁻⁵⁰ from the CRTS data alone and ∼4 × 10⁻¹⁰ and ∼1 × 10⁻⁹ for W1 and W2, where these estimates of F were calculated assuming evenly spaced data. The extremely small F values imply with high likelihood that the source is periodic in all three wavebands.

We next attempt to estimate the most reasonable time lag between the optical and mid-infrared light curves, τ in Equation (1), out of the infinite number of solutions allowed. First, we assume that the time lags are positive such that they represent the reverberation of the incident optical photons through dust re-processed emission, requiring $n\geqslant 0$ in Equation (1).

Further refinement of the value of n requires a physical model of the system. In the general picture of a toroidally shaped dusty torus, most of the mid-infrared emission at these wavelengths will come from the inner edge of the torus, with longer-wavelength emission coming from slightly larger radii, corresponding to cooler material (e.g., Nenkova et al. 2008). For an axisymmetric model of the torus simply responding to variations in the accretion disk luminosity, the time delay then depends on the light-travel time to the inner edge of the torus and the orientation of the torus relative to our line of sight, $\tau =R\;(1+\mathrm{sin}i)/c,$ where R is the radius of the inner edge of the torus, c is the speed of light, and i is the inclination angle ($i={0}^{\circ }$ corresponds to face-on; $i={90}^{\circ }$ corresponds to edge-on). Alternatively, D'Orazio et al. (2015b) present a compelling and physically motivated model for PG 1302-102 where the modulated UV/optical emission is due to relativistic Doppler boosting of emission from a "minidisk" in a steadily accreting, unequal mass binary. In 2D and 3D hydrodynamical simulations of supermassive black hole binaries, circumprimary and circumsecondary minidisk accretion flows form around the two black holes, with a larger circumbinary accretion disk on larger scales (e.g., Shi et al. 2012; D'Orazio et al. 2013). For an unequal mass system, the secondary carves out an annular region in the circumbinary disk, limiting the amount of material that reaches the primary. Most of the emission at optical through X-ray energies will thus be associated with the rapidly accreting secondary, which D'Orazio et al. (2015b) show is likely moving at ∼0.07c with an orbital inclination of $i={60}^{\circ }-{90}^{\circ }$ for the PG 1302-102 system. Relativistic Doppler boosting then produces a system akin to a rotating lighthouse, producing a sinusoidal light curve with a period T. Whereas the peaks of the UV and optical light occur when the "lighthouse" is pointed toward the observer's line of sight, peaks in the thermal mid-infrared emission, which is predominantly coming from the inner wall of the far side of the torus, correspond to when the "lighthouse" is pointed away from the observer, giving an additional T/2 to the infrared time lag.

We can now use these models to place limits on the dust lags based on the sublimation radius of hot dust. The physically motivated lighthouse model is effective when the inclination is large so that the relativistic boosting is working, and the near side of the torus is blocked from our line of sight. Thus, the variable circumbinary accretion disk model predicts $\tau =R\;(1+\mathrm{sin}i)/c$ and is appropriate for small i, while the relativistic Doppler boosting model predicts $\tau =R\;(1+\mathrm{sin}i)/c+T/2$ and requires large i. The inner radius of the torus corresponds to the dust sublimation radius R_sub, i.e., the closest distance to the central engine that dust is not destroyed, and we solve for n in the W1 time lag from Equation (1) where the two radii are closest to each other, or $R\simeq {R}_{{\rm{sub}}}.$ Graphite dust grains are better able to survive in the hotter, inner region of the torus than silicate, and they are thought to be responsible for the rest-frame near-infrared emission from quasars (e.g., Mor & Netzer 2012; Jun & Im 2013). We use the graphite dust sublimation radius ${R}_{{\rm{sub}}}\simeq 0.5{({L}_{{\rm{bol}}}/{10}^{46}\;{\rm{erg}}\;{{\rm{s}}}^{-1})}^{0.5}$ pc for a sublimation temperature ${T}_{{\rm{sub}}}=1800$ K (Mor & Netzer 2012). At z = 0.278 and a median $V\simeq R\simeq 15.0$ mag for PG 1302-102 (Marziani et al. 1996; Ojha et al. 2009; Prochaska et al. 2011; Graham et al. 2015a), ${L}_{{\rm{bol}}}\simeq 6.78\times {10}^{46}\;{\rm{erg}}\;{{\rm{s}}}^{-1}$ adopting a 5100 Å to bolometric correction of 10.33 (Richards et al. 2006). Thus, we have ${R}_{{\rm{sub}}}\simeq 1550$ rest-frame light days (1.30 pc), or ${R}_{{\rm{sub}}}\simeq 1980$ observed-frame light days. For a near face-on orientation of the torus ($i\simeq {0}^{\circ }$), $\tau =R\;(1+\mathrm{sin}i)/c\simeq R/c$, and the n = 1 solution from Equation (1), (${\tau }_{W1},{\tau }_{W2})=(2219\pm 153,2408\pm 148)$ days, best meets $R\simeq {R}_{{\rm{sub}}}$ from the W1 time lag, with $R=1.12{R}_{{\rm{sub}}}$ = 1.46 pc. Meanwhile, for a relativistic Doppler-boosted torus with i = 60°–90°, the n = 2 solution, (${\tau }_{W1},{\tau }_{W2})=(4103\pm 153,4292\pm 148)$ days, gives the best match between R and R_sub with $R=0.80{\rm{-}}0.86\;{R}_{{\rm{sub}}}=1.04{\rm{-}}1.11$ pc. We note that the uncertainty in the $R\simeq {R}_{{\rm{sub}}}$ approximation could be a few tens of percent (Mor & Netzer 2012), and that even with the same adopted T_sub = 1800 K, the value of R_sub varies across the literature by 20%–60% due to slightly different model assumptions (e.g., Barvainis 1987; Nenkova et al. 2008). Also, the uncertainty in the bolometric correction from the mean quasar template could be as large as 50% (Richards et al. 2006), translating into 25% uncertainty in R_sub.

We now consider the most probable WISE time lags in the context of the observed R_dust–L relation from Koshida et al. (2014). (Note that since Koshida et al. 2014 studied variable but non-periodic AGNs and the Doppler-boosted lighthouse model is not likely appropriate for their single black hole systems, their work does not suffer the same degeneracy with integer multiples of the variability period.) The WISE W1 and W2 filters correspond to rest-frame 2.6 and $3.6\;\mu {\rm{m}}$ at $z=0.278.$ Using the median optical magnitudes of $V\simeq R\simeq 15.0,$ the rest-frame V-band luminosity of PG 1302-102 is ${L}_{V}\simeq 6.56\times {10}^{45}\;{\rm{erg}}\;{{\rm{s}}}^{-1}.$ Based on the K-band ${R}_{{\rm{dust}}}-L$ relation (Koshida et al. 2014), the expected observed-frame time lag at z = 0.278 is $\tau \simeq 1590\pm 150$ days. Our W1 and W2 lags are above the Koshida et al. (2014) relation by 0.15 and 0.18 dex for n = 1 and 0.41 and 0.43 dex for n = 2. The former offsets are comparable to the observed scatter of the relation (0.14 dex), but the latter are not, consistent with the $T/2$ delay effect and the large inclination angle required by the lighthouse model. We also note that the K-band ${R}_{{\rm{dust}}}-L$ relation generally predicts ${R}_{{\rm{dust}}}$ smaller than R_sub by factors of a few, and Koshida et al. (2014) describe how R_sub depends both on dust grain size and the details of the dust grain absorption efficiency. Also, Nenkova et al. (2008) note that the largest grains could survive to closer radii than R_sub, presumably detected by the K-band time lags.

Previous near- and mid-infrared studies of local AGNs have shown dust lags are typically longer at longer wavelengths since cooler material (e.g., farther out material) will dominate the emission at progressively longer wavelengths (e.g., Kishimoto et al. 2011b; Vazquez et al. 2015). The W2 time lag is consistent with, or only marginally larger than, the W1 lag with a ratio of 1.09 ± 0.10 for the n = 1 and 1.05 ± 0.05 for the n = 2 solution in Equation (1). Following Vazquez et al. (2015), we quote the radial distribution of dust grains in radiative equilibrium (e.g., Barvainis 1987; Nenkova et al. 2008) scaling with temperature as ${R}_{{\rm{dust}}}\sim {T}^{-2.6}$ and find ${R}_{{\rm{dust,}}W2}/{R}_{{\rm{dust,}}W1}\simeq 2.3,$ much higher than the value from this study. Alternative explanations such as clumpy tori (Nenkova et al. 2008) are successful in emitting a broad infrared spectrum from a confined distribution of dust.

Graham et al. (2015a, 2015b) list several possible physical mechanisms to produce the periodic optical light curves in PG 1302-102: non-thermal contribution from a precessing jet in a supermassive black hole binary, periodic mass accretion in a circumbinary accretion disk, or a precessing warped disk (which might also be related to a binary supermassive black hole). Our results suggest that the observed mid-infrared variability is being driven by variations in the accretions disk(s). For example, in the precessing jet scenario, one would expect very little temporal offsets, if any, between the optical and infrared variability. Indeed, the observed periodic signal in the WISE data and the cross-correlation results support that dust is reverberating input luminosity from the accretion disk.

The elliptical host galaxy morphology of PG 1302-102 and the presence of close companion galaxies around it (Bahcall et al. 1995; Disney et al. 1995) are suggestive of a recent merger driving high-luminosity AGN activity (e.g., Treister et al. 2012; Hong et al. 2015). Indeed, the host galaxy of PG 1302-102 measured from H- and R-band imaging are classified as asymmetric and tidal tail morphologies, respectively (Guyon et al. 2006; Hong et al. 2015). The galaxy-scale morphology hints that this system may be in an ongoing merger.

The mid-infrared light curves of PG 1302-102 are indicative of periodic dust reverberation lagging behind the optically periodic light curve due to the light-travel time to the mid-infrared-emitting dusty torus, with a possible additional phase lag in the relativistic Doppler boosting lighthouse model. For both models, the time lags imply that the dust is echoing outside the accretion disk but from a small, possibly clumpy medium, and we predict that the narrow flourescent Fe Kα X-ray emission should vary on the same period and with the same phase lag as the mid-infrared emission. Time domain studies are starting to reveal events in AGNs occurring on physical scales that are not resolvable from the Earth, helping to constrain the source of variability and properties of the variable structure (e.g., Risaliti et al. 2013; LaMassa et al. 2015; Mehdipour et al. 2015; Stern et al. 2015). In particular, multiwavelength time domain surveys of AGNs combining optical data with mid-infrared data from WISE, as done here, or from the Spitzer Space Telescope (e.g., Gorjian et al. 2014; Vazquez et al. 2015), probe the physical structure and accretion mechanisms of the AGN accretion disk and surrounding dusty material.

We thank the anonymous referee for the comments that greatly improved the paper, as well as Daniel D'Orazio, Moshe Elitzur, Saavik Ford, Zoltan Haiman, Barry McKernan, and Robert Nikutta for helpful discussions. This research was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities through a contract with NASA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. CRTS was supported by the NSF grants AST-1313422 and AST-1413600. D.S. acknowledges support from NASA through ADAP award 12-ADAP12-0109.

Facilities: - WISE - Wide-field Infrared Survey Explorer - .