A PERIODICALLY VARYING LUMINOUS QUASAR AT z = 2 FROM THE PAN-STARRS1 MEDIUM DEEP SURVEY: A CANDIDATE SUPERMASSIVE BLACK HOLE BINARY IN THE GRAVITATIONAL WAVE-DRIVEN REGIME

The expectation for the existence of supermassive black hole binaries (SMBHBs) in galaxy nuclei is supported by two well-established properties of galaxies: (1) high spatial resolution observations of nearby galaxies have demonstrated that SMBHs are ubiquitous in the centers of galaxy bulges (Magorrian et al. 1998) with masses tightly correlated with the mass and structure of their host galaxies (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001), and (2) galaxies in a ΛCDM universe build up their structure hierarchically through mergers (e.g., Springel et al. 2005). When two galaxies merge, their SMBHs will sink to the center through dynamical friction and through three-body interactions with stars and viscous exchange of angular momentum with circumbinary gas form a gravitationally bound binary that eventually coalesces due to the radiation of gravitational waves (GWs; Begelman et al. 1980).

Recent progress has been made in the detection of "dual active galactic nuclei (AGNs)," double active nuclei in assumed merged galaxy systems with kiloparsec-scale separations (Komossa et al. 2003; Comerford et al. 2009). These dual AGNs, while a product of a galaxy merger, are not yet gravitationally bound, and thus are not necessarily fated to coalesce. A true SMBHB becomes gravitationally bound on the scale of parsecs, which, beyond our Local Group of galaxies, is well below the angular resolving power of the most powerful current, or even future, telescopes. However, several promising candidates have been identified indirectly via spectroscopy: quasars with offset and/or drifting broad-line peaks attributed to a broad-line region in orbit around an SMBH's binary companion (Boroson & Lauer 2009; Dotti et al. 2009; Barrows et al. 2011). However, alternative scenarios have been proposed that do not require an SMBHB, including double-peaked lines from a single accretion disk (Chornock et al. 2010).

A promising observational signature of SMBHBs is their variable accretion luminosity. One of the first sub-parsec SMBHB candidates, OJ287, was identified by its variability behavior (Lehto & Valtonen 1996). OJ287 is a quasar that demonstrates regular optical outbursts on a timescale of 12 yr that has been modeled as the result of a secondary SMBH companion passing through the primary SMBH's accretion disk (Valtonen et al. 2008). Such a configuration should be rare since the secondary BH's orbital axis must be highly misaligned with the primary BH's accretion disk axis in order for it to pass through its disk. A more generic signature of an SMBHB is likely to be related to accretion through its circumbinary disk.

