{'fraction_non_alphanumeric': 0.03855860265214665, 'fraction_numerical': 0.03270063168198712, 'mean_word_length': 4.623774298985033, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery of a Candidate Binary Supermassive Black Hole in a Periodic Quasar from Circumbinary Accretion Variability

Binary supermassive black holes are expected to be a generic byproduct from hierarchical galaxy formation. The final coalescence of BSBHs is thought to be the loudest gravitational wave siren in the universe, yet no confirmed BSBH is known in the GW-dominated regime. While periodic quasars have been proposed as BSBH candidates, the physical origin of the periodicity has been largely uncertain. Here we report discovery of a periodicity (P~1607 days) at 99.95% significance in the optical light curves of a quasar at z=1.53, SDSS J025214.67-002813.7. Combining archival Sloan Digital Sky Survey data with new, high signal to noise imaging from the Dark Energy Survey, the total ~20yr time baseline spans ~4.6 cycles of the observed 4.4-yr (restframe 1.7-yr) periodicity. The light curves are best characterized by a bursty model predicted by hydrodynamic simulations of circumbinary accretion disks. The periodicity is likely caused by accretion rate variability periodically modulated by a milli-parsec BSBH emitting gravitational waves while still being dynamically coupled to the circumbinary accretion disk. A characteristic bursty hydrodynamic variability model is statistically preferred over a smooth, sinusoidal model expected from relativistic Doppler boost, a kinematic effect proposed for PG1302-102. Furthermore, the frequency dependence of the variability amplitudes disfavors Doppler boost as the dominant mechanism, lending independent support to the circumbinary accretion variability hypothesis. Given our detection rate of one BSBH candidate from circumbinary accretion variability out of 625 quasars, it suggests that future large, sensitive synoptic surveys such as the Legacy Survey of Space and Time at Vera C. Rubin Observatory may be able to detect hundreds to thousands of candidate BSBHs from circumbinary accretion variability with direct implications for Laser Interferometer Space Antenna.

However, no confirmed case is known at sub-milliparsec scales, i.e., separations close enough to be in the GWdominated regime. While ∼150 periodic quasars have been suggested as close BSBH candidates (e.g., Valtonen et al. 2008;Graham et al. 2015;Charisi et al. 2016;Liu et al. 2019;Zhu & Thrane 2020), even the most promising candidates are subject to false positives due to quasar's stochastic, red noise variability, given the limited time baseline and relatively low sensitivity of existing surveys (e.g., see Vaughan et al. 2016 for evidence against any significant periodicity in PG 1302−102 and Goyal et al. 2018 in the case of OJ 287). The study of periodic quasars is important to the searches for close BSBHs in order to test theories of BSBH evolution to shed light on the expected rate of BSBH mergers as GW sources. The study is also important for understanding the physical origin of quasar periodicity, which is largely unknown.
Circumbinary accretion disks are generally expected around close BSBHs at the inferred binary separations of the candidate periodic quasars. Theory suggests that hydrodynamic variability in the circumbinary accretion disks may cause periodic light curves due to accretion rate modulation from the binary torque (e.g., Farris et al. 2014;Gold et al. 2014b;Shi & Krolik 2015;Duffell et al. 2019). This should be useful for finding BSBHs that are close enough to be emitting gravitational waves. However, no evidence has been found for the generic "sawtooth" pattern (i.e., with a sharp rise and a gradual decay, in contrast to a more smooth, sinusoidal modulation expected from Doppler beaming (e.g., D'Orazio et al. 2015a;Duffell et al. 2019), largely limited by the relatively low sensitivity of previous surveys.
In this paper, we present a significant periodicity discovered in the optical light curves of a redshift z=1.53 quasar, SDSS J025214.67−002813.7 (hereafter J0252 for short). Our systematic search combines new, highly sensitive light curves from the Dark Energy Survey (Dark Energy Survey Collaboration et al. 2016) Supernova (DES-SN) fields (2012-2019; Bernstein et al. 2012;Goldstein et al. 2015;Kessler et al. 2015;Tie et al. 2017) with archival data from the SDSS Stripe 82 (S82) survey (1998Ivezić et al. 2007). Unlike previous studies, which were based on few-cycle (e.g., ∼1.5) searches given the limited time baselines, the periodicity of J0252 was discovered based on ∼5 cycles enabled by a ∼20-yr long baseline. The long baseline and high sensitivity are instrumental in rejecting false positives and recovering false negatives caused by stochastic quasar variability (e.g., Vaughan et al. 2016). Furthermore, we show that the distinct "sawtooth" pattern (expected from hydrodynamic circumbinary accretion disk variability models) is favored over a smoother, sinusoidal expected from Doppler beaming (e.g., D'Orazio et al. 2015a;Duffell et al. 2019). In addition, the frequency dependence of the variability amplitudes disfavors Doppler beaming, lending further support to the circumbinary accretion model.
The paper is organized as follows. §2 describes the data and methods. §3 presents our results on the detection of a significant periodicity in J0252 and its relevant physical properties. We discuss the implications of our results in §4 and conclude in §5. A concordance ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 is assumed throughout. We use the AB magnitude system (Oke 1974) unless otherwise noted.

Sample Selection
We start with 763 spectroscopically confirmed quasars in the 4.6 deg 2 overlapping region between the SDSS Stripe 82 survey and the DES-SN fields (S1 and S2). They include 758 objects in the SDSS DR7/DR14 quasar catalogs Pâris et al. 2018) and/or the OzDES quasar catalog (Tie et al. 2017;Childress et al. 2017), as well as 5 objects supplemented from the Million Quasars Catalog (Flesch 2015) (v5.5, 14 November 2018. We focus on spectroscopically confirmed quasars to ensure a clean sample in this pilot study. Only point sources are included in the analysis to avoid systematics from host galaxy contamination. We request that the DES flag SPREAD MODEL <0.005, i.e., the difference between the source point spread function (PSF) and the local PSF model is smaller than 0.5%, or the source PSF is smaller than the local model PSF. We further require that a quasar has at least 30 >3σ SDSS epochs and 50 DES epochs in at least two bands. The final parent sample consists of 625 quasars. They have a median spectroscopic redshift of ∼1.8 and a median average i-band PSF magnitude of 21.0 mag (AB). The median epoch of observations is 80 from the SDSS and 135 from the DES.

Light Curve Data
We combine archival light curves from the SDSS Stripe 82 survey with new observations from the DES-SN fields (described in detail below). The time baseline of the combined light curves extends ∼20 yr (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) from SDSS Stripe 82 and 2012-2019 from DES-SN). For a typical quasar at z ∼1, the time baseline spans ∼ 10 yr in the quasar rest-frame to encompass 5 cycles for a period of 2 yr, which is the recommended number of cycles to minimize false periodicity (Vaughan et al. 2016). We have rejected > 5σ outliers from the running median in each band. We have binned the data within the same Julian date for a better S/N. We quote AB magnitudes throughout unless otherwise noted.
The DES is a wide-area 5000 deg 2 survey of the Southern Hemisphere in the grizY bands (Flaugher 2005 (Flaugher et al. 2015;Bernstein et al. 2017) with a 2.2 degree diameter field of view mounted at the prime focus of the Victor M. Blanco 4m telescope on Cerro Tololo in Chile. The typical single-epoch 5σ point source depths (Abbott et al. 2018) are g=24.3, r=24.1, i=23.5, z=22.9, and Y =21.4, much deeper than other surveys of larger area (e.g., SDSS and PanSTARRS1). The data quality varies due to seeing and weather variations. The DES absolute photometric calibration has been tied to the spectrophotometric Hubble CALSPEC standard star C26202 and has been placed on the AB system (Oke & Gunn 1983), with an estimated single-epoch photometric statistical precision of 7.3, 6.1, 5.9, 7.3, 7.8 mmag in grizY bands (Abbott et al. 2018). The DES contains a 30 deg 2 multi-epoch survey (DES-SN) to search for SNe Ia that has a mean cadence of ∼7 days in the griz bands. Two of the ten DES-SN fields (S1 and S2) are overlapped with the SDSS Stripe 82 (with an overlapping area of 4.6 deg 2 ). We adopt light curves generated from the Y6A1 Gold data (Morganson et al. 2018). We have also included the Science Verification data to maximize the time baseline.
The SDSS equatorial Stripe 82 region was observed from September 1998 to December 2007 with ∼70-90 total epochs of images in the ugriz bands obtained in yearly "seasons" about 2-3 months long (Adelman-McCarthy et al. 2007;Ivezić et al. 2007;Frieman et al. 2008). The typical single epoch 5σ point source depths are 22.2, 22.2, 21.3, and 20.5 in the griz bands (Adelman-McCarthy et al. 2007). The photometric calibration over the survey area is accurate to roughly 0.02 mag in the gri bands, and 0.03 mag in the uz bands (Ivezić et al. 2004). All SDSS magnitudes have been calibrated to be nearly on the AB system (Abazajian et al. 2009).
To stitch together the light curves for each quasar, we apply the appropriate corrections to convert the SDSS photometry to be on the DES system. The corrections are to compensate for the filter coverage and system throughput differences between the two surveys. We estimate the corrections empirically by calculating two sets of synthetic magnitudes by convolving each quasar spectrum with the SDSS and DES system transmission curves (including both instrument and atmosphere). For J0252, the corrections are where the errors are 1σ uncertainties estimated from 100 bootstrap resampling of the observed quasar spectrum randomly perturbed by the error spectrum.
J0252 is a spectroscopically confirmed quasar contained in the SDSS DR14 quasar catalog (Pâris et al. 2018). Figure  1 shows its griz multi-color optical light curves. The ∼20year observations combine archival SDSS data with new, higher signal to noise imaging from DES. The SDSS (DES) observations included 83, 83, 84, and 85 epochs (131, 140, 141, and 143 epochs) in the griz bands with a median separation of 4 days (7 days) between epochs and yearly seasonal gaps. The variability of J0252 is more coherent with a larger amplitude than typically observed for stochastic quasar variability (Morganson et al. 2014) generally believed to result from thermal fluctuations in the accretion disks driven by an underlying stochastic process such as a turbulent magnetic field. Figure 1 also shows archival photometry from the Palomar Transient Factory (PTF), Panoramic Survey Telescope and Rapid Response System (Pan-STARRS or PS1), and Zwicky Transient Facility (ZTF) survey, as well as our new observations from the Las Cumbres Observatory Global Telescope (LCOGT; DDT Program 2018B-004 and NOAO Program 2019A-0279; PI Liu). This is for the purposes of independent double checks only. We do not include them in our baseline analysis to: 1. ensure homogeneity for the analysis of the parent sample and 2. minimize possible uncertainties and/or caveats in the available data as we describe below. However, we have also verified that our results do not change qualitatively even when including them in the analysis. Further details of the archival photometry can be found in Appendix A, with photometry data listed in Table  A1 -A4.

Periodicity Detection
For any periodicity detection, we implement the following three selection criteria: (i) At least two bands have a 3σ detection of the same periodicity in the periodogram analysis.
(ii) The detected periodicity is the dominant component compared to the background noise.
(iii) The same periodicity is also identified in the autocorrelation function (ACF).
First, we adopt the generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009), which is appropriate for detecting periodicity in unevenly sampled data using the astroML package (Vanderplas et al. 2012). Comparing the observed power to that from the simulated light curves (see details below), we identify a significant periodicity candidate if at least two bands show a 3σ detection (i.e., the detected periodicity cannot be reproduced by >99.7% synthetic light curves) in the same periodicity grid window. To quantify the statistical significance of any periodogram peak, we adopt an approach similar to that of Charisi et al. (2016). Since the noise spectrum is frequency-dependent (due to the stochastic red noise quasar variability), it is more appropriate to quantify the statistical significance at a given frequency, i.e., as compared to the local background. Adopting a false-alarm probability that is flat over different frequencies instead would overestimate the true statistical significance of periodogram peaks (Liu et al. 2015). In addition, we reject any detection where fewer than three cycles are spanned by the observations or where the periodicity is shorter than 500 days. The former criterion is imposed to minimize false periodicity due to the stochastic red noise quasar variability (Vaughan et al. 2016), whereas the latter is to mitigate artifacts caused by seasonal gaps and low cadence sampling on short timescales (i.e., an aliasing effect, e.g., MacLeod et al. 2010).
Second, we fit a sine curve to the selected candidates, and reject any of them if the residue noise dominates over the periodicity signal, i.e., if σ 2 residue /A 2 sin > 1, where A sin is the periodicity amplitude and σ 2 residue is the variance of the residue light curve after subtracting the pe- The solid curves show the best-fit models from hydrodynamic circumbinary accretion disk simulations assuming a mass ratio q=0.11 (Farris et al. 2014), of which the thick solid denotes our baseline model assuming a background of random, red noise variability whereas the thin solid assumes white (flat spectrum) noise for comparison purposes only. Note that because we assume red noise in our baseline model for the background signal (from stochastic variability), the residual is not supposed to be zero, unlike the case of a white noise background. Also shown for comparison are a q=0.43 accretion model (dotted gray) and a sinusoidal model (dashed gray) expected from Doppler boost (D'Orazio et al. 2015b) both assuming red noise.
riodic signal. Finally, as a complementary test, we search for periodicity by fitting the ACF with the ZDCF package (Alexander 1997). For a periodicity on top of a stochastic background, ACF has a damped periodic oscillation with ACF(t)=cos (ωt) exp (−λt), where t is the lagging time, ω = 2π P , and λ is the decay rate of the stochastic background (Graham et al. 2015). We require that the GLS periodicity be consistent with that from the ACF test. Besides the above criteria, we have tested alternatives using the multiband GLS by VanderPlas & Ivezić (2015) and the modified GLS adopted by Zheng et al. (2016) and found that they provide no further constraint in our candidate selection, i.e., the candidates selected by the three criteria also have 3σ detections in these alternative methods.

Simulated Light Curves
To quantify the statistical significance of any periodicity, first we generate 50,000 evenly sampled mocked light curves assuming a damped random walk (DRW) model with variability parameters tailored to the observed properties of each quasar. A DRW model uses a self-correcting term added to a random walk model that acts to push any deviations back toward the mean. It captures the stochastic properties of quasar variability on a timescale 10 days (Kelly et al. 2009;MacLeod et al. 2010;Mushotzky et al. 2011;Koz lowski 2016;Smith et al. 2018). The DRW model is known as a red noise model, which has a higher spectrum power in lower frequencies. Models that failed to account for this red noise feature would likely identify false positives due to the generic power spectrum feature. The DRW model is governed by two parameters: σ 2 and τ, which describe the flux variance and correlation timescale of the variability.
To measure the DRW parameters and uncertainties for each quasar, we fit the light curve directly in the time domain by treating each data point as a state space with a Gaussian uncertainty due to both the stochastic process and measurement error, following Equations (6)-(12) of Kelly et al. (2009). For unevenly sampled data, fitting the light curve directly in the time domain is preferred over fitting the power spectrum density for better recovering the true DRW parameters. The structure function due to the observational cadence may induce an anomalous power in the power spectrum that could potentially bias the fitting result. We apply a Bayesian model using the Markov Chain Monte Carlo (MCMC) method with the emcee package adopting a non-informative prior (Foreman-Mackey et al. 2013). The fitting starts with 200 walkers and samples for 1500 steps. The first 750 steps are removed as a burn-in process. To test for the convergence, we repeat the above processes but with only half of the steps (750 steps), and the resulting parameter distribution is consistent. Figure 2 shows the parameter estimation. For J0252, we have also tested a different prior with a log-normal distribution centered at 0.08 mag and 200 days for σ and τ, respectively, with a standard deviation of 1.15 (Kelly et al. 2009;MacLeod et al. 2010). The analysis is consistent across different choices of prior. Figure 2 shows the best-fit DRW parameters for J0252. We have verified that the light curve baseline (∼7300 days) is more than 10 times larger than the correlation timescale (∼630 days in the observed frame), so that the correlation timescale recovers the true value (Koz lowski 2017). We then generate 50,000 mock light curves with parameters pairs (σ 2 and τ) randomly drawn from the posterior distribution in the DRW parameter fitting. We also down sample the mock light curves to match the cadence of the observations and add measurement errors.
We also consider a bending power law (BPL) model as an alternative to the DRW model assumption for the simulated light curves. This is motivated by results based on high-cadence Kepler observations that suggest deviations from the DRW model at the high frequency end f 1 10 day −1 (Mushotzky et al. 2011;Edelson et al. 2013Edelson et al. , 2014Smith et al. 2018). For the BPL model, we assume a −3 power spectrum index at the high frequency f > 1 10 day −1 and keep the DRW model at the low frequency f < 1 10 day −1 . We also test different power law indexes and different high frequency breaks. The result is not sensitive to these choices. The low-frequency breaks are drawn from the timescales τ in the DRW parameter fitting. We first generate 50,000 evenly sampled mocked light curves assuming the BPL model with the pyLCSIM 1 package. We then down sample the light curve and measure the power spectrum density using the GLS periodogram. Our result is consistent with that assuming a pure DRW model.
In addition to a DRW model, we have also considered the CAR(2,1) ) model, i.e., a damped harmonic oscillator, which is often used to describe a periodic signal (Graham et al. 2015;Moreno et al. 2019). The quality factor Q, defined as the ratio of the detected frequency to the corresponding frequency width, is used as a measure of the significance level. For J0252, we have Q 1 for the CAR(2,1) model, which could suggest a low significance level for the detected period or a higher order noise. We have tested that the significance of the periodic signal decreases when we assume the CAR(2, 1) model for the "stochastic" component instead of a DRW, although a > 99.74% detec-1 http://pabell.github.io/pylcsim/html/index.html tion remains in g-band. However, we also show in Table 1 that, using carma-pack 2 , CAR(1,0) has the lowest BIC value and is thus a proper noise model.

eBOSS Spectrum and Analysis
J0252 has an optical spectrum available from the SDSS DR14 data archive (Plate = 7820, Fiber ID = 470, MJD = 56984). It was taken by the BOSS spectrograph within the SDSS-IV/eBOSS survey (Dawson et al. 2016). The BOSS spectrum covers 3650-10400Å with a spectral resolution of R =1850-2200. Figure 3 shows its optical (rest-frame UV) spectrum, where multiple broad emission lines are detected including C IV λ1549, He II λ1640, C III] λ1909, and Mg II λ2800.
To estimate the viral black hole mass from the broad emission lines, we follow the procedure described in Shen & Liu (2012); Shen et al. (2019) by fitting the spectral models to the observed spectra. The spectral models contain a linear combination of power-law continuum, a pseudo continuum generated from Fe II emission templates, and single or multiple Gaussian components for the emission lines. Since the errors in the continuum model might change the fitting of the weak emission lines, we perform the a global fit to the mission-line free region first to construct the continuum model better. Then, we fit multiple Gaussian models to the emission lines around the Mg II λ2800 region locally. The Mg II λ2800 line is fitted by a combination of up to two Gaussians for the broad component and one Gaussian for the narrow component. For the FWHM of the narrow lines, we also impose an upper limit of 1200 km s −1 . Figure 3 shows our spectral models for J0252.

Discovery of A Significant Periodicity in J0252
Using the three criteria described in §2.3, we identify five significant periodic candidates out of the parent sample of 625 quasars in a 4.6 deg 2 field. J0252 was the most significant detection with >4 cycles spanned. We focus on J0252 here in this paper, whereas the other four candidates will be discussed in Chen et al. in prep. Figure 4 shows the generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009). An observed 4.4-yr (corresponding to restframe 1.7-yr at the redshift of 1.53) periodicity is detected at 99.95%, 99.43%, 99.78%, and 99.59% significance in the griz bands. The confidence level in each band was determined from 50,000 Monte Carlo simulations (described in §2.4) tailored to the observed variability flux variance and characteristic timescale assuming a damped random walk (DRW; Kelly et al. 2009) or a more general bending power-law (BPL) model. There is a ∼0.05% probability that the periodogram peak is produced by stochastic quasar variability (i.e., assuming a correlated red noise), but the fact that we have found five candidates at >99.74% significance in a parent sample of 625 (in which two cases are expected from a red noise by chance) suggests that we are not just seeing stochastic quasar variability. Similar to CAR(p, q) (1,0) (2,0) (2,1) (3,0) (3,1) (3,2) BIC -967 -927 -922 -917 -920 -906 Table 1. BIC values for CAR(p,q) model using g-band data with p ≤ 3 and q < p. With CAR(1,0) having the smallest BIC value, it suggests that DRW is the proper noise model for J0252.
2 . 4 2 . 8 3 . 2 3 . 6 log 10 ⌧ (day)  Figure 2. DRW model parameter estimates for J0252. The 2D contours show the 68% and 95% confidence levels estimated from the MCMC analysis. The histograms show the projected 1D probability density distributions for σ and τ (observed frame). Labeled on each panel are their best-fit value and the 1-σ (estimated from the 68% confidence levels denoted by the shaded histograms) uncertainties. The total light curve baseline (∼7300 days) is more than 10 times larger than the correlation timescale (∼630 days in observed frame), so that the correlation timescale recovers the true value (Koz lowski 2017).
Figure 4, Figure 5 shows the significance level assuming a CAR(2,1) noise. The detected periodicity is found to have 99.8%, 98.8%, 99.5% and 98.0% confidence level in griz bands under this more complex noise model. The false alarm probability of seeing such a significant peak in the periodograms is 10 −20 assuming a pure white (i.e., flat spectrum) noise instead (Zechmeister & Kürster 2009). We also want to caution that a rigorous estimation of the DRW parameter distribution over the whole quasar sample is required in order to accurately address the false detection probability, and is out of the scope of the current analysis. Although we do account for the DRW parameter distribution (Fig 2), the parameter distribution is the result of MCMC fitting on the light curve of J0252 and is designated for J0252 only.
Archival observations from the PTF (in the gR bands) and from Pan-STARRS (in the griz bands) provide independent verification of our baseline observations. They also partially filled the cadence gap between the SDSS and DES observations. New observations from the LCOGT and the ZTF provide independent support and verification to our baseline DES observations. Despite having significant gaps, the combined time baseline spans ∼4.6 cycles of the periodicity, approaching the number of observed cycles recommended for minimizing false positives from stochastic quasar variability (Vaughan et al. 2016).

Black Hole Mass Estimation
The black hole mass is estimated using the single-epoch spectrum by assuming virialized motion in the broad-line region clouds (Shen 2013). The broad-line region gas clouds would see the candidate BSBH as a single source. From the spectral fit to the eBOSS spectrum, the Mg II λ2800-based estimator gives a virial black hole mass of M = 10 8.4±0.1 M Figure 3. Optical spectrum and modeling of J0252. Shown are the data (black), the 1σ error (gray), the best fit model (orange), the Fe II pseudo-continuum (yellow), and the broken power-law model for the emission-line-and Fe II-subtracted continuum (with the griz bands plotted in blue, green, red, and magenta, respectively).
(1σ statistical error), using the parameters in Vestergaard & Osmer (2009). Shen (2013) suggests that Mg II λ2800-based masses are more reliable than C IV λ1549-based masses, given that C IV λ1549 is likely to suffer from non-virial motion like outflows and there is larger scatter between C IV λ1549 and Hβ masses for quasars at high redshift (Shen & Liu 2012).

Radio Loudness Upper Limit
J0252 was undetected by FIRST (Becker et al. 1995) with a 3σ flux density upper limit of <0.5 mJy at 1.4 GHz. It was covered by the VLA Sky Survey (Villarreal Hernández & Andernach 2018) (VLASS) footprint at 3 GHz to a sensitivity of 0.12 mJy RMS. It was also undetected by VLASS according to its quicklook image 3 , suggesting a 3σ upper limit of f obs 3 GHz < 0.36 mJy. Assuming that the radio flux follows a power law f ν ∝ ν α , this translates into f rest 6 cm < 0.18 mJy (6 cm corresponding to 5 GHz) for a spectral index α = −0.5 (Jiang et al. 2007), or f rest 6 cm < 0.20 mJy assuming α = −0.8 (Gibson et al. 2008). Combining the f 2500 measurement from the optical spectrum, the inferred radio loudness parameter (Kellermann et al. 1989) is R ≡ f 6 cm/ f 2500 < 34 assuming α = −0.5, or R < 39 assuming α = −0.8. While the VLASS upper limit cannot exclude the possibility of J0252 being radio loud (i.e., R > 10 according to the traditional definition based on PG quasars), it does rule out its optical emission being dominated by emission from a radio jet (i.e., R > 100 (Chiaberge & Marconi 2011)).
3 http://archive-new.nrao.edu/vlass/HiPS/VLASS1.1/Quicklook/ Figure 6 shows the SED of J0252. It is similar to a control sample of ordinary optically selected SDSS quasars that are matched in redshift and luminosity. The available SED observations include a radio flux density upper limit from the VLASS, MIR photometry from WISE (Wright et al. 2010), NIR photometry from UKIDSS (Lawrence et al. 2007), optical photometry from the SDSS (York et al. 2000) and an optical spectrum from eBOSS, UV photometry from GALEX (Martin et al. 2005) (including a detection in the NUV and an upper limit in the FUV), and an X-ray upper limit from ROSAT (Voges et al. 2000).

Spectral Energy Distribution
A generic prediction from circumbinary accretion disk simulations is a flux deficit in the optical/UV SED. The flux deficit may be a cutoff from a central cavity opened by the secondary black hole (Milosavljević & Phinney 2005) or a notch from minidisks formed around both black holes (Roedig et al. 2014;Farris et al. 2015). There is tentative evidence for an NUV deficit compared to the control sample from the existing optical spectroscopy and GALEX UV photometry, but the existing data are too uncertain to draw a firm conclusion. Future HST UV spectroscopy could confirm the potential UV deficit as a complementary test of circumbinary accretion disk models.

Physical Origins of the Periodicity
In addition to a pure stochastic quasar variability (i.e., the null hypothesis), we consider two common, competing models for the optical light curve periodicity. The first is a smooth, sinusoidal model, which is expected from Doppler boosting. It has been proposed to explain the periodic quasar candidate PG1302−102 (D'Orazio et al. 2015b). The highly relativistic motion of the secondary black hole drives an apparent periodicity in the light curve, assuming that the optical emission is dominated by contribution from a mini accretion disk fueling the secondary black hole.
The second is a more bursty, quasi-periodic variability model predicted by hydrodynamic simulations of circumbinary accretion disks. We adopt the bursty hydrodynamic circumbinary accretion disk variability model of (Farris et al. 2014). The model was generated from two-dimensional (2D) hydrodynamical simulations of circumbinary disk accretion using the finite-volume code DISCO (Duffell 2016). It solves the 2D viscous Navier-Stokes equations on a high-resolution moving mesh. The moving mesh shears with the fluid flow and thereby reduces the advection error in comparison to a fixed grid. Unlike previous simulations that have excised the innermost region surrounding the binary by imposing an inner boundary condition, and so potentially neglecting important dynamics occurring inside the excised region, the model was the first 2D study to include the inner cavity using shock-capturing Godunov-type methods. The simulations last longer than a viscous time such that the solutions represent a quasi-steady accretion state.
More specifically we consider two models, mass ratio q=0.11 and q=0.43. These two values are chosen because they represent two characteristic regimes in the light-curve behaviors (Figure 9 of (Farris et al. 2014)). In the simula-  Table 1) is detected at 99.95%, 99.43%, 99.78%, and 99.59% significance in the griz bands. The confidence levels are calculated from 50,000 tailored simulations assuming random, red noise variability. The gray curves show 200 examples drawn from the 50,000 for clarity. The cyan shaded region indicates the period uncertainty estimated using ranges above the >99.74% significance for the gi bands and above the >99.00% significance for the rz bands. The gray shaded regions mark the small timescales (<500 days) on which a periodicity may be subject to artifacts due to seasonal gaps and low cadence, and the large timescales (defined as total time baseline <3 cycles) where the data is more subject to false periodicity from stochastic quasar variability (Vaughan et al. 2016).
tions, significant periodicity in the accretion rates emerges only for q 0.1, where the binary torques are large enough to excite eccentricity in the inner cavity and create an overdense lump. The passing BHs interact with the overdense lump, producing periodicity in the accretion rate. There is a strong peak in the periodograms corresponding to the orbital frequency of the lump, which is also the binary frequency for q=0.11 but is ∼1/5 of the binary frequency for q=0.43. The quality of the existing light curves does not justify model comparison over an even finer parameter grid in mass ratio.
One caveat is that the 2D models only predict the accretion rate and miss 3D effects and radiative transfer processes. While more realistic simulations are still needed to capture the complex physics in the binary system in order to make reliable predictions, the dominant characteristic timescale, the orbital period, and harmonics that might arise, should emerge in the light curve. The gas has to be accelerated by the binary potential, and the emission of the gas has to reflect, at some level, this behavior. Whether or not we can get an accurate estimate of the mass ratio is indeed uncertain, but circumbinary accretion variability is still preferred over relativistic Doppler boosting both for the more bursty light curve characteristic and the frequency dependent variability amplitudes as discussed further below.

Light Curve Model Fitting and Model Comparison
We have shown that a periodic model is preferred over a correlated red noise (i.e., modeled with a DRW model) based on the periodogram analysis using tailored simulations (Figure 4). As an independent analysis, here we also fit the light curve with a covariance matrix that includes a correlated red noise between measurements. It allows us to test if the data favors an additional periodic signal on top of a background of pure random, red noise variability (i.e., from stochastic quasar variability), as well as to perform a comparison between a smooth, sinusoidal model and the more bursty accretion models. The covariance matrix between measurements is given by where σ i is the 1σ measurement error at the observation time t i . The nonzero off-diagonal terms come from a correlated red noise, where σ 2 and τ are the variance and the correlation time of the variability. The null hypothesis is a flat mean amplitude with a correlated red noise, equivalent to a pure DRW model. We consider three periodic models. These include a sinusoidal model as well as two fiducial bursty accretion models, assuming mass ratios q=0.11 and 0.43. We adopt a maximum likelihood approach for the parameter fit-  ting and model comparison. The likelihood function is give by where X i is the observed flux and M i the model flux at the observation time t i . First, we test if the q=0.11 model could explain our detected periodicity by maximizing the likelihood function without any limitation on the parameters. We use the emcee package to determine the best-fit parameters and their uncertainties. We initiate 100 individual chains to sample the maximum likelihood function for 500 steps. Then, we remove the first 250 steps as a burn-in process. The 1σ error is determined by the remaining 250 steps from 100 chains at the 84.14 and 15.86 percentiles. The best-fit q=0.11 bursty model period along with the 1σ error are listed in Table  1, consistent with the periodicity found in the periodogram analysis.
Then, we compare three models (sinusoidal + red noise, circumbinary accretion + red noise, and a pure stochastic red noise) using maximum likelihood estimation. All the calculations are done in flux units. In a single-band fit, the sinusoidal model has six free parameters: red noise amplitude, red noise correlation time, period, phase, amplitude and average magnitude. The more bursty, circumbinary disk accretion variability model also has six free parameters: red noise amplitude, red noise correlation time, period, phase shift, amplitude of variation, and the magnitude zero point. A DRW model has three free parameters: red noise amplitude, red noise correlation time, and mean magnitude.
We also do a combined fit making use of the light curves from all four bands. To help break parameter degeneracy, the periodicity, phase and red noise correlation timescale are fixed to be the same across different bands. In a combined fit with the periodic models, we have fifteen free parameters, including the mean flux, model amplitude and red noise am-plitude in each band, as well as the periodicity, the phase and the red noise correlation timescale which are the same across different bands. For the pure DRW model, there are nine model parameters, including the mean flux and the red noise amplitude in each band, and a red noise correlation time which is the same across different bands.
To compare different models, we adopt the Bayesian information criterion (BIC), which is defined as where k is the number of free model parameters and N the number of data points. A lower BIC value indicates the more preferred model. We adopt the proportional constant to be unity in Eq (3). Table 2 lists the BIC from our MCMC analysis for the model fitting and comparisons. Three periodic models are compared against the null hypothesis of a pure stochastic variability. A lower BIC value indicates the more preferred model, and a BIC difference of <−10 suggests strong evidence. In each band, the q = 0.11 accretion model always has a negative BIC difference (i.e., suggesting that it is more preferred than a pure stochastic variability), which is also the smallest among all the three periodic models considered. The observed light curves of J0252 statistically prefer the q=0.11 accretion model over the other models in all four bands. Taking the fit that combines all bands for example, the BIC difference between the q=0.11 accretion model and the pure stochastic quasar variability model translates to a likelihood ratio of, at least, exp[(−96.7)/(−2)]∼10 21 (Eq (4)). We thus conclude that the q = 0.11 accretion model to be the best model for the observed light curves.
We have tested that our qualitative conclusion (that the q = 0.11 accretion model is preferred over a sinusoidal model from having smaller BIC values) still holds assuming a background of pure white (flat spectrum) noise instead (i.e., with zero off-diagonal terms in the covariance matrix). We show the best-fit q = 0.11 accretion models under white noise (thin solid curves) in Figure 1 for illustration purposes only.

Relativistic Doppler Boost Modeling
The multi-band light curves enable us to conduct an independent, quantitative test of the relativistic Doppler boost hypothesis. The relativistic Doppler boost predicts unique and robust frequency-dependent variability amplitudes in different bands that can be tested with multi-color data (D'Orazio et al. 2015b). We adopt the total mass of the hypothesized binary in J0252 as M = 10 8.4±0.1 M (1σ statistical error) assuming the virial black hole mass estimated from Mg II λ2800. We measure the spectral indices of the continuum by fitting broken power law models over four wavelength windows corresponding to the griz bands. Table 2 lists the resulting broken power-law indices.
For a binary in a relativistic circular orbit, the observed frequency of the emitted photons from the secondary's accretion disk will change due to the relativistic motion (D'Orazio et al. 2015b;Charisi et al. 2018). The number of photons N, where N ∝ F ν /ν 3 , with F ν being the flux at a specific frequency ν, is Lorentz invariant. The photons are Doppler-boosted by a factor where v is the orbital velocity, v is the line-of-sight component, and γ = 1 − β 2 −1/2 . Assuming that the emitted radiation has a power-law spectrum F ν ∝ ν α ν , the observed flux is For a binary in a circular orbit, the Doppler-boost variability to first order in β is where v is the orbital velocity of the more luminous black hole (assumed to be the less massive secondary black hole, whereas the primary black hole is assumed to contribute negligible flux), i is the inclination of the binary orbit with respect to the line-of-sight (defined such that i=90 degrees for an edge-on view and 0 degrees for a face-on view), and 0 ≤ φ ≤ 2π is the phase of the orbit. We take the orbital separation to be effectively constant over the course of the observation since the time elapsed in the rest frame is much less than the coalescence timescale of the binary. Taking the g band for example, the amplitude of the variability is 0.229 mag (Table 1), corresponding to ∆F ν /F ν ∼ ±0.229. To explain this, a line-of-sight velocity amplitude of vsin(i) ≈0.069c would be needed, considering the g-band power-law index α g ∼−0.32 (Table 1). We calculate the frequency-dependent variability amplitude ratios expected from relativistic Doppler boost (i.e., relativistic beaming, or RB for short) to compare with the observations. Taking the gr bands for example, the RB model predicts A g,RB /A r,RB = (3 − α g )/(3 − α r ) = 0.78 ± 0.11 (1σ), where α ≡ dln(F ν )/dln(ν). The observed A g,obs /A r,obs is 1.41 ± 0.03. The RB hypothesis is therefore ruled out at 5σ. Figure 7 shows the observed variability amplitude ratio (A i /A j where i and j represent two bands) compared with the expected value inferred from relativistic beaming (RB) for each band pair, which is (3 − α i )/(3 − α j ). The RB model is being ruled out at 5σ considering the gr bands and at ∼2σ for the gi and ri bands. Figure 8 shows the parameter space that allows for a flux variability greater than a fiducial value of 16%-23% in order to explain the observed values ( Table 2). The parameters considered are the total black hole mass M, mass ratio q, orbital inclination i, and the fraction of the total emission coming from the secondary black hole f 2 (D' Orazio et al. 2015b). Our other fiducial model parameters are P orb = 1.7 yr, and α = −0.32, −1.25, −0.20, and 0.41 in the griz bands (Table 2). There is little to no parameter space for the RB hypothesis to work, because the required total black hole mass would be too large to reproduce the observed, large variability amplitudes in J0252, unless all the following three requirements are met: 1. the total black hole mass is significantly underestimated by the virial estimate, even when accounting for a 0.5-dex systematic uncertainty (Shen 2013), 2. >80% of the optical light is contributed by emission from  3/6/6/6 3/6/6/6 3/6/6/6 3/6/6/6 9/15/15/15 Table 2. Measurements of J0252. Line (1): Period and 1σ error (estimated from bootstrap re-sampling) from the generalized Lomb-Scargle (GLS) periodogram. Line (2): Best-fit period and 1σ error (statistical) from MCMC fitting the q=0.11 accretion model independently in different bands assuming a correlated red noise. Lines (3)-(5): Bayesian information criterion (BIC) differences between a periodic model and the null hypothsis, i.e., stochastic quasar variability characterized by a damped random walk (DRW) model. The periodic models considered include two bursty, circumbinary accretion models assuming q=0.11 and 0.43, and a sinusoidal model (expected for relativistic Doppler boost). A negative ∆BIC indicates that a periodic model is more preferred over a pure stochastic variability. ∆BIC<−10 suggests strong evidence. Line (6): Number of data points. Line (7): Power-law index of the continuum from spectral modeling. Errors represent 1σ uncertainties generated from 100 Monte Carlo simulations. Line (8) a mini accretion disk fueling the secondary black hole, and 3. the system is viewed close to being edge-on. Our estimates on the periodic variability amplitudes (Line 8 in Table 2) do not include contribution from a stochastic background of red noise; accounting for all the observed variability amplitudes instead would make the tension even stronger.

Gravitational-wave implications and
prospects.
The GW strain amplitude of a circular binary in the quadrupolar approximation is where µ = Mq/(1 + q) 2 is the reduced mass, M is the total mass, ω = 2π f orb , and D L is the luminosity distance to the source. From our parameter estimation, the inferred strain amplitude is h 0 ∼ 9.8×10 −19 , which makes this binary effectively undetectable by current Pulsar Timing Arrays (PTAs) as an individual source (Zhu et al. 2014). Recent PTA upper limits on the stochastic background have been used to constrain the ensemble properties of BSBH candidates Holgado et al. 2018). A growing census of millipc BSBH candidates will be further constrained as the PTA sensitivity improves over time. LISA would be sensitive to BSBH mergers at these given masses and mass ratios. We estimate the SNR from the latest LISA sensitivity curve (Cornish & Robson 2018). Figure 9 shows that the BSBH candidate would eventually merge in the LISA frequency band and a merging binary with the same mass and mass ratio would be detectable during a 5-yr observation with a nominal SNR ∼15 at redshift ∼1.5.

Alternative Interpretations
Unlike the two previously best known BSBH candidates OJ287 and PG1302, J0252 is not a blazar, nor is its optical emission dominated by contribution from a radio jet, and therefore jet precession cannot explain the periodicity. Precession of a warped accretion disk is unlikely because the amount of obscured continuum emission required would be too large to explain the observed variability amplitude in J0252 and that the effect is geometrical rather than bursty. The periodicity in J0252 (i.e., rest-frame 1.7 yr) is close to the expected value (∼200 days) inferred from its black hole mass assuming a scaled-up quasar version (King et al. 2013) of low-frequency accretion disk quasi-periodic oscillations (QPOs; e.g., from strong resonances in the accretion flow) as seen in the X-ray light curves of X-ray binaries (Vaughan & Uttley 2005). However, the characteristic bursty light curves would be difficult to explain with Lense-Thirring precession (Bardeen & Petterson 1975) of a geometrically thick accretion flow near the primary black hole, with which low-frequency QPOs in X-ray binaries are as- . Prospect for LISA detection of a source similar to the candidate BSBH in J0252 but 5 years before coalescence. The purple curve represents the expected LISA sensitivity limit assuming a 5-yr observation (Cornish & Robson 2018). The black curve denotes the gravitational-wave signal of a BSBH at z = 1.53 with mass 10 8.4 M and mass ratio q = 0.1 beginning at 5 yrs before coalescence, i.e., from the inspiral phase (low frequency) to the final merger and ringdown (high frequency). The blue, orange, green, and red shaded regions correspond to mergers with a primary mass of 10 6 M , 10 7 M , 10 8 M , and 10 9 M , respectively, at the same redshift with mass ratios ranging within 0.05 ≤ q ≤ 0.5. The blue, orange, green, and red lines correspond to q = 0.1.
sociated (Stella & Vietri 1998;Ingram et al. 2009). QPOs in X-ray binaries show drifts in period, phase, or amplitude (van der Klis 1989), and future continued monitoring observations can constrain the possibility of an optical QPO.

Implication for the Merger Hypothesis and Comparison with Theoretical Event Rates
Among a parent sample of 625 quasars we have detected one strong candidate BSBH, whose estimated gravitational-wave inspiral time (0.17 Myr) is about ∼10 2 -10 3 times shorter than estimates for quasar lifetimes (Martini & Weinberg 2001;Yu & Tremaine 2002). This implies that most quasars could be binary systems with a much larger binary separation that the circumbinary disk does not yet exist, which is unsurprising given the merger hypothesis (Volonteri et al. 2003;Hopkins et al. 2008;Shen 2009;Haiman et al. 2009a). Previous work has predicted the event rates of BSBHs that are detectable as periodic quasars (Haiman et al. 2009b). The most recent work by Kelley et al. (2019) combines cosmological hydrodynamic simulations, semi-analytic binary merger models, and analytic quasar spectra and variability prescriptions. Given DES sensitivity (assuming a typical single epoch 5σ point source depth of ∼23.5 AB mag), the expected number of detectable periodic quasars from circumbinary disk accretion variability at redshift z∼1.5 with observer-frame periods between 0.5 and 5.0 yr is ∼80 in an all-sky survey (∼30,000 deg 2 ), or ∼1 per 380 deg 2 (see right panel in Figure 6 of Kelley et al. 2019). This is ∼70 times lower than our detection rate 4 of ∼one strong BSBH candidate from circumbinary accretion variability per 5 deg 2 at face value. As further discussed by Chen et al. (2020), this apparent discrepancy is likely explained by the fact that our sample is dominated by less massive quasars at high redshift given our deep survey over a small area. As a result we are effectively measuring the differential detection rate (which is a function of redshift and BH mass) rather than the cumulative detection rate as quoted by Kelley et al. (2019). There are still significant uncertainties that prevent a fair comparison between our detection rate and theoretical predictions and PTA limits. First, theoretical event rates are still highly uncertain. The most significant uncertainty is on the inspiraling timescales, which could lead to highly uncertain estimates on the number of detectable binaries in the circumbinary accretion disk phase.
Second, the PTA upper limits are still subject to model uncertainties regarding the evolutionary history of a binary from large to small separations where GW emission dominates. PTA upper limit is model independent only for a particular binary separation range that corresponds to the PTA frequency. To extrapolate this to other separations (i.e., going from PTA frequency to the frequency relevant for periodic quasars), one needs to invoke assumptions on the evolutionary timescales. However, there is still no self-consistent model that can deal with the full evolution considering the effects of gas and stars, and so a high binary fraction at mili-parsec scales may not necessarily be in direct tension with the PTA upper limits. For example, if a binary stalls at large separations, or if it sweeps quickly through the PTA sensitivity range, there would be no PTA signal, even if the true binary fraction were high.
Finally, PTA is most sensitive to the most massive binaries at low redshift. However, our sample is most sensitive to the ∼ 10 8 M systems at intermediate redshift. A small binary fraction for the most massive black holes at z=0 from the PTA upper limits does not directly translate into the same binary fraction for the less massive black holes at z ∼ 1.5.

Comparison with Previous Work
As further shown in Chen et al. (2020), our detection rate of all candidate periodic quasars (not just BSBH candidates from circumbinary accretion variability), i.e., ∼0.8%, is ∼4-80 times of those from previous searches using other surveys (Graham et al. 2015;Charisi et al. 2016;Liu et al. 2019), even though our selection criteria are more stringent. For example, we request > three cycles in at least two bands, whereas only 1.5 cycles were adopted by Graham et al. (2015) and only one-band data were available. This is not a fair comparison, however, because our sample is probing less massive quasars at higher redshifts than those studies in previous shallower surveys over larger areas. As suggested in Chen et al. (2020), the significantly higher detection rate of periodic quasars found in our sample may be interpreted as the redshift evolution of the fraction of BSBHs, i.e., the binary fraction is larger at higher redshifts at a fixed BH mass.
In addition, previous datasets lacked the long time baseline and/or sensitivity to discover similar systems as J0252. Given shorter time baselines and lower sensitivities, false positives and/or false negatives would have been more likely to significantly bias the apparent detection rates because of stochastic background variability. In particular, Liu et al. (2019) has rejected most of the candidates found in their previous searches (Liu et al. 2015(Liu et al. , 2016 by continued monitoring of the "best candidates". While this demonstrates the importance of a long time baseline in rejecting false positives due to stochastic background variability, it does not address the question of possibly missing false negatives in those that have not been continuously monitored. A long time baseline for the full parent sample (i.e., not just the "best candidates" selected based on short-baseline light curves) is needed to robustly quantify the true detection rate.
In summary, the quality of the data (i.e., long time baseline, high sensitivity) is more important than the quantity of the data (i.e., size of the parent quasar sample) because the systematic error (e.g., bias from false positives and/or false negatives caused by stochastic background variability) is likely to be larger than the statistical error. Even though we have a much smaller sample of quasars in the parent sample, our detection rate is still likely to be more reliable than those from previous work based on shorter and shallower surveys of larger areas.

CONCLUSION AND FUTURE WORK
Our results on J0252 may provide the first, strong evidence for circumbinary accretion variability as the physical origin for periodic quasar (optical) light curves. Sensitive, longterm, multi-color light curves are key in disfavoring the competing relativistic-Doppler-boost hypothesis for J0252. Relativistic Doppler boost has been previously shown to best explain the characteristic periodic optical light curves and UV observations of PG1302−102 (D'Orazio et al. 2015b). We speculate that various mechanisms may be at work in different systems, such that the case for PG1302−102 may not necessarily apply to J0252 or other periodic quasars.
Recently, using a combination of cosmological, hydrodynamic simulations, comprehensive semi-analytic binary merger models, and analytic active galactic nucleus spectra and variability prescriptions, Kelley et al. (2019) suggests that hydrodynamic variability should be ∼5-25 times more common than relativistic Doppler boost in producing periodic quasar light curves in synoptic surveys. Our result suggests that hydrodynamic circumbinary accretion variability may indeed be a viable option to explain periodic light curves, at least for some, if not most, quasars as BSBH candidates, although we cannot draw large inferences from just a single detection. Alternatively, precession of a radio jet is likely ruled out, because unlike OJ287 (Valtonen et al. 2008), J0252 is not a blazar (with a 3σ radio flux density upper limit of <0.5 mJy at 1.4 GHz and <0.4 mJy at 3 GHz), nor is its optical emission dominated by contribution from a radio jet.
While we have adopted the simulated light curves of Farris et al. (2014) as the baseline model, our conclusion is not sensitive to this particular choice because similar characteristic bursty light curves are seen in other independent simulations of circumbinary accretion disks around BSBHs (MacFadyen & Milosavljević 2008;Shi et al. 2012;Roedig et al. 2012;D'Orazio et al. 2013;Gold et al. 2014a;Tang et al. 2018). While the archival SDSS data has been necessary in extending the time baseline for a statistically significant periodicity detection, the light curve was only well sampled by the new DES observations in terms of sensitivity and cadence, and there were significant observational gaps. The existing data cannot definitively discriminate between the q=0.11 and q=0.43 circumbinary accretion variability models, although q=0.11 is tentatively preferred ( Table 2). We consider these two q values as baseline examples because they represent two characteristic regimes in the lightcurve behaviors (Figure 9 of Farris et al. (2014)). In both regimes there is a strong peak in the periodograms of the simulation-predicted light curves corresponding to the orbital frequency of the overdense lump. Adopting a mass ratio of q=0.11, t orb = t period (whereas t orb ≈ 0.2t period for q=0.43 instead (Farris et al. 2014)), the inferred binary separation is d ∼ 4.4 milli-parsec (i.e., 5.1 light days, or ∼200 Schwarzschild radii), assuming a circular orbit. So the confirmation of this candidate would imply that the system has passed the "final-parsec" barrier (Begelman et al. 1980) at a redshift of z=1.53.
The inferred gravitational-wave inspiral time t gw with the preferred system parameters is ∼0.17 Myr. This implies that the candidate binary is efficiently emitting gravitational waves and will merge well within the age of the universe, even if environmental effects are neglected. BS-BHs with masses of ∼ 10 8 -10 9 M at redshift z 1 are generally expected around the time of pre-decoupling (Kocsis & Sesana 2011), i.e., when t gw > t visc , where t visc is the viscous timescale of the accretion disk. The gravitationalwave strain amplitude is ∼10 −18 at ∼37 nHz, which, as an individual source, is ∼10 5 below the current best sensitivity limit of pulsar timing arrays (Arzoumanian et al. 2018b) to continuous-wave sources, and will also be below the expected SKA sensitivity (Wang & Mohanty 2017). Laser Interferometer Space Antenna (Amaro-Seoane et al. 2017) would be able to detect a source similar to J0252 but ∼5 years before coalescence at 0.01 mHz with a signal-to-noise ratio of ∼15 at redshift 1.5 (Figure 9). Future sensitive, continued multi-band follow-up imaging is needed to further constrain the significance and nature of the optical light-curve periodicity observed in J0252. While the existing data spans 4.6 cycles, only ∼3 are well sampled in multiple bands. There is a 0.05% probability that the periodogram peak is caused by stochastic quasar variability. The significance of a real periodicity should increase as more cycles are covered . Continuous, sensitive follow up with the Blanco 4m/DECam is on-going to better characterize the light curve properties. Hydrodynamic simulations of circumbinary accretion disks predict additional, weaker peaks in the light curve periodograms at different characteristic frequencies depending on the mass ratio, with many associated harnomics for q>0.43 (Farris et al. 2014). Future observations may be able to better distinguish between the q=0.11 and the q=0.43 models (e.g., by searching for evidence for additional weaker peaks in the pe-riodogram and quantifying their characteristic relationships with the primary peak (Charisi et al. 2015)).
The observed SED of J0252 is similar to normal optically selected quasars that are matched in redshift and luminosity. Future more sensitive UV and/or X-ray observations are needed to put further independent constraints on any potentially characteristic SED features to compare with predictions from circumbinary accretion disk simulations (Roedig et al. 2012;Tang et al. 2018). While the broad-line region is expected to be well outside the radius of the binary, the circular velocity is about 0.05c, which is much greater than the width of the broad emission lines. Any emission lines originating from the disk could in principle show such shifts, but in practice the broad emission line profile becomes more complex and there are no expected coherent radial velocity drifts in the emission lines with time (Shen & Loeb 2010). There could be a shift in the Fe K-α line which probes the inner accretion disk and future sensitive X-ray spectroscopic monitoring is needed to test this (McKernan & Ford 2015).
Our detection of one strong BSBH candidate due to circumbinary accretion variability in a sample of 625 spectroscopically confirmed quasars from a 4.6 deg 2 survey implies a detection rate of ∼0.16%, or 1 per 5 deg 2 , which is ∼70 times higher than the expected event rate (Kelley et al. 2019) at face value, although the theoretical rate is still highly uncertain considering unconstrained model assumptions. Our detection rate of candidate periodic quasars is ∼4-80 times times of those from previous searches using other surveys (Liu et al. 2019), although this is not a fair comparison because previous datasets lacked the long time baseline and/or sensitivity to discover similar systems as J0252. Given shorter time baselines and lower sensitivities, false positives and/or false negatives would have been more likely to significantly bias the apparent detection rates because of stochastic background variability. We have demonstrated using J0252 that multi-band light curves with high sensitivity and a long time baseline is key to not only identifying periodicity but also sorting out its physical origin. Future large, sensitive synoptic surveys such as the Legacy Survey of Space and Time at Vera C. Rubin Observatory (Ivezić et al. 2019) may be able to detect hundreds to thousands of BSBH candidates from circumbinary accretion variability. and NSF grant AST-1715579. This work makes use of observations from the LCOGT network.
Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparoà Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.
The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.
Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
The DES data management system is supported by the National Science Foundation Ofek et al. 2012). We have further applied the SDSS-DES corrections listed in Equation 1 for the PTF photometry to be on the DES system. The PTF-to-SDSS correction depends on the (r-i) and (g-r) colors which are variable, however, on the timescales of a few yr. We have adopted the median colors averaged in the last year of the SDSS light curves and the first year of the DES observations that bracketed the PTF R-band observations. Furthermore, the current version of the PTF photometric pipeline uses MAG AUTO (not aperture or PSF magnitudes) which adjusts the aperture used to extract the source magnitude for each object. This introduces biases in the magnitudes for sources near the survey detection limit such as J0252.
The resulting effect on the color correction is a systematic bias toward larger negative values of r SDSS -R PTF/SDSS starting around r SDSS of magnitude 19.5 (Figure 2 of Ofek et al. 2012). We have empirically corrected for this systematic bias using the median value inferred for sources with similar luminosities of J0252 (i.e., at r SDSS ∼21 mag). We have further verified the empirical correction by comparing the four PTF R-band data points that overlapped with the DES Y1 observations (i.e., around MJD of 56,600), finding a general consistency. Nevertheless, given these significant uncertainties and caveats in the magnitude conversion, as well as the fact that J0252 is already at the PTF survey detection limit, we do not include the PTF photometry in our baseline analysis. The PS1 griz filters are similar to those of the SDSS. We apply the PS1-to-SDSS correction using a third order polynomial provided by Finkbeiner et al. (2016) that shifts the photometry to the SDSS system. The correction depends on the (g-i) color. The color is determined by averaging over PS1 light curve. After the correction to the SDSS system, Eq (1) is then applied for the conversion between the SDSS and DES systems. For ZTF, the photometry has been calibrated to the PS1 system (Masci et al. 2019). We thus follow the same steps in the PS correction and correct the light curves to be on the DES system. The LCOGT filters are similar to the SDSS. We convolve each quasar spectrum with the DES and LCOGT filter transmission curves to calculate the synthetic magnitude difference and correct the LCOGT to be on the DES system.  Table A3. Photometry data of J0252 in i-band. SDSS photometry data has been corrected, based on Eq 1, to match the DES system .  mjd  mag  err  survey  mjd  mag  err  survey  mjd  mag  err  survey  mjd  mag  err  survey  51075