# Discovery of a Candidate Binary Supermassive Black Hole in a Periodic Quasar
from Circumbinary Accretion Variability

Wei-Ting Liao,1,2 Yu-Ching Chen,1,2 Xin Liu,1,2 A. Miguel Holgado,1,2 Hengxiao
Guo,1,2 Robert Gruendl,1,2 Eric Morganson,2 Yue Shen,1,2 Tamara Davis,3
Richard Kessler,4,5 Paul Martini,6,7 Richard G. McMahon,8,9 Sahar Allam,10
James Annis,10 Santiago Avila,11 Manda Banerji,8,9 Keith Bechtol,12 Emmanuel
Bertin,13,14 David Brooks,15 Elizabeth Buckley-Geer,10 Aurelio Carnero
Rosell,16,17 Matias Carrasco Kind1,2, Jorge Carretero18, Francisco Javier
Castander19,20, Carlos Cunha21, Chris D’Andrea22, Luiz da Costa17,23,
Christopher Davis21, Juan De Vicente16, Shantanu Desai24, H. Thomas Diehl10,
Peter Doel15, Tim Eifler25,26, August Evrard27,28, Brenna Flaugher10, Pablo
Fosalba19,20, Josh Frieman10,4, Juan Garcia-Bellido29, Enrique Gaztanaga19,20,
Karl Glazebrook30, Daniel Gruen21,31, Julia Gschwend17,23, Gaston Gutierrez10,
Will Hartley15,32, Devon L. Hollowood33, Klaus Honscheid7,34, Ben Hoyle35,36,
David James37, Elisabeth Krause25, Kyler Kuehn38, Marcos Lima39,17, Marcio
Maia17,23, Jennifer Marshall40, Felipe Menanteau1,2, Ramon Miquel41,18, Andrés
Plazas Malagón26, Aaron Roodman21,30, Eusebio Sanchez16, Vic Scarpine10,
Michael Schubnell28, Santiago Serrano19,20, Mathew Smith42, R. Chris Smith43,
Marcelle Soares-Santos44, Flavia Sobreira45,17, Eric Suchyta46, Molly
Swanson2, Gregory Tarle28, Vinu Vikram47, Alistair Walker43, the DES
Collaboration E-mail: wliao10<EMAIL_ADDRESS>P.
Sloan Research Fellow

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Binary supermassive black holes (BSBHs) are expected to be a generic byproduct
from hierarchical galaxy formation. The final coalescence of BSBHs is thought
to be the loudest gravitational wave (GW) siren, 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$\pm 7$ days) at
99.95% significance (with a global p-value of ${\sim}10^{-3}$ accounting for
the look elsewhere effect) in the optical light curves of a redshift 1.53
quasar, SDSS J025214.67$-$002813.7. Combining archival Sloan Digital Sky
Survey data with new, sensitive imaging from the Dark Energy Survey, the total
$\sim$20-yr time baseline spans $\sim$4.6 cycles of the observed 4.4-yr
(restframe 1.7-yr) periodicity. The light curves are best fit by a bursty
model predicted by hydrodynamic simulations of circumbinary accretion disks.
The periodicity is likely caused by accretion rate modulation by a milli-
parsec BSBH emitting GWs, dynamically coupled to the circumbinary accretion
disk. A 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, 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 Vera C. Rubin Observatory Legacy Survey of Space
and Time may be able to detect hundreds to thousands of candidate BSBHs from
circumbinary accretion with direct implications for Laser Interferometer Space
Antenna.

###### keywords:

black hole physics – galaxies: active – galaxies: high-redshift – galaxies:
nuclei – quasars: general – surveys

††pubyear: 2020††pagerange: Discovery of a Candidate Binary Supermassive Black
Hole in a Periodic Quasar from Circumbinary Accretion Variability–affliations

##  1 Introduction

LIGO has detected gravitational waves (GWs) from stellar-mass binary black
hole mergers (Abbott et al., 2016), yet many GW sources are expected outside
the LIGO frequency (Sesana, 2017; Schutz, 2018). A binary supermassive black
hole (BSBH) consists of two black holes, each with a mass of $\sim
10^{6}$–$10^{9}$ M⊙. BSBHs are expected to frequently form in galaxy mergers
(Begelman et al., 1980; Haehnelt & Kauffmann, 2002; Volonteri et al., 2003),
given that most massive galaxies harbor SMBHs (Kormendy & Richstone, 1995;
Ferrarese & Ford, 2005). Their final coalescences should produce the loudest
GW sirens in the universe (Thorne & Braginskii, 1976; Haehnelt, 1994; Vecchio,
1997; Jaffe & Backer, 2003), which will be the primary source of low-frequency
GW experiments (Amaro-Seoane et al., 2017; Arzoumanian et al., 2018a; Sesana
et al., 2018). BSBHs are important for testing general relativity in the
strong field regime and for the studies of galaxy evolution and cosmology
(Centrella et al., 2010; Merritt, 2013; Colpi, 2014; Berti et al., 2015).

However, no confirmed case is known at sub-milliparsec scales, i.e.,
separations close enough to be in the GW-dominated regime. While ${\sim}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., 2020). 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., 2020), 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 (1998–2007; Ivezić et
al., 2007). Unlike previous studies, which were based on few-cycle (e.g.,
$\sim$1.5) searches given the limited time baselines, the periodicity of J0252
was discovered based on $\sim$5 cycles enabled by a $\sim$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; Barth & Stern, 2018). 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., 2020). 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 $\Lambda$CDM cosmology with
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=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.

##  2 Data and Methods

###  2.1 Sample Selection

