# Reliable Identification of Binary Supermassive Black Holes
from Rubin Observatory Time-Domain Monitoring

Megan C. Davis NSF Graduate Research Fellow  Department of Physics, 196A
Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA
Kaylee E. Grace Department of Physics & Astronomy, 217 Sharp Lab, Newark, DE
19716, USA  Jonathan R. Trump Department of Physics, 196A Auditorium Road,
Unit 3046, University of Connecticut, Storrs, CT 06269, USA  Jessie C. Runnoe
Department of Physics and Astronomy, Vanderbilt University, Nashville, TN
37235, USA  Department of Life and Physical Sciences, Fisk University, 1000
17th Avenue N, Nashville, TN, 37208, USA  Amelia Henkel Department of Physics,
University of New Hampshire, 9 Library Way, Durham, NH 03824, USA  Laura
Blecha Department of Physics, University of Florida, Gainesville, FL 32611,
USA  W. N. Brandt Department of Astronomy & Astrophysics, 525 Davey Lab, The
Pennsylvania State University, University Park, PA 16802, USA  Institute for
Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA  Department of Physics, 104 Davey Lab, The Pennsylvania
State University, University Park, PA 16802, USA  J. Andrew Casey-Clyde
Department of Physics, 196A Auditorium Road, Unit 3046, University of
Connecticut, Storrs, CT 06269, USA  Maria Charisi Department of Physics and
Astronomy, Vanderbilt University, Nashville, TN 37235, USA  Caitlin Witt
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, Evanston, IL 60208, USA  Adler Planetarium, 1300 S.
DuSable Lake Shore Dr., Chicago, IL 60605, USA

###### Abstract

Periodic signatures in time-domain observations of quasars have been used to
search for binary supermassive black holes. These searches, across existing
time-domain surveys, have produced several hundred candidates. The general
stochastic variability of quasars, however, can masquerade as a false-positive
periodic signal, especially when monitoring cadence and duration are limited.
In this work, we predict the detectability of binary supermassive black holes
in the upcoming Rubin Observatory Legacy Survey of Space and Time (LSST). We
apply computationally inexpensive sinusoidal curve fits to millions of
simulated LSST Deep Drilling Field light curves of both single, isolated
quasars and binary quasars. Period and phase of simulated binary signals can
generally be disentangled from quasar variability. Binary amplitude is
overestimated and poorly recovered for two-thirds of potential binaries due to
quasar accretion variability. Quasars with strong intrinsic variability can
obscure a binary signal too much for recovery. We also find that the most
luminous quasars mimic current binary candidate light curves and their
properties: false positive rates are 60% for these quasars. The reliable
recovery of binary period and phase for a wide range of input binary LSST
light curves is promising for multi-messenger characterization of binary
supermassive black holes. However, pure electromagnetic detections of binaries
using photometric periodicity with amplitude greater than 0.1 magnitude will
result in samples that are overwhelmed by false positives. This paper
represents an important and computationally inexpensive way forward for
understanding the true and false positive rates for binary candidates
identified by Rubin.

Gravitational Wave Sources – Surveys – Methods: Data Analysis – Techniques:
Photometric – Quasars: General – Quasars: Supermassive Black Holes

††journal: ApJ††software:  astroML (Vanderplas et al., 2012), SciPy (Virtanen
et al., 2020), NumPy (van der Walt et al., 2011), LMFits (Newville et al.,
2014)

##  1 Introduction

Binary supermassive black holes (SMBHs) are an expected consequence of galaxy
mergers (e.g., Kormendy & Ho, 2013). Close, sub-parsec binary SMBHs should be
the loudest gravitational wave sources in the Universe, likely producing the
recently announced low-frequency gravitational wave background detection from
the North American Nanohertz Observatory for Gravitational Waves (NANOGrav;
Agazie et al., 2023a, b) and other global Pulsar Timing Arrays (PTAs;
Antoniadis et al., 2023a; Reardon et al., 2023a). The binary SMBH population
traces galaxy merger rates, helping our understanding of both SMBH-galaxy co-
evolution and the hierarchical structure formation of galaxies (Comerford et
al., 2018). Constraints on inferred merger rates and the gravitational wave
background indicate false-positive contamination of the binary quasar
candidate population (see binary candidate discussion below; Sesana et al.,
2018). Identifying binary SMBHs via electromagnetic methods makes individual
gravitational wave sources easier to detect with pulsar timing arrays (e.g.,
Arzoumanian et al., 2020; Liu & Vigeland, 2021) and the upcoming,
revolutionary space-based gravitational wave observatory, LISA (Amaro-Seoane
et al., 2017, 2023). Beyond gravitational wave astrophysics, an
electromagnetic binary SMBH detection further opens the door to study black
hole inspiral and accretion environments, cosmology, and more (De Rosa et al.,
2019).

There are confirmed electromagnetic detections of dual quasars that are
closely separated (embedded within the same galaxy via merger, kiloparsec to
parsec scales) but not yet gravitationally bound (parsec to sub-parsec scales)
as binary SMBHs (e.g.,Van Wassenhove et al., 2012; Blumenthal & Barnes, 2018;
De Rosa et al., 2019; Chen et al., 2022; Koss et al., 2023). Candidate binary
SMBH systems are identified from observations spanning the electromagnetic
spectrum, such as notable radio interferometry candidates in Rodriguez et al.
(2006) and Bansal et al. (2017). In this work, we focus on photometric
searches of binary SMBHs in the optical regime.

Binary SMBH systems are expected to exhibit periodic modulation of photometric
light curves due to Doppler boosting or periodic accretion (more detailed
descriptions of these binary models can be found in Section 2.2; Charisi et
al., 2022). Equal-mass binaries (of mass ratio within a factor of a few) will
generally have the strongest periodic signatures in photometric observations,
in addition to being the loudest expected gravitational wave sources (D'Orazio
et al., 2013). Graham et al. (2015); Charisi et al. (2016); Liu et al. (2019);
Chen et al. (2020, 2022b) are all examples of optical photometric searches for
periodic behavior that have produced several hundred binary SMBH candidates
(see De Rosa et al., 2019; D’Orazio & Charisi, 2023 and references therein for
a review). Despite these efforts, we have yet to confirm these binary SMBH
candidates and true binaries remain elusive. A strongly varying single quasar
can appear indistinguishable from a binary SMBH system because intrinsic
optical quasar variability can masquerade as a periodic signature of a binary
SMBH. Limited survey cadence and duration can especially lead to false
positive detections (Vaughan et al., 2016; Barth & Stern, 2018; Witt et al.,
2022).

Upcoming massive time-domain surveys, like the Vera Rubin Observatory
(Rubin)’s 10-year Legacy Survey of Space and Time (LSST), will open a new
frontier for electromagnetic detection of binary SMBHs. LSST is expected to
observe between 20-100 million active galactic nuclei (AGN), with time-domain
monitoring that has a cadence and duration appropriate to detect potential
sub-parsec binary SMBHs (Ivezić, 2017; Kelley et al., 2019, 2021; Xin &
Haiman, 2021; Charisi et al., 2022). These AGN will be observed via the main
LSST survey, Wide-Fast-Deep (WFD; a footprint of roughly 18,000 deg2), that
uses 80%-90% of the available observation time and the Deep Drilling Fields
(DDFs; each field has a footprint of roughly 9.6 deg2), which will use a
portion of the remaining 10%-20% of the observation time (Bianco et al.,
2022). In this work we focus on photometric light curves in DDF observations,
representing the highest-cadence and deepest photometry for time-domain
observations of AGN by Rubin/LSST. The DDFs are expected to observe over
40,000 AGN (Brandt et al., 2018).

Rubin’s LSST, our best chance at identifying binary SMBHs with electromagnetic
observations, also pushes us further into the era of big data, as it is
predicted that it will produce over 20 terabytes of data per
night111https://www.lsst.org/about/dm (). In order to triage data as it
arrives, LSST will have to employ pipeline-style brokerage, or target-type
sorting, using various, computationally inexpensive analysis methods, from
curve fitting to machine learning (Saha et al., 2016; Narayan et al., 2018;
Matheson et al., 2021). In this work, we test a computationally inexpensive
method for identifying binary SMBHs using simulated periodic and single-quasar
light curves that resemble the expected observational properties of Rubin’s
LSST. This method, simple sinusoidal curve fits, could work as an initial
filter to identify binary candidates for further, more intense follow-up like
that of Witt et al. (2022). Witt et al. (2022) applies complex Bayesian
methods to simulated CRTS and LSST light curves that require several CPU hours
per source, which is computationally prohibitive to apply to the entire
database. Our work can predict LSST detection rates for binary SMBHs and give
insight into what types of quasars will produce false detections from millions
of light curves within a few hours.

In the next section, Section 2, we describe the construction and underlying
models of the Rubin LSST light curves. Section 3 details the simple sinusoidal
curve fitting. We present the results of the light curve fitting, including
predicted false positive rates and population statistics, in Section 4. We
discuss the results, with recommendations for future Rubin LSST observations
and comparisons with gravitational wave searches, in Section 5. We summarize
our conclusions and outline future prospects in Section 6.

We assume a flat, $\Lambda$CDM cosmology with $\Omega_{\Lambda}$ = 0.7,
$\Omega_{M}$ = 0.3, and $H_{0}$ = 67.4 $\pm$ 0.5 km s-1 Mpc-1 for this work
(Planck Collaboration et al., 2018).

##  2 Light Curve Generation

Our goal is to create realistic light curves of quasars, both isolated
(single-SMBH) and binary systems, for Rubin’s LSST Deep Drilling Fields
(DDFs). Light curve generation is accomplished in three steps: 1. generating
the quasar variability, 2. creating the binary signal (for binaries only), and
then 3. customizing for the observational noise, cadence, and monitoring
duration specific to the LSST DDFs in the $i$-band. Figure 1 summarizes these
steps, with the details in the sections that follow. Both isolated and binary
quasars are expected to have accretion disks that exhibit stochastic
variability, as close-separation binaries are predicted to have a circumbinary
accretion disk and mini-disks around each SMBH (Charisi et al., 2022;
Bogdanović et al., 2022; D’Orazio & Charisi, 2023). This is shown in the first
panel of Figure 1. The binary systems are expected to have an additional
periodic variability component, as described in Section 2.2 and seen in panel
two of Figure 1. The customization of the light curves included implementing
the expected cadence, duration, and noise properties of the LSST DDFs, as laid
out in Section 2.3 and seen in the last panel of Figure 1. The range of quasar
and binary parameters used to create the light curves are shown in Table 1.

Figure 1: A summary of the light curve assembly as outlined in Section 2. The top most panel is the variability model/base light curve creation for a $10^{8}$ $M_{\odot}$ quasar at z $\sim$ 1 that is accreting at 10% of its Eddington limit, with 10 additional realizations in faded colors behind the main black continuous line. The second panel reveals the periodicity model of the binary (some combination of sine and sawtooth waves). A black smooth sine is depicted here with 10 additional realizations behind it (in 10% increments of pure sine to sawtooth). A light curve for an isolated, single quasar would skip the second and third panel steps. We add the two models together in the third panel. We finally then add noise and apply a cadence and seasonal duration to match the expected Rubin LSST DDF observation in the final panel.  Parameter | [lower bound, upper bound]
---|---
log($M/M_{\odot}$) | [6.5, 9.5]
log($L/L_{edd}$) | [-2, 0]
z | [0.1, 2.5]
amplitude (mag) | [0.0, 1.0]
period (days) | [360, 1080]
sawtoothiness (%) | [0.0, 1.0 ]
phase | [0.0, 2$\pi$]
Table 1: The range of quasar properties of total binary mass, luminosity, and
redshift and binary signal parameters of amplitude, period, signal type, phase
used to generate the simulated light curves. A zero amplitude binary signal
represents the non-binary (single) quasar population. We sample each parameter
from a uniform distribution within its stated range.