In a gas-rich galaxy merger, strong gravitational torques drive gas inward, triggering both star formation and BH accretion (Hopkins et al. 2006). In particular, hydrodynamical simulations of circumbinary disks show that accretion via "hot streams" onto the BHs is strongly modulated by the binary's orbital motion for mass ratios of $0.05\lt q\leqslant 1$ (MacFadyen & Milosavljević 2008; Shi et al. 2012; Noble et al. 2012; D'Orazio et al. 2013; Gold et al. 2014). Simulations (D'Orazio et al. 2013; Farris et al. 2015) also detect a $t\sim 6$ ${{t}_{{\rm orb}}}$ timescale originating from a surface density "lump" just outside the central cavity of the circumbinary disk, which is a persistent but secularly evolving feature.

While working on our paper, we became aware of the report of PG 1302-102, a periodically variable quasar discovered by the Catalina Real-Time Transient Survey (CRTS; Graham et al. 2015). It is a 15 magnitude quasar at z = 0.2784, varying at the 0.14 mag level with a period of 5.2 ± 0.2 yr, with good sampling over 1.8 cycles and extended archival data going back 20 yr. Their physical interpretation for its variability is an SMBHB (${\rm log} (M/{{M}_{\odot }})\sim 8.5$, $a\sim 0.01$ pc), its luminosity being modulated due to either a precessing jet or an overdensity ("hot spot") in the inner edge of its circumbinary disk.

In this Letter, we present our most significant detection from the systematic search for periodically varying quasars in the Pan-STARRS1 Medium Deep Survey (PS1 MDS) field MD09. PSO J334.2028+01.4075 is a radio-loud quasar at z = 2.060 with archival spectroscopy from FQBS and extended baseline photometry from CRTS. The 8.5 yr baseline of the PS1+CRTS light curve is well described by a simple sinusoid, consistent with theoretical simulations for the modulated accretion rate in a $0.05\lt q\lt 0.25$ mass-ratio SMBHB. We use the rest-frame period and virial black hole mass estimate to infer an orbital separation of the binary that is in the GW-driven regime.

The dynamics of a gravitationally bound SMBHB system can be described by Kepler's third law:

$$\begin{eqnarray}&&{{t}_{{\rm orb}}}=0.88\;{\rm yr}\left( \frac{M}{{{10}^{7}}{{M}_{\odot }}} \right){{\left( \frac{a}{{{10}^{3}}{{R}_{s}}} \right)}^{3/2}},\end{eqnarray} \tag{ 1 }$$

where R_s is the Schwarzschild radius, ${{R}_{{\rm s}}}=\frac{2GM}{{{c}^{2}}}$; a is the separation between the BHs; and M is the total mass of the system. Haiman et al. (2009) calculate the minimum survey area to detect a statistically significant sample of quasars powered by SMBHBs as a function of variable magnitude depth, assuming reasonable values for the quasar luminosity function, quasar lifetime (${{t}_{Q}}={{10}^{7}}$ yr), the Eddington fraction (${{f}_{{\rm Edd}}}=0.3$), and the fractional variability amplitude (${\Delta}f/f=0.1$). The variability detection threshold we have achieved in MD09 (Section 4) corresponds to a ${\Delta}f/f\gt 0.1$ sensitivity for point sources brighter than $m\sim 21$ mag, and thus a variable magnitude of 23.5 mag. At this depth, Haiman et al. (2009) require an area of ∼100 deg² to yield a sample of over 100 SMBHBs; an excellent match to the area of PS1 MDS (80 deg²). Furthermore, the baseline (4.2 yr) and cadence (3 days) of PS1 MDS make us sensitive to timescales for which BHs (with $M\gt {{10}^{7}}{{M}_{\odot }}$) are in the GW-driven regime of orbital decay.

The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) is a wide-field imaging system designed for dedicated survey observations on a 1.8 m telescope on Haleakala, Hawaii, with a 1.4 Gigapixel camera and a 7 deg² field of view (Kaiser et al. 2010). The PS1 telescope is operated by the Institute for Astronomy (IfA) at the University of Hawaii and has just completed over 4 yr of operation in 2014 March. We present data from the Medium Deep Survey (MDS), a deep, multi-epoch survey of 10 circular fields distributed across the sky, each ∼8 deg² in size, whose daily observing cadence in five filters is excellent for studying persistent variable sources, including quasars. The PS1 MDS cadence of observation cycles through the ${{g}_{{\rm P}1}}$, ${{r}_{{\rm P}1}},$ ${{i}_{{\rm P}1}}$, and ${{z}_{{\rm P}1}}$ bands every three nights, with observations in the ${{y}_{{\rm P}1}}$ band close to the full Moon. Due to the poorer time sampling of the ${{y}_{{\rm P}1}}$ observations, we do not use them in this analysis.

We began our systematic search for SMBHB candidates among color-selected quasars in the PS1 MD09 field. This is the first MD field that was made available to the PS1 Science Collaboration in the Pan-STARRS Science Interface (PSI) online database. In order to maximize our sensitivity to intrinsic variability, we first applied the technique of differential ensemble photometry. This technique is able to correct for local systematic errors due to variable atmospheric conditions by comparing a target object with nearby non-variable stars (Honeycutt 1992; Bhatti et al. 2010). We created a color-selected reference star sample and quasar sample by cross-matching point sources ($m\lt 23$ mag) in the MD09 field with a custom catalog extracted from full-survey deep stacks from PS1 MDS in the ${{g}_{{\rm P}1}}$, ${{r}_{{\rm P}1}}$, ${{i}_{{\rm P}1}}$, ${{z}_{{\rm P}1}}$, and ${{y}_{{\rm P}1}}$ bands, as well as from observations with the Canada–France–Hawaii Telescope (CFHT) in the u band (S. Heinis et al. 2015, in preparation). We converted all magnitudes to the SDSS photometric system (Tonry et al. 2012) to take advantage of the SDSS color selection of stars and quasars already available in the literature (Schmidt et al. 2010; Sesar et al. 2007). Figure 1 shows the color–color diagrams of the point sources in MD09 selected as quasars and non-variable stars. This query resulted in 8158 reference stars and 316 quasars, each with an average of 350 detections in four filters.