We start with 763 spectroscopically confirmed quasars in the 4.6 deg2
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
(Schneider et al., 2007; 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$\sigma$ 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 $\sim$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.

###  2.2 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 $\sim$20 yr (1998–2007 from SDSS
Stripe 82 and 2012–2019 from DES-SN). For a typical quasar at $z\sim$1, the
time baseline spans $\sim 10$ yr in the quasar rest-frame to encompass
$\gtrsim$5 cycles for a period of $\lesssim$2 yr, which is the recommended
number of cycles to minimize false periodicity (Vaughan et al., 2016). We have
rejected $>5\sigma$ 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 deg2 survey of the Southern Hemisphere in the
$grizY$ bands (Flaugher, 2005; The Dark Energy Survey Collaboration, 2005;
Dark Energy Survey Collaboration et al., 2016). It uses the Dark Energy Camera
(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$\sigma$ 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 deg2 multi-epoch survey (DES-SN) to search for
SNe Ia that has a mean cadence of $\sim$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 deg2). 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 $\sim$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$\sigma$ 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

$\begin{split}g_{{\rm DES}}&=g_{{\rm SDSS}}-0.000\pm 0.002\\\ r_{{\rm
DES}}&=r_{{\rm SDSS}}-0.116\pm 0.005\\\ i_{{\rm DES}}&=i_{{\rm SDSS}}+0.053\pm
0.009\\\ z_{{\rm DES}}&=z_{{\rm SDSS}}+0.022\pm 0.016,\end{split}$ (1)

where the errors are 1$\sigma$ 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 $\sim$20-year 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 provided in the supplementary online data.

Figure 1:  SDSS and DES optical light curves of J0252. All observations have
been corrected to be on the DES system. Also shown are archival light curves
from the PTF, CRTS, ZTF, PS1, and new observations from the LCOGT. Error bars
represent 1$\sigma$ (statistical). 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.

###  2.3 Periodicity Detection

For any periodicity detection, we implement the following three selection
criteria:

  1. 1.

At least two bands have a 3$\sigma$ detection of the same periodicity in the
periodogram analysis.

  2. 2.

The detected periodicity is the dominant component compared to the background
noise.

  3. 3.

The same periodicity is also identified in the auto-correlation 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$\sigma$ 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
$\sigma_{\rm residue}^{2}/A_{\rm sin}^{2}>1$, where $A_{\rm sin}$ is the
periodicity amplitude and $\sigma_{\rm residue}^{2}$ is the variance of the
residue light curve after subtracting the periodic 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{(\omega
t)}\exp{(-\lambda t)}$, where $t$ is the lagging time,
$\omega=\frac{2\pi}{P}$, and $\lambda$ 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 multi-band 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$\sigma$ detections in these alternative
methods.

###  2.4 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 $\gtrsim$10 days
(Kelly et al., 2009; MacLeod et al., 2010; Mushotzky et al., 2011; Kozłowski,
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: $\sigma^{2}$ and $\tau$, 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 $\sigma$ and $\tau$, 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 ($\sim$7300 days) is more than 10 times larger than
the correlation timescale ($\sim$630 days in the observed frame), so that the
correlation timescale recovers the true value (Kozłowski, 2017). We then
generate 50,000 mock light curves with parameters pairs ($\sigma^{2}$ and
$\tau$) randomly drawn from the posterior distribution in the DRW parameter
fitting. We 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\gtrsim\frac{1}{10}$ day-1 (Mushotzky et
al., 2011; Edelson et al., 2013, 2014; Smith et al., 2018). For the BPL model,
we assume a $-3$ power spectrum index at the high frequency $f>\frac{1}{10}$
day-1 and keep the DRW model at the low frequency $f<\frac{1}{10}$ day-1. We
have tested different power law indices and different high frequency breaks.
Our result is insensitive to these choices. The low-frequency breaks are drawn
from the timescales $\tau$ in the DRW parameter fitting. We first generate
50,000 evenly sampled mocked light curves assuming the BPL model with the
pyLCSIM111http://pabell.github.io/pylcsim/html/index.html 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 the DRW model, we have also considered the CAR(2,1) model
(Kelly et al., 2014), i.e., a damped harmonic oscillator. CAR(2,1) 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\ll 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\%$ detection holds in the $g$ band. Table 1 shows that using carma-
pack222https://github.com/brandonckelly/carma_pack, CAR(1,0) has the lowest
BIC value and is thus a proper noise model.

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\leq 3$ and
$q<p$. CAR(1,0) has the smallest BIC value, suggesting that DRW is the proper
noise model for J0252.

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 $\sigma$ and
$\tau$ (observed frame). Labeled on each panel are their best-fit value and
the 1-$\sigma$ (estimated from the 68% confidence levels denoted by the shaded
histograms) uncertainties. The total light curve baseline ($\sim$7300 days) is
more than 10 times larger than the correlation timescale ($\sim$630 days in
observed frame), so that the correlation timescale recovers the true value
(Kozłowski, 2017).

###  2.5 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 $\lambda$1549, He II $\lambda$1640, C III]
$\lambda$1909, and Mg II $\lambda$2800.

Figure 3: Optical spectrum and modeling of J0252. Shown are the data (black),
the 1$\sigma$ 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).

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
$\lambda$2800 region locally. The Mg II $\lambda$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.

##  3 Results

###  3.1 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 deg2
field. J0252 was the most significant detection with $>$4 cycles spanned whose
light curves prefer a bursty circumbinary accretion model. We focus on J0252
in this paper, whereas the other four candidates are presented in Chen et al.
(2020).

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% single-peak significance in the $griz$ bands. The empirically estimated
global p-values (Barth & Stern, 2018) are 1.2$\times 10^{-3}$, 5.8$\times
10^{-3}$, 2.6$\times 10^{-3}$, and 8.4$\times 10^{-3}$, accounting for the
look elsewhere effect (e.g., Gross & Vitells, 2010) in the whole frequency
range being searched (see Chen et al., 2020, for details). 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 $\sim$0.1%
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% single-peak significance in a parent sample
of 625 (in which $\lesssim$two cases are expected from red noise; Chen et al.
2020) suggests that we are not just seeing stochastic quasar variability in
our small sample. Similar to Figure 4, Figure 5 shows the significance level
assuming a CAR(2,1) noise. The candidate periodicity is found at $99.8$%,
$98.8$%, $99.5$% and $98.0$% single-peak significance in the $griz$ bands
under the CAR(2,1) assumption for the stochastic noise. For context, the false
alarm probability of seeing such a significant peak in the periodograms is
${\ll}10^{-20}$ assuming a pure white (i.e., flat spectrum) noise instead
(Zechmeister & Kürster, 2009). Chen et al. (2020) has estimated the DRW
parameter distributions for the parent quasar sample to address the false
alarm probability for the candidates as a population.