In order to explore the full parameter range and unique combinations of their
respective uniform and continuously sampled ranges, we generated a total of
4,812,500 $i$-band LSST light curves. Over 3.6 million light curves met a
$i$-band magnitude limit of 23.5 that represents a light curve with sufficient
signal-to-noise for reliable fitting, and they were kept for the analysis
described in Section 3.

###  2.1 Quasar Properties

Our goal is to simulate light curves for a broad and representative range of
the quasar population to be observed by Rubin. We made a grid of quasar
properties spanning black hole mass (total mass of the binary, in the case of
a binary SMBH light curve), luminosity, and redshift, as seen as the first
three entries of Table 1. These ranges of quasar properties are representative
of the quasar population observed by the Sloan Digital Sky Survey and other
large surveys (e.g. Kollmeier et al., 2006; Vestergaard & Osmer, 2009; Kelly
et al., 2010; Shen et al., 2011; Wu & Shen, 2022). We chose this redshift
range, specifically, because a majority of $i<23.5$ of DDF quasars will be in
a redshift range of $z<2.5$ (see Table 10.2 of Chapter 10 in LSST Science
Collaboration et al., 2009). This restriction will also avoid Ly$\alpha$
forest absorption in the bluest LSST filter and is further motivated by GW
searches of binary candidates (Charisi et al., 2022). The parameters are
uniformly randomly sampled from within each bound to approximate continuous
distributions. Having such a broad range and thousands of unique combinations
of quasar parameters allows us to probe the properties of detectable binaries
and of false-positive isolated SMBHs with this method.

Our light curves do not include host-galaxy contamination, which may be
important for low luminosity AGN within the DDF footprints (LSST Science
Collaboration et al., 2009). Host-galaxy light is constant and will dilute the
fractional variability observed, equally for binary and variability
components. It would not change the relative variability of DRW compared to
binary signals. For this reason, we predict that the addition of host-galaxy
light will not change our general conclusions about binary identification and
false-positive rates relative to DRW variability contribution.

####  2.1.1 Variability Model

The UV/optical variability of a quasar is effectively described as a
stochastic process specifically in the form of a damped random walk (DRW)
(Kelly et al., 2009; MacLeod et al., 2010, 2012; Zu et al., 2013). This
variability is likely driven by thermal fluctuations in the accretion disk
(Kelly et al., 2009; MacLeod et al., 2010; Kozłowski et al., 2010; Zu et al.,
2013; Suberlak et al., 2021). The DRW structure function, SF, is defined as:

$\rm SF\left(\Delta t\right)=SF_{\infty}\left(1-e^{-|\Delta
t|/\tau}\right)^{1/2}$ (1)

The two parameters of the DRW are the characteristic damping timescale,
$\tau$, and DRW amplitude, $\rm SF_{\infty}$. The values of each parameter
depend on rest-frame wavelength, quasar luminosity, black hole mass, and
redshift (e.g., Uomoto et al., 1976; Hook et al., 1994; Giveon et al., 1999;
Hawkins, 2002; Vanden Berk et al., 2004; de Vries et al., 2005; Wilhite et
al., 2008; Bauer et al., 2009; Kelly et al., 2009; Kozłowski et al., 2010;
MacLeod et al., 2010, 2012). We chose to describe them via Equation 7 of
MacLeod et al. (2010):

$\displaystyle\log f=A+B\log\left(\frac{\lambda_{\rm
RF}}{4000\textup{\AA}}\right)+C\left(M_{i}+23\right)$ (2)
$\displaystyle+D\log\left(\frac{M_{\rm
BH}}{10^{9}M_{\odot}}\right)+E\log\left(1+z\right)$

where $\lambda_{\rm RF}$ is the rest-frame wavelength, $M_{i}$ is the absolute
$i$-band magnitude of the target, $M_{\rm BH}$ is the mass of the black hole,
and $z$ is redshift. We used the equation above with the best-fit coefficients
of MacLeod et al. (2010, 2012) summarized in Table 2, for quasars with black
hole mass, luminosity, and redshift drawn from the ranges described in Table
1. We use the generate_dampedRW function of the AstroML python package to
create each DRW from the coefficients calculated by Equation 2.

The DRW is an effective empirical model for quasar variability on timescales
of months to years (although it appears to be less effective on shorter and
longer timescales, Mushotzky et al., 2011; Kelly et al., 2014; Vaughan et al.,
2016; Kozłowski, 2017a; Moreno et al., 2019; Kozłowski, 2021; Yu et al.,
2022). We chose the DRW for its simplicity and accuracy in describing the
variability of quasar accretion disks with the parameter space described
above. We note, however, that the damping timescale is not well-constrained by
observations for much of the quasar parameter space (MacLeod et al., 2010,
2012; Kozłowski, 2017a, b). The damping timescale $\tau$ is likely to drive
the “recovered” periods of false-positive signals. There is also some debate
in the literature about the details of the correlations of the DRW parameters
and quasar properties (see Table 4 in Suberlak et al., 2021), although a
detailed analysis by Suberlak et al. (2021) found similar correlations to
MacLeod et al. (2010, 2012). We use the best-available model for $\tau$ in our
simulated light curves, while noting that future long-duration monitoring
(e.g., with Rubin/LSST) is likely to improve the quantitative description of
quasar variability on long timescales.

The general schema of this project is flexible enough to incorporate other
variability models in future work, like the continuous-time autoregressive
moving average (CARMA) models (Kelly et al., 2014) and the damped harmonic
oscillator (DHO) model (Yu et al., 2022). Future work can also include a more
detailed consideration of the differences between accretion disks around
single quasars with the circumbinary and mini-disks in binary quasar
(Artymowicz & Lubow, 1996; MacFadyen & Milosavljević, 2008; Westernacher-
Schneider et al., 2023).

f | A | B | C | D | E
---|---|---|---|---|---
$\rm SF_{\infty}$[mag] | -0.36 | -0.479 | 0.131 | 0.18 | 0.0
$\tau$ [days] | 2.7 | 0.17 | 0.03 | 0.21 | 0.0
Table 2: The coefficients of the best-fit DRW parameter formulae from MacLeod
et al. (2010), with updates on the A coefficient from MacLeod et al. (2012) to
account for long-term variability.

###  2.2 Binary Signal Parameters

We describe the electromagnetic signal of a binary quasar with a semi-
sinusoidal periodic function that is parameterized with amplitude, period,
phase, and “sawtoothiness.” The ranges of amplitude and period are shown in
Table 1 and are designed to include close-separation binaries that have values
detectable within the first half of the LSST DDF monitoring duration, but not
so short than the binary evolves within the duration of the light curve. We
sample phase from the full range of [0, 2$\pi$]: phase is generally treated as
a nuisance parameter in electromagnetic binary searches, although its recovery
is important for gravitational wave detection (Charisi et al., 2022). We
additionally allowed for an extra parameter (“sawtoothiness”/“burstiness”) to
account for non-sinusoidal, sudden variability and allow for combinations of
both functions.