Zoom In Zoom Out Reset image size

We modified the ensemble photometry software developed by Bhatti et al. (2010) for SDSS to the PS1 data format and ran it on the reference stars. In Figure 1 (bottom right panel), we plot the "corrected" magnitude error as a function of mean magnitude compared with the "raw" values before ensemble photometry. The ensemble photometry reduces the measured errors significantly, lowering the error floor from 0.045 to 0.025 mag, and resulting in a 2σ variability threshold of 0.05 mag on the bright end to 0.34 mag on the faint end.

We then applied ensemble photometry to the 316 color-selected quasars and flagged quasars as variable based on their magnitude error relative to their neighbors of a similar brightness; we set $2\sigma$ as our criterion for variability and required a variability flag in at least two filters. This selection yielded 168 variable quasars in MD09. Among these color-selected variable quasars, we searched for potential periodic signatures using the Lomb–Scargle (LS) periodogram, a Fourier analysis technique of unevenly spaced data with noise (Lomb 1976; Scargle 1982; Horne & Baliunas 1986). For N₀ data points in the time series spanning a total length of T in units of MJD, we sampled the periodogram at the number of recommended independent frequencies (N_i), from Horne & Baliunas (1986), from $1/T$ to ${{N}_{0}}/(2\;T)$ (which would be the Nyquist frequency if data were evenly sampled), resulting in a frequency resolution in the periodogram of ${\Delta}f=({{N}_{0}}/2-1)/(T{{N}_{i}})\sim 2\times {{10}^{-4}}$ d⁻¹. When identifying periodic sources, we took advantage of the redundancy of PS1 MDS monitoring in four filters (${{g}_{{\rm P}1}}$ ${{r}_{{\rm P}1}}$ ${{i}_{{\rm P}1}}$ ${{z}_{{\rm P}1}}$), each with a slightly different observing cadence due to weather and scheduling constraints to help filter out false detections from aliasing by requiring that periodogram peaks be coherent across multiple filters; 40 of the 168 variable quasars survived this test.

Among the candidate periodic quasars from our periodogram analysis, here we focus on our most significant detection, PSO J334.2028+01.4075 (J2000). In Figure 2, we show its periodogram in four PS1 filters, with the strongest peak marked with a dashed line. We fold each filter light curve on this period (Figure 3) and measure the scatter of the residuals from the best-fit sine curve (${{\sigma }_{r}}$). The error of the periodogram peak frequency ($\delta f$) and the signal-to-noise ratio of the peak power (ξ) can then be calculated from ${{\sigma }_{r}}$ and the amplitude of the signal A₀ (Horne & Baliunas 1986) as $\delta f=\frac{3{{\sigma }_{r}}}{4\sqrt{{{N}_{0}}}T{{A}_{0}}}$ (which gives us an error on the detected period of $\delta P=\delta f/{{f}^{2}}$) and $\xi \equiv A_{0}^{2}/(2\sigma _{r}^{2})$, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The resulting average period across all four filters is $P=541.8\pm 15.3$ days, with the highest signal-to-noise ratio in the ${{g}_{{\rm P}1}}$ filter with $\xi =3.19$, and periodogram peaks in all four filters well above a $1.5\times {{10}^{-23}}$ false-alarm probability (corresponding to 10σ) plotted with a dotted line in Figure 2. The PS1 data cover 2.6 cycles, just shy of the "rule of thumb" number of cycles (three) for a periodic variation to be apparent to the eye (Press 1978). From our Monte Carlo simulations of 1000 stochastic Damped Random Walk (Kelly et al. 2009) light curves, we find a false periodic detection rate of 6.3% using our selection criteria from Section 5. We further disfavor a false alarm from stochastic quasar variability since the 0.6% of the simulations that successfully mimic the periodic timescale of our candidate have short-timescale variances a factor of $\gtrsim 2$ larger than expected for the quasar's luminosity and inferred black hole mass (Section 7). Note that there is a secondary peak in the ${{g}_{{\rm P}1}}$ and ${{r}_{{\rm P}1}}$ periodograms that, if real and not an artifact from the PS1 data sampling, is at twice the primary peak frequency, a signature of $0.05\lt q\lt 0.25$ mass-ratio SMBHBs, which show an accretion rate modulation most closely described by a simple sinusoid (D'Orazio et al. 2013). The amplitude of PSO J334.2028+01.4075ʼs sinusoidal modulation increases with decreasing wavelength, consistent with the exponential dependence on wavelength found in previous quasar variability studies (e.g., Vanden Berk et al. 2004; MacLeod et al. 2010).