Figure 4: Generalized Lomb-Scargle periodogram showing the periodicity
detection of J0252. A periodicity (see $P_{{\rm GLS}}$ in Table 2) is detected
at 99.95%, 99.43%, 99.78%, and 99.59% single-peak significance (with global
p-values of 1.2$\times 10^{-3}$, 5.8$\times 10^{-3}$, 2.6$\times 10^{-3}$, and
8.4$\times 10^{-3}$ accounting for the look elsewhere effect in the whole
frequency range being searched; see Chen et al. 2020 for details) 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).

Figure 5: Similar to Figure 4, but with a CAR(2,1) noise model. The cyan
shaded regions are the same as those shown in Figure 4. The periodicity is
detected at single-peak significance of $99.8$%, $98.8$%, $99.5$% and $98.0$%
in the $griz$ bands.

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 ${\sim}4.6$ cycles of the
periodicity, approaching the number of observed cycles recommended for
minimizing false positives from stochastic quasar variability (e.g., Vaughan
et al., 2016).

###  3.2 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 $\lambda$2800-based estimator
gives a virial black hole mass of $M=10^{8.4\pm 0.1}M_{\odot}$ (1$\sigma$
statistical error), using the parameters in Vestergaard & Osmer (2009). Shen
(2013) suggests that Mg II $\lambda$2800-based masses are more reliable than C
IV $\lambda$1549-based masses, given that C IV $\lambda$1549 is likely to
suffer from non-virial motion like outflows and there is larger scatter
between C IV $\lambda$1549 and H$\beta$ masses for quasars at high redshift
(Shen & Liu, 2012).

###  3.3 Radio Loudness Upper Limit

J0252 was undetected by FIRST (Becker et al., 1995) with a 3$\sigma$ 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 image333http://archive-
new.nrao.edu/vlass/HiPS/VLASS1.1/Quicklook/, suggesting a 3$\sigma$ upper
limit of $f^{{\rm obs}}_{{3\,{\rm GHz}}}<0.36$ mJy. Assuming that the radio
flux follows a power law $f_{\nu}\propto\nu^{\alpha}$, this translates into
$f^{{\rm rest}}_{{6\,{\rm cm}}}<{0.18}$ mJy (6 cm corresponding to 5 GHz) for
a spectral index $\alpha=-0.5$ (Jiang et al., 2007), or $f^{{\rm
rest}}_{{6\,{\rm cm}}}<{0.20}$ mJy assuming $\alpha=-0.8$ (Gibson et al.,
2008). Combining the $f_{{\rm 2500}}$ measurement from the optical spectrum,
the inferred radio loudness parameter (Kellermann et al., 1989) is $R\equiv
f_{{6\,{\rm cm}}/f_{{\rm 2500}}}<{34}$ assuming $\alpha=-0.5$, or $R<{39}$
assuming $\alpha=-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.4 Spectral Energy Distribution

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).

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.

Figure 6: Spectral energy distribution (SED) of J0252. Also shown for
comparison are the mean, and 1$\sigma$, and 2$\sigma$ confidence levels of the
SEDs of a control quasar sample matched in redshift and luminosity with J0252,
the optically selected quasar SEDs from (Richards et al., 2006) (“R06-All” for
all quasars, “R06-OL” for optically luminous quasars, and “R06-OD” for
optically dim quasars), and the mean SED of (Hatziminaoglou et al., 2005)
(Hat05). Errorbars are 1$\sigma$ whereas upper limits are 3$\sigma$. Plotted
on top are the multi-wavelength postage stamps of J0252 with a FOV of 30′′
each. The green circles are 10′′ in diameter indicating the position of J0252.

##  4 Discussion

###  4.1 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
simulations, significant periodicity in the accretion rates emerges only for
$q\gtrsim$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
$\sim$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.

###  4.2 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

$C_{ij}={\sigma_{i}}^{2}\delta_{ij}+\sigma^{2}\exp{\bigg{[}\frac{-|t_{i}-t_{j}|}{\tau}\bigg{]}}\,,$
(2)

where $\sigma_{i}$ is the 1$\sigma$ measurement error at the observation time
$t_{i}$. The nonzero off-diagonal terms come from a correlated red noise,
where $\sigma^{2}$ and $\tau$ 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 fitting and model comparison. The likelihood
function is give by

$L\propto\det|C|^{-\frac{1}{2}}\exp\bigg{[}-\frac{1}{2}(X_{i}-M_{i})\left(C^{-1}\right)_{ij}(X_{j}-M_{j})\bigg{]}\,,$
(3)

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$\sigma$ 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$\sigma$ error are listed in Table 2, 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 amplitude 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

${\rm BIC}=-2\ln(L)+k\ln(N)\,,$ (4)

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).