Sinusoidal variability is predicted for binary SMBHs from relativistic Doppler
boosting (D'Orazio et al., 2015a; Charisi et al., 2018; De Rosa et al., 2019;
Charisi et al., 2022). Especially at the small separations associated with
gravitational wave emission, relativistic Doppler boosting will cause periodic
modulation of the photometric light curve in the form of a sine wave if the
orbits are circular (D'Orazio et al., 2015b; Tang et al., 2018). This is the
most ideal binary signature as all of the binary signal parameters (amplitude,
period, and phase) can be directly linked to the binary’s geometry and orbital
dynamics, and, thus, to the gravitational wave signal counterpart (Charisi et
al., 2022). The photometric binary amplitude is related to the orbital
velocity of the binary in the Doppler Boost model. The amplitude has uncertain
connections to binary parameters for other binary models (Duffell et al.,
2020). We sample periodic binary amplitude from a broad, generic range of
magnitudes to test general recovery with Rubin. Our model for the binary
signal is empirical rather than explicitly connected to binary mass or
luminosity ratio, although our underlying assumptions are generally associated
with equal-mass binaries with an active secondary SMBH, since these binary
systems are expected to produce the strongest periodic signatures in
photometry (D'Orazio et al., 2013; D’Orazio & Charisi, 2023).

The “bursty” sawtooth can be tied to eccentric orbits and to orbital periods
of hot spots (overdensities) on the edge of the internal cavity of the
circumbinary accretion disk and the sudden inflow of that built-up material
(periodic accretion; Duffell et al., 2020). This is the less ideal binary
signal as it is more difficult to tie the amplitude, period, and phase back to
the orbital dynamics of the binary (Charisi et al., 2022). Hot spots in the
inner edge of the circumbinary accretion disk can cause a photometric period
that is 5-8 times the actual orbital period for equal-mass binaries (D’Orazio
et al., 2013; Farris et al., 2014; Miranda et al., 2017; Charisi et al.,
2022). We do not explicitly model longer observed-frame period sawtooth waves
than sine waves nor do we impose a higher frequency sinusoid on top of
sawtooth (Duffell et al., 2020). Instead our range of binary parameters is
defined to include a representative range of binaries detectable by Rubin and
the DDF cadence. To model a sawtooth, we used the reverse of the Scipy.signal
sawtooth wave, such that there is a “burst”, or sudden increase, in
brightness, and then a gentle decrease. An example of the sawtooth signal can
be seen in the second panel of Figure 1.

We allow the light curves to have a combination of both signals through the
“sawtoothiness” parameter of Table 1. If “sawtoothiness” is 1.0, the binary
signal is a pure sawtooth signal. If it was 0.0, it is a smooth sine. Anything
in between is a convolved signal of the two models that takes the form:

$\displaystyle\Delta m(t)=(1-S)*A\sin\left(ft+\phi\right)$ (3)
$\displaystyle+~{}S*A\ \textrm{sawtooth}\left(ft+\phi\right),$

where $\Delta m$ is the change in magnitude, $m$, as a function of time $t$
with parameters $A$ (amplitude), $f$ (frequency, $f=2\pi/P$), $\phi$ (phase),
and $S$ for the sawtooth parameter percentage. Panel two of Figure 1 shows
examples of the convolved signal. The binary signal for each simulated light
curve is drawn from a uniform distribution of sawtoothiness between 0 and 1.

Figure 1 displays full light curves in the bottom panel, in addition to the
binary signals described in this section in the second panel. More examples of
full light curves, with and without binaries, over a range of quasar masses
and luminosities can be seen in Figure 2. This figure especially highlights
how much a DRW can hide or mimic a binary signal.

Figure 2: Examples of synthetic 5-year LSST quasar light curves with seven
month seasonal duration, featuring a single quasar (grey) and two potential
binary systems with varying amplitudes (red: weak, blue: strong). Periodic
binary signals can be difficult to distinguish from quasi-periodic quasar
variability. As seen in the upper left (low mass and low luminosity quasar
variability), a strong binary signal looks similar the underlying DRW. With a
high mass and high luminosity quasar, seen in the lower right, variability
amplitudes can be suppressed enough to directly mimic a “weak” binary signal.
Moving down the y-axis, the black hole mass increases, moving from left to
right across the x-axis, the luminosity increases. Representative values of
these two properties were plugged into the DRW best-fit coefficient formulae
of MacLeod et al. (2010) and MacLeod et al. (2012) to generate the overlaid
quasar variability. Reliable binary detection with sinusoidal curve fits will
be a strong function of quasar properties.

###  2.3  Rubin Deep Drilling Field Observational Parameters

The Deep Drilling Fields are the most interesting Rubin/LSST pointings for
identifying binary quasars. The DDFs are five observation fields (COSMOS, XMM-
LSS, ELAIS-S1, EDF-S, and CDF-S) that will receive high-cadence, multi-band
observations and also already have extraordinary multi-wavelength coverage
over the last two decades. They are expected to have an average cadence222In
detail, the Rubin DDF cadence is planned to cycle through each of the $ugrizY$
filters every $\sim$3 nights, with different filters observed each night
(Brandt et al., 2018; Bianco et al., 2022). The detailed DDF cadence may also
follow an accordion pattern rather than a constant cadence (Scolnic et al.,
2018). We simplify by assuming that observations in different filters can be
combined to create a higher-cadence light curve. of 1-3 night(s) over the full
10-year baseline of the survey (Brandt et al., 2018; Ivezić et al., 2019;
Bianco et al., 2022). Their long baselines and high cadences will provide the
best opportunity for electromagnetic detection of binary SMBHs. We can naively
estimate a binary population for the LSST DDFs by extrapolating from the
binary candidate population of current/previous time-domain monitoring
campaigns, like that of the Catalina Real-Time Transient Survey (CRTS, Drake
et al., 2009). Graham et al. (2015) found 111 CRTS candidate objects out of
243,500 spectroscopically confirmed quasars. For the DDFs, this implies about
20 binary candidates, although the higher depth and signal-to-noise of the
DDFs compared to CRTS means that binaries may be identified in an entirely
different parameter space such that this is an uncertain prediction.

We make our simulated light curves representative of Rubin/LSST DDF monitoring
by adding flux uncertainty (including the effects of weather and lunation,
described further below) and masking to match the anticipated cadence and
duration. We apply a seasonal monitoring duration of at least seven months out
of a year to simulate light curves for one of the equatorial DDFs (COSMOS or
XMM-LSS). The seasonal monitoring duration would increase for either of the
fields at lower declination (ELAIS-S1 or CDF-S).

The photometric error for each light curve was calculated with Equations 4 and
5 from Ivezić et al. (2019). From their Equation 4, the error of a single
visit, $\sigma_{1}$, is found from adding the random photometric error,
$\sigma_{rand}$, and the systematic photometric error, $\sigma_{sys}$ in
quadrature. The LSST team predicts that $\sigma_{sys}<0.005$ mag. We set it to
$0.003$ mag. Then $\sigma_{rand}$ can be found from their magnitude-dependent
Equation 5:

$\sigma^{2}_{rand}=\left(0.04-\gamma\right)x+\gamma x^{2}\left(mag^{2}\right)$
(4)

where $\gamma$ is a unitless filter-dependent sky brightness parameter and
$x\equiv 10^{0.4\left(m-m_{5}\right)}$ with $m_{5}$ being the $5\sigma$ depth
for a point source. For our $i$-band light curves, $\gamma=0.039$ and
$m_{5}=23.92$, as described in Table 2 of Ivezić et al. (2019). The $x$
parameter leads to sources with $i\gtrsim 23.5$ having exponentially
increasing flux uncertainties. We ultimately apply a $i<23.5$ quality cut on
our light curves as a representative range for LSST light curves of
sufficiently high quality for binary characterization.

We make these light curves more realistic by further including error
multipliers for the moon phase and environmental factors such as weather and
photometric conditions. The moon phase error multiplier is a sinusoidal
function based on the 28-day moon cycle, with values ranging from 1.0x (dark
time, new moon) to 2.5x (bright time, full moon). The weather multiplier was a
random, discrete choice between 0.5x (excellent photometric conditions), 1.0x
(average), and 1.5x (poor photometric conditions), randomly applied to the
light curves with weights of 0.15, 0.75, and 0.15, respectively. A weather
multiplier for variations in seeing on an observation-to-observation basis,
with the average resulting in a multiplier of 1x. We also simulate
observations lost due to weather or engineering time by randomly removing 30%
of the points in each light curve.

Once the variability model and binary signal were added together, the
individual points were then randomly resampled within their photometric error
and replaced (bootstrap re-sampling). The bootstrap re-sampling was done by
randomly selecting a new magnitude from a normal distribution of a single
observation’s magnitude and its corresponding error.

Summarizing, we apply the following treatment to the light curves to resemble
anticipated LSST DDF light curves:

  * •

1-day cadence

  * •

7-month seasonal duration

  * •

Flux uncertainty predicted for LSST (Ivezić et al., 2019)

  * •

Additional flux uncertainty associated with lunation and (random) seeing
conditions

  * •

30% (random) cadence losses due to weather

This final stage of the light curve generation could be easily adapted to the
cadence, duration, and noise properties of other surveys. In this work we
focus on the upcoming LSST DDF monitoring, as the most efficient near-future
survey for discovering binary quasars, but in future work we anticipate
applying our pipeline to both the LSST main survey and the lower cadence and
noisier monitoring observations of CRTS. If we assume a similar number of
binary candidates as identified in Graham et al. (2015), then the LSST WFD
main survey could have approximately 23,000 candidate binaries which is
promising for a photometric binary detection confirmation in the near future.

##  3 Light Curve Fitting

We chose to implement simple sinusoidal curve fits to the light curves as the
detection method in this work. This method is computationally inexpensive and
can provide insight into both binary parameter recovery rates and false
positive rates. A sine function fit is a common method used in finding binary
SMBHs (Graham et al., 2015; Witt et al., 2022). The goal of curve fitting the
data was to characterize the ability to recover the amplitude, phase, and
period of the binary signal from the simulated DRW + binary light curves
sampled with the expected observation pattern for DDFs. We fit the same sine
function to the isolated-quasar (pure-DRW) light curves to estimate the false
positive detection rate of periodic binary variability. Our computationally
inexpensive method could be used as a first “filter” to triage data before
applying multi-component and computationally expensive methods, such as the
Bayesian model selection analysis of Witt et al. (2022).

Figure 3: Examples of simulated LSST Deep Drilling Field light curves, with
and without binaries, and their respective fits. Both light curves have
underlying variability associated with a $10^{8}$ $M_{\odot}$ quasar at z
$\sim$ 1 that is accreting at 10% of its Eddington limit. Top: The best
sinusoidal fit to a non-binary, isolated quasar. Bottom: The best sinusoidal
fit to a binary quasar with a signal amplitude of 0.9 mag, period of 570 days,
and a phase of $\pi$.

We utilized the Python package LMFIT, a non-linear least-squares minimization
and curve-fitting tool for Python that utilizes the Levenberg–Marquardt
algorithm, for the curve fits (Newville et al., 2014). Examples of these fits
to simulated light curves can be seen in Figure 3. We used their built-in
SineModel to fit magnitude $m$ as a function of time $t$ with parameters $A$
(amplitude), $f$ (frequency, $f=2\pi/P$), and $\phi$ (phase):
$m(t)=A\sin\left(ft+\phi\right)$. We do not fit a sawtooth function to the
light curves, only the sinusoid. This is to investigate how much confusion
missing the sawtooth might introduce.

We evaluate goodness-of-fit using the reduced chi-squared statistic
$\chi_{\nu}^{2}\;$. We introduced an error floor of 0.5% of the light curve’s
mean magnitude, to be used as a fit weight, for the reduced chi-squared
evaluation to avoid unreasonable $\chi_{\nu}^{2}\;$values for bright quasars
(that have very small uncertainties from Equation 4). The relative
$\chi_{\nu}^{2}\;$ranges and visual inspection results can be seen in Figure
4. Fitting a sine to a DRW-only light curve gives a better numerical fit than
an DRW+sine because of the multi-component interaction of the two models
(i.e., different damping timescale and binary period values produce light
curves that are not well-fit by a single sine function). Larger amplitude
binary signals, especially sawtooth signals, worsen $\chi_{\nu}^{2}\;$values
overall as well.

Figure 4: Top: The normalized distribution of reduced chi-squared for the sine
fits to the simulated light curves, separated by light curve type. Shaded
region indicates that 50% of the total distribution has a
$\chi_{\nu}^{2}\lesssim 5.0$. Bottom: The distributions of reduced chi-squared
with the DRW amplitude (left) and damping timescale (right). Plotted are the
median y-axis value for the data binned along the x-axis. The error bars shown
are the 16th and 84th percentile values for the binned data. The fits are
worse when the DRW amplitude is strong, but there is not a simple monotonic
relationship with $\tau$. The best-fit sine function is generally a good fit
to the data, especially for light curves with strong binary and/or weak DRW
contributions. We note that the fits for strong binaries are generally “worse”
in terms of $\chi_{\nu}^{2}\;$, this is a consequence of large values in the
sawtoothiness parameter.

For the following analysis, we define a “weak” binary as one with an injected
binary signal amplitude $0.1<A\leq 0.5$ mag and a “strong” binary with
amplitude $0.5<A<1.0$ mag. We note that the “weak” label reflects binary
signals that are actually quite moderate when considering most candidate
binaries have amplitudes near 0.1 mag (Graham et al., 2015). This label is
only introduced for ease of the following analysis. The bottom panels of
Figure 4 demonstrate that the sine functions provide better fits (lower
$\chi_{\nu}^{2}$) to the simulated light curves with weaker DRW amplitudes and
stronger binary signals as expected. The sine function fits are also better
when the DRW damping timescale is longer that the underlying binary period,
likely because it makes the shorter binary period easier to distinguish from
the quasi-periodic DRW variability.

The next section discusses how the rates of binary detections and false
positives depend on quasar and binary properties. We emphasize that these
simple sinusoidal fits are computationally inexpensive: we achieved
$\sim$7,500 light curve fits per minute on a modest 4-core personal laptop.
This makes the fits well-suited to survey data applications like the LSST
event broker333https://www.lsst.org/scientists/alert-brokers ().

##  4 Results

###  4.1 Binary Signal Recovery

Figure 5: Left: Fractional binary parameter recovery performance for amplitude
($A$), period ($P$), and phase ($\phi$) for all binary types. For all
binaries, 28% can have all three parameters recovered within 10% (0.10) of
their input values. The phase spikes are a consequence of possible negative
amplitudes when fitting and phase folding. Middle: 2-D histogram of the period
and amplitude recovery to show any parameter dependence. Period and amplitude
recovery performance are generally uncorrelated. Only 30% of all binaries can
have both the binary signal amplitude and period recovered together. Right:
Period and phase recovery dependence. Where period is likely recovered, the
phase is recovered, especially within 10% of their input values. Over 61% of
all binaries can have their signal phases and periods recovered within 10%.

We first investigate the ability of our sinusoidal fits to recover the input
binary signal. For most parameters we characterize the accuracy of recovery as
a fractional error: (fit parameter output - parameter input) / (parameter
input), i.e. $\Delta{P}/{P_{in}}$. Phase recovery is instead quantified by the
fraction of a period, $\Delta\phi/(2\pi)$. Figure 5 shows the fit recovery of
the input binary period ($P$), amplitude ($A$), and phase ($\phi$). The left
panel shows the distribution of individual parameter recovery: binary period
and phase have similar recovery rates, while amplitude is more poorly
recovered. The right panels of Figure 5 demonstrate that the phase and period
are recovered in tandem: where the period is generally well-recovered, the
phase is as well. The accuracy of amplitude recovery is independent of the
period and phase recovery. This is due to the additional variability of the
DRW.

Period is generally a well-recovered parameter, with 71% of the fits to all
binary light curves recovering the input period with $\Delta{P}/P_{\rm in}<$
0.10 (or 10%). Binary signals with strong (0.5 $<A<$ 1.0) amplitudes can have
their period recovered 84% of the time. Figure 6 shows that binary period is
modestly overestimated by the fits for quasars with larger DRW amplitudes. It
also shows that all binary types are most affected by the DRW amplitude,
rather than the damping timescale.

Figure 6:  Period recovery, binned by DRW parameters (left: damping timescale
$\tau$, right: DRW amplitude $\rm SF_{\infty}$). Plotted here, and in all
similar plots that follow, are the mean y-axis values for the data binned
along the x-axis. The error bars shown are the standard deviations within
those bins. The period is overestimated for weaker binaries (by factor of
1.5-2) because the red noise variability acts to smooth out the periodic
signal, leading to best-fit sinusoids with longer periods. This is seen in
both parameters. The stronger the DRW amplitude, the more overestimated the
best-fit periods become, regardless of binary type.

The binary phase is recovered as well as the period, with a trend of
underestimation. Roughly 69% of all binary light curves had their phase
recovered $\Delta{\phi}/(2\pi)<$ 0.10, with strong binary phase recovery at
81%. Figure 7 shows that again the DRW amplitude has more influence on phase
recovery than the damping timescale. While period and phase recovery are
generally correlated, the phase is less affected by DRW parameters than the
period for larger values of $\rm SF_{\infty}$ and $\tau$. This is evident in
the lack of strong correlations and small standard deviations of
$\Delta\phi/(2\pi)$ with DRW parameters.

Figure 7: Fractional phase recovery means (points) and standard deviations
(error bars) as binned by DRW parameter (left: damping timescale $\tau$,
right: DRW amplitude $\rm SF_{\infty}$). Phase recovery is modestly
underestimated across both DRW parameters. It has the strongest dependence on
$\rm SF_{\infty}$, for all types of binaries.

Binary amplitude is the least effectively recovered parameter, with only 36%
of the fits recovering $\Delta{A}/A_{\rm in}<$ 0.10. Strong binaries see a
higher recovery rate at 44%. Figure 8 shows the binary amplitude recovery when
binned by DRW parameters. Amplitude recovery has the same general trends as
phase and period, with the DRW amplitude having the strongest effect on
recovery. The binary amplitudes are generally overestimated when the DRW
amplitude is larger. The weak binary behavior shows that as the DRW amplitude
becomes many more times the binary signal amplitude, the binary signal becomes
more difficult to recover. Only 28% of all binary light curves can have all
three binary signal parameters recovered at the same time to within 10% of
their input values. When the binary amplitude is near or larger than the DRW
amplitude (e.g., strong binaries with $A>0.5$ and $\rm SF_{\infty}$ $<0.5$),
recovery for all parameters is better, as similarly pointed out by Witt et al.
(2022). In general the period and phase can be recovered even when the
amplitude cannot.

Figure 8: Fractional recovery of binary amplitude, binned by the DRW
parameters (left: damping timescale $\tau$, right: DRW amplitude $\rm
SF_{\infty}$). The binary amplitude is effectively recovered for light curves
with strong binary signals (binary amplitude $A>0.5$ mag), with little
dependence on the DRW signal. Light curves with weak binary signals generally
have best-fit amplitudes that are significantly overestimated, especially for
large DRW amplitude.

Next, we investigate the effectiveness of the binary signal parameter recovery
as a function of binary signal structure. Figure 9 depicts the recovery
performance of each parameter as a function of sawtoothiness. The sinusoidal
fits still accurately recover the parameters even if the light curve is quite
“sawtoothy.” As discussed above with the fit results versus DRW parameters,
the amplitude and period are again overestimated and phase is underestimated,
with weaker binaries displaying larger deviations for all three parameters.
Larger values of sawtoothiness drive the best-fit values of amplitude and
phase lower compared to their smoother sine counterparts. However, these
differences are small at 10-20%. The DRW parameters have a stronger
correlation with binary parameter recovery than the binary signal structure
type.

Figure 9: Parameter recovery compared to binary signal “sawtoothiness” (0.0 is
a smooth sine, 1.0 is a pure sawtooth). From top to bottom, panels show the
parameter recovery for binary amplitude $A$, period $P$, and phase $\Phi$ for
weak (red points) and strong (blue) binaries. In all cases, the stronger the
binary signal, the better the recovery is. There is interesting behavior at
large values of sawtoothiness: the period and amplitude become more over- and
under- estimated, respectively, by 10-20 % compared to their smoother sine
counterparts. Despite this, DRW parameters have a stronger influence on
recovery than signal type/structure.

###  4.2 Non-Binary Quasar Fit Results

We next investigate how the single, isolated quasars performed when fit with
the periodic functions. Figure 4 shows that more than 80% of the single quasar
light curves have a best-fit model with $\chi_{\nu}^{2}\;$$<5.0$. We further
explore what the best-fit results are for these light curves. Figure 10
depicts the best-fit amplitude versus $\chi_{\nu}^{2}\;$, DRW amplitude, and
damping timescale for the single (non-binary) quasar light curves. The fit
amplitudes directly correlate to DRW amplitude, with slight underestimation.
The sinusoidal fits have a difficult time fitting larger $\rm SF_{\infty}$, as
most light curves with a $\chi_{\nu}^{2}\;$$\leq$ 5.0 have best-fit amplitudes
0.1 $<\rm A<$ 0.6. From Figure 2 and the results of the previous section, most
of the quasars with large structure functions would obscure a binary signal
and make recovery of the underlying parameters difficult.

Amplitude has a clear dependence with DRW parameters, as seen from Figure 10,
and while best-fit period has no relationship with the DRW parameters, it is a
non-uniform distribution (see contours of Figure 11). Phase was not
investigated further as its distribution was uniform and showed no dependence
across $\chi_{\nu}^{2}\;$, $\rm SF_{\infty}$, and $\tau$.

Investigating the best-fit binary models for non-binary quasars is important
because single quasars are likely to vastly outnumber binary quasars (Kelley
et al., 2019). If binary quasars are selected to have $\rm A>$ 0.1, as
reflected by the candidate populations in previous surveys (e.g. Graham et al.
2015), then putative samples of binary quasars will include (or be dominated
by) a large number of false positives. More discussion of false positive rates
and population demographics can be found in Section 5.2.

Figure 10: Best-fit amplitude results for single, isolated quasars. Left:
Single quasars that are well-fit by the sinusoidal binary light curve model
(with low $\chi_{\nu}^{2}\;$) generally have small best-fit amplitudes.
Middle: Best-fit amplitude and $\rm SF_{\infty}$have a positive, linear
correlation: the stronger the DRW amplitude, the larger the sinusoidal fit
amplitude. Right: The is no significant connection between fit amplitude and
damping timescale, $\tau$.  Figure 11: Best-fit amplitudes and periods of
sinusoidal binary model fits to non-binary quasars (left panel), shown in the
contours with false positive rates shown in the underlying 2-D histogram. The
best-fit amplitudes are small and the periods have a bimodal distribution with
peaks around 900 and 1800 days with preference for a longer period as a
consequence of the length of our light curves. The upper bound on the false
positive rate is also tied to light curve and survey duration. Right: the
best-fit amplitudes and periods of our binary light curves shown in the
contours, with recovery rates (using period recovery) of the full binary
population shown in 2-D histogram. Section 5.2 defines our false positive
threshold. The white dashed lines in both plots indicate the duration of our
light curves: five years. The stars are example binary candidates, discussed
in Section 5.3.

##  5 Discussion

In the following section, we discuss the general trends seen in the data. We
further define what kinds of binary SMBH detections can be trusted and
introduce false positive detection rates, provide recommendations for binary
SMBH searches from Rubin/LSST light curves, and discuss the multi-messenger
connections of these results.

###  5.1 Implications for True Binary Recovery

Our analysis of binary parameter recovery establishes that fitting simple (and
computationally inexpensive) sinusoidal binary light curve models usually
recovers the period and phase to within 10% of the input value (Figure 5). We
recover the period from 84% of our strong binary light curves and 71% of all
binary light curves. Figure 6 establishes that smaller values of damping
timescale and DRW amplitude make it easier to accurately recover the binary
signal. The phase can be recovered for 81% of all strong binaries and 69%
across all binaries. Phase is typically underestimated and most affected by
DRW amplitude. Binary signal amplitude is not effectively recovered for most
light curves, as only 44% of strong binaries and 36% of all binaries can have
their amplitudes accurately recovered within 10% of the input value. Both
phase and period can be recovered simultaneously for 62% of all binaries.
Sinusoidal curve fits can only recover all three binary parameters for 28% of
all binary light curves.

In general, the “stronger” (higher amplitude) the binary signal, the easier it
is to disentangle from a DRW. This remains true even if the DRW amplitude is
large. This is aligned with the results of Witt et al. (2022). Weaker binaries
are more affected by strong DRW parameters and best-fit amplitude especially
is overestimated because of the additional variability of a strong DRW signal.
The true binary amplitude will be smaller than measured for most binary
quasars by about a factor of 1.5-2 because of the contribution of DRW
variability. While strong binaries often have their binary signal parameters
recovered within 10% of the true value, we note that most binary candidates in
the literature (e.g., Graham et al., 2015; Charisi et al., 2016) have
amplitudes that fall within our “weak” binary category.

Simple, computationally inexpensive fits are ambivalent to sawtoothiness for a
wide range of binary signal parameters, as in the binary parameters can still
be recovered. The contribution from the DRW matters more than the details of
the binary signal. This is encouraging news as we do not know the true nature
of these binary signals (see Section 2.2 for description of the different
types and emission mechanisms of electromagnetic binary signals). The accuracy
of period and phase recovery is promising as they can be directly linked to
the gravitational wave signal parameters (Charisi et al., 2022). Section 5.4
offers more discussion on joint EM and GW endeavours.

###  5.2 False Positive Detection and Reliability Thresholds

Figure 12: 2-D histograms of the false-positive demographics, contours (white)
of the full non-binary population, and stars of five example binary candidates
from Table 3. Left: The luminosity and mass grid of false positives. We find
that the most luminous and massive quasars have the highest rates of false
positive detection. Right: False-positive rates in terms of DRW parameters,
damping timescale $\tau$ and DRW amplitude $\rm SF_{\infty}$. We having that
small $\tau$ and $\rm SF_{\infty}$drive the false positive rate, with the
light curves made of those values reporting upwards of 60% false positive
rates. Figure 13: Reliability thresholds: the best-fit amplitudes and periods
of false positive binaries and the “recovered” fraction of true positive
binaries (recovered, here, is defined as having the binary signal periods
recovered within 10% of their input value via the sinusoidal fits). Also
plotted are the Eddington ratio distributions for false positives and
recovered true positives. At small best-fit amplitudes, the false positive
detections dominate. The stronger the binary signal, the more likely we obtain
a true positive. False positive rates are highest for quasars with high
Eddington ratio.

We next address the rate of false-positive detections from the simulated
single-quasar (non-binary) light curves. For this purpose we define a binary
detection as:

  * •

good fit by the sinusoidal light curve model, $\chi_{\nu}^{2}\;$$<5.0$

  * •

best-fit binary amplitude $A>0.1$ mag

  * •

best-fit binary period $180<P<1100$ days

The $A>0.1$ requirement is motivated by the typical amplitudes of binary SMBH
candidates identified from previous CRTS observations by Graham et al. (2015).
The typical amplitudes of a binary candidate population in a survey are driven
by the photometric error associated with the survey. Compared to CRTS, LSST
will have lower photometric error and, thus, lower resolvable binary
amplitudes, but we leave this threshold as the requirement to be conservative.
The $\chi_{\nu}^{2}\;$$<5.0$ requirement indicates that the light curve is
well-fit by a sinusoidal model (this requirement is analogous to the auto-
correlation tests of Graham et al. 2015). The range of $180<P<1100$ days
requires that the period is well-detected (1.5 period cycles or more) within 5
years of LSST monitoring.

The left panel of Figure 11 presents the distribution of best-fit period and
amplitude for the sinusoidal fits to the non-binary quasars. For comparison,
the right panel of Figure 11 presents the period and amplitude parameter space
for reliable recovery (defined as $|\Delta P/P|<$ 0.10) of the simulated true-
binary light curves. Many of the non-binary light curves are best fit by
periods of $\sim$5 yr, comparable to the survey duration and rejected by our
period requirement for a binary detection. This is in contrast to the shorter
periods, by construction, of the simulated true binary light curves. Most of
the non-binary light curves are best fit by low-amplitude sinusoid models,
again in contrast to the true binary population having most accurate recovery
at high amplitude. The false-positive rate is very high for low to moderate
amplitudes. In other words, if binary candidates are selected to have well-fit
sinusoid models with amplitudes of $A>0.1$ mag, then false-positive non-binary
quasars will dominate the sample.

Figure 12 presents the false-positive rate as a function of quasar luminosity
and mass (left panel) and DRW parameters (right panel). The false-positive
rate is highest for luminous quasars accreting near their Eddington Limit,
which tend to have smaller DRW amplitudes and modestly variable light curves
that can be well-fit by periodic signals. Between 30 and 40% of quasars with
$\log(M_{\rm BH}/M_{\odot})>$ 9 and $\log(L/(erg/s))>$ 47 are identified as
false-positive binaries according to the selection criteria described above.
For luminous quasars near a mass of $\log(M_{\rm BH}/M_{\odot})\sim$ 7.5, the
false positive rate is over 50%. Witt et al. (2022) demonstrate that more
computationally intensive approaches can better recover binary parameters by
disentangling the DRW contribution, although that work still finds a
dependence of binary recovery on DRW amplitude.

Figure 13 compares the fraction of false-positive detections (black lines)
with the recovery fraction of true-positive binaries (dotted red lines), with
binary “recovery” defined as $|\Delta P/P|<$ 10% as before. In Figure 13, we
see that at best-fit amplitudes below 0.3 mag, the false positive binaries
values dominate over the true positive binaries. The higher the amplitude fit,
the more likely it is that the binary is real. Single, isolated quasars have a
difficult time producing a fitted amplitude above 0.4 mag and a decent
$\chi_{\nu}^{2}\;$statistic. We establish this as the reliability threshold
for trusting a binary quasar detection via sine fits, where we can begin to
trust strong binary signal amplitudes greater than 0.4 mag (see Section 5.4
for details of physical binary parameters for this value). We also note that a
sine fit to quasar light curves will always return a non-zero amplitude value
with the small predicted photometric errors of LSST. Fits with smaller periods
are more likely to be true-positive binaries, as shown in the middle panel and
also because almost half of the non-binary light curves are best fit by
periods comparable to the survey duration (left panel of Figure 11). Quasars
with lower Eddington ratio are also less likely to identified as false-
positive binaries, with a factor of $\sim$2 decrease in the false-positive
rate from $L_{\rm bol}/L_{\rm Edd}$ $\sim$0.1 to $\sim$0.01. In contrast the
rate of true-positive recovery decreases more modestly with Eddington ratio,
declining by a factor of only $\sim$1.5 over the same range.

Our overall conclusions are that false positives will dominate samples of
binary candidates found via sinusoidal curve fits for best-fit light curves of
small amplitude ($0.1<A<0.2$ mag), large period (especially if near the survey
length), and for luminous and massive quasars (with smaller stochastic
accretion-disk variability). Binary candidates with larger amplitude ($A>0.5$
mag), smaller period, and lower Eddington ratio are more likely to be genuine
binary quasars.

Figure 14: 2-D histograms of the true positive candidate demographics,
contours (white) of the full binary population, and stars of five example
binary candidates from Table 3. Left: The luminosity and mass grid of true
positives. We find that the most luminous and massive have the highest rates
of true positive detection. Right: True positives in terms of DRW parameters,
damping timescale $\tau$ and DRW amplitude $\rm SF_{\infty}$. Light curves
made of small values of the DRW parameters report 60% true positive detection
rates. We see similar trends in Figure 12 for the false positive population.

If we apply the same definition of binary detection to our true binary light
curves, we find that only 38% of all binary light curves are considered true
positives. The low true-positive rate is a consequence of requiring
$\chi_{\nu}^{2}\;$$<5.0$: about half of the binary light curves do not have
best-fit binary models below this threshold, as shown in Figure 4 and
discussed in Section 3. Figure 14 depicts the demographics of the true
positive population. The true positive detection rates peak at about 60% for
high-mass and high-luminosity quasars that have small DRW parameters. These
trends mirror that of the false positive detections for this method. We also
emphasize that a non-negligible number of binary light curves produce fit
periods near the light curve baseline of five years, as seen in the right
panel of Figure 11. Most of these are binaries with “weak” binary amplitudes
($A<0.5$), with small $\rm SF_{\infty}$and over the full range of potential
$\tau$ and binary periods.

###  5.3 Recommendations for Rubin

We recommend exercising caution when using sinusoidal fits for binary SMBH
detection. The very nature of pipeline-style sine fits and binary signals (in
combination with the natural noise properties of the quasars we observe), will
result in false positive detections. Even with analysis methods that model the
quasar variability, this point holds (Witt et al., 2022). If a sinusoidal fit
is to be applied to quasar light curves in search of periodic signals, we can
generally rely on the period and phase being accurate to within 10%. This is
especially true for binary signals that are strong, periods that are many
times the damping timescale of the variability, and for luminous quasars that
have smaller damping timescales (less massive) or structure functions (more
massive). The more luminous a quasar is, the more likely it is to be a false
positive detection due to its inherent variability properties. Up to 60% of
these quasars return false positive results across various black hole masses.

If the returned fit amplitude is larger than 0.4 mag, it’s more likely to be a
true binary than a false positive, as discussed in the previous subsection.
There is large uncertainty associated with amplitude parameter recovery and
the true underlying amplitude will usually be smaller than the best-fit value.
Especially for weaker amplitude binary signals, the DRW amplitude strongly
influences the fit and parameter recovery.

We see some correlation with DRW-only fit periods and length of light curves
as discussed in Section 4.2 and seen in the bimodal distribution of single
quasar fit periods in Figure 11 (left panel, the white dashed line denotes our
five-year baseline). If the best-fit period for a five-year light curve is
near five or two and a half years, then it should be treated with caution and
as a potential false positive that needs further analysis.

We can apply our recommendations to current binary SMBH candidates, with the
caveat that they are candidates from different surveys, done using analysis
methods that account for quasar variability. LSST is a higher-quality survey,
so it would produce the lower limit of false positives. An example of the
application of these recommendations is to attempt to apply them to long-
reigning binary candidates like PG1302-102 Graham et al. (2015). Assuming that
when observed with Rubin/LSST, its calculated mass of
log10($M_{BH}/M_{\odot}$) = 8.3-9.4, its luminosity near its Eddington limit,
and binary amplitude of 0.14 mag from its CRTS data would remain relatively
unchanged, this candidate would be flagged as a likely false-positive
candidate. Table 3 contains four more examples of candidates in addition to
PG1302-102 and our recommendation application, these candidates are plotted as
stars on all panels on Figures 11 and 12. It is notable that none of the best-
fit CRTS candidates have amplitudes above our threshold $A>0.4$ for lower
false-positive rates. We note that Graham et al. (2015) used auto-correlation
functions as well as sinusoid fits to identify the binary candidates shown in
Figures 11 and 12; our results demonstrate that fitting only by periodic
sinusoids will result very high false-positive rates for binary candidates
with similar properties.

Candidate (marker color) | Amplitude [mag] | Period [days] |  $log_{10}$ (M/M⊙) |  $log_{10}$ (Lbol/[ergs/s])
---|---|---|---|---
PG1302-102 (white) | 0.14 | 1, 884 $\pm$ 88 | 8.3-9.4 | 46.694
SDSS J143621+0727 (purple) | 0.08 | 1886 | 9.03 | 46.740
SNU J13120+0641 (cyan) | 0.09 | 1492 | 9.14 | 46.213
SDSS J080648+1840 (green) | 0.2 | 892 | 7.99 | 45.173
SDSS J114749+1631 (gray) | 0.15 | 1449 | 8.22 | 45.726
Table 3: Candidates from Graham et al. (2015) that approximately fall within
our parameter space. The candidate amplitudes are estimated values from their
respective light curves. The marker colors of the candidates are used in
Figures 11, 12, and 14.

###  5.4 Connections to Multi-Messenger Astrophysics

An exciting prospect of electromagnetic detections of binary SMBHs are their
connection to gravitational wave detection and multi-messenger astrophysics.
Our work shows that the binary phase and period are well-recovered in LSST DDF
light curves. Our work also shows a high false-positive rate for photometric
detections of binary SMBH, thus they are not sufficient on their own. Rubin
will be excellently poised for follow-up for binary candidates found via
gravitational waves and in other wavelengths across the electromagnetic
spectrum.

Our results are complementary to constraints on binary candidate populations.
Casey-Clyde et al. 2023 (submitted) has shown that the most luminous
candidates are also likely false positives. We also find that high luminosity
quasars can mimic binary SMBHs when fitting with sine waves. This places more
tension on the current lists of candidates, in which high luminosity sources
are well-represented. Despite this tension, we should continue to observe
these quasars as binary candidate population of N quasars goes as $N\approx
T_{obs}^{8/3}$.

We can further connect our recommended amplitude threshold of 0.4 mag to
examples of physical binary parameters and to gravitational wave sensitivity
with Pulsar Timing Arrays (PTAs). The change in apparent magnitude due to
Doppler boost variability is directly tied to the line-of-sight velocity of
the secondary black hole. By following the prescriptions laid out in Charisi
et al. (2022), we can determine that a secondary black hole would need to be
moving between 6.5% and 10% of the speed of light to produce that binary
amplitude. These lower and upper bounds are determined by the redshift of the
candidate. With a period detectable within a five-year LSST light curve (e.g.,
periods of 0.5 to 5.0 yr), these binary candidates could have a total masses
between $10^{9.0}$ and $10^{10.6}$ $M_{\odot}$. With a period, distance, and
total mass, we can determine if these candidates would lie within current and
future PTAs sensitivity ranges. Figure 6 of Charisi et al. (2022) is an
example of the overlap of EM and GW binary SMBH parameter space. Binaries with
an orbital period between 180 and 1,825 days and masses of $10^{9.0}$ and
$10^{10.6}$ $M_{\odot}$ would fall within potential PTA coverage. Keeping all
other parameters equal, a larger change in apparent magnitude corresponds to a
larger line-of-site velocity and a larger total binary mass, inferring higher
sensitivity to the GW emission.

##  6 Conclusion and Future Work

We have presented our work that included the generation and sine fitting of
over 3.6 million Rubin/LSST quasar light curves, many of which contain binary
SMBH signals. We conservatively estimate that over 40% of isolated, single
quasars will result in a false positive detection of a binary SMBH system with
a simple sinusoidal fit. The most luminous quasars in our study are the most
likely candidates for a false positive detection with this method, with the
false positive rate of these objects, specifically, being over 50%. They are
also, conversely, the binary types that produce the best recovery of our input
parameters, as the structure function of the variability is low for these
types of quasars. Less massive, less luminous quasars are more likely to
dilute a binary signal because of their variability properties (namely large
DRW amplitudes).

From our work, we provide the following recommendations for Rubin/LSST binary
SMBH searches and identification:

  * •

Periods and phases of binary signal can reliably be recovered within 10% of
the input binary signal value via sinusoidal curve fits roughly 70% of the
time, separately, or 62%, together.

  * •

Binary amplitude is the least-reliable parameter to recover, but best-fit
amplitudes larger than 0.4 magnitudes can be used to boost confidence in a
binary candidate.

  * •

Best-fit periods that are near the length or half-length of their respective
light curves should be treated with caution.

  * •

Sources that are massive and luminous are more likely to be a false positive
than a real binary if the DRW is not accounted for.

###  6.1 Future Work

The purpose of this paper was to explore a the detectability of binary SMBHs
for a representative quasar population in Rubin/LSST DDF observations. In the
future, the light curve generation could be adapted to use parameter
distributions that are more directly tied to observed quasar population and/or
limits on binary SMBHs inferred from current gravitational wave searches.

Beyond light curve generation, there are more analysis methods we plan to
apply to the light curves, like Lomb-Scargle periodograms and Auto-Correlation
Functions (Lomb, 1976; Scargle, 1982; VanderPlas, 2018; Edelson & Krolik,
1988; McQuillan et al., 2013) and machine learning algorithms. We hope to
revisit this analysis with the development and inclusion of computationally
inexpensive DRW fits, in addition to the sine fits. This could result in more
effective triage of false positives. There is more room to explore how these
results differ with different cadences (like main LSST survey WFD cadence
versus the DDFs) and longer baseline light curves (the full 10-year LSST
baseline and beyond).

We also plan to apply these analysis pipelines to existing large-scale
surveys, like the Catalina Real-time Transient Survey, as done in Graham et
al. (2015) and Witt et al. (2022). Our methods can also be applied to make
predictions for other upcoming time-domain surveys, like the High Latitude
Time Domain Survey of the Nancy Grace Roman Space Telescope (Roman) which is
expected to launch in 2027 (at the time of writing this paper Spergel et al.,
2013, 2015). Roman may also be able to detect massive binary systems through
its observation of $10^{4}$ \- $10^{6}$ quasars (Haiman et al., 2023; Shen et
al., 2023).

JRT, MCD, KEG, and AH acknowledge support from NSF grant CAREER-1945546. MCD
acknowledges support from a NSF Graduate Research Fellowship. LB acknowledges
support from Cottrell Scholar Award #27553 from the Research Corporation for
Science Advancement and from NSF award AST-1909933. WNB acknowledges financial
support from NSF grant AST-2106990. JACC is supported by the National Science
Foundation’s NANOGrav Physics Frontier Center, Award Number 2020265 and by the
National Science Foundation under Grants NSF PHY-2020265, and AST-2106552. MC
and JCR acknowledge support from NSF award AST-2007993. The authors wish to
thank London E Willson and Chiara MF Mingarelli for their helpful discussion
and insight. MCD wishes to especially thank Dani Lipman and Hugh Sharp for
their support.

## References

  * Agazie et al. (2023a) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. (2023b) —. 2023b, ApJ, 952, L37, doi: 10.3847/2041-8213/ace18b
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, Laser Interferometer Space Antenna, arXiv, doi: 10.48550/ARXIV.1702.00786
  * Amaro-Seoane et al. (2023) Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Living Reviews in Relativity, 26, 2, doi: 10.1007/s41114-022-00041-y
  * Antoniadis et al. (2023a) Antoniadis, J., Babak, S., Nielsen, A.-S. B., et al. 2023a, Astronomy & Astrophysics, 678, A48
  * Artymowicz & Lubow (1996) Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, L77, doi: 10.1086/310200
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2020, ApJ, 900, 102, doi: 10.3847/1538-4357/ababa1
  * Bansal et al. (2017) Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ, 843, 14, doi: 10.3847/1538-4357/aa74e1
  * Barth & Stern (2018) Barth, A. J., & Stern, D. 2018, The Astrophysical Journal, 859, 10, doi: 10.3847/1538-4357/aab3c5
  * Bauer et al. (2009) Bauer, A., Baltay, C., Coppi, P., et al. 2009, ApJ, 696, 1241, doi: 10.1088/0004-637X/696/2/1241
  * Bianco et al. (2022) Bianco, F. B., Ivezić, Ž., Jones, R. L., et al. 2022, ApJS, 258, 1, doi: 10.3847/1538-4365/ac3e72
  * Blumenthal & Barnes (2018) Blumenthal, K. A., & Barnes, J. E. 2018, MNRAS, 479, 3952, doi: 10.1093/mnras/sty1605
  * Bogdanović et al. (2022) Bogdanović, T., Miller, M. C., & Blecha, L. 2022, Living Reviews in Relativity, 25, 3, doi: 10.1007/s41114-022-00037-8
  * Brandt et al. (2018) Brandt, W. N., Ni, Q., Yang, G., et al. 2018, Active Galaxy Science in the LSST Deep-Drilling Fields: Footprints, Cadence Requirements, and Total-Depth Requirements.  https://arxiv.org/abs/1811.06542
  * Charisi et al. (2016) Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145, doi: 10.1093/mnras/stw1838
  * Charisi et al. (2018) Charisi, M., Haiman, Z., Schiminovich, D., & D'Orazio, D. J. 2018, Monthly Notices of the Royal Astronomical Society, 476, 4617, doi: 10.1093/mnras/sty516
  * Charisi et al. (2022) Charisi, M., Taylor, S. R., Runnoe, J., Bogdanovic, T., & Trump, J. R. 2022, MNRAS, 510, 5929, doi: 10.1093/mnras/stab3713
  * Chen et al. (2022) Chen, Y.-C., Hwang, H.-C., Shen, Y., et al. 2022, ApJ, 925, 162, doi: 10.3847/1538-4357/ac401b
  * Chen et al. (2020) Chen, Y.-C., Liu, X., Liao, W.-T., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 2245, doi: 10.1093/mnras/staa2957
  * Chen et al. (2022b) Chen, Y.-J., Zhai, S., Liu, J.-R., et al. 2022b, arXiv e-prints, arXiv:2206.11497, doi: 10.48550/arXiv.2206.11497
  * Comerford et al. (2018) Comerford, J., Nevin, R., Stemo, A., et al. 2018, The Astrophysical Journal, 867, 66, doi: 10.3847/1538-4357/aae2b4
  * De Rosa et al. (2019) De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New A Rev., 86, 101525, doi: 10.1016/j.newar.2020.101525
  * de Vries et al. (2005) de Vries, W. H., Becker, R. H., White, R. L., & Loomis, C. 2005, AJ, 129, 615, doi: 10.1086/427393
  * D’Orazio & Charisi (2023) D’Orazio, D. J., & Charisi, M. 2023, arXiv e-prints, arXiv:2310.16896, doi: 10.48550/arXiv.2310.16896
  * D'Orazio et al. (2013) D'Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2997, doi: 10.1093/mnras/stt1787
  * D’Orazio et al. (2013) D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2997, doi: 10.1093/mnras/stt1787
  * D'Orazio et al. (2015a) D'Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015a, Nature, 525, 351, doi: 10.1038/nature15262
  * D'Orazio et al. (2015b) —. 2015b, Nature, 525, 351, doi: 10.1038/nature15262
  * Drake et al. (2009) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870, doi: 10.1088/0004-637X/696/1/870
  * Duffell et al. (2020) Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020, Circumbinary Disks: Accretion and Torque as a Function of Mass Ratio and Disk Viscosity.  https://arxiv.org/abs/1911.05506
  * Edelson & Krolik (1988) Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646, doi: 10.1086/166773
  * Farris et al. (2014) Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014, The Astrophysical Journal, 783, 134, doi: 10.1088/0004-637X/783/2/134
  * Giveon et al. (1999) Giveon, U., Maoz, D., Kaspi, S., Netzer, H., & Smith, P. S. 1999, MNRAS, 306, 637, doi: 10.1046/j.1365-8711.1999.02556.x
  * Graham et al. (2015) Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Monthly Notices of the Royal Astronomical Society, 453, 1562–1576, doi: 10.1093/mnras/stv1726
  * Haiman et al. (2023) Haiman, Z., Xin, C., Bogdanović, T., et al. 2023, Massive Black Hole Binaries as LISA Precursors in the Roman High Latitude Time Domain Survey.  https://arxiv.org/abs/2306.14990
  * Hawkins (2002) Hawkins, M. R. S. 2002, MNRAS, 329, 76, doi: 10.1046/j.1365-8711.2002.04939.x
  * Hook et al. (1994) Hook, I. M., McMahon, R. G., Boyle, B. J., & Irwin, M. J. 1994, MNRAS, 268, 305, doi: 10.1093/mnras/268.2.305
  * Ivezić (2017) Ivezić, Ž. 2017, in New Frontiers in Black Hole Astrophysics, ed. A. Gomboc, Vol. 324, 330–337, doi: 10.1017/S1743921316012424
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, The Astrophysical Journal, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Kelley et al. (2021) Kelley, L. Z., D’Orazio, D. J., & Di Stefano, R. 2021, MNRAS, 508, 2524, doi: 10.1093/mnras/stab2776
  * Kelley et al. (2019) Kelley, L. Z., Haiman, Z., Sesana, A., & Hernquist, L. 2019, Monthly Notices of the Royal Astronomical Society, 485, 1579, doi: 10.1093/mnras/stz150
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kelly et al. (2014) Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., & Uttley, P. 2014, ApJ, 788, 33, doi: 10.1088/0004-637X/788/1/33
  * Kelly et al. (2010) Kelly, B. C., Vestergaard, M., Fan, X., et al. 2010, The Astrophysical Journal, 719, 1315, doi: 10.1088/0004-637x/719/2/1315
  * Kollmeier et al. (2006) Kollmeier, J. A., Onken, C. A., Kochanek, C. S., et al. 2006, The Astrophysical Journal, 648, 128, doi: 10.1086/505646
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Koss et al. (2023) Koss, M. J., Treister, E., Kakkad, D., et al. 2023, The Astrophysical Journal Letters, 942, L24, doi: 10.3847/2041-8213/aca8f0
  * Kozłowski (2017a) Kozłowski, S. 2017a, A&A, 597, A128, doi: 10.1051/0004-6361/201629890
  * Kozłowski (2017b) —. 2017b, ApJ, 835, 250, doi: 10.3847/1538-4357/aa56c0
  * Kozłowski et al. (2010) Kozłowski, S., Kochanek, C. S., Udalski, A., et al. 2010, ApJ, 708, 927, doi: 10.1088/0004-637X/708/2/927
  * Kozłowski (2021) Kozłowski, S. 2021, Acta Astronomica, 71, 103–112, doi: 10.32023/0001-5237/71.2.2
  * Liu & Vigeland (2021) Liu, T., & Vigeland, S. J. 2021, ApJ, 921, 178, doi: 10.3847/1538-4357/ac1da9
  * Liu et al. (2019) Liu, T., Gezari, S., Ayers, M., et al. 2019, ApJ, 884, 36, doi: 10.3847/1538-4357/ab40cb
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, LSST Science Book, Version 2.0.  https://arxiv.org/abs/0912.0201
  * MacFadyen & Milosavljević (2008) MacFadyen, A. I., & Milosavljević, M. 2008, The Astrophysical Journal, 672, 83, doi: 10.1086/523869
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž ., Kochanek, C. S., et al. 2010, The Astrophysical Journal, 721, 1014, doi: 10.1088/0004-637x/721/2/1014
  * MacLeod et al. (2012) MacLeod, C. L., Ivezić, Ž., Sesar, B., et al. 2012, ApJ, 753, 106, doi: 10.1088/0004-637X/753/2/106
  * Matheson et al. (2021) Matheson, T., Stubens, C., Wolf, N., et al. 2021, AJ, 161, 107, doi: 10.3847/1538-3881/abd703
  * McQuillan et al. (2013) McQuillan, A., Aigrain, S., & Mazeh, T. 2013, MNRAS, 432, 1203, doi: 10.1093/mnras/stt536
  * Miranda et al. (2017) Miranda, R., Muñoz, D., & Lai, D. 2017, Monthly Notices of the Royal Astronomical Society, 466, 1170, doi: 10.1093/mnras/stw3189
  * Moreno et al. (2019) Moreno, J., Vogeley, M. S., Richards, G. T., & Yu, W. 2019, PASP, 131, 063001, doi: 10.1088/1538-3873/ab1597
  * Mushotzky et al. (2011) Mushotzky, R. F., Edelson, R., Baumgartner, W., & Gandhi, P. 2011, The Astrophysical Journal, 743, L12, doi: 10.1088/2041-8205/743/1/l12
  * Narayan et al. (2018) Narayan, G., Zaidi, T., Soraisam, M. D., et al. 2018, ApJS, 236, 9, doi: 10.3847/1538-4365/aab781
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, 0.8.0, Zenodo, doi: 10.5281/zenodo.11813
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, Astronomy & Astrophysics, 641, A6, doi: 10.1051/0004-6361/201833910
  * Reardon et al. (2023a) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023a, The Astrophysical Journal Letters, 951, L6, doi: 10.3847/2041-8213/acdd02
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49, doi: 10.1086/504825
  * Saha et al. (2016) Saha, A., Wang, Z., Matheson, T., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, 99100F, doi: 10.1117/12.2232095
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Scolnic et al. (2018) Scolnic, D. M., Lochner, M., Gris, P., et al. 2018, arXiv e-prints, arXiv:1812.00516, doi: 10.48550/arXiv.1812.00516
  * Sesana et al. (2018) Sesana, A., Haiman, Z., Kocsis, B., & Kelley, L. Z. 2018, The Astrophysical Journal, 856, 42, doi: 10.3847/1538-4357/aaad0f
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, The Astrophysical Journal Supplement Series, 194, 45, doi: 10.1088/0067-0049/194/2/45
  * Shen et al. (2023) Shen, Y., Casey-Clyde, J. A., Chen, Y.-C., et al. 2023, Discovery and Characterization of Galactic-scale Dual Supermassive Black Holes Across Cosmic Time.  https://arxiv.org/abs/2306.15527
  * Spergel et al. (2013) Spergel, D., Gehrels, N., Breckinridge, J., et al. 2013, WFIRST-2.4: What Every Astronomer Should Know.  https://arxiv.org/abs/1305.5425
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report.  https://arxiv.org/abs/1503.03757
  * Suberlak et al. (2021) Suberlak, K. L., Ivezić, Ž., & MacLeod, C. 2021, ApJ, 907, 96, doi: 10.3847/1538-4357/abc698
  * Tang et al. (2018) Tang, Y., Haiman, Z., & MacFadyen, A. 2018, Monthly Notices of the Royal Astronomical Society, 476, 2249, doi: 10.1093/mnras/sty423
  * Uomoto et al. (1976) Uomoto, A. K., Wills, B. J., & Wills, D. 1976, AJ, 81, 905, doi: 10.1086/111969
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
  * Van Wassenhove et al. (2012) Van Wassenhove, S., Volonteri, M., Mayer, L., et al. 2012, ApJ, 748, L7, doi: 10.1088/2041-8205/748/1/L7
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692, doi: 10.1086/380563
  * Vanderplas et al. (2012) Vanderplas, J., Connolly, A., Ivezić, Ž., & Gray, A. 2012, in Conference on Intelligent Data Understanding (CIDU), 47 –54, doi: 10.1109/CIDU.2012.6382200
  * VanderPlas (2018) VanderPlas, J. T. 2018, The Astrophysical Journal Supplement Series, 236, 16, doi: 10.3847/1538-4365/aab766
  * Vaughan et al. (2016) Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145, doi: 10.1093/mnras/stw1412
  * Vestergaard & Osmer (2009) Vestergaard, M., & Osmer, P. S. 2009, The Astrophysical Journal, 699, 800, doi: 10.1088/0004-637x/699/1/800
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Westernacher-Schneider et al. (2023) Westernacher-Schneider, J. R., Zrake, J., MacFadyen, A., & Haiman, Z. 2023, arXiv e-prints, arXiv:2307.01154, doi: 10.48550/arXiv.2307.01154
  * Wilhite et al. (2008) Wilhite, B. C., Brunner, R. J., Grier, C. J., Schneider, D. P., & Vanden Berk, D. E. 2008, Monthly Notices of the Royal Astronomical Society, 383, 1232, doi: 10.1111/j.1365-2966.2007.12655.x
  * Witt et al. (2022) Witt, C. A., Charisi, M., Taylor, S. R., & Burke-Spolaor, S. 2022, The Astrophysical Journal, 936, 89, doi: 10.3847/1538-4357/ac8356
  * Wu & Shen (2022) Wu, Q., & Shen, Y. 2022, The Astrophysical Journal Supplement Series, 263, 42, doi: 10.3847/1538-4365/ac9ead
  * Xin & Haiman (2021) Xin, C., & Haiman, Z. 2021, Monthly Notices of the Royal Astronomical Society, 506, 2408, doi: 10.1093/mnras/stab1856
  * Yu et al. (2022) Yu, W., Richards, G. T., Vogeley, M. S., Moreno, J., & Graham, M. J. 2022, ApJ, 936, 132, doi: 10.3847/1538-4357/ac8351
  * Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, The Astrophysical Journal, 765, 106, doi: 10.1088/0004-637X/765/2/106


120 Tbit/s Transmission over Single Mode Fibre using Continuous 91 nm Hybrid Raman-EDFA Amplification

L Galdino1,*,
A Edwards2,
M Ionescu2,
J James2,
W Pelouch3,
E Sillekens1,
D Semrau1,
D Lavery1,
R I Killey1,
S Barnes2,
P Bayvel1,
and S Desbruslais2

1Optical Networks Group, Dept. of Electronic & Electrical Engineering, University College London, UK

2Xtera, Bates House, Church Road, Harold Wood, Essex, UK

3Xtera, 500 West Bethany Drive,
Allen,Texas, USA

Hybrid distributed Raman-EDFA amplifiers, with a continuous 91 nm gain bandwidth, are used to enable a record single mode fibre transmission capacity of 120.0 Tbit/s using 312$\times{}$35 GBd DP-256QAM over 9$\times{}$70 km spans.

(060.0060) fibre optics and optical communications; (060.2360) fibre optics links and subsystems

§ INTRODUCTION

xlabel=Transmitted Distance [km],
ylabel=SMF Throughput [Tbit/s]
[mark=square*,Dark2-A] coordinates (100,115.9);
[anchor=north west] at (axis cs:100,115.9) [3];

[mark=*,Dark2-B] coordinates (165,101.7);
[anchor=south east] at (axis cs:165,101.7) [5];

[mark=pentagon*,mark size=3pt,Dark2-C] coordinates (240,102.3);
[anchor=south west] at (axis cs:240,102.3) [4];

[mark=text,text mark=$\bigstar$,Dark2-D] coordinates (630,120);
[anchor=north east] at (axis cs:630,120) This work $^\nearrow$;

Capacity records in a single-mode fibre

Fundamental limits of an optical communication system are imposed by a combination of noise from the transceiver subsystem, optical amplifier noise and optical fibre nonlinearity. Advances in digital signal processing (DSP), modulation formats [1], high speed electronics, transmission fibres, and broadband optical amplification, have resulted in recent experimental single-mode fibre (SMF) capacity increases. Recent years have seen several landmark transmission results using single mode fibre, with demonstrations achieving throughputs exceeding  70 Tbit/s over transatlantic distances [2] and an overall throughput record of 115 Tbit/s over 100 km [3].

With the exception of [3] [4], previous record demonstrations have relied on C+L-band EDFA transmission, in which a spectral gap between the two bands was required in order to use EDFAs in each transmission window. In [3], a continuous-band 100 nm amplifier was developed based on a semiconductor optical amplifier (SOA), enabling an SMF capacity record of 115.9 Tbit/s over 100 km.
Although the bandwidth is impressive, SOAs have a relatively high noise figure compared with distributed Raman amplifiers, and so system performance decreases rapidly with distance.

In this work, we demonstrate record throughput of 120.0 Tbit/s transmitted over single mode fibre using a 91 nm continuous gain bandwidth hybrid distributed Raman-EDFA (HRE), designed as a prototype for wideband amplification. As highlighted in Fig. <ref>, the data throughput and distance are both increased over the previous capacity record, with an extra 4 Tbit/s capacity over more than 6 times longer transmission distance.

Schematic of the transmission experiment. Inset a) signal spectrum after the GEA (the modulated channels are centered at 1551 nm. Inset b) 256QAM constellation measured in back-to-back configuration.

§ EXPERIMENTAL CONFIGURATION

The experimental configuration used for this work is shown in Fig. <ref>.  Four carriers, spaced at 35.5 GHz, were connected to two independent dual-polarization IQ optical modulators, each driven by four 92 GS/s digital-to-analogue converters (DACs) to generate four odd/even channels. A digital root-raised cosine (RRC) filter with 0.01 roll-off was used to spectrally shape the signals and pre-emphasis was applied to overcome the electrical response of the transmitter components. The channels were generated at carrier frequencies which were tuned across the range 185.5405 to 196.5810 THz,
allowing the measurement of 312$\times{}$35 GBd dual-polarization 256-ary quadrature amplitude modulation (DP-256QAM) channels, covering the spectrum from 1524.9 to 1615.9 nm.

The modulated channels were amplified using a pair of 97 nm bandwidth discrete Raman amplifiers with 12.5 dB gain and combined with wideband amplified spontaneous emission (ASE) noise, which emulated co-propagating channels over the entire
transmitted bandwidth. The ASE noise, with continuity across the entire 97 nm bandwidth, was generated by a pair of discrete Raman amplifiers to achieve a total output power of 22 dBm within the required 91nm bandwidth, following which, a band stop filter (BSF) was used to create a notch in the ASE within which the modulated channels were positioned.
The validity of using ASE noise to emulate agressor channels was verified in [6], showing that this technique provides a conservative measure of system performance.
This technique has been used in [3] and [7] to estimate transmission system capacity.

A gain equalizing amplifier (GEA) with a continuous 91.04 nm bandwidth was used to amplify and spectrally shape (SS) the combined ASE and modulated channels. The combined SS-ASE and modulated channels occupied a total useable bandwidth of 11.0758 THz (91.04 nm) with a total output power of 22 dBm. The optical spectrum after the GEA is shown as inset (a) in Fig. <ref>. A power tilt of -2 dB across the bandwidth was applied to optimize channel performance, taking into account the wavelength-dependent noise figure of the transmission line amplifiers. The back-to-back 256QAM constellation is illustrated in the inset (b) of Fig. <ref>.

The transmission link comprised a straight-line of 9 spans, with 70 km of single-mode fibre and a hybrid distributed hybrid Raman-EDFA amplifier (HRE) in each span. Each HRE provided continuous gain from 1525 nm to 1616 nm and used two counter-propagating pumps at 1427 and 1495 nm with output powers of 300 mW and 310 mW into the transmission fibre, delivering a total power of 19.5 dBm to the EDFA stage. The single EDFA stage boosted the signal to a total output power of 22 dBm, and included a 91 nm gain flattening filter (GFF) to equalize the gain across the entire HRE bandwidth. Each 70 km fibre span comprised two fibre types; the first part of the hybrid span used 40 km of Sumitomo Z+150 fibre with an average attenuation of 0.148 dB/km and an effective core area of $149$ µm$^{2}$. The second part of the span used 30 km of Sumitomo Z fibre with an average attenuation of 0.16 dB/km and an effective core area $81$ µm$^{2}$. The second segment of the hybrid span used a lower effective area to increase the Raman gain in the second part of the span for a given pump power.

At the receiver, a bandpass filter (BPF) with a 40 GHz bandwidth was used to filter the channel under test. The BSF and the transmitter BPF were moved together in order to create a notch for the channel under test (channel 2 in the 4-channel group) and exclude the out-of-band noise from the receiver.  The coherent detection was carried out using a phase- and polarization-diverse coherent receiver incorporating 70 GHz bandwidth photodetectors, and the signal was digitized using a real-time oscilloscope with 63 GHz bandwidth, sampling at 160 GSa/s. Digital signal processing was performed as described in [6], which included matched filtering, single step chromatic dispersion compensation, a 21-tap blind adaptive equalizer, frequency offset compensation and decision directed carrier phase estimation. To quantify the system performance all 312 measured
channels were individually decoded with soft decision forward error correction (SDFEC) and were accepted if the post-SDFEC was below the assumed outer hard decision-FEC (HDFEC) threshold. For the FEC,
12 rate adapted low density parity check (LDCP) codes from the DVB-S2X standard were implemented.
The assumed outer HD-FEC with a bit error rate (BER) threshold of $3\times 10^{-4}$ allows for a BCH code with $0.5\%$ overhead. System performance was also quantified in terms of information rates.  Signal-to-noise ratio (SNR) was estimated for every channel and the generalised mutual information (GMI) was calculated on a per channel basis.

§ RESULTS

Fig.<ref> plots the information rate of all 312 x 35 GBd DP-256QAM channels, including the 0.5 GHz channel spacing. The green

diamonds show the Shannon limit for the received SNR given by $\log_2(1+\mathrm{SNR})$ summed over both polarisations. The total received SNR was evaluated as the ratio between the variance of the transmitted symbols $E[|X|^2]$ and the variance of the noise $\sigma^2$, where $\sigma^2=E[|X-Y|^2]$ and $Y$ represents the received symbols. The blue triangles depict the pre-FEC GMI calculated using received log-likelihood ratios. The post-FEC rate after inner and outer FEC is plotted by the red circles. A mean penalty of 0.94 bit/symbol between the GMI (11.62 bit/symbol) and the Shannon limit (12.56 bit/symbol) is due to the use of non-optimal finite constellation and bit labeling. For the net rate after FEC, the occupied spectrum yields a net bit rate of 10.99 bits/symbol providing a record single mode fibre capacity of 120.0 Tbit/s.

Fig. <ref> shows the pre- and post-SDFEC BER for all 312 channels. The channels were decoded using 12 rate adapted LDPC codes implemented from the DVB-S2x. An outer BCH HD-FEC code with a bit error rate (BER) threshold of $3\times 10^{-4}$  with $0.5\%$ overhead was assumed. All 312 channels were measured, confirming the total net throughput of 120.0 Tbit/s. Fig. <ref> shows the rates used for the SDFEC after which all channels were below the BER threshold for the outer HDFEC. The bulk of the codewords have a low spread, 293 out of 312 are between 40-50% FEC overhead. Reducing the number of 7 code rates applied was found to reduce the total net data throughput to 119.0 Tbit/s.

width = 1.05,
height = 6cm,
xlabel=Channels [nm],
ylabel=Throughput [bit/symbol],

legend reversed=true,
legend style=
legend cell align=left,

[Set1-A,only marks,mark=o,mark size=2pt] table[col sep=comma,x expr=299792458/Fc/1000,y expr=16/1.005*coderate]postfec.csv;
Net rate after FEC

[Set1-B,only marks,mark=triangle,mark size=2pt] table[col sep=comma,x expr=299792458/Fc/1000,y expr=GMI]postfec.csv;

[Set1-C,only marks,mark=diamond,mark size=2pt] table[col sep=comma,x expr=299792458/Fc/1000,y expr=2*log2(1+10^(SNR/10)]postfec.csv;
$\log_2(1+\mathrm{SNR}) $

Throughput per channel over 2 polarisations after 630 km.

width = 1,
height = 6cm,
xlabel=Carrier wavelength [nm],
legend reversed=true,
legend style=

ymin = 1e-8,
legend cell align=left,
legend columns=2,
log origin=infty,

only marks,

text mark=$\square$,
] table[col sep=comma,x expr=299792458/Fc/1000,y expr=postFECBER]postfec.csv;

only marks,
text mark=$\mathrlap{_\downarrow}\square$,
forget plot,
y filter/.code=ifthenelse(rawy==0,-16,nan),
] table[col sep=comma,x expr=299792458/Fc/1000,y=postFECBER]postfec.csv;