PSO J334.2028+01.4075 is a radio-loud quasar (FBQS J221648.7+012427) with an archival spectrum from the FIRST Bright Quasar Survey (Becker et al. 2001). We are also fortunate that this candidate has an archival V-band light curve from CRTS (Drake et al. 2009), which we use to test the persistence of the periodic variation over an extended baseline of 8.5 yr (corresponding to 5.7 cycles). To compare to the CRTS light curve, we convert the PS1 ${{g}_{{\rm P}1}}$-band light curves to the SDSS system (Tonry et al. 2012) and then to the Johnson V magnitude using the photometric transformation for quasars from Jester et al. (2005) and an average ${{g}_{{\rm P}1}}$−${{r}_{{\rm P}1}}\;=\;+0.10$ mag. We had to apply an additional offset of −0.17 mag to the pseudo-V PS1 magnitudes in order to match the average of the CRTS V-band data. Though the photometric errors are relatively large, the CRTS measurements are consistent with those of PS1 during their overlap (Figure 4) and have residuals over the entire CRTS baseline from the PS1-fitted sinusoidal model that are Gaussian with a $\sigma =0.17$ mag that is comparable to the mean photometric error of 0.18 mag.

Zoom In Zoom Out Reset image size

We use the width of the quasar's C IV line and its nearby continuum luminosity to make a virial estimate of the black hole mass from Vestergaard & Peterson (2006):

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\rm log} \left( \frac{{{M}_{{\rm BH}}}}{{{M}_{\odot }}} \right) & = & {\rm log} \left[ {{\left( \frac{{\rm FWHM}({\rm C}\ {\rm IV})}{1000\;{\rm km}\;{{{\rm s}}^{-1}}} \right)}^{2}}{{\left( \frac{\lambda {{L}_{\lambda }}}{{{10}^{44}}\;{\rm erg}\;{{{\rm s}}^{-1}}} \right)}^{0.53}} \right] \\ {} & {} & +6.66, \\ \end{array}\quad \end{eqnarray} \tag{ 2 }$$

where $\lambda =1350$ Å. The FBQS spectrum, though not publicly available in an electronic format, was measured with a ruler to determine ${{F}_{\lambda ,{\rm obs}}}$(1350 Å (1+z)) $\sim 8.5\times {{10}^{-17}}$ ergs s⁻¹ cm⁻² Å⁻¹ and FWHM (CIV $\lambda 1550$) ∼200 Å. The CIV line is symmetric in shape, and its width corresponds to a velocity in the rest-frame of ∼12,650 km s⁻¹. We correct for a Galactic extinction of $E(B-V)=0.0406$ mag (Schlafly & Finkbeiner 2011), using the extinction law from Cardelli et al. (1989), to find ${{L}_{\lambda ,{\rm em}}}(1350\;\overset{\circ}{\rm A} )=4\pi d_{L}^{2}{{F}_{\lambda ,{\rm obs}}}{{10}^{{{A}_{\lambda }}/2.5}}(1+z)\sim 9.5\times {{10}^{42}}$ erg s⁻¹Å⁻¹, where d_L is the luminosity distance assuming ${{H}_{0}}=70\;{\rm km}\;{{{\rm s}}^{-1}}\;{\rm Mp}{{{\rm c}}^{-1}}$, ${{{\Omega}}_{m}}=0.3$, ${{{\Omega}}_{{\Lambda}}}=0.7$, and k = 0, and ${{A}_{\lambda }}=0.1790$. This results in a black hole mass of ${\rm log} \left( \frac{{{M}_{{\rm BH}}}}{{{M}_{\odot }}} \right)=9.97$, with a scatter from the uncertainty in the relation of 0.5 dex. Applying a mean quasar bolometric correction at 1350 Å of BC=3.81 from Richards et al. (2006), one gets a bolometric luminosity of ${{L}_{{\rm bol}}}=\lambda {{L}_{\lambda }}{\rm BC}=4.9\times {{10}^{46}}$ erg s⁻¹. Note that this object is also radio loud, with a radio luminosity at the rest-frame frequency of 5 GHz of ${\rm log} ({{L}_{R}}$ (erg s⁻¹)) = 32.8 from Becker et al. (2001).