Parameter | g | r | i | z | griz
---|---|---|---|---|---
$P_{{\rm GLS}}$ (days)(1) | 1607$\pm$7 | 1615$\pm$9 | 1632$\pm$8 | 1607$\pm$10 | –
$P_{{\rm Acc,\,q=0.11}}$ (days)(2) | 1511${}^{+34}_{-55}$ | 1466${}^{+64}_{-12}$ | 1506${}^{+128}_{-61}$ | 1562${}^{+248}_{-99}$ | 1476${}^{+128}_{-5}$
BIC${}_{{\rm Acc,\,q=0.11}}$$-$BICDRW(3) |  $-$19.7 |  $-$23.7 |  $-$13.3 |  $-$13.7 |  $-$96.7
BIC${}_{{\rm Acc,\,q=0.43}}$$-$BICDRW(4) |  $-$9.7 |  $-$7.0 |  $-$2.5 |  $-$1.9 |  $-$78.3
BICsin$-$BICDRW(5) |  $-$5.5 |  $-$5.7 |  $+$0.9 |  $-$2.9 |  $-$47.2
$N$ (6) | 212 | 223 | 222 | 227 | 884
$\alpha\equiv dln(F_{\nu})/dln(\nu)$(7) |  $-$0.32$\pm$0.40 |  $-$1.25$\pm$0.34 |  $-$0.20$\pm$0.34 | 0.41$\pm$1.01 | –
$A_{{\rm obs}}$ (mag)(8) | 0.229$\pm$0.003 | 0.162$\pm$0.002 | 0.162$\pm$0.002 | 0.157$\pm$0.004 | –
$\tau_{\rm DRW}$ (days) (9) | 653 | 716 | 629 | 849 | 701
$k_{\rm DRW}$ / $k_{\rm q=0.11}$ / $k_{\rm q=0.43}$ / $k_{\rm sin}$ (10) | 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$\sigma$ error
(estimated from bootstrap re-sampling) from the generalized Lomb-Scargle (GLS)
periodogram. Line (2): Best-fit period and 1$\sigma$ 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 $\Delta$BIC indicates that a
periodic model is more preferred over a pure stochastic variability.
$\Delta$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$\sigma$ uncertainties generated from 100 Monte Carlo simulations.
Line (8): Variability amplitude from the best-fit sinusoidal model. Errors
represent 1$\sigma$ statistical uncertainties. Line (9): Best-fit correlation
time in the DRW model. Line (10): Number of free parameters for each of the
model.

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)]${\sim}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. We have also tested an eccentric
Doppler boost model, but the $q=0.11$ circumbinary accretion model still has
the lowest BIC.

###  4.3 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\pm 0.1}M_{\odot}$ (1$\sigma$ statistical error) assuming the virial
black hole mass estimated from Mg II $\lambda$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\propto{\cal F}_{\nu}/\nu^{3}$, with ${\cal F}_{\nu}$
being the flux at a specific frequency $\nu$, is Lorentz invariant. The
photons are Doppler-boosted by a factor

${\cal D}=\frac{1}{\gamma\left(1-v_{\parallel}/c\right)}\ ,$ (5)

where $v$ is the orbital velocity, $v_{\parallel}$ is the line-of-sight
component, and $\gamma=\left(1-\beta^{2}\right)^{-1/2}$. Assuming that the
emitted radiation has a power-law spectrum ${\cal
F}_{\nu}\propto\nu^{\alpha_{\nu}}$, the observed flux is

${\cal F}_{\nu}^{\rm obs}={\cal D}^{3-\alpha_{\nu}}{\cal F}_{\nu}^{\rm em}\ .$
(6)

For a binary in a circular orbit, the Doppler-boost variability to first order
in $\beta$ is

$\frac{\Delta{\cal F}_{\nu}}{{\cal
F}_{\nu}}=\left(3-\alpha_{\nu}\right)\beta\cos\phi\sin i\ ,$ (7)

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\leq\phi\leq 2\pi$ 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 2),
corresponding to $\Delta F_{\nu}/F_{\nu}\sim\pm$0.229. To explain this, a
line-of-sight velocity amplitude of $vsin(i)\approx$0.069$c$ would be needed,
considering the $g$-band power-law index $\alpha_{g}{\sim}$$-$0.32 (Table 2).

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,{\rm RB}}/A_{r,{\rm
RB}}=(3-\alpha_{g})/(3-\alpha_{r})=0.78\pm 0.11$ (1$\sigma$), where
$\alpha\equiv dln(F_{\nu})/dln(\nu)$. The observed $A_{g,{\rm obs}}/A_{r,{\rm
obs}}$ is $1.41\pm 0.03$. The RB hypothesis is therefore ruled out at
$\gtrsim$5$\sigma$.

Figure 7: Observed frequency dependent variability amplitude ratio for each
band combination compared with the expected values from relativistic Doppler
boost. The black line represents the 1 to 1 relation. Error bars denote
1$\sigma$ uncertainties.

Figure 8: Parameter space estimates for the relativistic Doppler boost model.
The four panels represent $griz$ bands. In each panel, the dashed contours
represent $f_{2}=1.0$ whereas the shaded contours denote $f_{2}=0.8$, where
$f_{2}$ is the fraction of the total emission coming from the secondary black
hole. Different colors show different mass ratios with $q$=0.0, 0.05, 0.1, and
0.2 for blue, green, orange, and red, respectively. The vertical solid line
with gray shades indicate our virial mass estimate and its 1$\sigma$
statistical error for the total black hole mass. The orbital inclination angle
$i$=90 degrees for an edge-on view, and 0 degree for a face-on view.

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-\alpha_{i})/(3-\alpha_{j})$. The RB model is being ruled out at
$\gtrsim$5$\sigma$ considering the $gr$ bands and at $\sim$2$\sigma$ 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_{{\rm orb}}$ = 1.7 yr, and $\alpha$ =
$-$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 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.

###  4.4 Gravitational-wave implications and prospects.

The GW strain amplitude of a circular binary in the quadrupolar approximation
is

$h_{0}=\frac{4G^{5/3}}{c^{4}}\frac{\mu M^{2/3}\omega^{2/3}}{D_{\rm L}}\ ,$ (8)

where $\mu=Mq/(1+q)^{2}$ is the reduced mass, $M$ is the total mass,
$\omega=2\pi f_{\rm orb}$, and $D_{\rm L}$ is the luminosity distance to the
source. From our parameter estimation, the inferred strain amplitude is
$h_{0}\sim 9.8{\times}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 (Sesana et
al., 2018; Holgado et al., 2018). A growing census of milli-pc 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 (Robson et al., 2019).
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 ${\sim}$15 at redshift
${\sim}$1.5.