[Set1-B,only marks,mark=+] table[col sep=comma,x expr=299792458/Fc/1000,y expr=preFECBER]postfec.csv;

[black,dashed,thick,forget plot,domain=1515:1630] 3e-4;

[anchor=south west,font=] at (axis cs:1520,3e-4) $0.5\%$ HD-FEC threshold;

[<-,thick] (axis cs:1570,0.7e-7) – (axis cs:1580,0.25e-7) node[anchor=west,align=left,style=font=] no errors

Pre-FEC and Post-SDFEC bit error rate for all 312 channels. Arrows indicate channels where no errors were observed after post-SDFEC.

xlabel=LDPC OH [%],
ylabel=Number of Channels,
log ticks with fixed point,
[ybar,Set1-B,fill,bar width=0.2pt] table[col sep=comma,x expr=(1/x-1)*100]

LDPC overhead usages.

§ CONCLUSION

A record capacity of 120 Tbit/s in a single mode optical fibre was realised through the development of a wide-bandwidth hybrid EDFA-Raman amplification scheme. In contrast to the most recent record capacity demonstration [3], this technique offers distributed gain, and so inherently maintains a high SNR over long transmission distances and the relatively long 70 km spans used in this work. In contrast to EDFA-based C+L band transmission systems, the continuous 91 nm gain bandwidth offered by hybrid amplification means that no spectrum is wasted.