Assuming the rest-frame period ${{P}_{{\rm rest}}}={{P}_{{\rm obs}}}/(1+z)$ is on the order of the orbital timescale of the SMBHB, with a caveat that in addition to a strong dependence on mass ratio, there are a range of theoretical predictions for translating ${{P}_{{\rm rest}}}$ to ${{t}_{{\rm orb}}}$ (e.g., Noble et al. 2012, D'Orazio et al. 2013), we then calculate the orbital separation of the binary to be $7_{-4}^{+8}{{R}_{s}}$ ($\sim 0.006_{-0.003}^{+0.007}$ pc), securely placing it in the gravitationally bound regime of a physically viable SMBHB system—a circumbinary accretion disk system capable of maintaining a central cavity, stable to gravitational fragmentation, and in the regime of orbital decay driven by GWs (Haiman et al. 2009; Kocsis et al. 2012; D'Orazio et al. 2013). Also note that since the viscous time scales as r², one could expect to be able to detect modulations in the accretion rate fed by the streams in the circumbinary disk cavity, without being washed out by viscous processes in the "mini disks" around each BH (Roedig et al. 2014). Remarkably, the rest-frame inspiral time for the binary is ${{t}_{{\rm insp}}}=\frac{5}{256}\frac{{{c}^{5}}}{{{G}^{3}}}\frac{{{a}^{4}}}{{{M}^{2}}\mu }=7.0\;{\rm yr}{{\left( a/7{{R}_{s}} \right)}^{4}}{{\left( M/{{10}^{9.97}}{{M}_{\odot }} \right)}^{-3}}$ for q = 0.25, where $\mu \equiv \frac{{{M}_{1}}{{M}_{2}}}{{{M}_{1}}+{{M}_{2}}}$, opening up the possibility for detecting the decay of the orbital period $(\dot{P})$ with future monitoring, as well as providing a promising target for direct GW detection for pulsar timing arrays (Sesana et al. 2009).

We present the most statistically significant periodically variable quasar candidate from our search in PS1 MD09, PSO J334.2028+01.4075, a radio-loud quasar at z = 2.060. We combine an estimate of its black hole mass with its variability timescale (assuming ${{P}_{{\rm rest}}}\sim {{t}_{{\rm orb}}}$) to find orbital parameters consistent with model predictions of a stable accreting SMBHB system with a $0.05\lt q\lt 0.25$ in the GW-driven regime (Haiman et al. 2009).

The redshift of this SMBHB candidate coincides with the peak epoch for SMBHB mergers (Volonteri et al. 2003), and its large mass ($M\approx {{10}^{10}}{{M}_{\odot }}$) makes it favorable for detection in the GW-driven regime due to the strong dependence on M of the residence time at a given orbital separation in units of R_s (Haiman et al. 2009). Like the CRTS SMBHB candidate PG 1302-102 reported by Graham et al. (2015), our SMBHB candidate is also a radio-loud quasar. However, given the shorter rest-frame period of our candidate of 0.5 yr (versus 4 yr in PG 1302-102), it is even more unlikely for its variability to be driven by jet precession, either originating from a single SMBH (Lu & Zhou 2005) or a binary SMBH (Lobanov & Roland 2005).

This pilot program in PS1 MD09 is a promising start to our systematic search for periodic variability signatures of SMBHBs among the expected $\approx 1000$ variable quasars in the full ∼80 deg² PS1 MDS. At the start of the next decade (∼2023), the Large Synoptic Survey Telescope (Ivezic et al. 2008) will probe a volume several thousand times larger than PS1 MDS, yielding tens of millions of quasars, and potentially thousands of SMBHBs periodically varying on the timescale of years, fated to coalesce.

S.G. thanks the Aspen Center for Physics, which is supported in part by the NSF under grant No. PHYS-1066293, for hosting the "Black Holes in Dense Star Clusters" Winter Workshop which facilitated stimulating discussion of SMBHB theory. T.L. thanks J. Krolik, R. Mushotzky, and the anonymous referee for helpful comments on the manuscript. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the IfA, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, MPIA, Heidelberg and MPE, Garching, JHU, Durham University, the University of Edinburgh, QUB, the Harvard-Smithsonian CfA, LCOGT Inc., the National Central University of Taiwan, STScI, NASA under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, NSF under grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University.