Figure 9: 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 (Robson et al.,
2019). The black curve denotes the gravitational-wave signal of a BSBH at
$z=1.53$ with mass $10^{8.4}M_{\odot}$ 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_{\odot}$,
$10^{7}M_{\odot}$, $10^{8}M_{\odot}$, and $10^{9}M_{\odot}$, respectively, at
the same redshift with mass ratios ranging within $0.05\leq q\leq 0.5$. The
blue, orange, green, and red lines correspond to $q=0.1$.

###  4.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 ($\sim$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 associated
(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.

###  4.6 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
$\sim$102–103 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$\sigma$ point source depth of $\sim$23.5 AB mag), the expected number of
detectable periodic quasars from circumbinary disk accretion variability at
redshift z$\sim$1.5 with observer-frame periods between 0.5 and 5.0 yr is
$\sim$80 in an all-sky survey ($\sim$30,000 deg2), or $\sim$1 per 380 deg2
(see right panel in Figure 6 of Kelley et al. 2019). This is $\sim$70 times
lower than our detection rate444Our estimated detection rate depends on the
depth of the parent spectroscopic quasar sample, which is incomplete. We do
not have a complete quasar sample down to DES depth. of $\sim$one strong BSBH
candidate from circumbinary accretion variability per 5 deg2 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 $\sim 10^{8}M_{\odot}$ 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\sim 1.5$.

###  4.7 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., $\sim$0.8%, is $\sim$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, 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.