Finally, we note that this demonstration has been made using mature signal technologies (standard square QAM modulation and binary LDPC codes) and is, thus, eminently and readily realisable for commercial deployment.

Support from UK EPSRC UNLOC programme and the Royal Academy of Engineering under the Research Fellowships programme is gratefully acknowleged. The authors thank Oclaro for the modulators and Nick Hepden for the heroic effort in constructing the hybrid amplifiers used in this work.

[1] S. Chandrasekhar  et al., “High-spectral-efficiency transmission of PDM-256QAM with Parallel Probabilistic Shaping at Record Rate-Reach Trade-offs,” Proc. ECOC, PDP Th3C.1 (2016).

[2] T-X. Cai  et al., “70.46 Tb/s over 7,600 km in in C+ L Band Using Coded Modulation with Hybrid Constellation Shaping and Nonlinearity Compensation,” Proc. OFC, PDP Th5B.2 (2017).

[3] J. Renaudier et al., “First 100-nm Continuous-Band WDM Transmission System with 115Tb/s Transport over 100km Using Novel Ultra-Wideband Semiconductor Optical Amplifiers,” Proc. ECOC, PDP (2017).

[4] A. Sano et al., “102.3-Tb/s C-band extended L-band all Raman transmission over 240 km using PDM-64QAM single carrier FDM with digital pilot tone,” Proc. OFC, PDPSC.3 (2012).

[5] D. Qian et al., “101.7-Tb/s (370x294-Gb/s) PDM-128QAM-OFDM Transmission over 3x55-km SSMF using Pilot-based Phase Noise Mitigation,” Proc. OFC, PDPB.5 (2011).

[6] D. J. Elson et al., “Investigation of Bandwidth Loading in Optical Fibre Transmission using Amplified Spontaneous Emission Noise,” Optics Express 25 (16), (2017).

[7] J. Cho  et al., “Trans-atlantic field trial using probabilistically shaped 64-QAM at high spectral efficiencies and single-carrier real-time 250-gb/s 16-QAM,” Proc. OFC, PDP Th5B.3 (2017).