##  5 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, long-term, 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 $\sim$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$\sigma$ 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 (e.g.,
MacFadyen & Milosavljević, 2008; Shi et al., 2012; Roedig et al., 2012;
D’Orazio et al., 2013; Gold et al., 2014a; Shi & Krolik, 2015; 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 light-curve
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_{{\rm orb}}=t_{{\rm period}}$ (whereas $t_{{\rm
orb}}\approx 0.2t_{{\rm period}}$ for $q$=0.43 instead; Farris et al. 2014),
the inferred binary separation is $d\sim 4.4$ milli-parsec (i.e., 5.1 light
days, or $\sim$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_{{\rm gw}}$ with the
preferred system parameters is $\sim$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. BSBHs
with masses of $\sim 10^{8}$–$10^{9}M_{\odot}$ at redshift $z\gtrsim 1$ are
generally expected around the time of pre-decoupling (Kocsis & Sesana, 2011),
i.e., when $t_{{\rm gw}}>t_{{\rm visc}}$, where $t_{{\rm visc}}$ is the
viscous timescale of the accretion disk. The gravitational-wave strain
amplitude is ${\sim}10^{-18}$ at $\sim$37 nHz, which, as an individual source,
is ${\sim}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 $\sim$5 years before coalescence at
$\gtrsim$0.01 mHz with a signal-to-noise ratio of $\sim$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 $\sim$3 are
well sampled in multiple bands. There is a $\sim$0.1% probability that the
periodogram peak is caused by stochastic quasar variability (i.e., red noise).
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 periodogram 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 (e.g., Foord et
al., 2017) on any potentially characteristic SED features to compare with
predictions from circumbinary accretion disk simulations (e.g., 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-$\alpha$ 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
deg2 survey implies a detection rate of $\sim$0.16%, or 1 per 5 deg2, which is
$\sim$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 in the parent sample is $\sim$4–80 times times of those from
previous searches using other surveys (Chen et al., 2020), 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 Vera C. Rubin
Observatory Legacy Survey of Space and Time (Ivezić et al., 2019) may be able
to detect hundreds to thousands of BSBH candidates from circumbinary accretion
variability.

## Acknowledgements

We thank the anonymous referee for helpful suggestions that improved the
manuscript. X.L. thanks A. Barth, S. Dodelson, S. Gezari, and A. Palmese for
comments, and B. Fields, C. Gammie, K. Gültekin, Z. Haiman, D. Lai, A. Loeb,
D. D’Orazio, S. Tremaine, and X.-J. Zhu for discussions, and LCO director Tod
Boroson for granting us DDT observation. W.-T. L. is supported in part by the
Gordon and Betty Moore Foundation’s Data-Driven Discovery Initiative through
Grant GBMF4561 to Matthew Turk and a government scholarship to study aboard
from the ministry of education of Taiwan. Y.C.C. and X.L. acknowledge a Center
for Advanced Study Beckman fellowship and support from the University of
Illinois campus research board. A.M.H. is supported by the DOE NNSA
Stewardship Science Graduate Fellowship under grant number DE-NA0003864. Y.S.
acknowledges support from the Alfred P. Sloan Foundation 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
under Grant Numbers AST-1138766 and AST-1536171. The DES participants from
Spanish institutions are partially supported by MINECO under grants
AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and
MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE
is partially funded by the CERCA program of the Generalitat de Catalunya.
Research leading to these results has received funding from the European
Research Council under the European Union’s Seventh Framework Program
(FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We
acknowledge support from the Australian Research Council Centre of Excellence
for All-sky Astrophysics (CAASTRO), through project number CE110001020, and
the Brazilian Instituto Nacional de Ciênciae Tecnologia (INCT) e-Universe
(CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics. The United States Government retains
and the publisher, by accepting the article for publication, acknowledges that
the United States Government retains a non-exclusive, paid-up, irrevocable,
world-wide license to publish or reproduce the published form of this
manuscript, or allow others to do so, for United States Government purposes.

We are grateful for the extraordinary contributions of our CTIO colleagues and
the DECam Construction, Commissioning and Science Verification teams in
achieving the excellent instrument and telescope conditions that have made
this work possible. The success of this project also relies critically on the
expertise and dedication of the DES Data Management group.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High-Performance Computing at the University of Utah. The SDSS
web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the Brazilian
Participation Group, the Carnegie Institution for Science, Carnegie Mellon
University, the Chilean Participation Group, the French Participation Group,
Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de
Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley
National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-
Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für
Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik
(MPE), National Astronomical Observatories of China, New Mexico State
University, New York University, University of Notre Dame, Observatário
Nacional / MCTI, The Ohio State University, Pennsylvania State University,
Shanghai Astronomical Observatory, United Kingdom Participation
Group,Universidad Nacional Autónoma de México, University of Arizona,
University of Colorado Boulder, University of Oxford, University of
Portsmouth, University of Utah, University of Virginia, University of
Washington, University of Wisconsin, Vanderbilt University, and Yale
University.

Facilities: DES, LCOGT, Sloan

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author. The light curves of J0252 are available in the article’s
online supplementary material.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2018) Abbott T. M. C., et al., 2018, ApJS, 239, 18
  * Adelman-McCarthy et al. (2007) Adelman-McCarthy J. K., et al., 2007, ApJS, 172, 634
  * Alexander (1997) Alexander T., 1997, in Maoz D., Sternberg A., Leibowitz E. M., eds, Astrophysics and Space Science Library Vol. 218, Astronomical Time Series. p. 163, doi:10.1007/978-94-015-8941-3_14
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, ArXiv e-prints 1702.00786,
  * Arzoumanian et al. (2018a) Arzoumanian Z., et al., 2018a, ApJS, 235, 37
  * Arzoumanian et al. (2018b) Arzoumanian Z., et al., 2018b, ApJ, 859, 47
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Barth & Stern (2018) Barth A. J., Stern D., 2018, ApJ, 859, 10
  * Becker et al. (1995) Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Bernstein et al. (2012) Bernstein J. P., et al., 2012, ApJ, 753, 152
  * Bernstein et al. (2017) Bernstein G. M., et al., 2017, PASP, 129, 114502
  * Berti et al. (2015) Berti E., et al., 2015, Classical and Quantum Gravity, 32, 243001
  * Centrella et al. (2010) Centrella J., Baker J. G., Kelly B. J., van Meter J. R., 2010, Reviews of Modern Physics, 82, 3069
  * Charisi et al. (2015) Charisi M., Bartos I., Haiman Z., Price-Whelan A. M., Márka S., 2015, MNRAS, 454, L21
  * Charisi et al. (2016) Charisi M., Bartos I., Haiman Z., Price-Whelan A. M., Graham M. J., Bellm E. C., Laher R. R., Márka S., 2016, MNRAS, 463, 2145
  * Charisi et al. (2018) Charisi M., Haiman Z., Schiminovich D., D’Orazio D. J., 2018, MNRAS, 476, 4617
  * Chen et al. (2020) Chen Y.-C., et al., 2020, arXiv e-prints, p. arXiv:2008.12329
  * Chiaberge & Marconi (2011) Chiaberge M., Marconi A., 2011, MNRAS, 416, 917
  * Childress et al. (2017) Childress M. J., et al., 2017, MNRAS, 472, 273
  * Colpi (2014) Colpi M., 2014, Space Sci. Rev., 183, 189
  * D’Orazio et al. (2013) D’Orazio D. J., Haiman Z., MacFadyen A., 2013, MNRAS, 436, 2997
  * D’Orazio et al. (2015a) D’Orazio D. J., Haiman Z., Duffell P., Farris B. D., MacFadyen A. I., 2015a, MNRAS, 452, 2540
  * D’Orazio et al. (2015b) D’Orazio D. J., Haiman Z., Schiminovich D., 2015b, Nature, 525, 351
  * Dark Energy Survey Collaboration et al. (2016) Dark Energy Survey Collaboration et al., 2016, MNRAS, 460, 1270
  * Dawson et al. (2016) Dawson K. S., et al., 2016, AJ, 151, 44
  * Duffell (2016) Duffell P. C., 2016, ApJS, 226, 2
  * Duffell et al. (2020) Duffell P. C., D’Orazio D., Derdzinski A., Haiman Z., MacFadyen A., Rosen A. L., Zrake J., 2020, ApJ, 901, 25
  * Edelson et al. (2013) Edelson R., Mushotzky R., Vaughan S., Scargle J., Gandhi P., Malkan M., Baumgartner W., 2013, ApJ, 766, 16
  * Edelson et al. (2014) Edelson R., Vaughan S., Malkan M., Kelly B. C., Smith K. L., Boyd P. T., Mushotzky R., 2014, ApJ, 795, 2
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Farris et al. (2015) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015, MNRAS, 447, L80
  * Ferrarese & Ford (2005) Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523
  * Finkbeiner et al. (2016) Finkbeiner D. P., et al., 2016, ApJ, 822, 66
  * Flaugher (2005) Flaugher B., 2005, International Journal of Modern Physics A, 20, 3121
  * Flaugher et al. (2015) Flaugher B., et al., 2015, AJ, 150, 150
  * Flesch (2015) Flesch E. W., 2015, Publ. Astron. Soc. Australia, 32, e010
  * Foord et al. (2017) Foord A., Gültekin K., Reynolds M., Ayers M., Liu T., Gezari S., Runnoe J., 2017, ApJ, 851, 106
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Frieman et al. (2008) Frieman J. A., et al., 2008, AJ, 135, 338
  * Gibson et al. (2008) Gibson R. R., Brandt W. N., Schneider D. P., 2008, ApJ, 685, 773
  * Gold et al. (2014a) Gold R., Paschalidis V., Etienne Z. B., Shapiro S. L., Pfeiffer H. P., 2014a, Phys. Rev. D, 89, 064060
  * Gold et al. (2014b) Gold R., Paschalidis V., Ruiz M., Shapiro S. L., Etienne Z. B., Pfeiffer H. P., 2014b, Phys. Rev. D, 90, 104030
  * Goldstein et al. (2015) Goldstein D. A., et al., 2015, AJ, 150, 82
  * Goyal et al. (2018) Goyal A., et al., 2018, ApJ, 863, 175
  * Graham et al. (2015) Graham M. J., et al., 2015, MNRAS, 453, 1562
  * Gross & Vitells (2010) Gross E., Vitells O., 2010, European Physical Journal C, 70, 525
  * Haehnelt (1994) Haehnelt M. G., 1994, MNRAS, 269, 199
  * Haehnelt & Kauffmann (2002) Haehnelt M. G., Kauffmann G., 2002, MNRAS, 336, L61
  * Haiman et al. (2009a) Haiman Z., Kocsis B., Menou K., Lippai Z., Frei Z., 2009a, Classical and Quantum Gravity, 26, 094032
  * Haiman et al. (2009b) Haiman Z., Kocsis B., Menou K., 2009b, ApJ, 700, 1952
  * Hatziminaoglou et al. (2005) Hatziminaoglou E., et al., 2005, AJ, 129, 1198
  * Holgado et al. (2018) Holgado A. M., Sesana A., Sandrinelli A., Covino S., Treves A., Liu X., Ricker P., 2018, MNRAS, 481, L74
  * Hopkins et al. (2008) Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Ivezić et al. (2004) Ivezić Ž., et al., 2004, Astronomische Nachrichten, 325, 583
  * Ivezić et al. (2007) Ivezić Ž., et al., 2007, AJ, 134, 973
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jaffe & Backer (2003) Jaffe A. H., Backer D. C., 2003, ApJ, 583, 616
  * Jiang et al. (2007) Jiang L., Fan X., Ivezić Ž., Richards G. T., Schneider D. P., Strauss M. A., Kelly B. C., 2007, ApJ, 656, 680
  * Kellermann et al. (1989) Kellermann K. I., Sramek R., Schmidt M., Shaffer D. B., Green R., 1989, AJ, 98, 1195
  * Kelley et al. (2019) Kelley L. Z., Haiman Z., Sesana A., Hernquist L., 2019, MNRAS, 485, 1579
  * Kelly et al. (2009) Kelly B. C., Bechtold J., Siemiginowska A., 2009, ApJ, 698, 895
  * Kelly et al. (2014) Kelly B. C., Becker A. C., Sobolewska M., Siemiginowska A., Uttley P., 2014, ApJ, 788, 33
  * Kessler et al. (2015) Kessler R., et al., 2015, AJ, 150, 172
  * King et al. (2013) King O. G., et al., 2013, MNRAS, 436, L114
  * Kocsis & Sesana (2011) Kocsis B., Sesana A., 2011, MNRAS, 411, 1467
  * Kormendy & Richstone (1995) Kormendy J., Richstone D., 1995, ARA&A, 33, 581
  * Kozłowski (2016) Kozłowski S., 2016, ApJ, 826, 118
  * Kozłowski (2017) Kozłowski S., 2017, ApJS, 228, 9
  * Lawrence et al. (2007) Lawrence A., et al., 2007, MNRAS, 379, 1599
  * Liu et al. (2015) Liu T., et al., 2015, ApJ, 803, L16
  * Liu et al. (2016) Liu T., et al., 2016, ApJ, 833, 6
  * Liu et al. (2019) Liu T., et al., 2019, ApJ, 884, 36
  * MacFadyen & Milosavljević (2008) MacFadyen A. I., Milosavljević M., 2008, ApJ, 672, 83
  * MacLeod et al. (2010) MacLeod C. L., et al., 2010, ApJ, 721, 1014
  * Martin et al. (2005) Martin D. C., et al., 2005, ApJ, 619, L1
  * Martini & Weinberg (2001) Martini P., Weinberg D. H., 2001, ApJ, 547, 12
  * Masci et al. (2019) Masci F. J., et al., 2019, PASP, 131, 018003
  * McKernan & Ford (2015) McKernan B., Ford K. E. S., 2015, MNRAS, 452, L1
  * Merritt (2013) Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei. Princeton Series in Astrophysics, Princeton University Press, http://books.google.com/books?id=cOa1ku640zAC
  * Milosavljević & Phinney (2005) Milosavljević M., Phinney E. S., 2005, ApJ, 622, L93
  * Moreno et al. (2019) Moreno J., Vogeley M. S., Richards G. T., Yu W., 2019, PASP, 131, 063001
  * Morganson et al. (2014) Morganson E., et al., 2014, ApJ, 784, 92
  * Morganson et al. (2018) Morganson E., et al., 2018, PASP, 130, 074501
  * Mushotzky et al. (2011) Mushotzky R. F., Edelson R., Baumgartner W., Gandhi P., 2011, ApJ, 743, L12
  * Ofek et al. (2012) Ofek E. O., et al., 2012, PASP, 124, 62
  * Oke (1974) Oke J. B., 1974, ApJS, 27, 21
  * Oke & Gunn (1983) Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
  * Pâris et al. (2018) Pâris I., et al., 2018, A&A, 613, A51
  * Richards et al. (2006) Richards G. T., et al., 2006, ApJS, 166, 470
  * Robson et al. (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Roedig et al. (2012) Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P., Haardt F., 2012, A&A, 545, A127
  * Roedig et al. (2014) Roedig C., Krolik J. H., Miller M. C., 2014, ApJ, 785, 115
  * Schneider et al. (2007) Schneider D. P., et al., 2007, AJ, 134, 102
  * Schutz (2018) Schutz B. F., 2018, Philosophical Transactions of the Royal Society of London Series A, 376, 20170279
  * Sesana (2017) Sesana A., 2017, in Journal of Physics Conference Series. p. 012018 (arXiv:1702.04356), doi:10.1088/1742-6596/840/1/012018
  * Sesana et al. (2018) Sesana A., Haiman Z., Kocsis B., Kelley L. Z., 2018, ApJ, 856, 42
  * Shen (2009) Shen Y., 2009, ApJ, 704, 89
  * Shen (2013) Shen Y., 2013, Bulletin of the Astronomical Society of India, 41, 61
  * Shen & Liu (2012) Shen Y., Liu X., 2012, ApJ, 753, 125
  * Shen & Loeb (2010) Shen Y., Loeb A., 2010, ApJ, 725, 249
  * Shen et al. (2019) Shen Y., et al., 2019, ApJS, 241, 34
  * Shi & Krolik (2015) Shi J.-M., Krolik J. H., 2015, ApJ, 807, 131
  * Shi et al. (2012) Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, ApJ, 749, 118
  * Smith et al. (2018) Smith K. L., Mushotzky R. F., Boyd P. T., Malkan M., Howell S. B., Gelino D. M., 2018, ApJ, 857, 141
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * Tang et al. (2018) Tang Y., Haiman Z., MacFadyen A., 2018, MNRAS, 476, 2249
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration 2005, arXiv e-prints, pp astro–ph/0510346
  * Thorne & Braginskii (1976) Thorne K. S., Braginskii V. B., 1976, ApJ, 204, L1
  * Tie et al. (2017) Tie S. S., et al., 2017, AJ, 153, 107
  * Valtonen et al. (2008) Valtonen M. J., et al., 2008, Nature, 452, 851
  * VanderPlas & Ivezić (2015) VanderPlas J. T., Ivezić Ž., 2015, ApJ, 812, 18
  * Vanderplas et al. (2012) Vanderplas J., Connolly A., Ivezić Ž., Gray A., 2012, in Conference on Intelligent Data Understanding (CIDU). pp 47 –54, doi:10.1109/CIDU.2012.6382200
  * Vaughan & Uttley (2005) Vaughan S., Uttley P., 2005, MNRAS, 362, 235
  * Vaughan et al. (2016) Vaughan S., Uttley P., Markowitz A. G., Huppenkothen D., Middleton M. J., Alston W. N., Scargle J. D., Farr W. M., 2016, MNRAS, 461, 3145
  * Vecchio (1997) Vecchio A., 1997, Classical and Quantum Gravity, 14, 1431
  * Vestergaard & Osmer (2009) Vestergaard M., Osmer P. S., 2009, ApJ, 699, 800
  * Villarreal Hernández & Andernach (2018) Villarreal Hernández A. C., Andernach H., 2018, ArXiv e-prints 1808.07178,
  * Voges et al. (2000) Voges W., et al., 2000, IAUC, 7432, 3
  * Volonteri et al. (2003) Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
  * Wang & Mohanty (2017) Wang Y., Mohanty S. D., 2017, Physical Review Letters, 118, 151104
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * York et al. (2000) York D. G., et al., 2000, AJ, 120, 1579
  * Yu & Tremaine (2002) Yu Q., Tremaine S., 2002, MNRAS, 335, 965
  * Zechmeister & Kürster (2009) Zechmeister M., Kürster M., 2009, A&A, 496, 577
  * Zheng et al. (2016) Zheng Z.-Y., Butler N. R., Shen Y., Jiang L., Wang J.-X., Chen X., Cuadra J., 2016, ApJ, 827, 56
  * Zhu & Thrane (2020) Zhu X.-J., Thrane E., 2020, ApJ, 900, 117
  * Zhu et al. (2014) Zhu X.-J., et al., 2014, MNRAS, 444, 3709
  * van der Klis (1989) van der Klis M., 1989, ARA&A, 27, 517

##  Appendix A Details on archival photometries

The publicly available PTF photometry was in PTF $g$ and $R$ bands in Vega
mags. For consistency we have converted them to the SDSS g and r bands in AB
mags following the empirically calibrated relations based on PTF stars
(Equations 4 and 5 of 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_{{\rm SDSS}}$-RPTF/SDSS starting around $r_{{\rm 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_{{\rm SDSS}}{\sim}$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.

## affliations

1Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana,
IL 61801, USA
2National Center for Supercomputing Applications, University of Illinois at
Urbana-Champaign, 605 East Springfield Avenue, Champaign, IL 61820, USA
3School of Mathematics and Physics, University of Queensland, QLD 4072
Australia
4Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA
5Department of Astronomy and Astrophysics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL 60637, USA
6Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus, OH 43210, USA
7Center for Cosmology and Astro-Particle Physics, The Ohio State University,
191 West Woodfuff Avenue, Columbus OH, 43210, USA
8Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
9Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK
10Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA
11Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth, PO1 3FX, UK
12LSST, 933 North Cherry Avenue, Tucson, AZ 85721, USA
13CNRS, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France
14Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d’Astrophysique
de Paris, F-75014, Paris, France
15 Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK
16 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Madrid, Spain
17 Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José
Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil
18 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain
19 Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
20 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain
21 Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450,
Stanford University, Stanford, CA 94305, USA
22 Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA
23 Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ -
20921-400, Brazil
24 Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India
25 Department of Astronomy/Steward Observatory, 933 North Cherry Avenue,
Tucson, AZ 85721-0065, USA
26 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Dr., Pasadena, CA 91109, USA
27 Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA
28 Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA
29 Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049
Madrid, Spain
30 Centre for Astrophysics & Supercomputing, Swinburne University of
Technology, Victoria 3122, Australia
31 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
32 Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093
Zurich, Switzerland
33 Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA
34 Department of Physics, The Ohio State University, Columbus, OH 43210, USA
35 Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse,
85748 Garching, Germany
36 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians
Universität München, Scheinerstr. 1, 81679 München, Germany
37 Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA
38 Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113,
Australia
39 Departamento de Física Matemática, Instituto de Física, Universidade de São
Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil
40 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and
Astronomy, and Department of Physics and Astronomy, Texas A&M University,
College Station, TX 77843, USA
41 Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain
42 School of Physics and Astronomy, University of Southampton, Southampton,
SO17 1BJ, UK
43 Cerro Tololo Inter-American Observatory, National Optical Astronomy
Observatory, Casilla 603, La Serena, Chile
44 Brandeis University, Physics Department, 415 South Street, Waltham MA 02453
45 Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas,
13083-859, Campinas, SP, Brazil
46 Computer Science and Mathematics Division, Oak Ridge National Laboratory,
Oak Ridge, TN 37831
47 Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA