ABSTRACT
Title of dissertation: THE SEARCH FOR
SUPERMASSIVE BLACK HOLE BINARIES
IN THE TIME DOMAIN
Tingting Liu
Doctor of Philosophy, 2018
Dissertation directed by: Professor Suvi Gezari
Department of Astronomy
Supermassive black hole binaries (SMBHBs) are expected due to galaxy merg-
ers and the ubiquity of central supermassive black holes (SMBHs) in galaxies, but
direct evidence for close-separation SMBHBs has been elusive. This thesis presents
my search for SMBHBs in the optical time domain, i.e. by searching for their opti-
cal variability signatures. It is a novel approach that can potentially yield SMBHBs
in close, sub-pc orbits, a population of SMBH pairs or binaries that can not be
directly imaged or resolved by current telescopes or techniques. Further, searches
in the optical time domain are sensitive to SMBHBs in the low-frequency gravita-
tional wave-emitting regime, opening up the possibility of studying them in the era
of multi-messenger astronomy.
I developed a custom pipeline to systematically search in the Pan-STARRS1
Medium Deep Survey (PS1 MDS) for periodically varying quasars, which have been
predicted as the manifestations of SMBHBs at close separations. I constructed a
spectroscopically-complete sample of SMBHB candidates using observations with
the Gemini Telescope or the Discovery Channel Telescope and measured their black
hole masses and redshifts. I also followed up the candidates with a dedicated mon-
itoring program on the Las Cumbres Observatory telescopes, in order to put their
periodicity to the test and identify false positives that are due to the stochastic
variability of regular quasars that do not host SMBHBs. I set up the expectations
for a true periodic signal by modeling normal quasar variability and showed that
evidence for a true signal should strengthen over a longer temporal baseline. I then
used the expectations as a guide and applied a range of statistical criteria to select
more robust candidates from PS1 MDS. From this down-selected sample, I was able
to determine an upper limit on the SMBHB rate. I also discussed the search for
SMBHBs in the era of the Large Synoptic Survey Telescope and SMBHB candidates
as possible gravitational wave sources for the pulsar timing arrays.
THE SEARCH FOR SUPERMASSIVE BLACK HOLE BINARIES
IN THE TIME DOMAIN
by
Tingting Liu
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2018
Advisory Committee:
Professor Suvi Gezari, Chair/Advisor
Professor M. Coleman Miller
Professor Richard Mushotzky
Professor Scott C. Noble
Professor Peter Shawhan, Dean’s Representative
©c Copyright by
Tingting Liu
2018

Preface
Chapter 2 is published as Liu, T. et al. 2015, The Astrophysical Journal
Letters, 803, L16;
Chapter 3 is published as Liu, T. et al. 2016, The Astrophysical Journal, 833,
6;
Chapter 4 is published as Liu, T., Gezari, S., & Miller, M.C. 2018, The Astro-
physical Journal Letters, 859, L12; and
Chapter 5 is to be submitted as “Supermassive Black Hole Binary Candidates
from the Pan-STARRS1 Medium Deep Survey” to The Astrophysical Journal.
ii
Acknowledgments
I owe my gratitude first and foremost to Dr. Suvi Gezari, who has supervised
my thesis research.
I have also had the honor and pleasure to work with Dr. Cole Miller on the
work presented in Chapter 4.
I am also grateful to the other members of my thesis committee, Dr. Richard
Mushotzky, Dr. Scott Noble, and Dr. Peter Shawhan; I could not have assembled
a better committee.
The work presented in this thesis made extensive use of the following observ-
ing facilities: the Panoramic Survey Telescope and Rapid Response System (Pan-
STARRS), operated by the Institute for Astronomy, the University of Hawaii; the
Gemini Telescope and Las Cumbres Observatory telescopes, through the National
Optical Astronomy Observatory; and the Discovery Channel Telescope at the Lowell
Observatory. This thesis research was partially supported by the National Science
Foundation.
iii
Table of Contents
Preface ii
Acknowledgements iii
List of Tables vii
List of Figures ix
List of Abbreviations xi
1 Introduction 1
1.1 Supermassive Black Hole Binaries . . . . . . . . . . . . . . . . . . . . 1
1.1.1 SMBHBs as the Products of Galaxy Mergers . . . . . . . . . . 1
1.1.2 Dual AGNs and the Spectroscopic Search for Close-Separation
SMBHBs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.3 The SMBHB Candidate OJ 287 . . . . . . . . . . . . . . . . . 8
1.1.4 SMBHBs as the Sirens of Low-frequency Gravitational Waves 10
1.2 Theoretically Predicted Variability Signatures of an SMBHB . . . . . 13
1.2.1 Modulated Mass Accretion as the Cause of Periodicity . . . . 13
1.2.2 Doppler Boosting as the Cause of Periodicity . . . . . . . . . . 17
1.3 The Search for SMBHBs . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.1 The Detectability of SMBHBs in a Time Domain Survey . . . 18
1.3.2 Systematically Searching for Periodically Varying Quasars . . 20
1.4 Quasar Variability and the Challenges of Searching for Periodic Quasars 24
1.4.1 Optical Variability of Quasars . . . . . . . . . . . . . . . . . . 24
1.4.2 Red Noise Contamination of the SMBHB Candidate Sample . 28
2 A Periodically Varying Luminous Quasar at z = 2 from the Pan-STARRS1
Medium Deep Survey: A Candidate Supermassive Black Hole Binary in the
Gravitational Wave-Driven Regime 33
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2 Theoretical Predictions . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3 Pan-STARRS1 Medium Deep Survey . . . . . . . . . . . . . . . . . . 37
iv
2.4 Ensemble Photometry . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.5 Selecting Periodic Quasar Candidates . . . . . . . . . . . . . . . . . . 39
2.6 Periodic Quasar Candidate PSO J334.2028+01.4075 . . . . . . . . . . 40
2.7 Physical Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.8 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 43
3 A Systematic Search for Periodically Varying Quasars in Pan-STARRS1: An
Extended Baseline Test in Medium Deep Survey Field MD09 49
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 The Pan-STARRS1 Medium Deep Survey . . . . . . . . . . . . . . . 53
3.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3.1 Sample Selection . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3.2 Variability Selection . . . . . . . . . . . . . . . . . . . . . . . 58
3.3.3 Selection Bias . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.3.4 Selecting Periodic Quasar Candidates . . . . . . . . . . . . . . 63
3.4 Extended Baseline Photometry . . . . . . . . . . . . . . . . . . . . . 65
3.4.1 Follow-up Imaging . . . . . . . . . . . . . . . . . . . . . . . . 67
3.4.2 Pre-PS1 Archival Photometry . . . . . . . . . . . . . . . . . . 69
3.4.3 Testing the Persistence of Periodicity with Extended Baseline
Photometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5 Black Hole Mass Estimates and Inferred Binary Parameters . . . . . 72
3.6 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 75
4 Did ASAS-SN Kill the Supermassive Black Hole Binary Candidate PG 1302−102? 93
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.2 Extended Light Curve from ASAS-SN . . . . . . . . . . . . . . . . . . 95
4.3 Expectations for a True Periodic Signal . . . . . . . . . . . . . . . . . 97
4.4 Extended Baseline Analyses of PG 1302 . . . . . . . . . . . . . . . . 100
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5 Supermassive Black Hole Binary Candidates from the Pan-STARRS1 Medium
Deep Survey 110
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.2 The Pan-STARRS1 Medium Deep Survey . . . . . . . . . . . . . . . 114
5.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.3.1 Quasar and Variability Selection . . . . . . . . . . . . . . . . . 115
5.3.2 Searching for Periodicity . . . . . . . . . . . . . . . . . . . . . 118
5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.4.1 Periodic Quasar Candidates from MDS . . . . . . . . . . . . . 118
5.4.2 Spectroscopy and Black Hole Mass . . . . . . . . . . . . . . . 120
5.5 Extended Baseline Photometry . . . . . . . . . . . . . . . . . . . . . 125
5.6 Reevaluating Periodic Candidates Using a Maximum Likelihood Method129
5.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
5.7.1 The Detection Rate of SMBHBs . . . . . . . . . . . . . . . . . 149
v
5.7.2 Implications for Gravitational Wave Astrophysics in the PTA
Regime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
5.7.3 Periodic Quasar Detections in the LSST Era . . . . . . . . . . 153
5.7.4 Probing the SED Properties of SMBHBs . . . . . . . . . . . . 157
5.8 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 158
6 Conclusions and Future Work 182
6.1 What Would be a Convincing Detection of Periodicity? . . . . . . . . 183
6.2 Can We Find Robust Periodic Candidates with Extended-baseline
Imaging? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
6.3 Can We Search for Evidence for Periodicity in the X-rays? . . . . . . 187
6.4 Can We Search for Complementary Evidence for an SMBHB? . . . . 189
Bibliography 198
vi
List of Tables
1.1 Summary of Systematic Searches for Periodically Varying Quasars . . 24
3.1 PS1 Quality Flags Used in the Query . . . . . . . . . . . . . . . . . . 58
3.2 QLF Model Values Used . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.3 MD09 by the Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.4 PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5 Extended Baseline Photometry of Periodic Quasar Candidates . . . . 87
3.6 SDSS to PS1 Filter Offsets from Synthetic Quasar Magnitudes . . . . 88
3.7 Detected Periods and Significance Factors of Periodogram-selected
Candidates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.8 Spectroscopic Information and Inferred Binary Parameters of Peri-
odic Quasar Candidates . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.1 Maximum Likelihoods for the Simulated DRW+periodic Light Curve 107
4.2 Damped Random Walk Parameter Ranges Sampled by MCMC . . . . 108
4.3 Broken Power-law Parameter Ranges Sampled by MCMC . . . . . . . 109
5.1 Medium Deep Field Centers . . . . . . . . . . . . . . . . . . . . . . . 116
5.2 Medium Deep Fields by the Numbers (MD01-MD02) . . . . . . . . . 131
5.3 Medium Deep Fields by the Numbers (MD03-MD04) . . . . . . . . . 132
5.4 Medium Deep Fields by the Numbers (MD05-MD06) . . . . . . . . . 133
5.5 Medium Deep Fields by the Numbers (MD07-MD08) . . . . . . . . . 134
5.6 Medium Deep Fields by the Numbers (MD09-MD10) . . . . . . . . . 135
5.7 Period, Significance Factors, and Number of Cycles of Periodic Can-
didates (MD01-MD04) . . . . . . . . . . . . . . . . . . . . . . . . . . 136
5.8 Period, Significance Factors, and Number of Cycles of Periodic Can-
didates (MD05-MD10) . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.9 PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates (MD01-05) . . . . . . . . . . . . . . . . . . . . . . . . . . 138
5.10 PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates (MD06-10) . . . . . . . . . . . . . . . . . . . . . . . . . . 139
5.11 Spectroscopic Follow-ups . . . . . . . . . . . . . . . . . . . . . . . . . 140
vii
5.12 Spectroscopic Measurements and Inferred Binary Parameters (MD01-
04) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.13 Spectroscopic Measurements and Inferred Binary Parameters (MD05-
10) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
5.14 Extended Baseline Monitoring of Candidates (MD01-MD04) . . . . . 145
5.15 Extended Baseline Monitoring of Candidates (MD05-MD10) . . . . . 146
5.16 Re-analyses Using the Maximum Likelihood Method under the DRW
Model (MD01-MD04) . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.17 Re-analyses Using the Maximum Likelihood Method under the DRW
Model (MD05-MD10) . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
5.18 Re-analyses Using the Maximum Likelihood Method under the BPL
Model (MD01-MD04) . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
5.19 Re-analyses Using the Maximum Likelihood Method under the BPL
Model (MD05-MD10) . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
5.20 Comparing PS1 MDS and LSST Capabilities . . . . . . . . . . . . . . 163
5.21 Parameters Used in Figure 5.15 . . . . . . . . . . . . . . . . . . . . . 165
5.22 The SED of PSO J185 . . . . . . . . . . . . . . . . . . . . . . . . . . 167
viii
List of Figures
1.1 The VLBA image of CSO 0402+379 . . . . . . . . . . . . . . . . . . 6
1.2 The V-band Light Curve and the Binary Model of OJ 287 . . . . . . 10
1.3 Simulations of a Binary-Disk System . . . . . . . . . . . . . . . . . . 15
1.4 The Ṁ Variation as a Function of Time . . . . . . . . . . . . . . . . 16
1.5 The Detectability of SMBHBs as Periodic Quasars in a Time Domain
Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 The Light Curve of the SMBHB Candidate PG 1302−102. . . . . . . 22
1.7 Light Curves Simulated From the Damped Random Walk Model . . . 27
2.1 Color–Color Diagrams Used to Select Point Sources in MD09 . . . . . 45
2.2 Sinusoidal Fit to the PS1 Light Curves of PSO J334.2028+01.4075 . . 46
2.3 LS Periodogram of PSO J334.2028+01.4075 . . . . . . . . . . . . . . 47
2.4 The Archival CRTS Light Curve of PSO J334.2028+01.4075 . . . . . 48
3.1 Comparing the Light Curves of the Periodic Quasar Candidate PSO
J334.2028+1.4075 From PV1 and PV2. . . . . . . . . . . . . . . . . . 78
3.2 Color Diagrams of Cross-matched Stars and Quasars . . . . . . . . . 79
3.3 The Observed Scatter of the Light Curve as a Function of Magnitude 80
3.4 Magnitude Error vs. Magnitude in the gP1 Filter . . . . . . . . . . . . 81
3.5 The Number Fraction of Our MD09 Quasars that Vary at the > σint
Level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.6 The gP1 Band Apparent Magnitude Distribution of Our MD09 Quasar
Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.7 Selection Bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.8 The LS Periodograms of PS1 MD09 Candidates . . . . . . . . . . . . 85
3.9 The Light Curves of PS1 MD09 Candidates . . . . . . . . . . . . . . 86
3.10 Spectra of PS1 MD9 Candidates . . . . . . . . . . . . . . . . . . . . . 91
3.11 Comparing with the Periodic Candidates from Graham et al. (2015b)
and Charisi et al. (2016) . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.1 Combined Light Curve of PG 1302-102 . . . . . . . . . . . . . . . . . 105
4.2 A Simulated DRW+periodic Light Curve . . . . . . . . . . . . . . . . 106
ix
5.1 The Number of Detections for a Random Sample of ∼ 27, 000 Objects
in MDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.2 The Color Diagram of Quasars in PS1 MDS . . . . . . . . . . . . . . 117
5.3 The Mean g Band Magnitude Distribution of Periodic Candidates . . 119
5.4 The Distribution of the Observed Periods . . . . . . . . . . . . . . . . 120
5.5 The Variability Amplitudes of Candidates . . . . . . . . . . . . . . . 121
5.6 The Distribution of the Rest-frame Periods . . . . . . . . . . . . . . . 124
5.7 The Black Hole Mass Distribution of Candidates . . . . . . . . . . . . 125
5.8 The Redshift Distribution of Candidates . . . . . . . . . . . . . . . . 126
5.9 The Parameter Space Probed by PS1 MDS and LSST . . . . . . . . . 143
5.10 The Number of Cycles . . . . . . . . . . . . . . . . . . . . . . . . . . 144
5.11 The Magnitude Distribution of Candidates . . . . . . . . . . . . . . . 150
5.12 The Cumulative Number of SMBHB Candidates . . . . . . . . . . . . 162
5.13 Gravitational Wave Strain Amplitudes of SMBHB Candidates . . . . 163
5.14 Detection Capabilities of the LSST . . . . . . . . . . . . . . . . . . . 164
5.15 GW Frequencies and Strain Amplitudes Explored by the LSST . . . . 166
5.16 The SED of PSO J185 . . . . . . . . . . . . . . . . . . . . . . . . . . 168
5.17 The PS1-only and Extended Light Curve of PSO J185 . . . . . . . . 169
5.18 The SDSS Spectrum of PSO J185 . . . . . . . . . . . . . . . . . . . . 170
5.19 The LS Periodogram of PSO J185 . . . . . . . . . . . . . . . . . . . . 170
5.20 The PS1-only and Extended Light Curves of Periodic Candidates
from PS1 MDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
5.21 Figure 5.20 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 172
5.22 Figure 5.21 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 173
5.23 Figure 5.22 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 174
5.24 Figure 5.23 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 175
5.25 Figure 5.24 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 176
5.26 Figure 5.25 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 177
5.27 SDSS Spectra of PS1 MDS Candidates . . . . . . . . . . . . . . . . . 178
5.28 Figure 5.27 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 179
5.29 Figure 5.28 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 180
5.30 Figure 5.29 Continued . . . . . . . . . . . . . . . . . . . . . . . . . . 181
x
List of Abbreviations
ACS the Advanced Camera for Surveys
AGN active galactic nucleus/nuclei
AIC Akaike information criterion
ASAS-SN All-sky Automated Survey for Supernovae
BAT Burst Alert Telescope
BPL broken power-law
CFHT Canada-France-Hawaii Telescope
CRTS Catalina Real-time Transient Survey
CSO compact symmetric object
CSS Catalina Sky Survey
CW continuous wave
DCT Discovery Channel Telescope
DRW damped random walk
EM electromagnetic
EPTA European Pulsar Timing Array
FOV field of view
FWHM full-width at half-maximum
GALEX Galaxy Evolution Explorer
GR general relativity
GW gravitational wave
GWB gravitational wave background
HST Hubble Space Telescope
IfA Institute for Astronomy
IPP Image Processing Pipeline
IPTA International Pulsar Timing Array
iPTF Intermediate Palomar Transient Factory
LCO Las Cumbres Observatory
LEDE luminosity evolution and density evolution
LINEAR Lincoln Near-Earth Asteroid Research
LISA Laser Interferometer Space Antenna
xi
LMI Large Monolithic Imager
LS Lomb-Scargle
LSST Large Synoptic Survey Telescope
MCMC Markov Chain Monte Carlo
MDS Medium Deep Survey
MHD magnetohydrodynamic
MJD Modified Julian Date
MSP millisecond pulsar
NANOGrav North American Nanohertz Observatory for Gravitational Waves
NuSTAR Nuclear Spectroscopic Telescope Array
Pan-STARRS Panoramic Survey Telescope and Rapid Response System
PLE pure luminosity evolution
PPTA Parkes Pulsar Timing Array
PSD power spectral density
PSF point spread function
PSI Pan-STARRS Science Interface
PTA pulsar timing array
PTF Palomar Transient Factory
PV1 Processing Version 1
QLF quasar luminosity function
SDSS Sloan Digital Sky Survey
SED spectral energy distribution
SF Structure Function
SKA Square Kilometre Array
SMBH supermassive black hole
SMBHB supermassive black hole binary
SVM support vector machine
TOA time of arrival
ULIRG ultraluminous infrared galaxy
UV ultraviolet
VLBA Very Long Baseline Array
VLBI very long baseline interferometry
xii
XRT X-ray Telescope
ZTF Zwicky Transient Facility
xiii
Chapter 1: Introduction
1.1 Supermassive Black Hole Binaries
1.1.1 SMBHBs as the Products of Galaxy Mergers
Supermassive black holes (SMBHs) with masses between 106− 109M appear
to be ubiquitous in the center of massive galaxies (e.g. reviews by Ferrarese & Ford
2005; Kormendy & Ho 2013; Kormendy & Richstone 1995). Their masses are tightly
correlated with host galaxy properties such as the stellar velocity dispersion of the
bulge, bulge luminosity, and bulge stellar mass, suggesting that the evolution of the
SMBHs is closely linked to that of the galaxies (e.g. Gültekin et al. 2009; McConnell
& Ma 2013).
In the paradigm of structure formation in the standard ΛCDM cosmology
(e.g. White & Rees 1978), galaxies merge with each other and grow into larger
units. During the merger of galaxies, the SMBHs are expected to settle to the
center of the newly-merged system, eventually forming a supermassive black hole
binary (SMBHB) (Begelman et al., 1980). This process is first driven by dynamical
friction, which efficiently drags the SMBHs to the center of the merger remnant on
a timescale of ∼ 108 years. After dynamical friction becomes inefficient at ∼ parsec
1
(pc) scales, the SMBHs can further lose energy and angular momentum through
three-body interaction with stars in the so-called “loss cone”: the star encountering
the SMBHs is “slingshot” out of the system, carrying an amount of energy with
it. The decrease in the semimajor axis of the SMBHs is a function of the velocity
dispersion and the density of the isotropic background of stars and a hardening rate,
and therefore this process is efficient for a fixed stellar background.
However, the loss cone is gradually depleted of stars, but the binary separa-
tion is still ∼ 100 times larger than the separation where gravitational wave (GW)
radiation can take over so that the binary coalesces within a Hubble time (the age
of the Universe). Without replenishment of the loss cone, or other mechanisms, the
binary will forever stall at a ∼ parsec separation. This gives rise to the “final parsec
problem” (Milosavljević & Merritt, 2003) and is most severe in spherical galaxies.
However, recent N-body simulations show that in non-spherical galaxies (especially
triaxial galaxies), there is a larger phase space of stellar orbits which can keep the
loss cone filled, so that the timescale of the stellar hardening stage can be comfort-
ably less than ∼ Gigayear (Gyr) (e.g. Vasiliev et al. 2015). The final parsec problem
as the barrier to SMBHB formation therefore can be alleviated.
After forming a close orbit (a < 0.1 pc), the SMBHs are now embedded in
a circumbinary gas disk. The evolution of the binary is in two phases: at larger
separations, the binary exchanges angular momentum with the disk through the
accretion of mass and gravitational torques (e.g. Tang et al. 2017); at close enough
separations (centi to milli-parsec scales), GW emission dominates the orbital decay,
which rapidly drives the binary to coalesce. It is at this stage that the electromag-
2
netic variability signatures have been theoretically predicted (they are summarized
in Section 1.2), signatures for which this thesis is dedicated to searching.
1.1.2 Dual AGNs and the Spectroscopic Search for Close-Separation
SMBHBs
Observational evidence for the galaxy merger scenario described above is pairs
of SMBHs not yet gravitationally bound to each other. These SMBHs at ∼ kpc
separations are called dual SMBHs and are progenitors of SMBHBs. If gas is fun-
neled to the center of the galaxy merger remnant in a gas-rich merger and actively
accreted by the black hole(s), the dual SMBHs would appear as offset or dual active
galactic nuclei (AGNs). In the local universe, the two active nuclei separated by
∼ kpc can be imaged directly. A prototype of such systems is the ultraluminous
infrared galaxy (ULIRG) NGC 6240. It shows two nuclei separated by 1′′.8 (640
pc) with an unusual, disturbed morphology (e.g. Fried & Schulz 1983), and high
spatial resolution X-ray imaging with Chandra indicates that both nuclei are AGNs
(rather than star-forming or starburst regions) (Komossa et al., 2003).
However, while more kpc-separation dual AGNs have been uncovered serendip-
itously (e.g. Comerford et al. 2009a) or via hard X-ray imaging (e.g. Koss et al.
2011), identifying dual AGNs in large numbers remained difficult, and a new and
systematic approach was proposed to identify more candidates: if both SMBHs are
active as AGNs, the galaxy spectrum should display two sets of narrow emission
lines (e.g. [OIII]) that are offset from the stellar absorption lines of the host galaxy
3
by a few hundred km s−1 (Comerford et al., 2009b). This technique has been applied
to search for dual AGNs in the DEEP2 Galaxy Redshift Survey (Comerford et al.,
2009b) and the Sloan Digital Sky Survey (SDSS) (Liu et al., 2010; Smith et al., 2010;
Wang et al., 2009), and it has provided a larger sample of candidates for follow-up
studies.
For example, integral-field spectroscopy of one of the Type-1 AGNs with
double-peaked [OIII] profiles from the Smith et al. (2010) sample, SDSS J150243.1+111557
(hereafter SDSS J1502), shows that the lines are separated by 657 km s−1, and high
spatial resolution adaptive optics observations in the near-IR indicate that it has
two nuclei separated by 7.4 kpc (Fu et al., 2012). The object was later resolved
by the Expanded Very Large Array (Fu et al., 2011) into two radio components: a
primary “J1502P” and a secondary “J1502S”; both components are steep-spectrum
sources which indicate synchrotron emission. The radio luminosity of SDSS J1502
is inconsistent with that of star-forming galaxies, indicating that it is instead pow-
ered by AGN activity. Its X-ray luminosity and radio-to-X-ray luminosity ratio
further support the discovery that SDSS J1502 is a dual AGN (Fu et al., 2011).
Interestingly, very long baseline interferometry (VLBI) observations of SDSS J1502
with the European VLBI Network further resolved the secondary J1502S into two
components (J1502SE and J1502SW) separated by 26 milliarcseconds (mas), or a
projected distance of ∼ 138 pc (z = 0.39) (Deane et al., 2014), making it one of the
closest-separation resolved binary or dual AGNs discovered to date and surpassed
only by the radio galaxy CSO 0402+3791.
1More recently, Kharb et al. (2017) reported an SMBHB candidate at the center of NGC 7674,
4
CSO 0402+379 was first identified by Maness et al. (2004) via Very Long
Baseline Array (VLBA; an interferometer that consists of ten radio telescopes) as
a compact symmetric object (CSO), a class of radio sources that are compact (< 1
kpc) in size and have emission on both sides of the central engine. However, several
features of 0402+379 are atypical of CSOs; in particular, the object has a large-scale
jet structure, and the two parsec-scale components (C1 and C2) are flat-spectrum
sources. Both cores also appear to be variable, and possible motion of C2 is also
detected (Maness et al., 2004). Maness et al. (2004) proposed several interpretations
of this puzzling source, one of them being that C1 and C2 are two members of an
SMBHB system. VLBA observations by Rodriguez et al. (2006) confirmed the jets
and the two compact components seen previously in Maness et al. (2004) and that
C2 is located between the jets, while C1 is offset from C2 by 7 pc (z = 0.055).
They also measured the motion and the flux variability over a 15-yr baseline: both
components increased in flux, while the lobes showed no clear variability; relative
position measurements show that C2 appears to be moving towards C1 (however,
the significance of the detection cannot be claimed), and both lobes are moving
away from C2 (Rodriguez et al., 2006).
Bansal et al. (2017) further constrain the orbits of C1 and C2 from proper
motion measurements and infer an orbital period of 3× 104 years and a black hole
mass of 1.5 × 1010M. However, given the separation and the constraint on the
black hole mass, the merger timescale of the binary due to gravitational radiation is
where they detected two cores separated by 0.65 mas, or 0.35 pc (z = 0.029). However, the cores
have not been confirmed as AGNs.
5
Figure 1.1: Left: the VLBA image of CSO 0402+379 at 8 GHz that shows the two
compact components and the lobes. Right: the motion of each component relative
to C1. (Both figures from Rodriguez et al. (2006))
8 M2
tmerger = 5.8× 106yr( a )4(10 M )3 1 ∼ 1012 years (Peters, 1964) — much0.01pc M1 M2(M1+M2)
longer than a Hubble time. Therefore, while it is already a gravitationally-bound
binary at < 10 pc, energy and angular momentum losses due to GW emission at
this separation are negligible.
At closer, sub-parsec separations, SMBHBs cannot be spatially resolved by
current telescopes, and indirect techniques are therefore required to search for them.
If one of the SMBHs is active, the broad line emitted by virialized gas bound to the
SMBH may be shifted from its normal wavelength due to the SMBH’s Keplerian
velocity, and therefore offset broad lines could be signs of SMBHBs (or recoiling
SMBHs). Motivated by this physical picture, Eracleous et al. (2012) have system-
6
atically searched for quasars in SDSS whose Hβ lines are offset from the rest-frame
of the quasar by > 1, 000 km s−1. For an observed (therefore, projected) velocity
offset u ∼ 103 km/s and a black hole mass ∼ 108M, the binary orbital period is
P = 2652M ( sin i sinφ)38 yr, and the binary separation is a = 0.432M (
sin i sinφ 2
(1+q)u 8
) pc,
3 (1+q)u3
where M = M /108M , u = u/1038 tot  3 kms
−1, i is the inclination of the orbit, and
φ is the orbital phase. Among the 15, 900 SDSS quasars at z < 0.9, 88 met their
selection criteria.
Further, the evidence for a binary would be stronger if the line is shifted
during the time span of several observations over a baseline of a few years. In the
study by Eracleous et al. (2012), the intervals between the original SDSS epochs and
the follow-up observations are 5 − 10 years, which allow them to probe a velocity
change of 100−200 km/s. Using a cross-correlation analysis, Eracleous et al. (2012)
compute the change in χ2 by shifting one of the spectra by small steps and are able
to measure significant (at the 99 % level) shifts of the offset Hβ lines in 14 quasars.
Runnoe et al. (2015) followed up the 88 quasars with Hβ broad line offsets
in Eracleous et al. (2012) so that many of them have three or more epochs of
spectroscopy over a ∼ 12-yr baseline. In Runnoe et al. (2017), they measure the
line velocity offset at each epoch relative to the first-epoch SDSS spectrum using the
χ2 cross-correlation method in Eracleous et al. (2012), and construct a radial velocity
curve that shows the velocity offset at different epochs. For the 29 objects with at
least three epochs of spectroscopy and in which they are able to reliably measure
the velocity shifts (i.e. no significant line profile shape changes), they are able to
set lower limits on the black hole masses by exploring the shortest periods of the
7
hypothesized SMBHBs that the velocity curves can tolerate. If the minimum black
hole mass inferred is unphysically large, then the original hypothesis — that the
system is an SMBHB – must have been false. Those minimum periods are decades
to centuries, and the minimum black hole masses inferred from their method range
from 4.7×104−3.8×108M, which are consistent with SMBH masses, and they are
therefore not yet able to rule out the binary hypothesis for any of the candidates in
their sample.
At yet closer separations (. 0.01 pc), the broad emission lines are no longer
distinctly associated with one black hole, line profile becomes more complex, and
velocity shifts are no longer good possible indicators of the orbital motion of the
binary (Shen & Loeb, 2010). However, these SMBHBs have periods P ∼ 10 years,
which opens up the possibility of detecting them via periodic variability in their
light curves. More excitingly, these SMBHBs emit gravitational waves (GWs) at
frequencies between f ∼ 10−7GW − 10−9 Hz and can be probed by a pulsar timing
array (Foster & Backer, 1990) in addition to via their electromagnetic signatures
(Section 1.1.4).
1.1.3 The SMBHB Candidate OJ 287
The strongest SMBHB candidate identified via its optical variability is the
blazar OJ 287. It has been monitored since the late-nineteenth century, and despite
the poorer temporal sampling at earlier times, it was realized that outbursts in its
light curve occur at a quasi-regular interval of Pobs ∼ 12 years. It was proposed
8
in 1988 that these quasi-periodic bursts suggest that OJ 287 is an SMBHB system
with an orbital period of Prest ≈ 9 years in the rest frame (z = 0.306, P = Pobsrest )1+z
(Sillanpaa et al., 1988) . In this early binary model, Sillanpaa et al. (1988) proposed
that as the secondary black hole crosses the accretion disk of the primary, it tidally
perturbs the disk and causes an inflow of matter towards the black hole.
However, there were features in the light curve that were not explained by
the 1988 model. In 1996, it was proposed that the secondary black hole orbits
around the primary on a ∼ 12-year period (observed frame) and crosses its accretion
disk twice per orbit, causing two impact flares (which are separated by one to
two years) per 12 yr-period (Lehto & Valtonen, 1996). The impact times give
a unique solution of the binary orbit, and parameters including the orbital period,
binary eccentricity, and precession of the orbit were determined. In 2008, there were
enough well-sampled and well-defined outbursts such that the orbital shrinkage due
to gravitational radiation could be measured: without energy loss due to GWs, the
September 2007 burst would have occurred 20 days later (Valtonen et al., 2008). It
was also predicted that the next outburst would occur in early January 2016, which
was indeed observed in December 2015 (Valtonen et al., 2016). The profile shape of
the 2015 flare was modeled as an expanding bubble of hot plasma (Pihajoki, 2016),
and coordinated polarimetric and X-ray observations also agreed with the thermal
nature of the burst (Valtonen et al., 2016).
While the 1996 binary model of OJ 287 remains somewhat controversial and
other models including the precession of the radio jet (e.g. Britzen et al. (2018))
have been proposed, OJ 287 is one of the best SMBHB candidates to date and a
9
Figure 1.2: Left: the V-band light curve of the SMBHB candidate OJ 287 from
∼ 1900 to ∼ 2010 (Figure from Valtonen et al. (2008)). Right: the orbit of OJ 287
from 2000− 2023 based on its binary model (Figure from Valtonen et al. (2016)).
potential target for the pulsar timing arrays (Arzoumanian et al., 2014; Babak et al.,
2016; Zhu et al., 2014).
1.1.4 SMBHBs as the Sirens of Low-frequency Gravitational Waves
Pulsars are rotating, highly-magnetized neutron stars characterized by rapid,
regular radio radiation. As the pulsars rotate, radio beams sweep across the direction
of the Earth like lighthouse beacons and are detected as short pulses at regular
intervals of milliseconds to seconds. The fastest-spinning pulsars are hence called
millisecond pulsars (MSPs). They are thought to belong in binary systems and have
been “recycled” through mass transfer from their binary companions, and they are
characterized by particularly stable pulse periods.
The times of arrival (ToAs) of these pulses measured at the observatory are
determined by factors including the pulsar spin and spin-down, radio waves propa-
10
gating in the interstellar medium, and a time standard. Measuring the precise ToAs
of the pulses and fitting them to a timing model is a technique known as pulsar
timing (e.g. Edwards et al. 2006). Deviations between the ToAs predicted by the
timing model and the actually-measured ToAs are called “timing residuals” and are
therefore due to effects that are not included in the timing model. In particular,
Sazhin (1978) and Detweiler (1979) first realized that GWs would induce such “un-
explained” timing residuals. So far, timing residuals of precision δt ∼ 1µs have
been achieved (Arzoumanian et al., 2018; Lentati et al., 2015; Shannon et al., 2015),
and GWs are expected to induce timing residuals δtGW < 100 ns (e.g. Hobbs et al.
2010).
Not only can timing residuals be induced by GWs, but they should also cor-
relate as a function of the angular separations between pulsars. This distinctive
correlation, known as the Hellings-Downs correlation (Hellings & Downs, 1983),
distinguishes GW-induced residuals from those caused by systematic effects, and it
requires a spatial distribution, or an “array”, of MSPs. There are currently three
pulsar timing arrays (PTAs) in operation: the North American Nanohertz Obser-
vatory for Gravitational Waves (NANOGrav; The NANOGrav Collaboration et al.
2015), the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), and the
Parkes Pulsar Timing Array (PPTA; Reardon et al. 2016), and they have been ob-
serving MSPs since ∼ 1990 at the cadence of once every few weeks. The lower end of
the GW frequency band to which they are sensitive is therefore given by the current
baseline of ∼ 10− 15 years, or fGW ∼ 10−9 Hz, and the upper bound is limited by
the observing cadence, hence fGW = a few ×10−7 Hz (e.g. Hobbs et al. 2010).
11
There are several sources of GWs that emit at these frequencies, and the ones
that are most relevant to this thesis are SMBHBs with orbital periods of months
to years. If individual SMBHBs are louder than the stochastic GW background —
the superposition of signals from a population of SMBHBs — they may be resolved
by the PTAs. Sesana et al. (2009) explored the detectability of individual SMBHBs
and investigated the distribution of the binary parameters including GW frequency,
1/5
chirp mass M = (M 3/5ch 1M2) /Mtot (where Mtot is the total mass of the binary),
and redshift. They conclude that the loudest SMBHBs are massive (Mch > 10
8M)
and at higher redshifts z > 0.2, and at least one resolvable source would induce
a timing residual of δtGW ∼ 10 ns, assuming a five-year observation time. This is
achievable with future PTA experiments but is below the upper limits of the cur-
rent ones. Nevertheless, current experiments have put interesting constraints on
SMBHBs (Arzoumanian et al., 2014; Babak et al., 2016; Zhu et al., 2014): the cur-
rent most stringent limit (Babak et al., 2016) has put constraints on SMBHBs with
M > 1010ch M out to a luminosity distance of ∼ 1 Gpc. The SMBHB candidates
OJ 287 (Section 1.1.3) and PG 1302−102 (Section 1.3.2) are both below this limit
and therefore cannot yet be ruled out. The most notable case of using PTA to
constrain an SMBHB candidate is 3C66B, which was proposed as an SMBHB at
z = 0.02 with a total mass M = 5.4 × 1010M, a period of 1.05 years, and a mass
ratio of 0.1 (Sudou et al., 2003). The timing residual expected from these binary
parameters is δtGW ∼ 5µs, and therefore the signal would have been detected in the
timing residual of PSR B1855+09 (δt = 1.5µs) if 3C66B were a true SMBHB (with
those parameters) (Jenet et al., 2004).
12
1.2 Theoretically Predicted Variability Signatures of an SMBHB
1.2.1 Modulated Mass Accretion as the Cause of Periodicity
As previously mentioned in Section 1.1.1, a circumbinary disk is expected
to form around an SMBHB. Such binary-disk systems have been explored in two-
dimensional (2D) or 3D hydrodynamical or magnetohydrodynamical (MHD) simu-
lations for equal-mass binaries or binaries with a range of mass ratios (e.g. Bowen
et al. 2017, 2018; D’Orazio et al. 2013; Farris et al. 2014; MacFadyen & Milosavl-
jević 2008; Noble et al. 2012; Shi & Krolik 2015; Shi et al. 2012; Tang et al. 2018,
2017). These studies find that the tidal torque of the binary clears and maintains
a low-density gap (“cavity”) of a radius ∼ 2a (where a is the binary separation)
between the binary and the circumbinary disk with material flowing into the cavity
in two narrow streams and onto the black holes. The mass accretion rate measured
for an individual black hole or the total accretion rate of the binary is strongly
modulated by the binary orbital motion, resulting in periodic variation (or multiple
periodicities) in the accretion rate (Ṁ). The different simulations and their findings
of the periodic behavior of the binary are briefly summarized as follows2:
In the 2D hydrodynamical simulation of an equal-mass binary with a circumbi-
nary disk presented by MacFadyen & Milosavljević (2008), the mass accretion rate
as measured at the inner boundary of the disk is modulated on the most prominent
frequency f = 2forb/9 (where forb is the binary orbital frequency), accompanied by
2These simulations have adopted circular orbits.
13
a few harmonic frequencies. D’Orazio et al. (2013) and Farris et al. (2014) followed
up on this work by exploring a range of mass ratios q = M2/M1 (the black hole mass
ratio between the secondary and the primary) and investigated the dependence of
mass accretion (and variability) on q. The choice of mass ratios is motivated by the
wide distribution (q ∼ 0.1− 1) often assumed in galaxy merger models (e.g. Volon-
teri et al. 2003). Both studies find that significant modulation of mass accretion
only occurs for less extreme mass ratios q & 0.1; this can be understood by taking
the limit of q → 0: it therefore approaches a single black hole case and does not
result in modulation by the binary orbital motion. Farris et al. (2014) also measured
the mass accretion onto each BH (instead of through the inner boundary) and found
that the accretion rate onto the secondary is larger than that onto the primary for
all mass ratio cases: (Ṁ2/M2)/(Ṁ1/M1) > 1.
The first 3D MHD simulation of a binary-disk system was performed by Shi
et al. (2012); it is therefore a more realistic description of the angular momentum
exchange between the binary and the disk than previous hydro simulations. Similar
to the MacFadyen & Milosavljevic simulation, which also includes an equal-mass
binary, Shi et al. (2012) find two dominant frequencies that modulate Ṁ : one at
twice the binary orbital frequency, and the other at the orbital frequency of the
“lump”, an overdensity of matter at the inner edge of the disk. However, because of
the stronger internal stress due to MHD, it delivers more material into the cavity,
and the accretion rate at the inner boundary is ∼ 40 times larger than the previous
hydro simulation.
Those simulations mentioned above keep the binary at a fixed orbit and adopt
14
Figure 1.3: Left: a 2D hydrodynamical simulation of an SMBHB system with two
different mass ratios. The time-averaged surface density profiles show the accretion
disks around the individual (primary and secondary) black holes (“minidisks”), as
well as the circumbinary disk around both black holes (Figure from Farris et al.
(2014)). Center: a log-scale surface density profile from a MHD simulation of a
binary. It shows the accretion streams flowing towards the binary at four different
times in the simulation. Right: a linear-scale surface density profile from the same
simulation. It emphasizes the time-evolution of the overdensity (“lump”) in the
inner edge of the circumbinary disk. (Both figures from Noble et al. (2012))
Newtonian dynamics. Those assumptions are appropriate when the binary sepa-
ration is farther apart, but will no longer suffice when the binary has evolved to a
much closer separation (e.g. a few tens of gravitational radii) and the black holes are
rapidly inspiraling due to gravitational wave radiation — this would require post-
Newtonian (Noble et al., 2012) or a full general relativistic (GR) treatment (Gold
et al., 2014). Those simulations find a similar disk structure and mass accretion rate
15
Figure 1.4: Left: the Ṁ variation as a function of time, which depends on different
mass ratios and details of the simulations. Right: the Fourier transform of each of
the Ṁ curves on the left. (Figure from Gold et al. (2014))
as the Newtonian case.
However, the interpretation of an observed variability period (or periods) may
not be straightforward, as those simulations of different mass ratios have shown
(D’Orazio et al., 2013; Gold et al., 2014). It has also been pointed out that emission
from the minidisks or circumbinary disk could also smear out the periodic signatures
(Gold et al., 2014; Noble et al., 2012), making it difficult to distinguish an SMBHB
system from a conventional AGN. Despite these caveats, if Ṁ is directly translated to
EM radiation, these SMBHBs would manifest as luminous AGNs that show periodic
variability in their photometric light curves, which have periods that are within a
factor of a few of the orbital periods.
16
1.2.2 Doppler Boosting as the Cause of Periodicity
Another possible cause of periodic variability is the relativistic Doppler boost-
ing of the emission from the minidisk, which is a model that was proposed by
D’Orazio et al. (2015) to interpret the periodicity of the SMBHB candidate PG
1302−102 (Section 1.3.2). The velocity of the secondary black hole is v2 = (2π/1 +
q)(GMtot/4π
2P )1/3, and for the binary parameters inferred for PG 1302−102, the
velocity is moderately relativistic. Even if the minidisk radiates steadily, the radi-
ation will be beamed towards the forward-moving direction of the secondary black
hole. From the viewpoint of the observer on Earth, the apparent flux F is therefore
modulated sinusoidally on the orbital frequency: ∆Fν/Fν = ±(3−α)vcos(φ/c)sin(i),
where α is the power-law slope of the quasar spectrum Fν ∝ να, vsin(i) is the line-of-
sight velocity, and φ is the orbital phase angle. As can be seen from the equations
above, this model requires that the black holes have unequal masses (q is small)
and that the binary orbit is viewed not far from edge-on (sin(i) is large), so that
the beaming effect is strong. However, one caveat is the presence of a circumbinary
disk, which may obscure the binary due to the high inclination angle and make it
difficult to observe this variability signature (d’Ascoli et al., 2018).
A corollary of the Doppler boost model is that variability amplitudes A in
different bandpasses are directly related to the spectral slopes in the respective
bands, i.e. AFUV/Aopt = (3− αFUV)/(3− αopt) and ANUV/Aopt = (3− αNUV)/(3−
αopt). If the slopes can be measured from the far-UV and near-UV spectra of the
source, the amplitude ratio would be useful supporting evidence for the Doppler
17
boost model of the putative binary (Charisi et al., 2018; D’Orazio et al., 2013).
1.3 The Search for SMBHBs
1.3.1 The Detectability of SMBHBs in a Time Domain Survey
If SMBHBs spend a significant fraction of their lifetimes at orbital separations
that correspond to periods of months to years, they could be detected as a population
of periodic quasars or AGNs in time domain surveys that repeatedly visit a large
number of objects over a period of a few years. This was explored by Haiman et al.
(2009). For a binary whose orbital decay is driven by GW emission, the decay
8/3
timescale is given by t = −r/ṙ ∝ r4GW ∝ torb (for a fixed black hole mass; r and torb
indicate the orbital radius and period, respectively), and the residence time that the
8/3
binary spends at the orbit is then tres = tGW ∝ torb. For a population of GW-driven
SMBHBs, their number fraction should therefore be proportional to their residence
times: fvar = tres/tQ, where tQ is the “lifetime” that the binary is “active” as a
quasar and is a few ×107 years (Martini, 2004). Given the power-law dependence
of tres on torb, the prediction is that a statistical sample of periodic quasars should
also follow a power-law relation:
f 7
−5/3
var = Norb/Ntot = (10 yr/tQ)[tvar/50.2(1 + z)weeks]
8/3M −17 qs , (1.1)
where M7 = MBH/10
7M, qs = 4q/(1 + q)
2, and tvar = torb(1 + z). On the other
hand, a population of SMBHBs that are at farther separations and are still in the
18
Figure 1.5: Left: The residence times of SMBHBs of different masses as a function
of orbital periods. Right: the sky coverage and the limiting magnitude required to
find a population of periodic quasars. The black line is calculated for the fiducial
values in the text. (Both figures from Haiman et al. (2009))
gas-driven stage of orbital decay should have a flatter power-law scaling relation
f αvar ∝ torb, where α < 8/3.
One may realize from Equation 1.1 that in a very large sample of quasars,
SMBHBs that manifest as periodic sources may not be so rare; it is therefore possible
to search for them in a large time domain survey. Haiman et al. (2009) calculated
the required survey area and magnitude depth in order to detect a statistical sample
of periodic sources by adopting a set of fiducial values, including a quasar lifetime
of tQ = 10
7 years, a mass ratio q = 1, an Eddington fraction fEdd = 0.3, and a
variability fraction ηvar = 0.1, and predicted e.g. 20 sources varying at tvar = 35
weeks (the redshift of the sources is fixed at z = 2) in a survey with a sky coverage
19
of 100 deg2 and a limiting magnitude of i = 24 mag. According to Equation 1.1,
the same survey volume would find e.g. 28/3× 20 = 127 sources varying at tvar = 70
weeks. There is also a trade-off between the survey area and the limiting magnitude
of the survey: a shallower magnitude depth would require a larger survey area
in order to probe the same sky volume. For a fixed black hole mass, this survey
requirement also depends on the values of q, ηvar, and fEdd: a smaller q results in
a longer tres and therefore requires a smaller volume to find the same number of
periodic quasars; increasing either ηvar or fEdd (or both parameters) would also lead
to a smaller required survey volume.
More excitingly, an optical time domain survey with a baseline of a few years
would be able to detect periods of a few months to a few years. For typical SMBH
masses of 106 − 109M, these timescales correspond to binaries in the GW-driven
regime of orbital evolution (a ∼ 10−3 pc or ∼ 102 gravitational radii rg). A system-
atic search in a large-field optical time domain survey can therefore potentially find
a population of close-separation SMBHBs in the GW-regime.
1.3.2 Systematically Searching for Periodically Varying Quasars
The predicted survey requirement by Haiman et al. (2009) is an excellent match
to the sky volume probed by existing large surveys such as the Catalina Real-time
Transient Survey (CRTS; Drake et al. 2009), the Palomar Transient Factory (PTF;
Law et al. 2009), and the Pan-STARRS1 survey (PS1; Kaiser et al. 2010), and
systematic searches for periodic quasars in those survey datasets have been made in
20
the past few years.
Graham et al. (2015b) performed a systematic search for periodic quasars in
the CRTS, which operates three telescopes from both hemispheres, covering a sky
area of ∼ 30, 000 deg2 to a magnitude depth of V ∼ 19−21.5 mag (depending on the
telescope) over a nine-year baseline. Their initial sample of quasars was first selected
by cross-matching a catalog of spectroscopically-confirmed quasars with the CRTS
data set. After applying the minimum sampling criterion and removing spurious
detections, the final sample consists of 243, 486 quasars.
Graham et al. (2015b) then applied a combination of wavelet transform and
autocorrelation techniques to select candidates that display periodic behavior; they
also required that the periodic variation persists for at least 1.5 cycles over the
observational baseline. 111 candidates passed their selection criteria. To calculate
the number of regular variable quasars that show periodicity by chance, they also
simulated a comparable number of mock light curves characterized by the Damped
Random Walk model (DRW; Kelly et al. 2009; see Section 1.4.1 for a review of
AGN variability and the DRW model); none of the simulated light curves passed
their selection criteria.
The detection rate of Graham et al. (2015b) is roughly compatible with the
predicted number of SMBHBs at low redshifts (Volonteri et al., 2009) and the frac-
tion of binaries as a function of orbital period (Haiman et al., 2009). Assuming
all periodic candidates are true SMBHBs, the inferred binary orbital separations
would be in the milli- to centi-parsec regime, and their GW strain amplitudes are
individually consistent with the current upper limit from the pulsar timing arrays.
21
Figure 1.6: The combined light curve of the SMBHB candidate PG 1302−102 from
LINEAR, CRTS, and the archival data. (Figure from Graham et al. (2015a))
The most significant candidate from the Graham et al. (2015b) sample of
periodic candidates is PG 1302−102 (Graham et al., 2015a). It has a median mag-
nitude of V = 15 mag and a variability amplitude of ∼ 0.14 mag. It was covered
by the CRTS and the Lincoln Near-Earth Asteroid Research (LINEAR) for a total
of ∼ 2 well-sampled cycles. Its light curve can be well-fitted to a sinusoid of a pe-
riod P = 1, 884 ± 88 days and modeled by the relativistic Doppler boosting of the
minidisk of the secondary black hole (D’Orazio et al. 2015; Section 1.2.2).
Charisi et al. (2016) performed their systematic search in the PTF survey,
which covered a total sky area of ∼ 8, 000 deg2 down to a magnitude depth R ∼ 20.5
mag over a ∼ four-year period. They cross-matched PTF sources with a catalog of
spectroscopically-confirmed quasars and pre-processed the light curves by removing
22
outlier data points. They further binned the measurements on the same night and
only selected light curves with dense temporal sampling (N > 50 in the R-band).
The number of quasars in their sample is therefore 35, 383.
They then searched for periodic signals using a generalized version of the
Lomb-Scargle (LS) periodogram (Lomb, 1976; Scargle, 1982). The LS periodogram
is a standard Fourier method for detecting periods in light curves with uneven sam-
pling, and for a range of trial frequencies, the one where the periodogram power
is maximum indicates a possible periodic signal at that frequency. As the null hy-
pothesis of the LS periodogram is pure white noise, they also generated mock light
curves modeled by the DRW process and compared their periodograms with those of
selected candidates; a periodic quasar candidate from their sample is considered sig-
nificant if its p-value < 4×10−5. 50 candidates satisfy this threshold and have more
than 1.5 cycles of variation. They then obtained light curves of the candidates from
the intermediate-PTF (iPTF, the extension of the PTF survey) and CRTS surveys
and recalculate the p-values of the extended baseline light curves. 33 candidates
have p-values < 1× 10−3 and are considered significant by their analysis.
Periodic quasar candidates discovered in those surveys are summarized in Ta-
ble 1.1. In comparison with the candidates selected by Graham et al. (2015b), the
Charisi et al. (2016) search is more sensitive to shorter periods and fainter sources.
In both searches, the distribution of binary residence times is most consistent with
the expected distribution (Haiman et al. 2009; Section 1.3.1) if all candidates have
the mass ratio q = 0.01, assuming all periodic quasar candidates from Charisi et al.
(2016) host SMBHBs and 25% of the candidates from Graham et al. (2015b) host
23
Table 1.1: Summary of Systematic Searches for Periodically Varying Quasars
Systematic Search Survey Initial Sample of Quasars Number of (Significant)
Periodic Candidates
Graham et al. (2015b) CRTS 243,486 111
Charisi et al. (2016) PTF 35,383 50 (33)
SMBHBs (Charisi et al., 2016).
However, AGN variability is stochastic (see Section 1.4.1) and can easily mimic
a periodic variation over a short baseline (see Section 1.4.2). Therefore, normal
AGNs that do not host SMBHBs may contaminate a periodic sample and present a
challenge for the search for robust SMBHB candidates via their optical variability.
A significant part of this thesis presents the extended baseline monitoring of the
periodic candidates from PS1 and a more rigorous method that selects a more robust
sample of binary candidates.
1.4 Quasar Variability and the Challenges of Searching for Periodic
Quasars
1.4.1 Optical Variability of Quasars
Active galactic nuclei (AGNs) and quasars are long known to be variable on
timescales from hours to years and in wavelengths from X-ray to radio (e.g. review
24
by Peterson (2001) and review on variability of blazars by Ulrich et al. (1997)).
While AGN variability has been measured since the age of photographic plates (e.g.
Matthews & Sandage 1963; Sillanpaa et al. 1988), modern, digital sky surveys have
produced a large amount of photometric data and have transformed the study of the
variable sky. Some examples of those surveys include the Sloan Digital Sky Survey
(SDSS; York et al. 2000), the Catalina Real-time Transient Survey (CRTS; Drake
et al. 2009), and Pan-STARRS1 (PS1; Kaiser et al. 2010), and some of the data
from those surveys will appear in this thesis.
Vanden Berk et al. (2004) used a sample of∼ 25, 000 spectroscopically-confirmed
quasars from SDSS with two epochs of photometry in three bands and measured the
dependence of variability on parameters such as luminosity, rest-frame wavelength,
and the time lag between the epochs. It was the largest sample of variable quasars
at that time. By measuring the structure function (SF), they find that the vari-
ability amplitude increases with time lag; they also find that quasars vary more at
shorter rest wavelengths, and more luminous quasars are less variable. The Stripe
82 survey of SDSS, which repeatedly visited a narrow stripe of the sky, allowed one
to move beyond few-epoch photometry and produced a large number of light curves,
the measurements of the brightness as a function of time. This has enabled one to
study individual variable quasars in greater quantitative detail. For example, Sesar
et al. (2007) found that ∼ 90 % of the quasars in their Stripe 82 sample vary at
the 0.03 mag level and 30 % vary at the > 0.1 mag level; the highly variable nature
of quasars also make it as good a selection algorithm as the color-magnitude or
color-color diagram.
25
The light curve variability of an AGN is stochastic and aperiodic (e.g. Peterson
2001), i.e. the variation appears to be random and does not repeat itself regularly.
Motivated by the potential to study AGN physics via its variability, Kelly et al.
(2009) modeled the AGN variability as a stochastic process of two parameters,
the characteristic timescale and the variability amplitude. This is known as the
Damped Random Walk (DRW) process, or an Ornstein-Uhlenbeck (O-U) process
(Gillespie, 1996). The power spectral density (PSD), or how much power is at a
certain frequency, of DRW is:
2σ2τ 2
P (f) = , (1.2)
1 + (2πτf)2
where τ is the characteristic timescale and σ is the variability amplitude on short
timescales. As one can see from Equation 1.2, the variability power decreases as f−2
at frequencies f > (2πτ)−1 and flattens towards white noise at lower frequencies. By
modeling quasar light curves as the DRW process, Kelly et al. (2009) find that the
characteristic timescale τ is correlated with black hole mass and luminosity of the
quasar but the short-timescale variability amplitude σ is anti-correlated with those
two physical parameters. The characteristic timescales, ranging between ∼ 10 days
to ∼ 10 years, are consistent with the thermal or orbital timescale, and the success
of DRW in fitting light curves also suggests that quasar variability is driven by
stochastic processes, such as fluctuations in the disk caused by a turbulent magnetic
field (Kelly et al., 2009).
The SF ensemble analysis of Vanden Berk et al. (2004) was based on the as-
26
Figure 1.7: Left: light curves simulated from the Damped Random Walk model
with different characteristic timescales. Right: the power spectra of the light curves
on the left. (Both figures from Kelly et al. (2009))
sumption that a large statistical sample of quasars (each was observed at two epochs)
would reveal the same properties as well-sampled light curves of individual quasars.
In order to understand the individual and ensemble properties of quasar variabil-
ity, MacLeod et al. (2010) applied the DRW model to ∼ 9, 000 spectroscopically-
confirmed quasars in the SDSS Stripe 82 survey and confirmed previous results that
DRW is a good model for quasar variability (Kelly et al., 2009; Koz lowski et al.,
2010). Since the structure function SF∝ ∆tβ (where ∆t is the time lag) is equiv-
alent to P (f) ∝ fα, where α = −2β − 1 (Bauer et al., 2009), the power-law slope
β = 0.5 is therefore steeper than what was previously measured by Vanden Berk
27
et al. (2004) (β = 0.3).
Despite the success and popularity of the DRW model in describing AGN vari-
ability, recent studies have found deviations from the DRW model in their analyses
of AGN light curves from ground-based surveys and high photometric precision light
curves produced by Kepler, the exoplanet hunting mission. Two results from these
studies (e.g. Aranzana et al. 2018; Caplar et al. 2017; Mushotzky et al. 2011; Smith
et al. 2018a) are: (1) there is a variation of power-law PSD slopes; and (2) those
slopes are in general steeper than α = −2, the high frequency slope of the DRW
model, and range between −2 and −3 (Aranzana et al., 2018; Smith et al., 2018a).
Understanding the variability of normal quasars (i.e. single black hole-powered)
therefore has farther-reaching implications than studying quasar accretion physics;
it also affects our search for periodic quasars that could host SMBHBs and our
estimation of their statistical significance, as I will discuss in Section 1.4.2 and
throughout this thesis.
1.4.2 Red Noise Contamination of the SMBHB Candidate Sample
There have been ∼ 200 SMBHB candidates reported based on their period-
icity (see Section 1.3.2). However, normal AGN variability is stochastic, and the
fluctuation of “red noise” can easily mimic a periodic variation over a few cycles,
especially when the photometric uncertainty is large and the sampling is uneven. To
demonstrate this, Vaughan et al. (2016) generated mock light curves under two dif-
ferent red noise models, a DRW model and a broken power law (BPL) PSD model;
28
the choice of the power-law slopes of the BPL model is based on the power spec-
trum analysis of Zw 229-15 (Edelson et al., 2014), and the break frequency fbr is
scaled down from that of Zw 229-15 based on a MBH−fbr relation. They simulated
100, 000 light curves under each model and looked for periodic candidates that have
comparable goodness-of-fit as PG 1302−102 and a similar period range as the Gra-
ham et al. (2015b) periodic sample (Section 1.3.2). They find “periodic candidates”
at a not-insignificant rate of one to two per 1, 000 simulations under both models,
and the period distribution of those false positives is similar to that of the Graham
et al. (2015b) candidate (i.e. more candidates at longer periods), suggesting that
one tends to find false periodicities with long periods. They also find that the rate
of false periods is highly dependent on the choice of the power law slopes and the
break frequency; for example, if the high-frequency slope is α = −3 instead of −2,
the number of false positives increases by a factor of a few (Vaughan et al., 2016).
Therefore, in a large time domain survey of ∼ 105 quasars (e.g. the Graham et al.
2015b CRTS sample), the effect of red noise contamination cannot be ignored, and
estimating the statistical significance of selected periodic candidates also requires a
close inspection of normal quasar variability.
Fortunately, as Vaughan et al. (2016) suggested, the detected periodic can-
didates can be easily tested with further observations; in the long term (at least
∼ 5 cycles), a true periodic signal should persist, while stochastic variability does
not. I will put the periodic candidates from PS1 to the test with extended baseline
imaging in this thesis (Chapters 3 - 5).
Sesana et al. (2017) tested the binary hypothesis of the periodic candidates
29
from Graham et al. (2015b) and Charisi et al. (2016) from the perspective of the
pulsar timing arrays (PTAs). In order to build a gravitational wave background
(GWB) from a population of SMBHBs, it first requires the knowledge of the chirp
mass of each candidate. Since the black hole mass measured from the quasar spec-
trum is taken to be the total mass Mtot of the system, they assigned mass ratios
q to candidates using three different mass ratio models: (1) the distribution of q is
lognormal, peaking at log q = 0 (but only for q < 1, by definition); this is motivated
by the expected q distribution of SMBHBs formed in major galaxy mergers; (2) the
distribution of q ranges from < 0.01 to 1; this distribution peaks at ∼ 0.1 but has a
minimum cut-off at q = 0.05; this is representative of the modulated mass accretion
model, which requires a larger mass ratio (Section 1.2.1); (3) same as Model (2),
but without the cut-off; this model is representative of the Doppler boosting case,
which requires a small mass ratio (Section 1.2.2). Since the black hole mass MBH
estimates also have a systematic uncertainty of ∼ 0.3 dex (e.g. Shen et al. 2008),
they also consider two MBH scenarios: one assumes the measurements are exact,
and the other includes the uncertainties. They therefore built a GWB under six
different models.
They then considered a number of effects due to incompleteness of the search,
including: (1) only consider binary candidates at z < 1.3, at which the Graham
et al. (2015b) periodic sample should be less incomplete; (2) include candidates at all
redshifts; (3) correct for the incompleteness of sky coverage. Since interactions with
stars and gas can add eccentricity to the binary and therefore reduce the GW signal
at low frequencies, they have also considered this effect. After those considerations,
30
they find that the GWB inferred from the Graham et al. (2015b) sample is in tension
with the PTA upper limits, i.e. the inferred amplitude is higher by at least a factor
of a few than the current upper limits reported in Shannon et al. (2015), Lentati
et al. (2015), and Arzoumanian et al. (2016), and the periodic sample from Charisi
et al. (2016) similarly infers a much larger GWB amplitude. While the tension
can be alleviated by e.g. fixing all mass ratios at q = 0.01 and assuming all black
hole masses are overestimated, this study places a hard, independent constraint on
the SMBHB samples, suggesting that many of the periodic candidates may be false
positives. As they also find by calculating individual contributions from the binary
candidates, only misidentifying a few loudest sources would have a large effect on
the inferred GWB amplitude. Down-selecting a more robust sample of SMBHB
candidates is therefore one of the main objectives of this thesis.
This thesis is organized as follows: Chapter 2 presents the most significant
periodic candidate from our pilot study in MD09, one of the fields in the Pan-
STARRS1 Medium Deep Survey (PS1 MDS). Chapter 3 presents a full analysis of
the PS1 field MD09 and follows up the candidates using extended imaging with the
Discovery Channel Telescope and archival data, in order to identify false positives
that are due to the stochastic variability (“red noise”) of normal quasars that do
not host SMBHBs. To further address the issue of red noise contamination of
periodic candidate samples, I set up the expectations for a true periodic signal (in
the presence of red noise) in Chapter 4 and re-analyze the well-known SMBHB
candidate PG 1302−102 using an extended light curve and a more rigorous method.
Those pieces in Chapter 2 - 4 are put together in Chapter 5, where I expand the
31
systematic search to all ten fields of PS1 MDS. Our more robust subsample of
periodic quasar and SMBHB candidates from PS1 MDS can then be followed up
with multi-wavelength observations in order to search for other electromagnetic
signatures that could independently verify the binary hypothesis (Chapter 6).
32
Chapter 2: A Periodically Varying Luminous Quasar at z = 2 from
the Pan-STARRS1 Medium Deep Survey: A Candidate
Supermassive Black Hole Binary in the Gravitational Wave-
Driven Regime
2.1 Introduction
The expectation for the existence of supermassive black hole binaries (SMB-
HBs) in galaxy nuclei is supported by two well-established properties of galaxies:
(1) high spatial resolution observations of nearby galaxies have demonstrated that
SMBHs are ubiquitous in the centers of galaxy bulges (Magorrian et al., 1998) with
masses tightly correlated with the mass and structure of their host galaxies (Fer-
rarese & Merritt, 2000; Gebhardt et al., 2000; Graham et al., 2001), and (2) galaxies
in a ΛCDM Universe build up their structure hierarchically through mergers (e.g.
Springel et al. 2005). When two galaxies merge, their SMBHs will sink to the center
through dynamical friction and through three-body interactions with stars and vis-
cous exchange of angular momentum with circumbinary gas form a gravitationally
bound binary that eventually coalesces due to the radiation of gravitational waves
(GWs; Begelman et al. 1980).
33
Recent progress has been made in the detection of “dual active galactic nuclei
(AGNs),” double active nuclei in assumed merged galaxy systems with kiloparsec-
scale separations (Comerford et al., 2009b; Komossa et al., 2003). These dual AGNs,
while a product of a galaxy merger, are not yet gravitationally bound, and thus are
not necessarily fated to coalesce. A true SMBHB becomes gravitationally bound on
the scale of parsecs, which, beyond our Local Group of galaxies, is well below the an-
gular resolving power of the most powerful current, or even future, telescopes. How-
ever, several promising candidates have been identified indirectly via spectroscopy:
quasars with offset and/or drifting broad-line peaks attributed to a broad-line re-
gion in orbit around an SMBH’s binary companion (Barrows et al., 2011; Boroson
& Lauer, 2009; Dotti et al., 2009). However, alternative scenarios have been pro-
posed that do not require an SMBHB, including double-peaked lines from a single
accretion disk (Chornock et al., 2010).
A promising observational signature of SMBHBs is their variable accretion
luminosity. One of the first sub-parsec SMBHB candidates, OJ 287, was identified
by its variability behavior (Lehto & Valtonen, 1996). OJ 287 is a quasar that
demonstrates regular optical outbursts on a timescale of 12 yr that has been modeled
as the result of a secondary SMBH companion passing through the primary SMBH’s
accretion disk (Valtonen et al., 2008). Such a configuration should be rare since
the secondary BH’s orbital axis must be highly misaligned with the primary BH’s
accretion disk axis in order for it to pass through its disk. A more generic signature
of an SMBHB is likely to be related to accretion through its circumbinary disk.
In a gas-rich galaxy merger, strong gravitational torques drive gas inward,
34
triggering both star formation and BH accretion (Hopkins et al., 2006). In par-
ticular, hydrodynamical simulations of circumbinary disks show that accretion via
“hot streams” onto the BHs is strongly modulated by the binary’s orbital motion
for mass ratios of 0.05 < q ≤ 1 (D’Orazio et al., 2013; Gold et al., 2014; MacFadyen
& Milosavljević, 2008; Noble et al., 2012; Shi et al., 2012). Simulations (D’Orazio
et al., 2013; Farris et al., 2015) also detect a t ∼ 6 torb timescale originating from
a surface density “lump” just outside the central cavity of the circumbinary disk,
which is a persistent but secularly evolving feature.
While working on our paper, we became aware of the report of PG 1302−102, a
periodically variable quasar discovered by the Catalina Real-Time Transient Survey
(CRTS; Graham et al. 2015b). It is a 15 magnitude quasar at z = 0.2784, varying
at the 0.14 mag level with a period of 5.2 ± 0.2 yr, with good sampling over 1.8
cycles and extended archival data going back 20 yr. Their physical interpretation
for its variability is an SMBHB (log(M/M) ∼ 8.5, a ∼ 0.01 pc), its luminosity
being modulated due to either a precessing jet or an overdensity (“hot spot”) in the
inner edge of its circumbinary disk.
In this Letter, we present our most significant detection from the systematic
search for periodically varying quasars in the Pan-STARRS1 Medium Deep Survey
(PS1 MDS) field MD09. PSO J334.2028+01.4075 is a radio-loud quasar at z =
2.060 with archival spectroscopy from FQBS and extended baseline photometry
from CRTS. The 8.5 yr baseline of the PS1+CRTS light curve is well described by a
simple sinusoid, consistent with theoretical simulations for the modulated accretion
rate in a 0.05 < q < 0.25 mass-ratio SMBHB. We use the rest-frame period and
35
virial black hole mass estimate to infer an orbital separation of the binary that is in
the GW-driven regime.
2.2 Theoretical Predictions
The dynamics of a gravitationally bound SMBHB system can be described by
Kepler’s third law:
( )( )3/2
M a
torb = 0.88 yr , (2.1)
107M 103Rs
where Rs is the Schwarzschild radius, R =
2GM
s 2 ; a is the separation between thec
BHs; and M is the total mass of the system. Haiman et al. calculate the minimum
survey area to detect a statistically significant sample of quasars powered by SMB-
HBs as a function of variable magnitude depth, assuming reasonable values for the
quasar luminosity function, quasar lifetime (tQ = 10
7 yr), the Eddington fraction
(fEdd = 0.3), and the fractional variability amplitude (∆f/f = 0.1) (Haiman et al.,
2009). The variability detection threshold we have achieved in MD09 (Section 2.4)
corresponds to a ∆f/f > 0.1 sensitivity for point sources brighter than m ∼ 21 mag,
and thus a variable magnitude of 23.5 mag. At this depth, Haiman et al. (2009)
require an area of ∼100 deg2 to yield a sample of over 100 SMBHBs (Haiman et al.,
2009); an excellent match to the area of PS1 MDS (80 deg2). Furthermore, the
baseline (4.2 yr) and cadence (3 days) of PS1 MDS make us sensitive to timescales
for which BHs (with M > 107M) are in the GW-driven regime of orbital decay.
36
2.3 Pan-STARRS1 Medium Deep Survey
The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS)
is a wide-field imaging system designed for dedicated survey observations on a 1.8m
telescope on Haleakala, Hawaii, with a 1.4 Gigapixel camera and a 7 deg2 field
of view (Kaiser et al., 2010). The PS1 telescope is operated by the Institute for
Astronomy (IfA) at the University of Hawaii and has just completed over 4 yr of
operation in 2014 March. We present data from the Medium Deep Survey (MDS),
a deep, multi-epoch survey of 10 circular fields distributed across the sky, each ∼ 8
deg2 in size, whose daily observing cadence in five filters is excellent for studying
persistent variable sources, including quasars. The PS1 MDS cadence of observation
cycles through the gP1, rP1, iP1, and zP1 bands every three nights, with observations
in the yP1 band close to the full Moon. Due to the poorer time sampling of the yP1
observations, we do not use them in this analysis.
2.4 Ensemble Photometry
We began our systematic search for SMBHB candidates among color-selected
quasars in the PS1 MD09 field. This is the first MD field that was made available to
the PS1 Science Collaboration in the Pan-STARRS Science Interface (PSI) online
database. In order to maximize our sensitivity to intrinsic variability, we first applied
the technique of differential ensemble photometry. This technique is able to correct
for local systematic errors due to variable atmospheric conditions by comparing a
37
target object with nearby non-variable stars (Bhatti et al., 2010; Honeycutt, 1992).
We created a color-selected reference sample and quasar sample by cross-matching
point sources (m < 23 mag) in the MD09 field with a custom catalog extracted from
full-survey deep stacks from PS1 MDS in the gP1, rP1, iP1, zP1, and yP1 bands, as
well as from observations with the Canada–France–Hawaii Telescope (CFHT) in the
u band (S. Heinis et al. 2015, in preparation). We converted all magnitudes to the
SDSS photometric system (Tonry et al., 2012b) to take advantage of the SDSS color
selection of stars and quasars already available in the literature (Schmidt et al.,
2010; Sesar et al., 2007). Figure 2.1 shows the color-color diagrams of the point
sources in MD09 selected as quasars and non-variable stars. This query resulted in
8158 reference stars and 316 quasars, each with an average of 350 detections in four
filters.
We modified the ensemble photometry software developed by Bhatti et al.
(2010) for SDSS to the PS1 data format and ran it on the reference stars. In Figure
2.1 (bottom right panel), we plot the “corrected” magnitude error as a function of
mean magnitude compared with the “raw” values before ensemble photometry. The
ensemble photometry reduces the measured errors significantly, lowering the error
floor from 0.045 to 0.025 mag, and resulting in a 2σ variability threshold of 0.05
mag on the bright end to 0.34 mag on the faint end.
38
2.5 Selecting Periodic Quasar Candidates
We then applied ensemble photometry to the 316 color-selected quasars and
flagged quasars as variable based on their magnitude error relative to their neighbors
of a similar brightness; we set 2σ as our criterion for variability and required a
variability flag in at least two filters. This selection yielded 168 variable quasars
in MD09. Among these color-selected variable quasars, we searched for potential
periodic signatures using the Lomb–Scargle (LS) periodogram, a Fourier analysis
technique of unevenly spaced data with noise (Horne & Baliunas, 1986; Lomb, 1976;
Scargle, 1982). For N0 data points in the time series spanning a total length of
T in units of MJD, we sampled the periodogram at the number of recommended
independent frequencies (Ni) from Horne & Baliunas (1986), from 1/T to N0/(2T )
(which would be the Nyquist frequency if data were evenly sampled), resulting in
a frequency resolution in the periodogram of ∆f = (N0/2 − 1)/(TNi) ∼ 2 × 10−4
d−1. When identifying periodic sources, we took advantage of the redundancy of
PS1 MDS monitoring in four filters (gP1 rP1 iP1 zP1), each with a slightly different
observing cadence due to weather and scheduling constraints to help filter out false
detections from aliasing by requiring that periodogram peaks be coherent across
multiple filters; 40 of the 168 variable quasars survived this test.
39
2.6 Periodic Quasar Candidate PSO J334.2028+01.4075
Among the candidate periodic quasars from our periodogram analysis, here we
focus on our most significant detection, PSO J334.2028+01.4075 (J2000). In Figure
2.3, we show its periodogram in four PS1 filters, with the strongest peak marked
with a dashed line. We fold each filter light curve on this period (Figure 2.2), and
measure the scatter of the residuals from the best-fit sine curve (σr). The error of
the periodogram peak frequency (δf) and the signal-to-noise ratio of the peak power
(ξ) can then be calculated from σr and the amplitude of the signal A0 (Horne &
Baliunas, 1986) as δf = √ 3σr (which gives us an error on the detected period of
4 N0TA0
δP = δf/f 2) and ξ ≡ A2 20/(2σr), respectively.
The resulting average period across all four filters is P = 541.8 ± 15.3 days,
with the highest signal-to-noise ratio in the gP1 filter with ξ = 3.19, and periodogram
peaks in all four filters well above a 1.5×10−23 false-alarm probability (corresponding
to 10σ) plotted with a dotted line in Figure 2.3. The PS1 data cover 2.6 cycles,
just shy of the “rule of thumb” number of cycles (three) for a periodic variation
to be apparent to the eye (Press, 1978). From our Monte Carlo simulations of
1000 stochastic Damped Random Walk (Kelly et al., 2009) light curves, we find a
false periodic detection rate of 6.3% using our selection criteria from Section 2.5.
We further disfavor a false alarm from stochastic quasar variability since the 0.6%
of the simulations that successfully mimic the periodic timescale of our candidate
have short-timescale variances a factor of ∼ 2 larger than expected for the quasar’s
luminosity and inferred black hole mass (Section 2.7). Note that there is a secondary
40
peak in the gP1 and rP1 periodograms that, if real and not an artifact from the PS1
data sampling, is at twice the primary peak frequency, a signature of 0.05 < q <
0.25 mass-ratio SMBHBs, which show an accretion rate modulation most closely
described by a simple sinusoid (D’Orazio et al., 2013). The amplitude of PSO
J334.2028+01.4075’s sinusoidal modulation increases with decreasing wavelength,
consistent with the exponential dependence on wavelength found in previous quasar
variability studies (e.g. MacLeod et al. 2010; Vanden Berk et al. 2004).
PSO J334.2028+01.4075 is a radio-loud quasar (FBQS J221648.7+012427)
with an archival spectrum from the FIRST Bright Quasar Survey (Becker et al.,
2001). We are also fortunate that this candidate has an archival V-band light curve
from CRTS (Drake et al., 2009), which we use to test the persistence of the periodic
variation over an extended baseline of 8.5 yr (corresponding to 5.7 cycles). To
compare to the CRTS light curve, we convert the PS1 gP1-band light curves to the
SDSS system (Tonry et al., 2012b) and then to the Johnson V magnitude using
the photometric transformation for quasars from Jester et al. (2005) and an average
gP1-rP1 = +0.10 mag. We had to apply an additional offset of −0.17 mag to the
pseudo-V PS1 magnitudes in order to match the average of the CRTS V -band data.
Though the photometric errors are relatively large, the CRTS measurements are
consistent with those of PS1 during their overlap (Figure 2.4), and have residuals
over the entire CRTS baseline from the PS1-fitted sinusoidal model that are Gaussian
with a σ = 0.17 mag that is comparable to the mean photometric error of 0.18 mag.
41
2.7 Physical Interpretation
We use the width of the quasar’s C IV line and its nearby continuum luminosity
to make a virial estimate of the black hole mass from Vestergaard & Peterson (2006):
(M ) [( )2( )0.53]BH FWHM(CIV ) λLλ
log = log + 6.66 (2.2)
M 1000km/s 1044erg/s
where λ = 1350 Å . The FBQS spectrum, though not publicly available in an
electronic format, was measured with a ruler to determine Fλ,obs(1350Å (1+z)) ∼
8.5 × 10−17 ergs s−1cm−2 Å−1, and FWHM (CIV λ1550) ∼ 200 Å. The CIV line
is symmetric in shape, and its width corresponds to a velocity in the rest-frame of
∼ 12, 650 km s−1. We correct for a Galactic extinction of E(B − V ) = 0.0406 mag
(Schlafly & Finkbeiner, 2011), using the extinction law from Cardelli et al. (1989),
to find Lλ,em(1350Å) = 4πd
2F 10Aλ/2.5L λ,obs (1 + z) ∼ 9.5× 1042 erg s−1Å−1, where dL
is the luminosity distance assuming H0 = 70 km/s/Mpc, Ωm =(0.3, )ΩΛ = 0.7, and
k = 0, and Aλ = 0.1790. This results in a black hole mass of log
MBH = 9.97, with
M
a scatter from the uncertainty in the relation of 0.5dex. Applying a mean quasar
bolometric correction at 1350Å of BC=3.81 from Richards et al. (2006a), one gets a
bolometric luminosity of Lbol = λLλBC = 4.9× 1046 erg s−1. Note that this object
is also radio loud, with a radio luminosity at the rest-frame frequency of 5 GHz of
log(LR(erg s
−1)) = 32.8 from Becker et al. (2001).
Assuming the rest-frame period Prest = Pobs/(1+z) is on the order of the orbital
timescale of the SMBHB, with a caveat that in addition to a strong dependence
on mass ratio, there are a range of theoretical predictions for translating Prest to
42
torb (e.g. MacLeod et al. 2010; Vanden Berk et al. 2004), we then calculate the
orbital separation of the binary to be 7+8−4Rs (∼ 0.006+0.007−0.003 pc), securely placing
it in the gravitationally bound regime of a physically viable SMBHB system — a
circumbinary accretion disk system capable of maintaining a central cavity, stable
to gravitational fragmentation, and in the regime of orbital decay driven by GWs
(e.g. D’Orazio et al. 2013; Haiman et al. 2009; Kocsis et al. 2012). Also note
that since the viscous time scales as r2, one could expect to be able to detect
modulations in the accretion rate fed by the streams in the circumbinary disk cavity,
without being washed out by viscous processes in the “minidisks” around each BH
(Roedig et al., 2014). Remarkably, the rest-frame inspiral time for the binary is
t = 5 c
5 a4 = 7.0 yr (a/7R )4(M/109.97M )−3 for q = 0.25, where µ ≡ M1M2insp 256 G3 M2µ s  ,M1+M2
opening up the possibility for detecting the decay of the orbital period (Ṗ ) with
future monitoring, as well as providing a promising target for direct GW detection
for pulsar timing arrays (Sesana et al., 2009).
2.8 Discussion and Conclusions
We present the most statistically significant periodically variable quasar can-
didate from our search in PS1 MD09, PSO J334.2028+01.4075, a radio-loud quasar
at z = 2.060. We combine an estimate of its black hole mass with its variability
timescale (assuming Prest ∼ torb) to find orbital parameters consistent with model
predictions of a stable accreting SMBHB system with a 0.05 < q < 0.25 in the
GW-driven regime (Haiman et al., 2009).
43
The redshift of this SMBHB candidate coincides with the peak epoch for
SMBHB mergers (Volonteri et al., 2003), and its large mass (M ≈ 1010M) makes
it favorable for detection in the GW-driven regime, due to the strong dependence on
M of the residence time at a given orbital separation in units of Rs (Haiman et al.,
2009). Like the CRTS SMBHB candidate PG 1302−102 reported by Graham et al.
(2015a), our SMBHB candidate is also a radio-loud quasar. However, given the
shorter rest-frame period of our candidate of 0.5 yr (versus 4 yr in PG 1302−102),
it is even more unlikely for its variability to be driven by jet precession, either orig-
inating from a single SMBH (Lu & Zhou, 2005) or a binary SMBH (Lobanov &
Roland, 2005).
This pilot program in PS1 MD09 is a promising start to our systematic search
for periodic variability signatures of SMBHBs amongst the expected ≈ 1000 variable
quasars in the full ∼ 80 deg2 PS1 MDS. At the start of the next decade (∼ 2023), the
Large Synoptic Survey Telescope (LSST) (Ivezic et al., 2008) will probe a volume
several thousand times larger than PS1 MDS, yielding tens of millions of quasars,
and potentially thousands of SMBHBs periodically varying on the timescale of years,
fated to coalesce.
44
Figure 2.1: Color–color diagrams used to select point sources in MD09 for the
8158 point sources in the reference star sample (red) and 316 point sources in the
quasar sample (blue). Photometry is measured from the CFHT+PS1 catalog and
converted to the SDSS system. Dashed lines show the color selection boundaries.
The stellar color-color selection box was chosen to avoid RR Lyrae stars, which
are intrinsically variable. Bottom right panel: observed standard deviation of the
reference sample of non-variable stars before and after applying the technique of
ensemble photometry (open circles and filled squares, respectively), compared to
the Poisson error expected from the reported PS1 flux errors (black dashed lines).
45
Figure 2.2: Upper panels: sinusoidal fit to the folded PS1 light curve of PSO
J334.2028+01.4075 in four filters, with the error bar indicating the typical pho-
tometric error for an object of similar brightness in that filter. The period corre-
sponding to the peak of the periodogram and its error bar, the amplitude of the
fitted sine wave, and the signal-to-noise ratio of the peak power are each labeled.
Lower panel: sinusoidal fit plotted over the complete PS1 light curves in the gP1 rP1
iP1 and zP1 bands. The data used to create this figure are available.
46
Figure 2.3: Our automated LS periodogram routine selects periodically variable
candidates by requiring that the strongest peak is detected at the same frequency in
at least three filters. This quasar candidate, PSO J334.2028+01.4075, was selected
through this method and had the periodogram peak with the highest signal-to-noise
ratio of all of our candidates. The dashed lines mark the strongest peak in each
filter. The dotted line corresponds to a false-alarm probability of 1.5× 10−23, or 10
σ.
47
Figure 2.4: PS1 light curve (asterisks) converted to the V band to compare to the
archival CRTS light curve (dots with error bars). The CRTS data points are binned
in one-day intervals, with the error bars measured from the standard deviation in
the bin and not including data points with a photometric error greater than 0.25
mag or nights with less than three measurements. This results in 34/113 nights of
data being thrown out. The CRTS measurements are overall consistent with the
PS1 light curve and the sine fit to the PS1 light curve (dashed curve).
48
Chapter 3: A Systematic Search for Periodically Varying Quasars in
Pan-STARRS1: An Extended Baseline Test in Medium
Deep Survey Field MD09
3.1 Introduction
Supermassive black holes (SMBHs) appear to be at the centers of most, per-
haps all, massive galaxies (e.g. Kormendy & Richstone 1995). Thus, when two
massive galaxies merge in the ΛCDM Universe, it is expected that their nuclei will
form a supermassive black hole binary (SMBHB; e.g. Springel et al. 2005). As the
binary coalesces, the early stage of its orbital decay is driven by exchanging angular
momentum with the circumbinary gas disk through viscosity; at smaller separations
(a < 1 pc), its orbital decay becomes more dominated by gravitational wave (GW)
radiation (e.g. Begelman et al. 1980).
However, sub-parsec separation SMBHBs at cosmological distances are too
compact to resolve with current, or even future, telescopes. Indirect searches so
far, therefore, have been focused on spectroscopy, looking for offset broad lines that
suggest two broad line emission regions, each likely associated with each black hole
in the binary system (Boroson & Lauer, 2009), or offset or shifted peak of the broad
49
line region (e.g. Dotti et al. 2009; Eracleous et al. 2012).
Another observational aspect of SMBHBs, however, was much under-exploited
until recently — their potential optical variability. One of the first sub-parsec
SMBHB candidates identified via its variability was OJ 287 (Sillanpaa et al., 1988),
which showed quasi-periodic optical outbursts at intervals of 12 years, with the
physical interpretation of the burst being the secondary black hole passing through
the accretion disk of the primary (e.g. Lehto & Valtonen 1996; Valtonen et al.
2011, 2008). More recently, another sub-parsec SMBHB candidate, PG 1302−102
(Graham et al., 2015a), was discovered by the Catalina Real-time Transient Survey
(CRTS; Drake et al. 2009). Its V -band light curve can be fitted to a sinusoidal func-
tion with period of 1,884 days and amplitude of 0.14 mag. A physical interpretation
of PG 1302−102’s periodic variability is relativistic Doppler boosting (D’Orazio
et al., 2015): in this scenario, where the luminosity is dominated by the steadily
accreting secondary black hole and the system is viewed at a high inclination angle,
emission from the minidisk of the secondary is Doppler-boosted as the black hole
orbits at a moderately relativistic speed (along the line of sight).
Another possible scenario that could give rise to periodic variability is mod-
ulated mass accretion in the system. Simulations of an SMBHB embedded in a
circumbinary disk show that although the binary tidal torque clears and maintains
a low gas density cavity at radius < 2a (where a is the binary separation), mate-
rials can penetrate the cavity through a pair of streams and be accreted onto the
binary. These simulations have the similar results that for a mass ratio 0.01 < q < 1
— as expected in the merger of two massive galaxies — mass accretion through
50
the circumbinary disk is strongly modulated as a result of the binary’s orbital mo-
tion within the circumbinary disk, including two-dimensional (2D) hydrodynamical
(MacFadyen & Milosavljević, 2008), 3D Newtonian magnetohydrodynamical (MHD;
Shi et al. 2012) and Post-Newtonian MHD (Noble et al., 2012) for an equal mass bi-
nary, and general relativistic (GR) MHD (Gold et al., 2014) and 2D hydrodynamical
simulations (D’Orazio et al., 2013) for various mass ratios. In these simulations, the
accretion rate varies on a time scale that is on the order of the binary orbital time
scale, which is in turn a function of the total black hole mass and orbital separation
by virtue of Kepler’s law. Assuming that luminosity tracks mass accretion of the
circumbinary disk, the former should then vary as the latter varies. For a typical
black hole mass of 107M and typical separation 10
3Rs, the orbital period is on
the order of ∼ year, an(observa)t(ionally)feasible time scale for current time-domain3/2
surveys: torb = 0.88 yr
M a
107 3
(where Rs is the Schwarzschild radius:M 10 Rs
Rs = 2GM/c
2).
These theoretically explored variability signatures of an SMBHB, as well as
encouraging predictions for the detection rates of periodically varying quasars from
SMBHBs in a cosmological context (Haiman et al., 2009), motivated several recent
systematic searches in large optical time-domain surveys with a temporal baseline of
several years — Graham et al. (2015a,b), with the CRTS; Charisi et al. (2016), with
Palomar Transient Factory (PTF) and additional data from intermediate-PTF and
CRTS; Zheng et al. (2016), with the Sloan Digital Sky Survey (SDSS) and CRTS;
and Liu et al. (2015), with the Pan-STARRS1 Medium Deep Survey (MDS).
In our pilot study (Liu et al. 2015, hereafter L15), we performed a systematic
51
search for SMBHB candidates in MDS’s MD09 field and reported our first significant
detection of such a candidate, PSO J334.2028+1.4075. As reported in L15, PSO
J334.2028+1.4075 has a coherent period of P = 542±15 days in gP1 rP1 iP1 zP1
filters, corresponding to almost 3 cycles of variation that is well fitted to a sinusoidal
function. It also has an archival V -band light curve from CRTS (Drake et al., 2009).
Even though the photometric precisions are not comparable, the CRTS light curve
is consistent (in the residual sense) with the PS1 only (PV1) sinusoidal fit over ∼ 9
years, or ∼ 6 cycles. It is also a radio loud quasar (R = log (f5GHz/f ) = 2.30;
2500Å
Becker et al. 2001) from the VLA FIRST catalog (FIRST J221648.6+012427; White
et al. 1997).
Since then, we have repeated our analysis of MD09 with data Processing Ver-
sion 2 (PV2) which was made available late-2014 and includes extra data from the
final phase of the PS1 survey (Figure 3.1). We find three periodic quasar candidates
that satisfy our selection criteria: a coherent period in at least three filters, an S/N
for a sinusoidal fit of > 3 in at least one filter, and a variation over at least 1.5 cy-
cles. In addition, we use extended baseline data (from archival and new monitoring
observations) to test the persistence of our periodic candidates over 5 − 12 cycles.
Recently, it has been pointed out by Vaughan et al. (2016) that the intrinsic red
noise (increasing power at lower frequencies) characteristic of quasar variability can
easily mimic periodic variability over a small number of cycles, and they emphasize
the importance of demonstrating persistence of periodicity over > 5 cycles.
This paper thus presents our detailed analysis with MD09 PV2 and is organized
as follows. In Section 3.2 we introduce the time domain data set used in this study:
52
MD09 from the Pan-STARRS1 MDS In Section 3.3 we describe our methods of
variability selection and periodicity search; we also discuss our biases in selecting
variable active galactic nuclei (AGNs) in a flux-limited survey like PS1 MDS. In
Section 3.4, we test the persistence of the candidates’ periodicity with archival light
curves and follow-up imaging. In Section 3.5, we measure the black hole mass of
binary candidates and calculate their inferred binary parameters. Finally, in Section
3.6, we conclude with implications for searches for periodic quasars in a large time-
domain survey. Throughout this paper, we adopt cosmological parameters for a flat
universe: Ωm = 0.3, Ωλ = 0.7, H0 = 70 km/s/Mpc.
3.2 The Pan-STARRS1 Medium Deep Survey
The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS;
Kaiser et al. 2010) is a multi- filter imaging system designed for sky surveys on a
1.8 m telescope on the summit of Haleakala in Maui, Hawaii, with a 1.4 gigapixel
camera and a 7 deg2 field of view. The Pan-STARRS1 (PS1) telescope is operated
by the Institute for Astronomy (IfA) at the University of Hawaii and completed its
4.2 years of operation in the spring of 2014. ∼ 25% of the PS1 telescope time was
spent on the MDS, a deep, time domain survey of 10 circular fields distributed across
the sky, totaling ∼ 80 deg2, chosen for their overlap with extragalactic legacy survey
fields that have multi-wavelength corollary data. The PS1 MDS cadence typically
cycles through the gP1, rP1, iP1, and zP1 bands every three nights during the 6–8
months when the field is visible, observing in gP1 and rP1 on the same night and in
53
the yP1 band close to the full Moon (though y band observations were not used in
our study due to the poorer sampling). Nightly observations consist of eight 113
s (gP1 rP1) or 240 s (iP1 zP1 yP1) exposures (Tonry et al., 2012a); over the course
of MDS, each object is observed ∼ 300 times to a 5σ (i.e. where Σ = 0.217 mag
in Figure 3.3) limiting magnitude of ∼ 22.5 mag in gP1 rP1 iP1 and ∼ 22.0 mag in
zP1 in a single exposure. Individual exposures can be combined into nightly stacks
or full-survey-depth “deep” stacks, to reach much deeper limits of ∼ 23.5 mag and
∼ 25 mag, respectively (in the gP1, rP1, and iP1 bands).
The PS1 photometric calibration includes a combination of “absolute” cal-
ibration, which translates the number of photons detected to the physical unit of
magnitude, and “relative” calibration, which removes variations due to the telescope
system and atmosphere over the course of the survey. The PS1 absolute photomet-
ric calibration is accomplished by observing photometric standard stars from HST’s
Calspec catalog (Tonry et al., 2012b), as part of the Image Processing Pipeline (IPP;
Magnier 2006). The relative calibration is based on the algorithm of Padmanab-
han et al. (2008) which is known as “Ubercalibration” (Ubercal). The PS1 Ubercal
(Schlafly et al., 2012) uses multiple observations of the same non-intrinsically vari-
able sources on photometric nights and demands that the observed magnitude does
not change over time and thereby minimizes variations in the zero point. PS1 data
are further calibrated through “Relphot” (Magnier et al., 2013), which solves for
an additional zero point offset for each exposure, using the Ubercal solutions as a
starting point.
In L15, we employed a similar technique adapted from Bhatti et al. (2010)
54
(“ensemble photometry”; Honeycutt 1992) which implements the ENSEMBLE1 soft-
ware package in our attempt to achieve precision photometry with MDS. We con-
structed an “ensemble” of point sources near each target object in a ∼0.1 deg×∼0.1
deg field and ran the algorithm iteratively to obtain a least-squares solution that
locally reduces the scatter for all observations of each source over the course of the
survey. However, since the PS1 data products had already been Ubercaled, the
overall improvement in our control sample of stars in L15 or our re-analysis with
PV2 was not significant enough to justify this time-intensive procedure; thus we
do not apply the method of ensemble photometry in the analysis presented in this
paper.
3.3 Methods
3.3.1 Sample Selection
We first extracted from the PS1 Science Archive all sources from MD09 that
matched the following criteria: 1) they are point sources selected as deep stack
magpsf−mag 2Kron< 0 that have a good point spread function (PSF) quality factor
from the IPP (psfQF > 0.85), 2) they have stackPSFMag < 23 mag, 3) they have
at least five detections in each filter, and 4) masks were applied to exclude bad and
poor detections (Table 3.1). The query resulted in ∼ 40,000 point sources, for which
we get PSF magnitudes, each with an average of ∼ 300 detections in each of the
1http://spiff.rit.edu/ensemble/
2Since the Kron radius captures more flux from an extended source than the PSF profile, while
for a point source its Kron magnitude should be close to its PSF magnitude.
55
four filters.
For our color-selection of quasars and stars, we use a catalog of Kron magni-
tudes extracted from deep stack images in the gP1, rP1, iP1, zP1, yP1 bands from
PS1 MDS as well as in the uCFHT band from the Canada-France-Hawaii Tele-
scope (CFHT) (Heinis et al. 2016b, hereafter the PS1×CFHT catalog). We use
the PS1×CFHT catalog star/galaxy classification, which was determined using a
machine learning method of support vector machine (SVM) that was trained on an
HST/ACS sample of stars and galaxies and has a completeness of 88.5% for stars
iP1 < 21 mag, or 97.4% of all objects down to iP1 = 24.5 mag (Heinis et al., 2016a).
We then cross-match (using a 1′′ radius) our PS1 Science Archive point sources with
point sources in the PS1×CFHT catalog with rP1 < 23 mag (where the star/galaxy
separation is the most reliable), and within a 1.59 deg radius from the center of the
MD09 field (to avoid edge effects).
For the ∼ 15,000 cross-matched point sources, we then converted their CFHT
u3 and PS1 g r i z (Tonry et al., 2012b) band magnitude to the SDSS magnitude
system
uSDSS = (uCFHT − 0.241 gSDSS)/0.759
gSDSS = 0.014 + 0.162 (gP1 − rP1) + gP1
rSDSS = −0.001 + 0.011 (gP1 − rP1) + rP1
iSDSS = −0.004 + 0.020 (gP1 − rP1) + iP1
3http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/megapipe/docs/filt.html
56
zSDSS = 0.013− 0.050 (gP1 − rP1) + zP1 ,
and selected those that have the following quasar colors for their clean separation
from stars (Sesar et al., 2007):
uSDSS − gSDSS < 0.7
−0.2 < gSDSS − rSDSS < 0.5 .
We also select stars for our control sample, carefully avoiding the region oc-
cupied by RR Lyrae variables (uSDSS−gSDSS ∼ 1.15; Sesar et al. 2007). The color
diagrams of selected quasars and stars are shown in Figure 3.2. In order to obtain
our variability detection limit (Section 3.3.2), we plot their error vs. magnitude
relation for our star sample in Figure 3.3 and fit the binned relation in each filter
to a parabola:
σ(g) = 2.64372− 0.293112 g + 0.00818841 g2 (3.1)
σ(r) = 2.39328− 0.267030 r + 0.00749830 r2 (3.2)
σ(i) = 2.13028− 0.237271 i+ 0.00666299 i2 (3.3)
σ(z) = 2.77188− 0.309921 z + 0.00874017 z2 . (3.4)
57
Table 3.1: PS1 Quality Flags Used in the Query
Flag Name Flag Name
FAIL POOR
PAIR SATSTAR
BLEND BADPSF
DEFECT SATURATED
CR LIMIT EXT LIMIT
MOMENTS FAILURE SKY FAILURE
SKYVAR FAILURE BELOW MOMENTS SN
BLEND FIT SIZE SKIPPED
PEAK ON SPIKE PEAK ON GHOST
PEAK OFF CHIP
3.3.2 Variability Selection
To select intrinsic variables from our quasar sample, we perform a variability
selection in such way that systematic effects local to the field are minimized: we plot
the magnitude error in terms of standard deviation of the light curve as a function of
magnitude for each object within ∆R.A. = 0.5 deg and ∆decl. = 0.5 deg from each
color-selected quasar. Each of these “ensembles” contains ∼ 1000 point sources. We
calculate the median value for each magnitude bin, while avoiding the bins with less
than five stars, and interpolate linearly between the bin centers. The intrinsic vari-
58
ables have a significantly higher magnitude scatter than stars of the same brightness
and thus appear as outliers that deviate from the error vs. magnitude trend estab-
lished by the majority of objects. We iteratively remove variables from the linear
interpolation, and after three iterations, those that passed the final 2σ detection
threshold (have at least twice the magnitude error than the linear interpolation) are
tagged as variables (Figure 3.4). Our piece-wise interpolation method is adapted
from the variability selection procedure in ENSEMBLE and gives better results than
the parabolic fitting method in the previous version that was applied in our analysis
in L15.
Of the 670 quasars processed through this stage, we flag variables indepen-
dently in each filter, and further require a variability flag in at least two filters.
To compare our quasar sample with previous studies, we calculate their intrinsic
variability σint by putting in quadrature the standard deviation of the light curve
√
Σ and the photometric error σ: σ 2 2int = Σ − σ (Sesar et al., 2007) for Σ > σ
and σint = 0 otherwise, where σ is the magnitude-dependent photometric error from
Equation 3.1–3.4. We find the number fraction of quasars varying at the > σint level
qualitatively agrees with the results from SDSS Stripe 82 (S82) quasars in Sesar
et al. (2007) at σint > 0.06 mag, where ∼ 60% of quasars vary at or above that
level, compared to a control sample of stars for which the fraction is < 5% (Figure
3.5). The lower quasar variability fraction that we find for smaller variability am-
plitudes is likely due to our factor of ∼ 2 larger photometric errors compared to S82
(σ(g)>0.04 mag vs. σ(g)>0.018 mag).
59
Table 3.2: QLF Model Values Used in Equation 3.5 – 3.8
QLF Model α β M∗i k
∗
1 k2 log(Φ )
PLE (0.3 < z < 2.2) -1.16 -3.37 -22.85 1.241 -0.249 -5.96
LEDE (2.2 < z < 3.0) -1.29 -3.51 -26.57 -0.689 -0.809 -5.93
3.3.3 Selection Bias
To investigate the possible biases in our variability selection in a flux-limited
survey, we simulated ∼9000 quasars whose population is derived from the quasar
luminosity function (QLF). The QLF is defined as the number of quasars per co-
moving volume per unit magnitude and is described by a broken power law:
Φ∗
Φ(M, z) = , (3.5)
100.4(α+1)[M−M∗(z)] + 100.4(β+1)[M−M∗(z)]
where M∗ is the characteristic break absolute magnitude, and α and β are the
slopes of the QLF at the faint end and bright end, respectively. At lower redshifts
(0.3 < z < 2.2), the QLF is described by a pure luminosity evolution (PLE) model,
where the characteristic number density Φ∗ remains constant while M∗ evolves with
redshift quadratically (in i band):
M∗i (z) = M
∗
i (z = 0)− 2.5(k1z + k2z2) . (3.6)
At higher redshifts (z > 2.2), it is necessary to model the QLF as the results
of luminosity evolution and density evolution (LEDE), where Φ∗ and M∗ evolve
60
independently with redshift:
log Φ∗(z) = log Φ∗(z = 2.2)− c1(z − 2.2) (3.7)
M∗i (z) = M
∗
i (z = 2.2)− c2(z − 2.2) . (3.8)
We adopt the values for the constants (Table 3.2) given in Ross et al. 2013, which
expanded the redshift range in previous QLF studies (e.g. Richards et al. 2006b) to
0.3 < z < 3.5.
In each redshift bin of size 0.4 from z = 0.3−3.0, we integrate over the absolute
magnitude range (−14 < Mi < −30 mag, ∆Mi = 1 mag) and use the cosmology
calculator by Wright (2006) to calculate the co-moving volume of each shell from z.
From this quasar redshift distribution, we then populate each absolute magnitude
bin for each redshift according to the Φ(M, z). To convert absolute magnitudes at
different redshifts to the observed frame, it is necessary to apply the K correction:
m = Mz=2 + distance modulus + K(z), where we adopt the values for K(z) from
Richards et al. (2006b). We also converted the absolute magnitude at z = 2 in the
QLF to z = 0: Mz=0 = Mz=2 + 0.596, assuming a constant quasar spectral power
law index of α = −0.5 (Richards et al., 2006b).
With a distribution of quasars in redshifts and absolute magnitudes, we then
adopt the empirical relations from Heinis et al. (2016b) to calculate the expected
variability amplitude given the quasar’s (rest frame) absolute magnitude. Using
difference imaging of ∼1000 variability selected AGNs from PS1 MDS, Heinis et al.
(2016b) measured the anti-correlation between fractional flux variability and the
61
AGN luminosity in the gP1 rP1 iP1 zP1 bands, and here we convert their relations to
the magnitude space:
(∆f )
log ( ) = 0.17Mg + 3.40 (3.9)f g∆f
log ( = 0.18Mr + 3.64 (3.10)f r∆f )
log( ) = 0.20Mi + 4.02 (3.11)f i∆f
log = 0.20Mz + 4.07 , (3.12)
f z
where log(∆f/f) was measured from the maximum flux on the difference images
over the course of PS1 MDS, and the absolute magnitude M of the AGN was derived
from SED fitting (Heinis et al., 2016b).
We calculate the expected fractional flux variability log (∆f/f) from Equations
3.9 – 3.12 and convert to ∆m: ∆m = 2.5 log [1 + 10log (∆f/f)]. Finally, we scale
up to the expected variability amplitude in its rest frame wavelength using the
relation from Vanden Berk et al. (2004): v(λ) = 0.616 exp(−λrest/988Å)+0.164
(where v is variability amplitude measured in terms of “structure function,” also in
units of magnitude). To estimate the variability detection threshold for simulated
quasars, we calculate the expected photometric errors σ for a given quasar’s apparent
magnitudes m in the four PS1 filters from Equations 3.1 – 3.4. We adopt average
quasar colors of gP1 − rP1 = 0.14 mag, rP1 − iP1 = 0.15 mag, and iP1 − zP1 = 0.08
mag from our sample of PS1×CFHT quasars in Figure 3.2 as a proxy for our color
selection. Note that this color box is a valid assumption for quasars at z < 3
(Richards et al., 2002). Next, we apply the same magnitude cut as our PS1 quasar
62
sample (r < 23 mag). And finally, we assume that Σ = 0.023+0.27 ∆m, the average
empirical relation we found for our MD09 variable quasar sample, and using the same
variability selection criteria as described in Section 3.3.2, a quasar varying at Σ > 2σ
level in any and at least two filters is flagged in our simulation as a variable.
The redshift, absolute magnitude, and apparent magnitude distributions of
the “variable” simulated quasars are compared with the “visible” simulated quasars
(m < 23 mag) in Figure 3.7. Among the 924 “visible” quasars (m < 23 mag), and
assuming all obey the Heinis et al. (2016b) relation, 106 (or 11.5%) are selected
as “variables,” comparable to the observed variable quasar fraction of 15.5% we
find in the MD09 sample (see Table 3.3), while both the “visible” and “variable”
samples have similar apparent magnitude distributions with MD09 quasars (Figure
3.6). We find that in a sample of normal quasars or AGNs, our 2σ variability
selection is biased toward brighter quasars (m < 21 mag) at lower redshifts (z <
1), with a relatively flat distribution in luminosity (−21 mag > M > −27 mag).
Understanding this selection bias will be important in calculating the volume density
of SMBHB candidates in our final sample.
3.3.4 Selecting Periodic Quasar Candidates
We then began to search for potential periodic signatures using the Lomb-
Scargle (LS) periodogram, a Fourier analysis technique for unevenly spaced data
with noise (Horne & Baliunas, 1986; Lomb, 1976; Scargle, 1982). For N0 data
points in the time series spanning a total length of T in unit of days, we sampled
63
the periodogram at Ni independent frequencies (Horne & Baliunas, 1986) from 1/T
to N0/(2T ) (which would be the Nyquist frequency if data were evenly sampled);
the resolution of the periodogram is thus ∆f = (N0/2−1)/(TNi). Plotting power as
a function of f for all test frequencies, the dominant peak at frequency f or period
P = 1/f then signals a significant variation at that frequency or period.
When identifying periodic sources from their periodogram peaks, we also took
advantage of the redundancy of PS1 MDS monitoring in four filters (g, r, i, z), each
with a slightly different observing cadence due to weather and technical downtime,
to help filter out false detections by requiring periodogram peaks are coherently
detected (within a 10% error) in at least three filters. In each filter, the error of
√
the peak due to noise can further be calculated as ∆f = 3σr/(4 N0TA0) (Horne
& Baliunas, 1986; Kovacs, 1981) — where we calculate A0 as the best-fit sinusoidal
amplitude of the light curve phase-folded on the averaged period P̄ and σr as the
standard deviation of residuals after subtracting the signal from the light curve —
which gives us an error on the detected period: δP = δf/f 2. The total uncertainty
of the detected period is ca√lc∑ulated by put∑ting the theoretical and measured errors
in quadrature: (∆P )2 = ( δP 2i /4)
2 + (P − P̄ )2i /(N − 1) where i = 1...N is
the index of the coherent filter.
We calculate the S/N ratio of the sinusoidal fit as ξ = A20/(2σ
2
r), where σr
is the standard deviation of the model-fit residuals (Horne & Baliunas, 1986). We
mask any outliers that deviate from the mean by more than 4.5 σ and require
that candidates have ξ > 3 in at least one filter. Finally we require at least 1.5
cycles of variation, in accordance with similar studies (Charisi et al., 2016; Graham
64
et al., 2015b); this limit on the maximum allowed period is also justified since
spurious periods are oftentimes found on a timescale close to the total data length
(MacLeod et al., 2010). Our selection leaves three candidate periodic quasars in
MD09 (Table 3.3). Their periodograms in gP1 rP1 iP1 zP1 are shown in Figure 3.8,
and their complete PS1 light curves are presented in Figure 3.9. Note that in all
three periodograms, variability power increases with lower frequency, which is a
characteristic of red noise and a cause for concern in searching for periodicity. It is
thus important to understand the false-alarm rate due to red noise and further test
the sinusoidal model with extended baseline data (see discussion in Section 3.4.3).
3.4 Extended Baseline Photometry
Historically, there have been claims of (quasi-) periodicity on a number of
AGNs, but they failed to withstand re-analyses or follow-up observations (see e.g.
review by Vaughan & Uttley 2006). In the case of searching for light curve period-
icity with a Fourier method, a finite temporal baseline makes the observer highly
susceptible to “red noise leak” (see e.g. review by Press 1978 on red noise), where
low frequency variations are transferred to the sampled high frequencies for objects
with “red” power spectra of increasing power at low frequencies, such as AGNs and
X-ray binaries.
Fortunately, all three of our candidates have extended baseline photometry, ei-
ther from the archival database or our ongoing imaging campaign (Table 3.5), giving
us an advantage of testing the persistence of their periodic behavior by extending
65
the baseline to ∼ 2− 3 times the length of PS1 MDS with comparable photometric
precision.
Table 3.3: MD09 by the Numbers
Category Number
PS1 point sources 40,488
PS1×CFHT quasars 670
PS1×CFHT variable quasars 104
Variable quasars with coherent periodogram peaks 77
Candidates with ξ >3.0 in at least one filter 6
Candidates with Ncycle > 1.5 3
Table 3.4: PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates
PS1 Designation m (g,r,i,z) A0 (g,r,i,z)
PSO J333.0298+0.9687 (21.42, 20.94, 20.96, 20.95) (0.68, 0.51, 0.53, 0.39)
PSO J333.9832+1.0242 (18.97, 18.85, 18.79, 18.57) (0.11, 0.10, 0.09, 0.07)
PSO J334.2028+1.4075 (19.38, 19.28, 19.14, 18.94) (0.13, 0.11, 0.08, 0.06)
66
3.4.1 Follow-up Imaging
We have an on-going observing program at the Discovery Channel Telescope
(DCT) in Happy Jack, Arizona, to further monitor candidate PSO J334.2028+1.4075
with its Large Monolithic Imager (LMI) in gSDSS rSDSS iSDSS zSDSS filters (Table 3.5).
Here we present data from four observing runs on UT 2015 May 28, 2015 September
17 and 19, 2016 May 15, and 2016 July 10.
Each observation had five exposures (taken in a dither pattern) in each filter
on UT 2015 May 28 (5×50s), UT 2015 September 17 (gSDSS rSDSS iSDSS, 5×50s)
and 19 (zSDSS, 5×50s), UT 2016 May 15 (5×100s), and UT 2016 July 10 (5×100s).
The images were reduced with standard IRAF routines, astrometry-corrected with
SCAMP (Bertin, 2006), and co-added with Swarp (Bertin et al., 2002). For zSDSS band
images which are affected by fringe patterns, we constructed a master fringe map
from all z band images (with different telescope pointings) taken on one night using
the IDL function create fringes (Snodgrass & Carry, 2013). Combining with a
series of “control pairs” which mark the positions of adjacent bright and dark fringes
in the map, we then subtracted a scaled fringe map from the image using the IDL
function remove fringes (Snodgrass & Carry, 2013).
Using SExtractor (Bertin & Arnouts, 1996), we performed aperture photom-
etry on the co-added image, with the aperture radius used in each filter being the
typical full-width at half-maximum (FWHM) of the image and produced a catalog
of detections in the LMI’s 12’.3×12’.3 field of view. We then cross-matched the
catalog using a 1′′ radius with all the point sources (type = “star”) that are within
67
a 6′ radius from the target quasar with clean photometry (clean = 1) from the
SDSS catalog. We excluded very bright objects (m < 16 mag) to avoid saturated
detections on the LMI images, and on the cloudy night (UT 20150917) and on all
z band images, we also constrained the fitting to the locus where m < 21 mag.
We iteratively removed outliers that systematically deviate from the residual fit by
more than 0.2 mag (for m < 22 mag only) and fitted a linear function to the PSF
magnitude psfMag vs. the SExtractor instrumental magnitude mag aper relation
(we exclude the target quasar, which is variable, from fitting). Each residual plot
was also visually inspected to confirm a good fit.
The magnitude error was calculated by taking the standard deviation of the
residuals in the ∆m = 0.5 mag vicinity of the target quasar. Finally, we converted
the LMI photometric data to the PS1 magnitude for direct comparison with the
light curves from MDS.
In their quasar variability study, Morganson et al. (2014) pointed out there
are non-zero, albeit small, offsets for quasars after converting to PS1 magnitudes
from the SDSS system. They adopt a third-order polynomial (derived for main
sequence stars) to convert from SDSS to PS1 and add an additional average offset
to correct for the color-dependent difference between the magnitudes. Since the
Tonry et al. (2012b) filter transformations were also derived for stars (from synthetic
magnitudes of stellar SEDs), we have the following options to correct between the
SDSS and PS1 magnitudes in our light curves: (1) adopt the Tonry et al. (2012b)
relation without any additional offset or correction; (2) adopt the Morganson et al.
(2014) filter offsets; (3) calculate redshift-dependent synthetic magnitudes (and thus
68
offsets) from a composite quasar spectrum, where we redshift the composite quasar
spectrum from Vanden Berk et al. (2001) to the respective redshift of the candidate
quasar and convolve it with the SDSS and PS1 filter sensitivity curves (airmass
= 1.3 and 1.2, respectively) to calculate the synthetic magnitude in the respective
bandpass and therefore the mP1-mSDSS filter offset for each target.
Even though we eventually adopted our redshift-dependent synthetic quasar
correction as the most generic method, we note that the difference between the
conversion equations are small (∼ 0.01 mag), and, for quasars varying at the > 0.1
mag level, as our candidates are, the different choices of filter conversion are unlikely
to significantly change our conclusions with regard to the persistence of the variation.
3.4.2 Pre-PS1 Archival Photometry
We retrieved pre-PS1 archival SDSS S82 PSF light curves in gSDSS rSDSS iSDSS
zSDSS from SDSS-III DR12 (Alam et al., 2015). The S82 magnitudes were converted
to the PS1 system (Section 3.4.1) before being “stitched” to the PS1 light curves.
The resulting PS1+SDSS light curves are shown in Figure 3.9.
Candidate PSO J334.2028+1.4075 also has a Galaxy Evolution Explorer (GALEX)
Time Domain Survey (Gezari et al., 2013) light curve available in the NUV band
(λeff = 2316Å) ≈ 1 year before the start of PS1 MDS. We superimpose on the NUV
light curve a sinusoid of the same period and phase as in the PS1 light curves and
scale up the sinusoidal amplitude of the gP1 band (λeff = 4810Å) by the observed
exponential relation of variability amplitude as a function of (rest-frame) wavelength
69
for quasars from Vanden Berk et al. (2004) (see Section 3.3.3). The model is visu-
ally consistent with the larger variability amplitude of the NUV light curve (Figure
3.9). (The ordinate offset of the sinusoid is chosen such that it matches the mean
magnitude of the NUV light curve.)
In explaining the observed periodic variability of the CRTS candidate PG
1302−102, D’Orazio et al. (2015) derived the expected variability amplitude ratio
between the GALEX FUV and NUV and the CRTS V -band from spectral slopes,
a corollary of their Doppler boosting model. However, we have shown, in addition
to numerous previous studies (e.g. Gezari et al. 2013; Heinis et al. 2016b; MacLeod
et al. 2010; Vanden Berk et al. 2004) that a larger variability amplitude at shorter
wavelengths is commonly observed in quasars and AGNs, and Doppler boosting is
not unique in explaining the phenomenon.
3.4.3 Testing the Persistence of Periodicity with Extended Baseline
Photometry
We recalculated the S/N parameter ξ for the extended light curves by forcing
the same P̄ detected in PS1 only light curves (Table 3.7). Though all candidates
still have high significance values (ξ(g) ∼ 3) and the extended data have variation
amplitudes similar to their model sinusoidal amplitudes (Table 3.4), the extended
light curves do not agree with the extrapolation of their respective “PS1 only”
sinusoidal models and the periodic oscillations are not persistent.
The three candidates were selected by first calculating their significance with
70
respect to the null hypothesis of white noise (i.e. constant power over frequencies)
(Section 3.3.4). Previous systematic searches also assumed the null hypothesis of
damped random walk noise (DRW; Kelly et al. 2009) to calculate the false-alarm
rate and thus statistical significance of their selected binary candidates (Graham
et al. 2015b; L15; Charisi et al. 2016). (However, we note that the extended data
in Charisi et al. (2016) show that their DRW simulations underestimate the false-
alarm rate) The DRW null hypothesis is motivated by results from quasar light curve
analyses which demonstrate that the DRW model is a good description of normal
quasar variability (Kelly et al., 2009; MacLeod et al., 2010). The power spectrum of
the DRW process is P (f) = 2σ2τ 2/[1 + (2πτf)2] – where σ2 is the short-timescale
variance and τ corresponds to the characteristic timescale — it has a power law
slope of −2 at high frequencies (f > (2πτ)−1) and flattens to 0 at low frequencies,
analogous to the X-ray power spectrum of AGNs and X-ray binaries in the “low
hard” state (e.g. review by McHardy 2010, pp 203-32).
However, regardless of the model chosen for the power spectrum of quasar
variability to evaluate the significance of the period detection, we are fundamen-
tally limited by the several-year temporal baseline of current time domain surveys.
Vaughan et al. (2016) show that mock light curves generated from both DRW and
a broken power law power spectrum cannot be distinguished from a periodic signal
over ∼ 2 cycles, especially when adding photometric noise and down-sampling the
light curve to the actual observing cadence. Fortunately, a periodic candidate can
be favored or disfavored by observing it for a longer period of time (for a total of >5
cycles, ideally with better sampling and photometric precision), as Vaughan et al.
71
(2016) have suggested and as we have demonstrated in this paper.
3.5 Black Hole Mass Estimates and Inferred Binary Parameters
In order to measure the total black hole mass of the system (MBH) and derive
parameters under the binary model, we extracted the archival SDSS spectrum of
candidate PSO J333.9833+1.0242. We were also able to acquire spectroscopic ob-
servations of the other two candidates from DCT or the Gemini-South Telescope:
the spectrum of PSO J333.0298+0.9687 was obtained in 2015 Quarter 3 with DCT’s
DeVeny spectrograph with 300 g/mm grating and 1” slit for an exposure time of
1400 s. The data were reduced with the standard IRAF routines. A Gemini GMOS-S
long-slit spectrum was obtained for PSO J334.2028+1.4075 in the 2015A Semester
(Program ID: GS-2015A-Q-17. PI: T. Liu) with R400 grating and 0”.75 slit for a
total exposure time of 720 s. The Gemini spectrum was reduced with the Gemini
IRAF package.
In both spectra acquired, we clearly captured the broad Mg II line, allowing
us to use a combination of the broad line velocity and luminosity of the nearby
continuum to estimate black hole mass of the system from McLure & Dunlop (2004):
(M ) ( λL )0.62 (BH λ FWHM(MgII))2
log = 3.2 − − . (3.13)M 1044 ergs s 1 km s 1
In order to measure the Mg II broad line width for the black hole mass esti-
mate, we based our line fitting process on the prescription given by Vestergaard &
Wilkes (2001), in order to subtract the iron pseudo-continuum emission that contam-
72
inates in the vicinity of the Mg II 2800Å line: we first broadened the iron template
presented in Vestergaard & Wilkes (2001) by convolving with a series of Gaussian
kernels in an incremental step of 250 km/s, such that 1000 km/s < FWHMQSO <
12000 km/s. Then, in a fitting window of [2250,2650]Å where iron emission is con-
spicuous (Forster et al., 2001), we compare the FWHM=2000 km/s template with
the spectrum (from which a power law continuum was already subtracted) and it-
eratively determined a scale factor. We then compared the series of scaled and
broadened templates with the spectrum to determine the best-fit FWHMQSO.
After fitting for the iron emission and subtracting from the spectrum, a Gaus-
sian was then fitted to the Mg II broad line for the fitting range [2700,2900]Å,
whose FWHM is subsequently substituted into Equation 3.13. Any uncleaned sky
lines were excluded from the fitting process, and the final continuum and iron-fitted
spectrum was visually inspected to ensure the fitting is satisfactory (Figure 3.10).
Unfortunately, part of the Gemini spectrum (PSO J334.2028+1.4075) was affected
by a misbehaving amplifier over the wavelength range where iron pseudo-continuum
emission is strong. We had to mask the affected region and were not able to obtain
a good fit of the iron emission; instead we only fit a power law continuum to the
spectrum.
In the spectrum, we measured the continuum flux density at λ = 3000Å
(fλ(3000 Å)) from the continuum fitting and corrected for Galactic extinction us-
ing the dust map from Schlafly & Finkbeiner (2011) and the extinction curve from
Cardelli et al. (1989). With the redshift measured from the spectrum, we were then
able to translate flux into luminosity λLλ = λ4πD
2
Lfλ(1 + z) and calculate the total
73
black hole mass.
Though candidate PSO J333.9833+1.0242 has a measured black hole mass
from Shen et al. (2008), we applied the same line fitting and mass measurement
routine to its SDSS spectrum to obtain a self-consistent measurement. We estimate
log (MBH/M) = 9.5± 0.4, consistent with the black hole mass of log (MBH/M) ≈
9.8 from Shen et al. (2008).
For PSO J334.2028+1.4075, we measured black hole mass log (MBH/M) = 9.1
(with an error of 0.4 dex associated with the black hole mass estimator Mg II;
McLure & Jarvis 2002); it is lower than log (MBH/M) = 9.97± 0.5 quoted in L15
— which was estimated from C IV (Vestergaard & Peterson, 2006) and was not
measured from an electronic spectrum — but consistent with the previous black
hole mass, considering the large scatter between the Mg II and C IV-based methods
(log (MMgII/MCIV) = −0.06 dex, with a dispersion of 0.34 dex; Shen et al. 2008).
Having obtained the black hole masses, we convert the observed variability
period to the rest frame of the presumed binary: tvar = torb(1 + z) and directly
calculate binary orbital separation from Kepler’s third law for a circular binary
orbit:
t2 4π2orb = , (3.14)
a3 GM
where torb is the rest-frame binary orbital period, a is the orbital separation, and
M is the total mass of the system. The candidates’ continuum flux density, Mg II
line width, black hole mass, redshift, rest frame variability period trest, and inferred
74
binary orbital separation a are tabulated in Table 3.8.
Even though the three periodic quasar candidates from our search in PS1
MD09 have been disfavored by our extended baseline analysis, we compare their
observed period and inferred binary separation with search results from two other
time domain surveys in Figure 3.11: CRTS (Graham et al., 2015b) and PTF (Charisi
et al., 2016). Assuming the typical CRTS baseline of 9 years, all but seven of the 111
candidates claimed by Graham et al. (2015b) have variations of less than 3 cycles
(they require a minimum number of 1.5 cycles in their search), an insufficient data
length for a robust periodicity detection according to Vaughan et al. (2016) (see
our discussion in Section 3.4.3). As for the 50 candidates from PTF (Charisi et al.,
2016), although the majority (82%) of the candidates have more than 3 cycles of
variation, a large fraction of them have observed periods clustered around one year
(42% of their candidates have periods between 300− 400 days), a potential sign of
the aliasing effect of periodograms due to seasonal sampling (MacLeod et al., 2010),
even though their DRW simulations were down-sampled to the observing cadence
to account for this effect.
3.6 Summary and Conclusions
Periodic variability in quasars on the timescales of months to years has been
theoretically predicted as a signature of an SMBHB. Recent simulations show that
in triaxial galaxies (e.g. Vasiliev et al. 2015), the “final parsec problem” (e.g. review
by Milosavljević & Merritt 2003) is no longer an insurmountable problem that stalls
75
binary evolution at a > 1 pc separations and that binaries can evolve into the
GW-dominated regime (a < 10−3 pc) within a few Gyrs. A systematic search for
periodic quasars in a large synoptic survey thus provides a novel method to search
for SMBHBs in the final phase of their evolution and can potentially yield GW
sources in the nano-Hz frequency regime which is accessible to pulsar timing arrays
(PTAs) including NANOGrav (McLaughlin, 2013) and the Parkes Pulsar Timing
Array (Hobbs, 2013).
Our systematic search in the Pan-STARRS1 (PS1) MD09 field resulted in three
periodic quasar candidates, from an initial sample of ∼ 700 color-selected quasars,
that are apparently periodic over the PS1 baseline of 4 years. We further tested
the persistence of their periodicity with archival light curves from SDSS Stripe 82
and followed up with imaging with the DCT. Archival GALEX photometry also
confirms a larger amplitude of variation at shorter wavelengths, consistent with
previous quasar variability studies. These extended-baseline data with photometric
precision comparable to that of PS1 disfavor a simple sinusoidal model for the three
candidates over an extended baseline of ∼ 5 – 12 cycles. This corresponds to a
detection rate of < 1 out of 670 quasars (< 1.5×10−3), which is still compatible with
the theoretically predicted sub-parsec binary quasar fraction of < 10−3 out to z = 1
from cosmological SMBH merger simulations (Volonteri et al., 2009). The detection
rate per area (< 1 in 5 deg2) is also in agreement with the theoretical prediction
of 100 quasars per 1000 deg2 of search area (or 0.5 periodic quasars in 5 deg2)
from Haiman et al. (2009) for a flux-limited survey of quasars with mi < 22.5 mag.
Our ongoing search over all 10 PS1 MDS fields, together with using nightly stacked
76
images in the future which are ∼ 1 mag deeper, should increase our sensitivity to
true SMBHBs by a factor of 100, and yield tens of promising SMBHB candidates
for extended baseline monitoring and multiwavelength studies.
In comparison to other SMBHB searches, we note that there are two binary
candidates with double broad-line features from a sample of ∼ 17,500 SDSS quasars
(or a detection rate of ≈ 10−4) for z < 0.7 (Boroson & Lauer, 2009), consistent
with the predicted SMBHB rate of ∼10−4 (z < 0.7) from Volonteri et al. (2009).
We also note that Graham et al. (2015b) imply a similar detection rate to our
study of 68/∼75,000 ∼ 0.9×10−3 (for quasars z < 1), and Charisi et al. (2016) find
a detection rate of ≈ 1.4×10−3 for z < 3 (or 0.9×10−3 for the sub-sample that
remained significant after their re-analysis with extended data); however, see our
discussions in Section 3.4.3 and the relevant parts in Section 3.5 on the robustness
of those claimed candidates.
We have demonstrated the power of an extended baseline in testing periodic
quasar candidates in surveys whose temporal baselines (covering only 1.5−4 cycles)
are susceptible to false detections from red noise characteristic of normal quasar
variability. Fortunately, for most of the periodic quasar candidates discovered in
recent optical time domain surveys, continued monitoring over the next few years
can robustly test the persistence of the periodicity over a necessary number of cycles
(> 5) to filter out false alarms, and verify strong SMBHB candidates for direct
detection in GWs by PTAs.
77
Figure 3.1: In L15, we analyzed the periodic quasar candidate PSO
J334.2028+1.4075 based on its light curves in PV1 (upper panel), while its analysis
in this paper is based on its light curves from PV2 (lower panel). We note the extra
data from the last phase of PS1 MDS are included in PV2 (dashed box), while our
conclusions from our new analysis on its significance as a periodic quasar candidate
did not change. 4.5 σ outliers in g and z filters in both versions have been clipped.
The dashed lines are a sinusoid of P = 558 days (see text for details).
78
Figure 3.2: First three panels: cross-matched stars and quasars are selected by their
SDSS colors (converted from uCFHT gP1 rP1 iP1 zP1). In the upper left panel, the
regions occupied by quasars (blue) and stars (red) are represented by dashed boxes,
and the stellar region does not include RR Lyrae variables. Bottom right panel:
spatial map of all stars and quasars (red dots and blue crosses, respectively) in MD09
that have cross-matches in the PS1×CFHT catalog. The deep stack photometry
(Heinis et al., 2016b) was performed with each PS1 “sky cell” as the smallest unit
(each MD field is divided into 10×10 such rectangular regions), hence the rectangular
shape. The actual search area is smaller than the total area of MD09 field and is
about 5 deg2.
79
Figure 3.3: For our sample of stars, we plot the observed scatter (standard deviation
of the light curve) as a function of magnitude and fit the binned relation to a parabola
(Equations 3.1–3.4). We have masked outliers at the 4.5 σ level in our scatter plot,
and they are excluded from the binned scatter vs. magnitude relation. The relations
are similar in the four filters (∼ 0.03 mag at the bright end and rising to ∼ 0.15 mag
at the faint end), and the size of the error bars is less than 20 mmag, reflecting the
stability of the PS1 system over the course of the survey and zero-point variations
mainly due to the atmosphere (Schlafly et al., 2012).
80
Figure 3.4: Magnitude error (light curve standard deviation) vs. magnitude in the
gP1 filter for all the objects in one “ensemble” (crosses). The majority of the 2821
sources in the ensemble are non-variable, and their error vs magnitude relation (black
crosses) can be represented as a piece-wise linear function (blue solid line), and any
objects varying above the 2σ level (red dashed line) are excluded (red crosses). After
three iterations, the target quasar of this ensemble (marked with an additional red
square) is selected as a variable from this ensemble.
81
Figure 3.5: The number fraction of our MD09 quasars that vary more than σint
(bin size = 0.01 mag) in the gP1 band (solid histogram) decreases with increasing
intrinsic variability, and is in agreement with results from SDSS Stripe 82 quasars
(Sesar et al., 2007) for σint > 0.06 mag. The variability fraction of a control sample
of stars is shown in the dashed histogram. Plotted with a dotted line are number
fractions estimated from Sesar et al. (2007) which were derived from their sample
of spectroscopically confirmed quasars.
82
Figure 3.6: The gP1 band apparent magnitude distribution of our MD09 quasar
sample. The full PS1×CFHT quasar sample (dashed histogram; NQSO = 670) is
similar to the distribution derived from quasar luminosity function (Section 3.3.3),
and our variability selection (solid histogram; Nvar = 104).
83
Figure 3.7: To simulate the detectability of quasars by our selection criteria whose
variability amplitude obeys the empirical relations of normal AGN variability found
in Heinis et al. (2016b), we draw ∼9000 quasars from the redshift and absolute
magnitude distributions derived from the quasar luminosity function in Ross et al.
(2013) (full distribution not included). Among the “visible” quasars (m < 23 mag;
dashed histogram), our variability selection is biased toward lower luminosity (and
thus in general more variable) quasars at lower redshifts (solid histogram).
84
Figure 3.8: As part of our periodic quasar selection, we ran the Lomb-Scargle peri-
odogram on PS1 light curves and selected the sources that have a coherent period
detected in all four filters with high significant factors. In each set of panels, the
coherent peak was marked with a dashed line in each filter.
85
Figure 3.9: Left panels: PS1-only light curves in the gP1 rP1 iP1 zP1 filters. Light
curves are offset for clarity. The light curves in different filters are fitted to sinu-
soidal functions of the same period (P̄ ) and phase and of their respective best-fit
amplitudes (dashed lines). The PS1 photometric error bars are omitted for clarity;
instead, the typical photometric error is indicated in the “extended” panel. Middle
panels: In the extended baseline light curves are fitted to sinusoidal functions of
the same “PS1 only” period with the phase and amplitude being free parameters.
S82 light curves and LMI data (taken in SDSS filters) have been converted to the
PS1 photometric system. Light curves are also offset for clarity. For candidate PSO
J334.2028+1.4075, its Galaxy Evolution Explorer (GALEX) UV light curve is also
included, and we superimpose on it sinusoids of the scaled-up amplitude (purple)
(see text for details). Right panels: the extended baseline light curves are folded on
the same period as the first two panels and fitted to sinusoidal functions.
86
87
Table 3.5: Extended Baseline Photometry of Periodic Quasar Candidates
PS1 Designation Archival Follow-up UT Date of Follow-up Observations
PSO J333.0298+0.9687 SDSS – –
PSO J333.9832+1.0242 SDSS – –
PSO J334.2028+1.4075 GALEX DCT15Q2/Q3/16Q2/Q3 2015 May 28, 2015 September 17, 2015 September 19
2016 May 15, 2016 July 10
The column is empty if no follow-up imaging program is presented in the analysis.
88
Table 3.6: SDSS to PS1 Filter Offsets from Synthetic Quasar Magnitudes
PS1 Designation z gSDSS-rSDSS gP1-gSDSS rP1-rSDSS iP1-iSDSS zP1-zSDSS
PSO J333.0298+0.9687 1.284 0.2842 -0.0197 0.0087 -0.0005 0.0087
PSO J333.9832+1.0242 2.234 0.1033 0.0052 0.0090 -0.0108 -0.0032
PSO J334.2028+1.4075 2.070 0.1170 -0.0127 0.0094 -0.0060 -0.0002
z is the spectroscopic redshift of the candidate. The SDSS colors g-r and SDSS to PS1 magnitude
offsets mP1-mSDSS are calculated from synthetic magnitudes by convolving the composite quasar
spectrum with the respective filter. The offset is then added to each LMI magnitude.
89
Table 3.7: Detected Period and Significance Factors of Periodogram-Selected Candidates
PS1 Designation P̄ ±∆P (day) ξ (g,r,i,z) Ncycle ξ (g,r,i,z) Ncycle
PSO J333.0298+0.9687 428±12 (3.54, 2.82, 2.75, 1.09) 3.8 (2.96, 2.48, 2.42, 0.98) 12.1
PSO J333.9833+1.0242 465±11 (3.93, 2.58, 2.16, 1.34) 3.5 (2.84, 1.99, 1.61, 1.07) 11.1
PSO J334.2028+1.4075 558±19 (3.84, 2.74, 1.80, 0.91) 2.8 (3.33, 2.42, 1.66, 0.90) 4.6
The PS1+LMI light curve of PSO J334.2028+1.4075 totals a time span that corresponds to 4.6 cycles of sinusoidal
variation (P = 558 days). The extended light curves span 5.2 cycles if the GALEX UV data are also taken into
account.
90
Table 3.8: Spectroscopic Information and Inferred Binary Parameters of Periodic Quasar Candidates
PS1 Designation Spectroscopy fλ (3000 Å) FWHM(Mg II) log (MBH/M) z trest a a
[erg s−1cm−2Å−1] [km s−1] [day] [pc] [Rs]
PSO J333.0298+0.9687 DCT15Q3 2.4× 10−17 8851 9.2±0.4 1.284 244.3 0.004 28
PSO J333.9833+1.0242 SDSS 4.2× 10−17 6157 9.5±0.4 2.234 144.0 0.003 13
PSO J334.2028+1.4075 GS15A 1.9× 10−17 5492 9.1±0.4 2.070 181.8 0.003 28
Since the three candidates are no longer significant periodic candidates with the extended baseline, the binary model is
disfavored, and therefore we only show their binary parameters for comparison with other systematic searches.
60 10
14 Lyα
PSO J333.0298+0.9687 50 PSO J333.9833+1.0242 C III] PSO J334.2028+1.4075
12
z = 1.284 z = 2.234 8C III] z = 2.070±0.005
10 40 Mg II
Mg II Atmospheric
6 Atmospheric
A−band A−band8 C IV
30
6
C III] 4
20
4 Mg II
2
10
2
0 0 0
2000 2500 3000 3500 1500 2000 2500 3000 1800 2000 2200 2400 2600 2800 3000 3200
rest wavelength [Å] rest wavelength (Å) rest wavelength [Å]
6 10 6
4 4
5
2
2
0
0 0
−2
−2
−4 −5
−4
2400 2600 2800 3000 2400 2600 2800 3000 2200 2400 2600 2800 3000
rest wavelength (Å) rest wavelength (Å) rest wavelength (Å)
Figure 3.10: Top panels: MD09 candidate spectra from DCT/DeVeny, SDSS, and
Gemini-South. Prominent emission features and the telluric absorption line at
7600 − 7630 Å (observed wavelength) are marked with red tick marks. Bottom
panels: we show the procedure by which we measure the Mg II line width: we fit
and subtract the power law continuum (red dashed line) and iron emission in a
spectral window (orange lines) and fit the Mg II line to a Gaussian (blue line).
91
−17 2
flux (arbitrary units) flux [10  erg/s/cm /Å]
−17 2
flux (arbitrary units) flux (10  erg/s/cm /Å)
−17 2
flux [10  erg/s/cm /Å]
flux (arbitrary units)
Figure 3.11: The blue plus signs mark our three periodic quasar candidates in PS1
MD09, which we classify as false alarms after failing the test of persistence over
an extended baseline. The blue lines mark, from left to right, the length that
corresponds to 3 cycles of variation over the MDS baseline and the 4.2-year MDS
baseline. The 111 candidates from a systematic search in CRTS (Graham et al.,
2015b) are in red circles. Red lines represent 3 CRTS cycles, 1.5 CRTS cycles (the
minimum number of cycles required in their search), and the CRTS survey baseline
(assuming the typical length of 9 years). The 50 candidates from PTF (Charisi
et al., 2016) are in orange dots. Orange lines mark 3 PTF cycles and the PTF
baseline (assuming the typical length of 3.8 years).
92
Chapter 4: Did ASAS-SN Kill the Supermassive Black Hole Binary
Candidate PG 1302−102?
4.1 Introduction
Periodic light curve variability of quasars has been predicted as an observa-
tional signature of supermassive black hole binaries (SMBHBs) at sub-parsec sep-
arations, due to modulated mass accretion onto the binary (e.g. D’Orazio et al.
2013; Farris et al. 2014; Gold et al. 2014), or relativistic Doppler boosting of the
emission of the secondary black hole minidisk (D’Orazio et al., 2015). This predicted
signature has motivated several systematic searches for periodically varying quasars
in large time domain surveys, including Graham et al. (2015a,b) (hereafter G15),
Liu et al. (2015, 2016), and Charisi et al. (2016), and spurred a number of recent
claims of (quasi-)periodicity (and binarity) that were discovered serendipitously or
in previously well-known active galactic nuclei (AGN)1 (e.g. Dorn-Wallenstein et al.
2017; Kovačević et al. 2018). G15 reported a periodic quasar and SMBHB candidate
PG 1302−102 (hereafter PG 1302), with a light curve from the Catalina Real-Time
Transient Survey (CRTS) that can be fitted to a sinusoid of an observed period of
1However, some of these claims have already been challenged: for example, Barth & Stern
(2018) pointed out some issues that affect the Dorn-Wallenstein et al. (2017) analysis.
93
P = 1884 ± 88 days over the ∼ 9 year CRTS baseline. Its light curve including
the Lincoln Near-Earth Asteroid Research (LINEAR; Sesar et al. 2011) data, which
extends ∼ 0.5 cycles before the CRTS data, is consistent with the sinusoidal fit, and
archival photometry data from various telescopes are largely consistent with the
extrapolation of the sinusoid ∼ 10 years before LINEAR, although their sampling
is sporadic.
While there have been multi-wavelength analyses of PG 1302 in the UV
(D’Orazio et al., 2015), IR (Jun et al., 2015), and radio (Kun et al., 2015), which
can provide key complementary clues about the true nature of a variability-selected
SMBHB candidate, the periodicity of PG 1302 remains unconvincing due to the
small number of cycles (Ncycle ∼ 2 over a combined LINEAR+CRTS baseline).
Vaughan et al. (2016) have cautioned against claiming periodicity over such a small
number of cycles, as the stochastic variability (“red noise”) of normal quasars and
AGN (i.e., those that do not host SMBHBs) can easily mimic periodic variation.
Indeed, Vaughan et al. (2016) showed that aperiodic light curves simulated using
the Damped Random Walk model (DRW; Kelly et al. 2009) or a broken power-law
(BPL) power spectrum can also be fitted to few-cycle data after down-sampling
and adding photometric noise. Moreover, an extended baseline analysis using new
monitoring data disfavors the persistence of the periodic quasar candidates from the
Pan-STARRS1 Medium Deep Survey (PS1 MDS) MD09 field (Liu et al., 2016).
Three years after G15 and five since its last published CRTS data, we revisit
the periodicity of PG 1302 in this Letter, by adding the publicly available light
curve from the All-sky Automated Survey for Supernovae (ASAS-SN). We describe
94
the ASAS-SN light curve in Section 4.2 and the maximum likelihood method that
we use in the analysis in Section 4.3. In Section 4.4 we describe and simulate the
expectations in the case where a genuine periodicity is present, and then compare
those expectations with our reanalysis of PG 1302. We conclude in Section 4.5.
4.2 Extended Light Curve from ASAS-SN
The ASAS-SN survey (Kochanek et al., 2017; Shappee et al., 2014) is regularly
monitoring the variable sky down to V ∼ 17 mag using multiple telescopes hosted
by the Las Cumbres Observatory. We retrieved the ASAS-SN light curve of PG
1302 (J2000 R.A. = 196.3875, decl. = −10.5553) from 2012 February 15 to 2018
March 1 (MJD = 55, 972− 58, 178) from the Sky Patrol2. For calibration purposes,
we choose the length of the ASAS-SN light curve (≈ 2200 days) to overlap with
the CRTS light curve by ∼ 400 days. Due to the dense sampling and the large
photometric uncertainty of the ASAS-SN light curve, we have binned the light curve
using a width of ∼ 100 days (such that there are 20 bins over ∼ 2000 days with an
average of 46 measurements per bin) using the arithmetic mean, and the uncertainty
of each bin is given by the standard deviation of the measurements.
The CRTS (Drake et al., 2009) light curve of PG 1302 was retrieved from
the Second Data Release of the Catalina Sky Survey (CSS). While VCSS is based
largely on the Johnson V magnitude system used in ASAS-SN, there are some
differences. Instead of calculating a color-dependent correction to convert between
2https://asas-sn.osu.edu
95
the V magnitudes of the two surveys, we simply apply a constant offset to the ASAS-
SN light curve before it was “stitched” to the CRTS light curve: after binning
the CRTS data via the same method described above (15 bins each of width of
∼ 180 days), we calculate the difference between the (binned) CRTS and ASAS-SN
magnitudes in each of the two overlapping seasons, i.e., MJD ≈ 55, 900 − 56, 100
and MJD ≈ 56, 200 − 56, 500, and offset the ASAS-SN light curve by the average
difference (0.17 mag) in order to match to CRTS. The LINEAR light curve of PG
1302 has also been offset and binned in the same way. Although early-time data
from Garcia et al. (1999), Eggers et al. (2000) and Pojmanski (1997) are generally
consistent with the extrapolated sinusoidal fit to LINEAR+CRTS data, we do not
include them in our analysis due to their much sparser sampling and less reliable
photometry.
The full baseline in our analysis is therefore given by LINEAR+CRTS+ASASSN.
We present both the binned and unbinned light curves in Figure 4.1. Although the
ASAS-SN light curve does undulate, the periodic fluctuation detected in the CRTS
light curve is not consistent with the ASAS-SN data. In particular, the extended
ASAS-SN light curve fluctuation is clearly out of phase with the sinusoid fitted
to the LINEAR+CRTS light curve, and the full data set favors a longer apparent
period and larger amplitude.
96
4.3 Expectations for a True Periodic Signal
Because the LINEAR+CRTS+ASASSN light curve is inconsistent with a sinu-
soid of the best-fit period and phase from G15, we now analyze the combined data
by considering a possible periodic signal in the presence of red noise. The basic
picture is that fluctuations in the accretion disk can produce a red noise component
in the power spectrum, whereas the binary is expected to produce a periodic signal.
We adopt the maximum likelihood method introduced by Bond et al. (1998),
which has been applied in a number of previous studies, including Miller et al.
(2010), Zoghbi et al. (2013), and Foreman-Mackey et al. (2017). The observed light
curve is the combination of signal and noise x = s + n, or in terms of a correlation
matrix
Cx = Cs + Cn , (4.1)
where Cs = 〈sisj〉 and Cn = 〈ninj〉, and the indices i and j indicate elements of the
light curve, which has a total of N elements. The noise terms are assumed to be
Gaussian (which is usually true in optical astronomy); further assuming that they
are uncorrelated, Cn is simply a diagonal matrix with elements nini. Each element
of the signal matrix Cs can be expressed using the autocorrelation function
∫ +∞
〈sisj〉 = A(∆t) = P (f) cos(2πf∆t)df , (4.2)
−∞
where P (f) is the power spectral density (PSD) of the signal, and ∆t is the time lag
97
between si and sj. Having calculated the signal matrix Cx for a set of parameters
p, we can then construct a likelihood function L(p) under the model P (f):
L 1(p) = (2π)−N/2|C |−1/2exp(− xTC−1x x x) , (4.3)2
where |C | and C−1x x are the determinant and inverse of the matrix Cx, respectively,
and xT is the transpose of the time series x. To calculate the likelihood under
the DRW (Kelly et al., 2009), which has been successful in characterizing quasar
variability (e.g. Kelly et al. 2009; MacLeod et al. 2010), P (f) in Equation 4.2 would
take the following form:
2σ2τ 2
P (f) = , (4.4)
1 + (2πτf)2
where σ2 is the short-timescale variance, and τ is the characteristic timescale. To
search for a periodic component of frequency f0 in addition to DRW noise (here-
after “DRW+periodic”), we can introduce a delta function δ(f − f0), so that the
autocorrelation function in Equation 4.2 becomes
[ ∫ +∞ ]
A(∆t) = P (f) cos(2πf∆t)df + A0 cos(2πf0∆t) . (4.5)
−∞
where A0 is the amplitude of the periodic signal.
To test our implementation of the method, we simulated 10 light curves under
the DRW model using the Timmer & Koenig (1995) method, uniformly sampling σ
from 0.00224 mag day−1/2, which is the minimum value from the Kelly et al. (2009)
98
quasar sample, to 0.0187 mag day−1/2, which corresponds to the value at 3σGaussian
after fitting the Kelly et al. (2009) σDRW distribution to a Gaussian; the input τ
ranges from ≈ 30 to 970 days3. We then add sinusoidal functions with amplitudes
measured from the periodic candidates from PS1 MDS (T. Liu et al. 2018, in
preparation) so that A0 ≈ 0.1−0.3 mag. The input periods range from P ≈ 50−970
days; the maximum period corresponds to two-thirds of the length of the baseline,
which is the requirement in previous work including Graham et al. (2015b) and
Charisi et al. (2016). We then down-sample the light curve to the observing cadence
of PS1 MDS and add typical PS1 photometric noise. We then use a C implementation
of an affine-invariant Markov Chain Monte Carlo (MCMC) sampler (Goodman &
Weare, 2010) to sample the parameter space. Our implementation is successful
in recovering the input period: the best-fit periods generally follow a one-to-one
correlation with the input values. To select those by which the DRW+periodic model
is at least moderately preferred, we further impose the cut AICDRW+periodic−AICDRW
< −2, where the Akaike information criterion AIC = 2n − 2 lnL when there are
n parameters in the model. The AIC imposes a penalty on the more complex
model, and between two models the model with the lower AIC value is therefore
the preferred one. Those best-fit periods that meet this criterion correspond to > 3
cycles, and they follow a yet tighter correlation.
Next, we apply the method to a simulated DRW +periodic light curve to
demonstrate the expected decrease in the p-value (and therefore increase in signifi-
cance) if the periodic signal is real. We down-sampled the simulated light curve to
3All temporal parameters explored in this analysis are in the observed frame.
99
the sampling of the LINEAR+CRTS+ASAS-SN light curve and added photomet-
ric uncertainties that are typical of the three different surveys (Figure 4.2). The
light curve is then binned using the same bin sizes as Figure 4.1. The relative
amplitudes of the sinusoid and DRW noise are such that the significance level at
which the DRW+periodic model is preferred is comparable between the (binned)
LINEAR+CRTS-sampled light curve from the simulation and that from PG 1302.
The input period of P = 2012 days is chosen to be the same as the best-fit pe-
riod from our reanalysis of the LINEAR+CRTS+ASASSN light curve of PG 1302
(Section 4.4), and the phase of the simulated light curve also mimics that of PG
1302. As Table 4.1 shows, the method consistently recovered the input period in
the LINEAR+CRTS and LINEAR+CRTS+ASASSN-sampled light curves, and the
longer baseline produced a best-fit period that is closer to the true value with a
smaller uncertainty. Furthermore, the p-value (for a chi-squared distribution with
two degrees of freedom) has decreased significantly (by a factor of ∼ 100) when
the mock ASAS-SN data are included, even though they have a larger photometric
uncertainty than the simulated CRTS data.
4.4 Extended Baseline Analyses of PG 1302
We now apply this method to PG 1302, and the ranges of the sampled pa-
rameters are summarized in Table 4.2; in particular, the ranges of τ and P are
sampled from 200 days to 3000 days (recall that the putative period is P = 1884
days). As calculating the inverse and determinant of a large N ×N matrix is com-
100
putationally intensive (Equation 4.3; both are typically O(N3) operations4), where
N ∼ 1000 for the unbinned full-baseline light curve, we apply the method only to
the binned light curve, where N = 19 for LINEAR+CRTS and N = 35 for LIN-
EAR+CRTS+ASASSN. When we first applied the method to the CRTS-only and
LINEAR+CRTS light curves (Table 4.2), the DRW+periodic model is preferred over
the DRW-only model at the 98.4% and 99.9% levels, respectively. If PG 1302 were
the only quasar analyzed, this would be intriguing evidence for periodicity. However,
given that it was selected from an initial sample of ∼ 200, 000 CRTS quasars, its
periodicity can easily be produced by chance alone; to demonstrate strong evidence
for periodicity, the candidate should instead have a p-value < 5× 10−6.
As we showed in Section 4.3, for a genuinely periodic source we expect that ad-
ditional data should strengthen the evidence. However, the p-value of the DRW+periodic
model has increased from p = 1.39 × 10−3 on the LINEAR+CRTS baseline to
p = 4.70 × 10−3 after including ASAS-SN data (Table 4.2). The decrease in sig-
nificance after adding new data is inconsistent with our expectation when a true
periodic signal is present, which suggests that the periodic signal is not persistent.
The decrease in significance after including extended data was also seen for the
sources in Charisi et al. (2016). Their initial systematic search in the Palomar Tran-
sient Factory (PTF) identified 50 periodic quasar candidates from ∼ 35, 000 spec-
troscopically confirmed quasars. They analyzed those candidates using additional
data from CRTS and/or the intermediate Palomar Transient Factory (iPTF). Of
4However, we note that the algorithm celerite (Foreman-Mackey et al., 2017) is able to
compute Equation 4.3 at a cost of O(N) for some classes of PSD models, which include DRWs.
101
the 47 candidates that have additional data, all but two had significantly increased
p-values. Although the CRTS measurements have larger photometric uncertainties
than PTF or iPTF and are in a different filter, the increase in the p-value may still
be an indication that the additional data are inconsistent with the claimed period-
icity. A similar phenomenon from the statistical perspective is also seen in a large
sample of SDSS Stripe 82 quasars by Andrae et al. (2013): although a small number
of quasars are better described by the DRW+periodic model than the DRW-only
model, more quasars are preferred by DRW-only as the number of observations in-
creases. The failure of PG 1302 and the many periodic candidates from Charisi
et al. (2016) to demonstrate persistent periodicity therefore seems typical of the
stochastic variability that is ubiquitous in normal (single black hole) quasars and
AGN.
While quasar variability can be characterized by the DRW process, high-
frequency power-law slopes that deviate from DRW have been found in a number
of studies, including those using large samples from ground-based surveys (Caplar
et al., 2017; Koz lowski, 2016; Simm et al., 2016) and the ones using high-quality
Kepler AGN light curves (Aranzana et al., 2018; Edelson et al., 2014; Mushotzky
et al., 2011; Smith et al., 2018a). Because assuming the incorrect PSD form would
result in an overestimate of the significance of the periodic signal, we have also an-
alyzed PG 1302 under the more general BPL model and take the PSD in Equation
4.2 to be
102
Af−αlo
P (f) = , (4.6)
1 + (f/f −α +αbr) lo hi
where A is the normalization, fbr is the break frequency, and αlo and αhi are the low-
and high- frequency slopes, respectively. The parameter ranges sampled are listed in
Table 4.3; as also shown in the table, while the BPL+periodic model is moderately
preferred over the BPL only model and the best-fit period is consistent with that
in the DRW+periodic model, evidence for the periodic signal also becomes weaker
when ASAS-SN data are included.
4.5 Conclusions
PG 1302 has been reported as an SMBHB candidate, having shown apparent
periodic variation over∼ 2 cycles on a LINEAR+CRTS baseline of∼ 10 years (G15).
Its variability has been modeled as the relativistic Doppler boosting of the secondary
minidisk (D’Orazio et al., 2015), and it has an inferred binary separation of ∼ 0.01
pc. If verified, PG 1302 would be one of the most compact SMBHB candidates dis-
covered to date, and searches using similar techniques can potentially uncover more
candidates in the gravitational wave-emitting regime for multi-messenger studies
with the pulsar timing arrays.
In this Letter, we have included the recent ASAS-SN data for this source, which
has regular and dense sampling spanning ∼ 5 years since CRTS and thus extends the
total baseline to ∼ 15 years. We have also applied a maximum likelihood analysis
to search for a periodic component in addition to red noise, which is modeled as the
103
DRW process or a BPL PSD. While we find that DRW+periodic or BPL+ periodic
is the preferred model for the LINEAR+CRTS light curve, evidence for either model
becomes weaker after adding ASAS-SN data. As Doppler boost from a binary should
produce persistent periodicity, and more data should only strengthen the signal, our
reanalysis suggests that the variability of PG 1302 may be inconsistent with this
proposed model.
In this Letter, we have highlighted the importance of the long-term monitoring
of SMBHB candidates that have been selected for their periodicity; it is also neces-
sary to evaluate the significance of the periodic signal in the presence of stochastic
variability. Any robust periodic quasar and SMBHB candidate should be able to
withstand those two tests.
104
13.5
LINEAR −0.09 mag binned LINEAR −0.09 mag
CRTS binned CRTS
14.0
ASAS−SN −0.17 mag binned ASAS−SN −0.17 mag
14.5
15.0
15.5
16.0
3000 4000 5000 6000 7000 8000
MJD−50000
Figure 4.1: Combined light curve of PG 1302−102 from LINEAR (pink), CRTS
(black) and ASAS-SN (blue). LINEAR and ASAS-SN have been offset to match
CRTS (see the text). Adopting the best-fit period and its uncertainties from G15,
sinusoids with periods of P = 1884 days (cyan dashed line) and P = 1884 ± 88
days (cyan dotted lines) have been fitted to the LINEAR+CRTS light curve and ex-
trapolated to guide the eye. Additionally, we have superimposed a best-fit sinusoid
of the period P = 2012 days (black dashed line), the best-fit period of the LIN-
EAR+CRTS+ASASSN light curve that we determined under the DRW+periodic
model. The binned light curve is also shown (LINEAR: green; CRTS: orange; ASAS-
SN: magenta). Individual images of the light curve have not been checked for spu-
rious photometric detections.
105
V mag
−1.5
4
2
0
−1.0 simulated DRW+periodic −2
LINEAR+CRTS+ASASSN −4
−6
binned LINEAR+CRTS+ASASSN
−8
−10
−3 −2 −1
−0.5 −1log(frequency[day ])
0.0
0.5
0 2000 4000 6000 8000
normalized MJD
Figure 4.2: We generate a light curve under the DRW model and inject a periodic
function. The light curve is initially nightly sampled (black line). We then down-
sampled the perfect light curve and added typical photometric noise of the LINEAR,
CRTS, and ASAS-SN data (gray circles with error bars). The resampled light curve
is then binned (blue squares with error bars). The inset shows the periodogram of
the evenly sampled light curve without photometric noise. The DRW model that
generates the light curve is superimposed (red line), and the input period is indicated
with a red tick mark. We find that despite the significant photometric uncertainties
in the simulated ASAS-SN data, its addition to the analysis strongly improves the
evidence for periodicity when a periodic signal is actually present.
106
normalized magnitude
log(power)
107
Table 4.1: Maximum Likelihoods for the Simulated DRW+periodic Light Curve
DRW DRW+periodic DRW DRW+periodic
(LINEAR+CRTS) (LINEAR+CRTS) (LINEAR+CRTS+ASASSN) (LINEAR+CRTS+ASASSN)
lnLmax 22.98 30.04 46.42 58.17
p-value ... 8.59×10−4 ... 7.87×10−6
Pbestfit (day) ... 2060.75
+229.75 +59.42
−430.24 ... 2026.83−70.57
108
Table 4.2: Damped Random Walk Parameter Ranges Sampled by MCMC
DRW DRW+periodic DRW DRW+periodic
(LINEAR+CRTS) (LINEAR+CRTS) (LINEAR+CRTS+ASASSN) (LINEAR+CRTS+ASASSN)
σ (mag day−1/2min ) 0.00224 0.00224 0.00224 0.00224
σ (mag day−1/2max ) 0.0187 0.0187 0.0187 0.0187
τmin (day) 200 200 200 200
τmax (day) 3000 3000 3000 3000
lnA0 min ... −10 ... −10
lnA0 max ... 5 ... 5
Pmin (day) ... 200 ... 200
Pmax (day) ... 3000 ... 3000
lnLmax 20.49 27.07 33.19 38.55
p-value ... 1.39×10−3 ... 4.70×10−3
AIC −36.98 −46.15 −62.39 −69.10
P +434.87 +280.03bestfit (day) ... 1773.60−125.12 ... 2012.62−219.96
109
Table 4.3: Broken Power-law Parameter Ranges Sampled by MCMC
BPL BPL+periodic BPL BPL+periodic
(LINEAR+CRTS) (LINEAR+CRTS) (LINEAR+CRTS+ASASSN) (LINEAR+CRTS+ASASSN)
A 10−3 10−3 −3 −3min 10 10
A −2 −2 −2 −2max 10 10 10 10
f −1br min (day ) 0.00033 0.00033 0.00033 0.00033
fbr max (day
−1) 0.005 0.005 0.005 0.005
αlo min 0 0 0 0
αlo max 2 2 2 2
αhi min 2 2 2 2
αhi max 4 4 4 4
lnA0 min ... −10 ... −10
lnA0 max ... 5 ... 5
Pmin (day) ... 200 ... 200
Pmax (day) ... 3000 ... 3000
lnLmax 19.89 26.82 32.90 38.37
p-value ... 9.78×10−4 ... 4.21×10−3
AIC −31.79 −41.65 −57.80 −64.75
P (day) ... 1713.78+133.90 +114.11bestfit −116.09 ... 1988.13−845.88
Chapter 5: Supermassive Black Hole Binary Candidates from the
Pan-STARRS1 Medium Deep Survey
5.1 Introduction
Supermassive black hole binaries (SMBHBs) are expected as a result of galaxy
mergers occurring the Universe (e.g. Begelman et al. 1980). As the supermassive
black holes (SMBHs) in the centers of massive galaxies sink to the center of the
merged system via dynamical friction, the pair of active SMBHs on ∼ a few-kpc
scale can be observable as a dual active galactic nucleus (AGN; e.g. Comerford
et al. 2015). As its separation continues to shrink by ejecting stars in the “loss
cone”, the pair becomes a gravitationally-bound SMBHB at a sub-parsec separation.
While spatially resolving close-separation SMBHBs has been achieved with very long
baseline interferometry (e.g. Rodriguez et al. 2006), the direct imaging of SMBHBs
at farther distances is beyond the capabilities of current, or even future, telescopes.
An indirect method to search for SMBHBs is via spectroscopy, where the broad
emission line from one black hole is shifted due to its radial velocity (Eracleous
et al., 2012; Runnoe et al., 2017), or there is the presence of a double broad line
feature that is due to the broad line region associated with each black hole (Boroson
110
& Lauer, 2009).
Another indirect technique to search for SMBHBs is via their temporal vari-
ability signatures. (Magneto-) hydrodynamical simulations of an SMBHB system
(e.g. D’Orazio et al. 2013; Farris et al. 2014; Gold et al. 2014; MacFadyen & Milosavl-
jević 2008; Noble et al. 2012; Shi et al. 2012) show that the binary tidal torque clears
and maintains a low-gas-density cavity of a radius∼ 2a (where a is the binary separa-
tion) in the circumbinary disk, and material is ushered in through a pair of accretion
streams. This distinct accretion pattern of a binary-disk system causes the accretion
rate to strongly modulate on the order of the orbital frequency. Therefore, assuming
the accretion rate directly translates to electromagnetic luminosity, these SMBHBs
would manifest as AGNs or quasars that periodically vary on a timescale of months
to years. More recently, D’Orazio et al. (2015) also proposed a relativistic Doppler
boosting model: the SMBHB system is viewed at a high inclination angle, and the
emission dominated by the minidisk of the secondary black hole is Doppler-boosted
as the secondary travels at a relativistic speed along the line of sight. In addition to
optical variability, periodicity in the X-ray bands has also been predicted for SMB-
HBs at the inspiral stage, due to gas being flung towards and hitting the cavity wall
(Tang et al., 2018).
Observationally, there have been a number of systematic searches for periodi-
cally varying quasars in large optical time domain surveys: Graham et al. (2015a,b),
using the Catalina Real-time Transient Survey (CRTS); Charisi et al. (2016), using
the Palomar Transient Factory (PTF); and Liu et al. (2015, 2016) (hereafter L15 and
L16, respectively), using the Pan-STARRS1 Medium Deep Survey (PS1 MDS). Gra-
111
ham et al. (2015b) claimed 111 SMBHB candidates from a search among ∼ 200, 000
spectroscopically confirmed quasars in the CRTS footprint, and Charisi et al. (2016)
found 50 SMBHB candidates from a sample of ∼ 35, 000 spectroscopic quasars in
PTF, 33 of which remained significant after their re-analysis with extended data.
However, due to the stochastic nature of normal (i.e. single black hole) quasar
variability, the search for a periodic signal is highly susceptible to red noise (i.e.
increasing variability power on longer timescales) masquerading as periodicity over
a small number of cycles (Vaughan et al., 2016) and thus could produce a large
number of false positive detections in a systematic search. In fact, from the candi-
dates reported by Graham et al. (2015b) and Charisi et al. (2016) (assuming they
are all genuine SMBHBs with their claimed binary parameters), Sesana et al. (2017)
concluded that the expected stochastic gravitational wave background would exceed
the current pulsar timing array (PTA) upper limit by a factor of a few to an order
of magnitude. We addressed this issue of false positives due to red noise contamina-
tion in L16, where we tested the persistence of the periodic candidates with archival
SDSS Stripe 82 light curves or new monitoring data taken at the Discovery Channel
Telescope (DCT) since 2015, extending the total length of the baseline to Ncycle >
5. We find 3 periodic candidates from ∼ 1, 000 color-selected quasars in one PS1
Medium Deep field MD09, though none of them appear to be persistent over an
extended baseline. Further, we have reanalyzed the best candidate from the CRTS
SMBHB sample, PG 1302−102 (Graham et al., 2015a), by including new photomet-
ric data from the All-sky Automated Survey for Supernovae (ASAS-SN; Kochanek
et al. 2017; Shappee et al. 2014), and we have shown that the detected periodicity
112
does not persist as expected for a true SMBHB (Liu et al., 2018).
Here, we expand our analysis in L16 to all ten fields in the Pan-STARRS1
Medium Deep Survey (Section 5.2), where we have searched for periodic quasars
(Section 5.3) and identified 26 candidates (Section 5.4). In addition to monitoring
the candidates on a longer baseline with DCT and the Las Cumbres Observatory
network telescopes (Section 5.5), we also adopt a maximum likelihood method to
reanalyze the significance of their periodicity and take into account normal quasar
variability (Section 5.6). We also discuss the cumulative SMBHB rate from our
search and compare with previous work and theoretical predictions (Section 5.7.1).
A sample of SMBHB candidates from a systematic search are also exciting in the era
of multi-messenger astronomy, as they are possible gravitational wave sources for the
PTAs (Section 5.7.2). We also look ahead to the era of the Large Synoptic Survey
Telescope (LSST), which will transform the study of the variable sky over its 10 years
of operation and can potentially discover tens of thousands of SMBHBs (Section
5.7.3). Finally, we explore the multi-wavelength properties of the best SMBHB
candidate from our sample, properties which have been explored both theoretically
and observationally and are complementary to searching for SMBHBs in the time
domain (Section 5.7.4). We summarize our results in Section 5.8. We adopt the
following cosmological parameters throughout this Paper: Ωm = 0.3, Ωλ = 0.7, H0
= 70 km/s/Mpc.
113
5.2 The Pan-STARRS1 Medium Deep Survey
The Pan-STARRS1 survey (PS1; Kaiser et al. 2010) operated from 2009 –
2014 on the 1.8-meter PS1 telescope at the summit of Haleakala on Maui, Hawaii.
∼ 25% of the survey time was dedicated to the Medium Deep Survey (MDS), a
multi-filter, high-cadence time domain survey of 10 circular fields (Table 5.1), each
of which is ∼ 8 deg2 in size. MDS observed in the g 1P1 rP1 iP1 zP1 yP1 filters on
the AB photometric system (Tonry et al., 2012b) and can reach a 5σ magnitude
depth of 22.5 mag in gP1 rP1 iP1 and 22.0 mag in the zP1 filter in a single exposure
of 113 sec (gP1 rP1) or 240 sec (iP1 zP1). The data were processed by the PS1 Image
Processing Pipeline (IPP; Magnier 2006) and were made available to members of
the PS1 Science Consortium through the PS1 Science Archive.
Each nightly observation consisted of eight single exposures; though the sub-
exposures can be combined to produce “nightly stacks”, we have used the single-
exposure detections in this work, as well as in our previous work presented in L15
and L16. The telescope visited the field during the 6− 8 months that it was visible
and rotated through the gP1 rP1 iP1 zP1 filters every three nights (observations in
gP1 rP1 were carried out on the same night). Therefore, in the full MDS data set,
most objects were observed ∼ 400 times over the ∼ 4-year baseline (Figure 5.1).
1Although the yP1 filter was not used in our work.
114
Figure 5.1: The number of detections for a random sample of ∼ 27, 000 objects in
MDS.
5.3 Methods
5.3.1 Quasar and Variability Selection
We first extract sources from the catalog from the PS1 Science Archive that
meet the same criteria in L16 for MD09 data: (1) they are point sources (defined as
deep stack magpsf−magKron< 0) with good PSF quality factors (psfQF > 0.85); (2)
they have at least five detections; and (3) the same quality flags in L16 were applied
to exclude bad or poor detections. The query returns ∼ 30, 000 sources from each
MD field.
We then cross-match the PS1 sources with a catalog of deep stacked images in
the Canada–France–Hawaii Telescope (CFHT) u band and the PS1 g r i z y bands
(hereafter the PS1×CFHT catalog; Heinis et al. 2016b) using a one-arcsec radius.
115
Table 5.1: Medium Deep Field Centers
MD Field RA (J2000) Dec (J2000)
MD01 02:24:50 −04:35:00
MD02 03:32:24 −28:08:00
MD03 08:42:22 +44:19:00
MD04 10:00:00 +02:12:00
MD05 10:47:40 +58:05:00
MD06 12:20:00 +47:07:00
MD07 14:14:49 +53:05:00
MD08 16:11:09 +54:57:00
MD09 22:16:45 +00:17:00
MD10 23:29:15 −00:26:00
To extract point sources from the PS1×CFHT catalog, we used the star/galaxy clas-
sification in the catalog that has been trained on an HST/ACS sample of stars and
galaxies (Heinis et al., 2016a). We then convert the uCFHT gP1 rP1 band magnitudes
to the SDSS system, so that the quasar selection box in SDSS colors from Sesar
et al. (2007) can be directly applied (Figure 5.2). This results in ∼ 9,000 quasars in
∼ 50 deg2 of the total cross-matched sky area.
We then follow the method in L16 to select variable quasars: we construct an
ensemble of objects within ∆RA = 0.5 deg and ∆Dec = 0.5 deg from each color-
selected quasar. Then, in each filter, we compute the standard deviation σ of the
116
Figure 5.2: The CFHT u and PS1 g i z magnitudes were first converted to the
SDSS system, and quasars (blue dots) and stars (red dots) were selected by their
uSDSS − gSDSS and gSDSS − rSDSS colors (dashed lines represent the color boxes).
∼9,000 quasars were selected in the full PS1 MDS.
light curve for each object in the ensemble and iteratively exclude outliers by fitting
a piece-wise linear function to the σ–m relation: σ = σ(m). While most objects
in the ensemble are stars and follow a tight σ–m trend, intrinsic variable objects
such as quasars have significantly larger σ than stars of similar brightness and thus
would appear as outliers from the trend. We identify the quasars with standard
deviation > 2σ(m) in at least two filters as variables, and ∼1,400 out of the ∼9,000
color-selected quasars are identified as variable quasars.
117
5.3.2 Searching for Periodicity
To search for periodicity among the variable quasars, we compute the Lomb-
Scargle (LS) periodogram (Horne & Baliunas, 1986; Lomb, 1976; Scargle, 1982)
and take advantage of the multi-filter observations and their different sampling to
determine a coherent periodic signal by a “majority vote”. We also adopt the signal-
to-noise (ξ) definition in Horne & Baliunas (1986), where ξ = A2 20/(2σr) (A0 is the
sinusoidal amplitude, and σr is the standard deviation of the residual) and select
periodic candidates with high significance by requiring ξ > 3 in at least one filter.
Finally, we require that the periodic variation has at least 1.5 cycles over the 4-year
PS1 baseline. The search results in 26 periodic candidates from 10 MD fields, and in
Tables 5.2-5.5, we break down the number from each step of the selection pipeline
by the MD field. The mean gP1 magnitude distribution in Figure 5.3 shows our
detection limit of periodic candidates at gP1 ∼ 20 mag.
5.4 Results
5.4.1 Periodic Quasar Candidates from MDS
The light curves in the gP1 rP1 iP1 zP1 filters are then fitted to sinusoids of the
same period, while the amplitude and phase in each filter are left as free parameters.
In Table 5.7, we tabulat∑e the ξ calculated for each filter and the best-fit period P̄N
of each candidate: P̄ = (Pi)/N , where i = 1...N is index of the filter in which a
i
coherent period has been detected, and the uncertainty of P̄ is determined from the
118
Figure 5.3: The mean g band magnitude distribution of the 26 periodic candidates.
The magnitude limit of our candidate selection pipeline is ∼ 20 mag.
√∑ ∑
uncertainty in each filter: (∆P )2 = ( δP 2/N)2i + (P
2
i − P̄ ) /(N − 1), where
√
the δP in each filter is given by uncertainty in the frequency δω = 3πσ/(2 N0TA)
(Horne & Baliunas, 1986). We also list the number of cycles in Tables 5.7 and 5.8,
which is simply Ncycle = [max(MJD) − min(MJD)]/P̄. In Figure 5.4, we show the
distribution of the detected periods in the observed frame.
In Tables 5.9 and 5.10, we list the amplitude of the best-fit sinusoid and
the mean PS1 magnitude in each filter for each candidate. To compare with the
relation of variability amplitude vs. rest-frame wavelength with the previous study
of normal AGNs by Vanden Berk et al. (2004), we calculate the rest wavelength
for each PS1 filter at the redshift of each quasar (λeff(g) = 4810 Å, λeff(r) = 6170
Å, λeff√(i) = 7520 Å, λeff(z) = 8660 Å) and define an intrinsic variability amplitude
V = π(∆m)2/2− σ2, where here ∆m is the amplitude A0 obtained from our
119
Figure 5.4: The distribution of Pobs, the observed period determined by the peri-
odogram.
sinusoidal fit, and the magnitude-dependent observed scatter from stars is used
as a proxy for σ (see L16). The intrinsic variability amplitude V of our candidates
decreases with longer rest wavelength, which is consistent with the empirical relation
from Vanden Berk et al. (2004) and has no apparent deviation from regular AGNs
(Figure 5.5). We note, however, the exception of PSO J334.0298+0.9687, which
shows much larger variability amplitudes in all filters and an apparently steeper
amplitude-wavelength trend; a visual inspection of its light curves (Figure 5.25)
also shows a large variation (∼ 0.8 mag in the g band).
5.4.2 Spectroscopy and Black Hole Mass
We retrieved archival spectra of 16 candidates from the SDSS Science Archive
Server. The remaining candidates with no archival spectra were observed at the
120
Figure 5.5: We measure the variability amplitude of each candidate in each filter
after subtracting the measurement uncertainty in quadrature (gP1: blue circles;
rP1: green squares; iP1: orange diamonds; zP1: red pentagons). The amplitude V
decrease with longer rest wavelength, consistent with the exponential relation from
Vanden Berk et al. (2004) (black curve).
Gemini-South Telescope (PI: Liu) or the Discovery Channel Telescope (DCT). The
Gemini spectra were obtained with the R400 slit with GMOS, while the DCT spectra
were obtained with the DeVeny spectrograph with a 300 g/mm grating; we summa-
rize the details of the observations in Table 5.11. The Gemini/GMOS spectra were
reduced with the Gemini IRAF package, and the DCT/DeVeny data were reduced
with standard IRAF procedures.
Due to the variable weather conditions under which the spectra were taken,
a standard star may not accurately calibrate the science object’s flux. Therefore,
121
in addition to the standard procedures to reduce the spectroscopic data, we also
calibrate the object’s flux to its latest photometric measurement. We first convolve
the DeVeny spectrum with the SDSS r filter sensitivity curve to calculate a syn-
thetic magnitude r′SDSS; if it differs from the latest photometric measurement rSDSS
by more than the variability amplitude of the object (Table 5.9) — where rSDSS is
either observed with DCT/LMI (see Section 5.5) or, in the absence of new obser-
vations, obtained from the SDSS Science Archive Server — we then re-normalize
the spectrum to match its synthetic magnitude to rSDSS. The procedure is repeated
iteratively until |r′SDSS− rSDSS| < 0.05 mag. We note that this re-normalization pro-
cedure is unlikely to significantly bias our black hole mass estimates: a ∆m ∼ 0.8
mag intrinsic variability (which is on the order of the maximum variability ampli-
tude in our sample of candidates) translates to a factor of ∼ 2 difference in the
continuum luminosity (assuming z = 1), which in turn corresponds to a ∼ 0.2 dex
error on the black hole mass – much smaller than the systematic uncertainty of
Equation 5.1 or 5.2. Spectra of all candidates (including the renormalized DeVeny
spectra) are presented in Figures 5.27-5.30.
To measure a virial black hole mass from the spectrum, we first use the follow-
ing procedure to measure the broad line width of MgII: we fit a power-law continuum
in the range [2200, 2675] and [2925, 3090] Å and subtract it from the spectrum; we
then broaden and scale the iron emission template from Vestergaard & Wilkes (2001)
by fitting it to the range [2250, 2650] Å where iron emission is strong, which is then
subtracted from the spectrum. In those spectra where S/N is low, we do not fit
the iron emission to avoid over-fitting and subtracting. Next, we fit a single Gaus-
122
sian to the emission line in the range [2700, 2900] Å and measure a Full Width
at Half Maximum (FWHM). Although McLure & Dunlop (2004) fit two compo-
nents (broad and narrrow) to the MgII line and adopt the broad component in the
black hole mass estimate, we do not find a clear presence of a narrow component
in every spectrum and thus only fit a single Gaussian. Then, we measure the flux
density fλ at 3000 Å in the fitted continuum and convert to a continuum luminosity:
λLλ = λ4πD
2
Lfλ(1 + z). We also correct for Galactic extinction using the dust map
by Schlafly & Finkbeiner (2011) and the extinction curve of Cardelli et al. (1989).
Finally, we substitute the FWHM and λLλ into the following equation from McLure
& Dunlop (2004) to calculate the black hole mass:
M (BH FWHM(MgII))2(λL (3000Å))0.62λ
= 3.2 − − , (5.1)M km s 1 1044 ergs s 1
In a spectrum where CIV is the black hole mass estimator, we fit the continuum
in the range [1445, 1465] and [1700, 1705] Å, and after subtracting the continuum,
we adopt the procedure in Shen et al. (2008) and use a three-component fit to fully
characterize the CIV line profile: a narrow component with FWHM < 1200 km s−1,
a broad component with FWHM > 1200 km s−1, and a broader hump component.
We then measure the FWHM from the fitted profile. The corresponding continuum
luminosity is calculated from the mean flux density in the range [1340, 1360] Å, and
the black hole mass estimate is adopted from Vestergaard & Peterson (2006):
(M ) [(BH FWHM(CIV))2(λLλ(1350Å))0.53]
log = log + 6.66. (5.2)
M 1000 km s−1 1044 ergs s−1
123
Figure 5.6: The distribution of Prest ranges from ∼ 100−500 days. These rest-frame
periods are also tabulated in Tables 5.12 and 5.13.
The fitting procedures are demonstrated in the last two panels of Figure 5.30.
The measurements of z, fλ, FWHM, and MBH are listed in Tables 5.12 and 5.13,
and the distribution of rest-frame periods Prest = Pobs/(1 + z) is shown in Figure
5.6. We show the distribution of the measured black hole masses in Figure 5.7.
We also compare the redshifts of the candidates from PS1 MDS with those
from CRTS and PTF in Figure 5.8: our search with PS1 MDS appears to be more
sensitive to candidates at higher redshifts (〈z̄〉 ∼ 2) than CRTS or PTF (〈z̄〉 ∼ 1). In
fact, the redshifts of MDS candidates follow an opposite trend to those of the variable
quasars that our selection pipeline can detect (see L16), suggesting a selection bias
towards high redshifts.
From the black hole mass and redshift, we can infer a binary separation a via
Kepler’s law, assuming the variation is on the rest-frame orbital period timescale:
124
Figure 5.7: The black hole masses of the candidates from PS1 MDS are 108 - 1010 M
(solid histogram). This distribution is similar to that of the candidates from Graham
et al. (2015b) (dashed histogram) and Charisi et al. (2016) (dash dot histogram).
a3/t2 2orb = GM/4π , where torb = Pobs/(1 + z) is the rest-frame orbital period. In
Figure 5.9, we show the a-MBH parameter space occupied by the SMBHB candidates
with more than three cycles from Graham et al. (2015b), Charisi et al. (2016), and
this work. While those candidates with more than three cycles may already be in the
GW-dominated regime of orbital decay, LSST will explore much larger parameter
space than any of the three surveys.
5.5 Extended Baseline Photometry
As been pointed out by Vaughan et al. (2016), red noise can easily mimic a pe-
riodic variation over a small number of cycles (∼ 3), especially when the sampling
is sparse and uneven and the photometric uncertainty is large. Therefore, cur-
125
Figure 5.8: The redshift distribution of our candidates from PS1 MDS (solid his-
togram) shows that our selection is sensitive out to z ∼ 2, while the redshift distri-
butions of the periodic candidates from CRTS and PTF peak at z ∼ 1 (dashed and
dash dot histograms, respectively).
rent efforts to systematically search for periodically varying quasars (e.g. Charisi
et al. 2016; Graham et al. 2015b) are limited by the several-years-long baseline of
the survey, and it is essential to test the persistence of periodicity with long-term
monitoring. Our extended baseline analysis of the periodic candidates from MD09
presented in L16 and of PG 1302−102 in Liu et al. (2018) further demonstrated this
effect, and here we present our ongoing campaign to monitor candidates from the
full PS1 MDS.
New imaging data presented in this work include those taken with the Large
Monolithic Imager (LMI) in gSDSS rSDSS iSDSS zSDSS filters at DCT from 2015 May to
2017 November. In Tables 5.14 and 5.15, we list the Modified Julian Dates (MJDs)
126
on which the observations were carried out, as well as the filters that were used.
The images were reduced using standard IRAF routines and corrected for astrometry
with SCAMP (Bertin, 2006). For the zSDSS band images which are affected by fringe
patterns, we also subtract a scaled master fringe pattern created via create fringes
(Snodgrass & Carry, 2013) from all zSDSS band images taken on the same night and
remove the fringes using the routine remove fringes (Snodgrass & Carry, 2013).
We then co-add five sub-exposures in each filter (taken in a dither pattern
to avoid bad pixels) with SWARP (Bertin et al., 2002) before performing aperture
photometry using SExtractor (Bertin & Arnouts, 1996). Following the method
described in L16, we cross-match SExtractor detections with an SDSS catalog of
point sources from DR12 (Alam et al., 2015), resulting in ∼ 200 cross-matched pairs
in LMI’s 12′.3× 12′.3 field of view (FOV). We exclude bright, saturated detections
(m < 16 mag), faint objects (m > 22 mag), outliers, and the target itself (which
is variable) and obtain a linear transformation from the SExtractor instrumental
magnitude to an SDSS magnitude. We then apply the transformation to the target
and obtain a measurement of its magnitude on the SDSS photometric system. To
convert the SDSS magnitudes to the PS1 system and thus directly compare with
PS1 data, we adopt the same customized method in L16 that is suitable for quasar
colors. We first calculate synthetic PS1 and SDSS magnitudes by convolving the
(redshifted) composite quasar spectrum from Vanden Berk et al. (2001) with the
respective filter sensitivity. We then apply the PS1-SDSS magnitude offset to the
LMI measurements to obtain their magnitudes on the PS1 system. Finally, we
construct an extended light curve of the object by “stitching together” its PS1 light
127
curves and new LMI data. Those extended baseline light curves are presented in
Figures 5.20 - 5.26.
We have also included data from our monitoring program with the Las Cum-
bres Observatory (LCO), a global network of telescopes on both hemispheres. The
observations were carried out with the Spectral imager on the 2-meter class tele-
scopes at the Haleakala Observatory on Maui, Hawaii and the Siding Spring Observa-
tory in Australia between 2017 April and 2017 November (Project ID: NOAO2017AB-
013; PI: Liu) in the gSDSS and rSDSS filters. The MJDs of those observations are
also included in Tables 5.14 and 5.15. The LCO images have been reduced by the
BANZAI pipeline and are retrieved from the LCO Science Archive. Coadding of the
sub-exposures and photometry on the coadded image are then run on the same
custom-developed pipeline that we apply to LMI data. However, due to the smaller
FOV of the 2-m class LCO telescope (10′ × 10′) and shallower magnitude depth
(∼ 22 mag), we instead obtain ∼ 50 − 100 SDSS-cross-matched point sources on
each image, and we avoid faint detections and potential saturated detections by
excluding objects with m > 21 mag or m < 15 mag when performing photometry.
We tabulate the number of cycles for each candidate on its PS1-only and
extended baselines in Tables 5.7 and 5.8. While most PS1-only light curves only
have 1.5 − 3 cycles, the LMI and LCO monitoring data extended the baseline to
3 − 8 cycles, and, in the cases where archival SDSS Stripe 82 light curves are also
available, as long as ≈ 14 cycles (Figure 5.10). We have also re-calculated the S/N
factor ξ on the extended baseline, by assuming a pure sinusoid with the same period
as detected in the PS1-only light curve. Three candidates appears to be persistent,
128
having the same or improved ξ in all filters (Tables 5.7 and 5.8).
5.6 Reevaluating Periodic Candidates Using a Maximum Likelihood
Method
In Section 5.3, we have assumed the null hypothesis of white noise when search-
ing for a periodic signal using the LS periodogram. However, quasar variability is
known to be stochastic and has the characteristic of “red noise”, where variability
power increases on longer timescales. Here we reevaluate the significance of our
periodic candidates using a maximum likelihood method and investigate whether a
periodic component is justified if a red noise background is also present. We have
recently applied this method to the SMBHB candidate PG 1302−102 (reported by
Graham et al. 2015a) in Liu et al. (2018), and we refer the reader to that paper for
details. Guided by the expectations when a genuine periodic signal is present in ad-
dition to DRW noise, we down-select candidates whose maximum likelihoods under
the DRW model (LDRW) and the DRW model plus a periodic signal (LDRW+periodic)
meet the following criteria:
1. lnLDRW+periodic > lnLDRW for both PS1-only and extended light curves;
2. (lnLDRW+periodic − lnLDRW)Extended >
(lnLDRW+periodic − lnLDRW)PS1−only;
3. Pbestfit = PLS; where Pbestfit is determined from the maximum likelihood method,
and PLS is measured from the LS periodogram;
129
4. p < 1 ; where 9314 is the size of the initial sample of quasars, and the p-value
9314
is calculated for a χ2 distribution of two degrees of freedom.
While 12 candidates passed Criteria (1)–(3) (Tables 5.16 and 5.17), only one
(PSO J185.8689+46.9752, hereafter PSO J185) met all four criteria with such a
p-value that its periodicity is unlikely due to chance alone. The PS1 and extended
light curves of PSO J185 are presented in Figure 5.17, its SDSS spectrum is shown
in Figure 5.18, and its LS periodogram is shown in Figure 5.19. We note that
PSO J334.2028+1.4075 (L15) only passed Criteria (1), which is consistent with our
previous analysis in L16 that while it has been selected as a periodic candidate over
its PS1 baseline, its periodicity does not persist.
If a BPL red noise background is assumed instead, then none of the candidates
are statistically significant (Tables 5.18 and 5.19). However, we note that PSO
J185 is one of the eight candidates that passed Criteria (1)–(3), i.e. its p-value
has decreased when extended data are included, and therefore could still become a
significant candidate as our monitoring program continues.
130
Table 5.2: Medium Deep Fields by the Numbers (MD01-MD02)
Category Number
MD01
PS1 point sources 30,109
PS1×CFHT quasars 983
PS1×CFHT variable quasars 109
Variable quasars with coherent periodogram peaks 88
Candidates with ξ >3.0 in at least one filter 5
Candidates with Ncycle > 1.5 2
MD02
PS1 point sources 28,845
PS1×CFHT quasars 1,147
PS1×CFHT variable quasars 112
Variable quasars with coherent periodogram peaks 97
Candidates with ξ >3.0 in at least one filter 3
Candidates with Ncycle > 1.5 1
131
Table 5.3: Medium Deep Fields by the Numbers (MD03-MD04)
Category Number
MD03
PS1 point sources 31,350
PS1×CFHT quasars 942
PS1×CFHT variable quasars 202
Variable quasars with coherent periodogram peaks 134
Candidates with ξ >3.0 in at least one filter 7
Candidates with Ncycle > 1.5 4
MD04
PS1 point sources 32,661
PS1×CFHT quasars 1,030
PS1×CFHT variable quasars 200
Variable quasars with coherent periodogram peaks 158
Candidates with ξ >3.0 in at least one filter 11
Candidates with Ncycle > 1.5 6
132
Table 5.4: Medium Deep Fields by the Numbers (MD05-MD06)
Category Number
MD05
PS1 point sources 29,517
PS1×CFHT quasars 1,083
PS1×CFHT variable quasars 163
Variable quasars with coherent periodogram peaks 102
Candidates with ξ >3.0 in at least one filter 3
Candidates with Ncycle > 1.5 3
MD06
PS1 point sources 34,112
PS1×CFHT quasars 854
PS1×CFHT variable quasars 115
Variable quasars with coherent periodogram peaks 77
Candidates with ξ >3.0 in at least one filter 1
Candidates with Ncycle > 1.5 1
133
Table 5.5: Medium Deep Fields by the Numbers (MD07-08)
Category Number
MD07
PS1 point sources 29,031
PS1×CFHT quasars 815
PS1×CFHT variable quasars 120
Variable quasars with coherent periodogram peaks 84
Candidates with ξ >3.0 in at least one filter 3
Candidates with Ncycle > 1.5 2
MD08
PS1 point sources 38,194
PS1×CFHT quasars 1,013
PS1×CFHT variable quasars 138
Variable quasars with coherent periodogram peaks 98
Candidates with ξ >3.0 in at least one filter 5
Candidates with Ncycle > 1.5 3
134
Table 5.6: Medium Deep Fields by the Numbers (MD09-MD10)
Category Number
MD09
PS1 point sources 40,488
PS1×CFHT quasars 670
PS1×CFHT variable quasars 104
Variable quasars with coherent periodogram peaks 77
Candidates with ξ >3.0 in at least one filter 6
Candidates with Ncycle > 1.5 3
MD10
PS1 point sources 28,455
PS1×CFHT quasars 777
PS1×CFHT variable quasars 106
Variable quasars with coherent periodogram peaks 68
Candidates with ξ >3.0 in at least one filter 3
Candidates with Ncycle > 1.5 1
Full MDS
Ntot (quasar) 9,314
Ntot (variable) 1,369
Ntot (periodic) 26
135
136
Table 5.7: Period, Significance Factors, and Number of Cycles of Periodic Candidates (MD01-MD04)
PS1 Designation P̄ ±∆P (day) ξ (g r i z) Ncycle ξ (g r i z) Ncycle
(PS1 Only) (PS1 Only) (PS1 Only) (Extended) (Extended)
PSO J35.7068−4.2314 427±4 (3.6 3.1 3.6 2.2) 3.6 (3.4 3.1 3.4 2.2) 6.8
PSO J35.8704−4.0263 829±23 (3.5 3.8 3.5 2.0) 1.9 (3.4 3.7 3.5 2.0) 3.5
PSO J52.6172−27.6268 992±33 (5.0 5.6 4.9 2.9) 1.6 ... ...
PSO J129.4288+43.8234 313±5 (2.6 3.2 1.9 1.8) 4.9 (2.5 3.1 1.9 1.8) 8.8
PSO J130.9953+43.7685 717±18 (2.9 3.1 3.0 2.7) 2.2 (2.9 3.1 3.0 2.7) 3.6
PSO J131.1273+44.8582 843±31 (3.5 3.5 3.0 2.1) 1.8 (3.4 3.3 3.0 2.1) 3.0
PSO J131.7789+45.0939 697±18 (3.2 3.0 2.0 1.0) 2.2 (3.0 2.8 1.9 1.0) 4.2
PSO J148.8485+1.8124 816±5 (3.7 4.0 2.9 1.4) 1.9 (3.7 3.7 2.8 1.4) 3.2
PSO J149.4989+2.7827 960±8 (2.0 2.7 3.1 2.2) 1.6 (1.8 2.4 2.8 2.2) 2.7
PSO J149.2447+3.1393 810±8 (4.0 3.1 2.0 1.2) 1.9 (3.8 2.8 2.0 1.2) 3.2
PSO J149.9400+1.5090 417±5 (2.8 3.3 2.9 1.6) 3.7 (2.8 3.3 2.9 1.6) 6.3
PSO J149.6873+1.7192 820±5 (2.8 4.3 4.5 3.3) 1.9 (2.7 4.0 4.5 3.3) 3.6
PSO J150.9191+3.3880 741±9 (1.9 2.7 3.8 2.6) 2.1 (2.0 2.7 3.8 2.6) 3.6
137
Table 5.8: Period, Significance Factors, and Number of Cycles of Periodic Candidates (MD05-MD10)
PS1 Designation P̄ ±∆P (day) ξ (g r i z) Ncycle ξ (g r i z) Ncycle
(PS1 Only) (PS1 Only) (PS1 Only) (Extended) (Extended)
PSO J160.6037+56.9160 988±17 (3.0 2.0 1.6 1.2) 1.6 (2.9 1.9 1.6 1.2) 2.7
PSO J161.2980+57.4038 982±10 (3.7 3.2 2.9 1.6) 1.6 (3.5 3.0 2.8 1.5) 2.6
PSO J163.2331+58.8626 1000±13 (2.1 3.2 3.3 2.1) 1.5 (2.1 3.1 3.3 2.1) 2.7
PSO J185.8689+46.9752 958±19 (3.3 2.9 2.1 1.6) 1.6 (3.3 2.8 2.1 1.6) 3.3
PSO J213.9985+52.7527 727±22 (5.2 5.0 3.7 2.5) 2.2 (3.4 3.1 2.8 2.0) 4.0
PSO J214.9172+53.8166 1003±21 (4.0 4.4 4.0 2.4) 1.6 (3.5 3.8 3.7 2.3) 3.0
PSO J242.5040+55.4391 862±24 (2.9 3.5 2.8 2.0) 1.8 (2.8 3.3 2.7 2.0) 3.6
PSO J242.8039+54.0585 735±22 (3.2 2.8 2.1 1.4) 2.1 (3.1 2.8 2.0 1.4) 4.2
PSO J243.5676+54.9741 984±17 (3.2 2.6 1.2 0.4) 1.6 (3.0 2.5 1.2 0.4) 3.1
PSO J333.0298+0.9687 428±12 (3.5 2.8 2.8 1.1) 3.8 (1.9 1.9 2.0 0.9) 15.3
PSO J333.9833+1.0242 466±11 (3.9 2.6 2.2 1.3) 3.5 (2.4 1.8 1.6 1.1) 14.1
PSO J334.2028+1.4075 556±17 (3.8 2.7 1.8 0.9) 2.8 (3.2 2.3 1.6 0.9) 5.4
PSO J351.5679−1.6795 805±6 (1.9 2.0 3.2 2.5) 1.9 (1.8 1.9 3.1 2.5) 3.7
Table 5.9: PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates (MD01-05)
PS1 Designation m (g,r,i,z) A0 (g,r,i,z)
PSO J35.7068−4.23144 (19.69, 19.64, 19.69, 19.53) (0.23, 0.18, 0.23, 0.14)
PSO J35.8704−4.0263 (19.52, 19.46, 19.52, 19.23) (0.24, 0.21, 0.24, 0.13)
PSO J52.6172−27.6268 (20.37, 20.20, 20.14, 19.93) (0.34, 0.29, 0.22, 0.16)
PSO J129.4288+43.8234 (19.53, 19.37, 19.50, 19.48) (0.17, 0.16, 0.15, 0.15)
PSO J130.9953+43.7685 (19.88, 19.65, 19.81, 19.88) (0.21, 0.17, 0.18, 0.18)
PSO J131.1273+44.8582 (20.57, 20.42, 20.12, 19.87) (0.21, 0.20, 0.19, 0.15)
PSO J131.7789+45.0939 (20.62, 20.29, 20.29, 20.37) (0.22, 0.15, 0.14, 0.12)
PSO J148.8485+1.8124 (20.43, 20.17, 20.10, 19.88) (0.25, 0.20, 0.16, 0.11)
PSO J149.4989+2.7827 (20.34, 20.25, 20.24, 20.04) (0.19, 0.17, 0.15, 0.13)
PSO J149.2447+3.1393 (20.72, 20.72, 20.48, 20.45) (0.31, 0.27, 0.18, 0.17)
PSO J149.9400+1.5090 (20.17, 19.91, 20.00, 20.09) (0.18, 0.15, 0.15, 0.14)
PSO J149.6873+1.7192 (20.42, 20.12, 20.08, 20.08) (0.19, 0.15, 0.14, 0.14)
PSO J150.9191+3.3880 (19.63, 19.49, 19.39, 19.20) (0.20, 0.20, 0.21, 0.15)
PSO J160.6037+56.9160 (19.52, 19.33, 19.28, 19.33) (0.19, 0.13, 0.11, 0.11)
PSO J161.2980+57.4038 (20.45, 20.44, 20.18, 20.22) (0.28, 0.22, 0.15, 0.15)
PSO J163.2331+58.8626 (19.59, 19.48, 19.43, 19.19) (0.17, 0.15, 0.13, 0.09)
138
Table 5.10: PS1 Mean Magnitudes and Variability Amplitudes of Periodic Quasar
Candidates (MD06-10)
PS1 Designation m (g,r,i,z) A0 (g,r,i,z)
PSO J185.8689+46.9752 (20.54, 20.50, 20.23, 20.28) (0.30, 0.21, 0.17, 0.18)
PSO J213.9985+52.7527 (19.94, 20.13, 19.90, 19.89) (0.22, 0.22, 0.16, 0.16)
PSO J214.9172+53.8166 (20.53, 20.32, 20.39, 20.44) (0.28, 0.23, 0.21, 0.20)
PSO J242.5040+55.4391 (20.17, 20.17, 19.91, 19.95) (0.22, 0.24, 0.18, 0.18)
PSO J242.8039+54.05853 (19.72, 19.64, 19.87, 19.89) (0.27, 0.22, 0.18, 0.19)
PSO J243.5676+54.9741 (19.97, 19.64, 19.58, 19.61) (0.18, 0.15, 0.11, 0.07)
PSO J333.0298+0.9687 (21.42, 20.94, 20.96, 20.95) (0.68, 0.51, 0.53, 0.39)
PSO J333.9832+1.0242 (18.97, 18.85, 18.79, 18.57) (0.11, 0.10, 0.09, 0.07)
PSO J334.2028+1.4075 (19.38, 19.28, 19.14, 18.94) (0.13, 0.11, 0.08, 0.06)
PSO J351.5679−1.6795 (18.91, 18.56, 18.54, 18.67) (0.15, 0.12, 0.13, 0.12)
139
140
Table 5.11: Spectroscopic Follow-ups
PS1 Designation Telescope/Instrument Semester or Quarter Grating Slit width Exposure time
(arcsec) (sec)
PSO J52.6172−27.6268 Gemini/GMOS 16B (Gemini ID: GS-2016B-Q-50) R400 0.75 2×1000
PSO J149.2447+3.1393 Gemini/GMOS 15B (Gemini ID: GS-2015B-Q-42) R400 0.75 2×1000
PSO J149.6873+1.7192 DCT/DeVeny 17Q1 300 g/mm 1.5 2×2000
PSO J161.2980+57.4038 DCT/DeVeny 17Q1 300 g/mm 1.5 2×1700
PSO J163.2331+58.8626 DCT/DeVeny 17Q1 300 g/mm 1.5 2×1800
PSO J242.5040+55.4391 DCT/DeVeny 17Q1 300g/mm 1.5 2100
PSO J243.5676+54.9741 DCT/DeVeny 16Q3 300g/mm 1.5 2×900
PSO J333.0298+0.9687 DCT/DeVeny 15Q3 300g/mm 1.5 1400
PSO J334.2028+1.4075 Gemini/GMOS 15A (Gemini ID: GS-2015A-Q-17) R400 0.75 720
PSO J351.5679−1.6795 DCT/DeVeny 17Q2 300g/mm 1.5 1200
141
Table 5.12: Spectroscopic Measurements and Inferred Binary Parameters (MD01-04)
PS1 Designation Spectroscopy MBH Estimator fλ FWHM log (MBH) z Prest a a
[erg s−1cm−2Å−1] [km s−1] [M] [day] [pc] [Rs]
PSO J35.7068−4.23144 SDSS MgII 1.4×10−17 5185 8.7 1.564 167 0.002 47
PSO J35.8704−4.0263 SDSS MgII 3.3×10−17 3810 8.8 1.916 284 0.004 55
PSO J52.6172−27.6268 GS16B MgII 1.3× 10−17 7384 9.2 2.134 317 0.005 32
PSO J129.4288+43.8234 SDSS MgII 4.5×10−17 3744 8.3 0.959 160 0.002 80
PSO J130.9953+43.7685 SDSS MgII 4.1×10−17 3850 8.4 0.986 361 0.003 133
PSO J131.1273+44.8582 SDSS MgII 1.6×10−17 2450 8.3 2.011 280 0.002 126
PSO J131.7789+45.0939 SDSS MgII 2.0×10−17 6773 8.8 1.233 312 0.004 58
PSO J148.8485+1.8124 SDSS MgII 7×10−18 5402 8.9 2.378 242 0.003 45
PSO J149.4989+2.7827 SDSS CIV 3.4×10−17 5173 9.1 2.376 284 0.004 38
PSO J149.2447+3.1393 GS15B MgII 8.6× 10−17 1955 8.5 1.859 283 0.003 94
PSO J149.9400+1.5090 SDSS MgII 2.4×10−17 3715 8.3 1.106 198 0.002 102
PSO J149.6873+1.7192 DCT17Q1 MgII 1.3×10−17 (n) 5755 8.6 1.354 348 0.004 85
PSO J150.9191+3.3880 SDSS MgII 6.9×10−17 1995 7.7 0.719 431 0.002 426
142
Table 5.13: Spectroscopic Measurements and Inferred Binary Parameters (MD05-10)
PS1 Designation Spectroscopy MBH Estimator fλ FWHM log (MBH) z Prest a a
[erg s−1cm−2Å−1] [km s−1] [M] [day] [pc] [Rs]
PSO J160.6037+56.9160 SDSS MgII 3.7×10−17 3251 8.5 1.445 404 0.004 119
PSO J161.2980+57.4038 DCT17Q1 MgII 2.0×10−17(n) 3043 8.5 1.798 351 0.003 114
PSO J163.2331+58.8626 DCT17Q1 CIV 6.7×10−17(n) 5611 9.2 2.165 316 0.005 33
PSO J185.8689+46.9752 SDSS MgII 1.3×10−17 6070 8.9 1.681 357 0.004 59
PSO J213.9985+52.7527 SDSS MgII 1.5×10−17 4123 8.7 1.867 253 0.003 67
PSO J214.9172+53.8166 SDSS MgII 1.5×10−17 4907 8.4 1.169 462 0.004 142
PSO J242.5040+55.4391 DCT17Q1 MgII 1.9× 10−17 (n) 5547 8.9 1.780 310 0.004 53
PSO J242.8039+54.0585 SDSS MgII 3.6×10−17 6581 8.8 0.960 375 0.004 70
PSO J243.5676+54.9741 DCT16Q3 MgII 3.5× 10−17(n) 2041 8.0 1.268 434 0.002 280
PSO J333.0298+0.9687 DCT15Q3 Mg II 2.4× 10−17 8851 9.2 1.284 244 0.004 28
PSO J333.9833+1.0242 SDSS MgII 4.2×10−17 6157 9.5 2.234 144 0.003 13
PSO J334.2028+1.4075 GS15A Mg II 1.9× 10−17 5492 9.1 2.070 182 0.003 28
PSO J351.5679−1.6795 DCT17Q2 MgII 10.7× 10−17 4702 8.9 1.156 373 0.005 59
Figure 5.9: The blue solid curves represent the MBH – a detection limit if we require
3 cycles of periodic variation over a 4-year baseline (i.e. sources located below the
curves have more than 3 cycles over the baseline); curves from dark to light shades
correspond to z = 0, 1, 2. The purple solid curves correspond to 3 cycles over 10
years (i.e. LSST); from dark to light shades: z = 0, 1, 2. The black dashed curve
represents a binary separation of 500 RS, which we use as fiducial value within which
the binary is in the gravitational wave driven regime (i.e. candidates that are below
the black dashed curve). Only a small number of candidates from the three surveys
satisfy the 3-cycle requirement (assuming 4-year MDS and PTF baselines and a 9-
year CRTS baseline). LSST, however, will be able to detect periodic candidates in a
larger part of the a−MBH parameter space, many of them will be in the GW-driven
regime.
143
Figure 5.10: While most candidates have only ∼ 2 cycles in their PS1-only light
curves, we have extended the baseline to > 3 cycles with new imaging data from
DCT/LMI and LCO/Spectral and archival SDSS data.
144
145
Table 5.14: Extended Baseline Monitoring of Candidates (MD01-MD04)
PS1 Designation Telescope/Instrument MJD(s) of follow-up observation(s) Filters
PSO J35.7068−4.23144 DCT/LMI, LCO/Spectral 57641, 57682, 57940, 57993 g r i z
PSO J35.8704−4.0263 DCT/LMI, LCO/Spectral 57641, 57682, 57939 g r i z
PSO J52.6172−27.6268 ... ... ...
PSO J129.4288+43.8234 DCT/LMI, LCO/Spectral 57682, 57901 g r i z
PSO J130.9953+43.7685 DCT/LMI 57682, 57741 g r i z
PSO J131.1273+44.8582 DCT/LMI 57682 g r i z
PSO J131.7789+45.0939 DCT/LMI 57682, 57741, 58075 g r i z
PSO J148.8485+1.8124 DCT/LMI 57787 r i
PSO J149.4989+2.7827 DCT/LMI 57787 g r i
PSO J149.2447+3.1393 DCT/LMI 57369, 57788 g r i z
PSO J149.9400+1.5090 DCT/LMI 57788 g r i z
PSO J149.6873+1.7192 DCT/LMI 57788, 58075 g r
PSO J150.9191+3.3880 DCT/LMI, LCO/Spectral 57833, 57845 g r i z
146
Table 5.15: Extended Baseline Monitoring of Candidates (MD05-MD10)
PS1 Designation Telescope/Instrument MJD(s) of follow-up observation(s) Filters
PSO J160.6037+56.9160 DCT/LMI, LCO/Spectral 57741, 57852 g r i z
PSO J161.2980+57.4038 DCT/LMI 57741 g r i z
PSO J163.2331+58.8626 DCT/LMI, LCO/Spectral 57741, 57851 g r i z
PSO J185.8689+46.9752 DCT/LMI, LCO/Spectral 57833, 57858, 58075 g r i z
PSO J213.9985+52.7527 DCT/LMI 57833 g r i z
PSO J214.9172+53.8166 DCT/LMI, LCO/Spectral 57170, 57284, 57369, 57522, 57977 g r i z
PSO J242.5040+55.4391 DCT/LMI, LCO/Spectral 57522, 57578, 57642, 57851, 57977, 58012 g r i z
PSO J242.8039+54.0585 DCT/LMI, LCO/Spectral 57522, 57579, 57641, 57851, 57977, 58012 g r i z
PSO J243.5676+54.9741 DCT/LMI, LCO/Spectral 57522, 57579, 57851, 57901, 57976, 58012 g r i z
PSO J333.0298+0.9687 DCT/LMI, LCO/Spectral 57579, 57641, 57940, 57990, 58012, 58016 g r i z
PSO J333.9832+1.0242 DCT/LMI 57579, 57641, 57935, 58016 g r i z
PSO J334.2028+1.4075 DCT/LMI 57170, 57282, 57284, 57523, g r i z
57579, 57641, 57682, 57990, 57935, 58016
PSO J351.5679−1.6795 DCT/LMI, LCO/Spectral 57578, 57641, 57682, 57940, 57985, 58016 g r i z
147
Table 5.16: Re-analyses Using the Maximum Likelihood Method under the DRW Model (MD01-MD04)
PS1 Designation lnLDRW lnLDRW+p P [day] p-value lnLDRW lnLDRW+p P [day] p-value
(PS1 Only) (PS1 Only) (PS1 Only) (PS1 Only) (Ext.) (Ext.) (Ext.) (Ext.)
PSO J35.7068−4.2314 123.199 125.139 436.10+350.50−169.49 0.143 144.519 150.933 426.22+213.63−336.36 0.00163
PSO J35.8704−4.0263 131.743 137.596 953.43+9.75−240.24 0.00287 150.884 157.592 878.89+74.18−15.81 0.00122
PSO J52.6172−27.6268 88.442 95.024 953.20+15.16−34.83 0.00138 ... ... ... ...
PSO J129.4288+43.8234 146.559 146.866 412.99+192.45 0.735 144.646 146.082 230.54+218.78−277.54 −211.21 0.237
PSO J130.9953+43.7685 139.766 142.989 658.15+310.40−349.59 0.0398 159.584 166.840 721.75
+190.72
−29.27 7.06×10−4
PSO J131.1273+44.8582 121.928 125.289 845.20+113.61−566.38 0.0347 136.706 141.580 967.15
+4.20
−55.79 0.00764
PSO J131.7789+45.0939 113.227 122.755 694.46+259.30−40.69 7.27×10−5 140.334 145.664 801.40+49.08−140.911 0.00484
PSO J148.8485+1.8124 115.019 120.755 829.88+50.10−59.89 0.00322 134.159 138.174 885.46
+80.96
−279.03 0.0180
PSO J149.4989+2.7827 113.043 121.036 873.75+35.77−44.22 3.37×10−4 127.629 131.386 966.54+0−80.76 0.0233
PSO J149.2447+3.1393 101.544 108.689 929.64+41.27 7.88×10−4 99.573 103.032 916.98+51.74−28.72 −118.25 0.0314
PSO J149.9400+1.5090 131.136 131.962 332.23+340.78−329.21 0.437 149.546 151.097 408.28
+238.52
−381.47 0.212
PSO J149.6873+1.7192 105.086 111.844 812.11+93.47−56.52 0.00116 136.715 145.201 864.87
+32.49
−17.50 2.06×10−4
PSO J150.9191+3.3880 94.701 95.448 768.34+134.88−145.11 0.4737 109.721 111.391 708.60
+241.80
−158.19 0.188
148
Table 5.17: Re-analyses Using the Maximum Likelihood Method under the DRW Model (MD05-MD10)
PS1 Designation lnLDRW lnLDRW+p P [day] p-value lnLDRW lnLDRW+p P [day] p-value
(PS1 Only) (PS1 Only) (PS1 Only) (PS1 Only) (Ext.) (Ext.) (Ext.) (Ext.)
PSO J160.6037+56.9160 160.342 162.407 840.95+124.35−135.64 0.126 159.561 162.939 885.10
+78.26
−341.73 0.0341
PSO J161.2980+57.4038 120.002 127.036 850.52+79.49 −4−180.50 8.8×10 104.117 106.059 759.86+212.38−197.61 0.1434
PSO J163.2331+58.8626 157.293 158.086 948.16+22.23−177.76 0.452 158.228 159.103 294.14
+181.23
−268.76 0.416
PSO J185.8689+46.9752? 112.525 114.019 279.05+603.95−166.04 0.224 112.272 122.465 954.32
+17.41 −5
−22.58 3.74×10
PSO J213.9985+52.7527 146.229 154.131 768.64+159.95 −4 +20.72−410.04 3.7×10 142.168 148.269 686.37−29.27 0.00224
PSO J214.9172+53.8166 118.971 120.498 89.91+315.98−44.01 0.217 120.820 122.185 892.77
+70.37
−129.62 0.255
PSO J242.5040+55.4391 150.588 155.394 950.25+20.73−39.26 0.00818 172.481 174.968 881.85
+78.69
−91.30 0.0831
PSO J242.8039+54.0585 140.902 146.990 720.49+40.56−39.43 0.00226 170.788 176.619 725.43
+95.77
−34.22 0.00293
PSO J243.5676+54.9741 152.054 152.959 113.69+295.66 0.404 184.189 186.841 189.08+155.22−84.33 −164.77 0.0705
PSO J333.0298+0.9687 62.335 63.660 408.09+365.78−254.21 0.265 106.265 107.822 971.06
+0
−230.31 0.210
PSO J333.9832+1.0242 181.015 184.466 466.40+298.54−441.45 0.0317 242.674 248.747 442.39
+464.69
−15.30 0.00230
PSO J334.2028+1.4075 164.979 167.303 578.32+108.16−441.83 0.0978 214.130 215.850 695.27
+262.37
−127.62 0.179
PSO J351.5679−1.6795 117.924 124.945 814.33+140.21 8.92×10−4 81.783 88.879 826.42+67.75 −4−89.78 −22.24 8.28×10
5.7 Discussion
In Section 5.6, we have applied a more rigorous method and selected one
statistically-significant candidate which is preferred by the DRW+periodic model
and 12 candidates that are consistent with the expectations for a true periodic
signal (and eight if BPL noise is assumed instead). In this Section, we discuss
the astrophysical implications: how does our detection rate of SMBHB candidates
compare with previous work? Are these variability-selected SMBHB candidates from
PS1 MDS within the reach of current and future pulsar timing array experiments?
Given the capabilities of the Large Synoptic Survey Telescope, how many periodic
quasars can it detect?
5.7.1 The Detection Rate of SMBHBs
Boroson & Lauer (2009) (hereafter BL09) searched for SDSS quasars that
have multiple redshift systems, which could indicate the presence of a binary, and
there are two candidates that show such features from ∼ 17, 500 SDSS quasars at
z < 0.7. This rate (∼ 0.01%) is consistent with the results from Volonteri et al.
(2009) (hereafter VMD09), who predicted an upper limit of ∼ 0.1% per quasar for
z < 0.7 or ∼ 1% for z < 1.
To compare with the results from BL09, we calculate the cumulative number
of SMBHB candidates (N(< z)) per 1, 000 quasars included in the search in PS1
MDS. We also compare with previous work by Graham et al. (2015b) (hereafter
G15) and Charisi et al. (2016) (hereafter C16): G15 searched among ≈ 243, 000
149
Figure 5.11: Dotted histogram: V band magnitude distribution of the candidates
from CRTS (G15). Dashed histogram: the R magnitude distribution of the candi-
dates from PTF (C16). Solid histogram: the gPS1 magnitude distribution of candi-
dates from this work (also shown in Figure 5.3).
spectroscopically confirmed quasars and claimed 111 candidates; 50 candidates from
C16 were selected among ≈ 35, 000 spectroscopic quasars (33 after re-analysis with
extended data).
We first calibrate the completeness of G15 and C16 in detecting periodic
quasars relative to this work: our candidates have a magnitude cut-off at m ∼ 20
mag (Figure 5.3), which results in our sensitivity out to z ∼ 2 (Figure 5.8). As-
suming this work is complete out to z ∼ 2 and the candidates from G15 and C16
are relatively complete down to m ∼ 18 mag and m ∼ 19 mag, respectively (Figure
5.11), that translates to a redshift limit at z ∼ 1.0 for G15 and z ∼ 1.4 for C16. We
then calculate N(< z) out to the redshift that each search is complete (relative to
150
this work) by counting the total number of < z candidates that is in the respective
sample.
We note that from our down-selection criteria outlined in Section 5.6, six
candidates passed Criteria (1)-(3) under both DRW and BPL models. They are still
consistent with the expectation that new data should strengthen the evidence when
a genuine periodic signal is present (Liu et al., 2018), and with more extended data
from our ongoing monitoring program, the p-value is expected to further decrease if
the periodicity persists. Since at most six candidates can potentially be statistically
significant, we use them to calculate an upper limit, which corresponds to an SMBHB
rate of 0.6 per 1, 000 quasars. The cumulative rate inferred from C16 has a higher
value out to a lower redshift and is tension with our rate (Figure 5.12, left panel).
We also compare the number of SMBHB candidates per deg2 of sky area
searched (Figure 5.12, right panel). We performed our search in the cross-matched
area between the PS1 MDS and PS1×CFHT catalogs, which covers an area of ∼ 50
deg2. This corresponds to a rate of < 0.12 SMBHBs per deg2 (out to z ∼ 2). To
compare with the predicted observability of periodic sources, we adopt the fiducial
values in Haiman et al. (2009) (hereafter HKM09), for which we expect 20 sources
varying at 245 days in a 104 deg2 sky area for a survey magnitude depth of 22
mag (which is the magnitude limit of our candidate selection). Since most of our
candidates vary on the timescale of ∼ 800 days (Figure 5.4), we then apply the
scaling relation for a population of purely GW-driven SMBHBs to calculate the
expected number of periodic quasars, i.e. (800/245)8/3 × 0.002 = 0.05, which is
consistent with our upper limit (Figure 5.12, right panel). We note the caveat,
151
however, that the redshift of the sources from HKM09 is fixed at z = 2, while we
have measured a cumulative rate out to z = 2. We have also compared with the
predicted upper limit from VMD09 of ∼ 0.1 per deg2 out to z = 1, and our measured
rate is still consistent with this upper limit2.
5.7.2 Implications for Gravitational Wave Astrophysics in the PTA
Regime
SMBHBs are potential sources of gravitational waves (GW) in the nanohertz
frequencies (10−9 − 10−7 Hz). In the late 70s, Sazhin (1978) and Detweiler (1979)
realized that GW signals would appear in the times of arrival (TOAs) of pulsar sig-
nals as “unexplained” residuals, after subtracting a timing model from the observed
TOA. For a number of pulsars distributed across the sky, their timing residuals
should correlate with their angular separations on the sky (Hellings & Downs, 1983).
There are currently three pulsar timing array (PTA) experiments dedicated to us-
ing well-timed pulsars to detect this effect induced by GWs — the North American
Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013),
the European Pulsar Timing Array (EPTA; Kramer & Champion 2013), and the
Parkes Pulsar Timing Array (PPTA; Hobbs 2013) — and they combine to form the
International Pulsar Timing Array (IPTA; Hobbs et al. 2010). In the PTA regime,
SMBHBs could appear as individual continuous wave (CW) sources, or contribute
to a stochastic gravitational wave background (GWB). Those PTA experiments,
2However, we note that the VMD09 prediction is motivated by SMBHBs with broad emission
line features and not optical periodicity.
152
having operated for ∼ 10 years, have started to put meaningful constraints on the
strain amplitudes of CW sources (Arzoumanian et al., 2014; Babak et al., 2016; Zhu
et al., 2014) and the GWB (Arzoumanian et al., 2016; Lentati et al., 2015; Shannon
et al., 2015).
Recently, Perera et al. (2018) have placed upper limits on the SMBHB candi-
dates from G15, C16, and L163, and here we calculate the GW strain amplitudes
of the 26 candidates and compare them with this upper limit. For an SMBHB in a
circular orbit, the strain is given by:
2π2/3G5/3(f )2/3(M )5/3GW ch
h = , (5.3)
c4dL
where the observed gravitational wave frequency fGW = 2forb, Mch is the chirp mass:
M = (M M )3/5(M + M )−1/5ch 1 2 1 2 , and dL is the luminosity distance of the source.
As shown in Figure 5.13, the strain amplitudes of our SMBHB candidates from PS1
MDS are lower than this upper limit by ∼ 2− 4 orders of magnitude and therefore
cannot be detected as individual CW sources by the current PTAs.
5.7.3 Periodic Quasar Detections in the LSST Era
Expected to start its operation in ∼2022, the Large Synoptic Survey Telescope
(LSST; Ivezic et al. 2008) will be thousands of times more powerful than PS1 MDS,
thanks to its magnitude depth, photometric precision, and large survey area (Table
5.20). Here, we explore its capabilities to detect periodic quasars by using our results
3We note that PSO J334.2028+1.4075, which is included in both L15 and L16, has a lower BH
mass in L16 which was measured from its Gemini spectrum.
153
from PS1 MDS as a benchmark. The notation Ñ represents the number of quasars
from a simulated population, while N is the observed or expected number from a
survey.
Following the method in L16, we first simulate a population of quasars from
0.3 < z < 3.1 given the quasar luminosity function. We then apply the magni-
tude cut at m < 25 mag (Table 5.20); from Ñtot,LSST = 8, 996 simulated quasars,
Ñsel,LSST = 1, 700 quasars can be “visible” in the survey (Figure 5.14). Next, we
assign a variability amplitude to each quasar based on the same amplitude–absolute
magnitude relation from Heinis et al. (2016b). To determine the variability detec-
tion threshold, we adopt the expected photometric error as a function of magnitude
from Ivezic et al. (2008):
σ2 = σ2sys + (0.04− γ)100.4(m−m5) + γ100.8(m−m5) , (5.4)
From the same simulation performed for MDS in L16, the Ñsel,MDS = 924
quasar are selected from an initial sample of Ñtot,MDS = 8, 996. To estimate the
total number of quasars in the LSST footprint, we simply scale up the number of
quasars selected in MDS (Nsel,MDS = 9, 314) by the survey area A (Table 5.20):
Ñtot,MDS
Ntot,LSST = Nsel,MDS × ×
A(LSST)
= 3.63× 107 (5.5)
Ñsel,MDS A(MDS)
Since Ñvar,LSST = 1, 199 quasars are selected as variables from Ñtot,LSST =
8, 996 quasars from our simulation, the number of variable quasars that can be
detected by LSST is:
154
× Ñvar,LSSTNvar,LSST = N 6tot,LSST = 4.84× 10 . (5.6)
Ñtot,LSST
Assuming the same periodic candidate selection method (which selected 26
candidates from Nvar,MDS = 1, 369 quasars, out of which Ncand,MDS = 3 is persistent
over a baseline spanning 3 cycles, Table 5.7) is applied to LSST variable quasars,
the number of periodic candidates it could yield is:
Ncand,MDS
Ncand,LSST = Nvar,LSST × ≈ 10, 600 , (5.7)
Nvar,MDS
a factor of ∼ 50 more than the number of SMBHB candidates from G15, C16, and
this work combined.
We also estimated the part of the h − fGW parameter space that LSST can
explore. For simplicity, we adopt the mean redshift of the variable quasars from our
simulation of LSST quasars: zLSST = 1.6. The distribution of absolute magnitudes
to which LSST is sensitive (−17 < M < −27) corresponds to a black hole mass
range 6.5 < log(MBH(M)) < 9.6 (which is calculated using the quasar sample
from Peng et al. 2006). To determine the range of gravitational wave frequencies to
which it is sensitive, we assume at least 3 cycles over its 10-year baseline and adopt
a ∼ 20–day cadence for LSST.
To compare them with parameters that CRTS is sensitive to, we take MBH to
range from the the black hole mass that corresponds to 3σ below the mean to the
maximum mass of the G15 sample, i.e. 7.5 < log(MBH) < 10.2, and adopt its mean
redshift z = 1.0. We also require at least 3 cycles over its 9–year baseline and adopt
155
a 11–day cadence for CRTS. By comparison, the MDS parameter space is calculated
for the mean redshift of the simulated variable quasar sample (zMDS = 1.0), at least
3 cycles over a 4-year baseline, and a typical 3–day cadence. We assumed q = 1 in
calculating the chirp mass Mch from MBH, and those parameters of z, Mch, and fGW
are listed in Table 5.21.
We then compare those parameters with those of the first most likely CW
sources calculated for the IPTA and simulated SKA arrays from Rosado et al. (2015):
for each of the parameters log z, logMch, and log fGW, we approximate the proba-
bility density of the first likely sources as a Gaussian and draw ∼ 300−400 SMBHB
realizations from the distribution. The minimum and maximum values of each pa-
rameter are listed in Table 5.21, as well as the number of realizations of SMBHB
sources for each PTA array. As shown in Figure 5.15, the SMBHB candidates
from PS1 MDS with more than 3 cycles are unlikely to be detected as individual
sources by IPTA or the SKA arrays due to their higher GW frequencies; neither are
the current N > 3 cycles CRTS candidates4 likely to be detected by IPTA. How-
ever, candidates occupying the LSST parameter space might be detectable by the
SKA array: if LSST yields SMBHBs candidates that are nearby (z < 0.5), massive
(log(MBH) ∼ 9), and at long orbital periods (fGW ∼ 10 nHz), they are the best
contestants for the first individual CW sources detectable by the SKA.
4Only candidates with black hole masses in G15 are included.
156
5.7.4 Probing the SED Properties of SMBHBs
While a long baseline is essential to break false signals due to red noise and
help to verify the variability behavior of SMBHB candidates (Sections 5.5 and 5.6),
analyses of these systems based on optical variability alone may not suffice to identify
robust SMBHB candidates, and follow-up multi-wavelength studies are needed to
independently verify an SMBHB candidate.
For example, Roedig et al. (2014) have predicted a deficit in the spectrum
(“notch”), due to missing radiation from the cavity in the circumbinary disk. For
typical SMBHB parameters, the notch is expected in the UV band5. However, a
multi-wavelength study of the SMBHB candidate PSO J334.2028+1.4075 (hereafter
PSO J334) by Foord et al. (2017) shows that its spectral energy distribution (SED)
is consistent with that of a radio-quiet quasar6 and does not show evidence for the
predicted notch. Since PSO J334 has been disfavored as a periodic source by L16
and this work, we now explore any possible notch for the best candidate from our
PS1 MDS sample, PSO J185.
We query the archival photometry data from the AllWISE (Cutri & et al.,
2013), SDSS, and GALEX (Bianchi et al., 2011) catalogs; in the radio and X-
ray bands, where no detections are reported, we instead use their respective up-
per limits. We summarize the calculated rest-frame ν and νLν in Table 5.22.
We then calculate the temperature range kT0 − 15kT0 where the spectral notch
5However, see Farris et al. (2015), who do not predict the presence of a notch.
6With R∼17 (Foord et al., 2017), PSO J334 is technically classified as radio-loud.
157
is expected, where T0 is the characteristic temperature of the notch: T0 = 3.3 ×
104[ṁ(η/0.1)−1M−1 −3 1/48 (a/100Rg) ] K (we have assumed a radiative efficiency η =
0.1) and the largest deficit is expected at ∼ 4kT0 (Roedig et al., 2014). As we show
in Figure 5.16, the SED of PSO J185 is consistent with that of a radio quiet quasar
and does not show evidence for a spectral deficit. However, due to the low spectral
resolution, the presence of a notch cannot be ruled out either.
5.8 Summary and Conclusions
We have conducted a systematic search for periodically varying quasars in the
PS1 Medium Deep Survey, following our previous work in L16. Periodic variability
has been predicted as a signature of a supermassive black hole binary (SMBHB)
system as the mass accretion is modulated by the binary’s orbital motion; in an
SMBHB viewed at a high inclination angle, periodic variation can also produced by
relativistic Doppler boosting. SMBHBs at sub-pc separations should be a natural
product of galaxy mergers; however, compelling observational evidence for their
existence has been elusive. A systematic search for periodic quasars in the time
domain is therefore a novel approach to identify SMBHB candidates that are not
resolvable via direct imaging.
One challenge to the SMBHB candidates identified in systematic searches (e.g.
G15; C16; L16) is a robust detection of periodicity, since stochastic, normal quasar
variability can easily mimic periodic variation over a small number of cycles. We
monitor the variability of our periodic candidates using the Discovery Channel Tele-
158
scope and the Las Cumbres Observatory network and are able to extend the total
baseline of observations to 3-15 cycles. We have also searched for a periodic signal in
the presence of red noise, which is modeled by the Damped Random Walk (DRW)
process, or a broken power-law power spectrum. Only one candidate is considered
statistically significant from this analysis (assuming DRW noise), but six candidates
have increased significance when extended data are included (assuming either noise
model) and are therefore still consistent with the expectations for a true periodic
signal. This sample of SMBHB candidates down-selected from our more rigorous
method translates to a rate of < 0.6 per 1000 quasars, or < 0.12 per deg2, which
is consistent with theoretical predications but in possible tension with the rates
inferred by previous searches.
159
160
Table 5.18: Re-analyses Using the Maximum Likelihood Method under the BPL Model (MD01-MD04)
PS1 Designation lnLBPL lnLBPL+p P [day] p-value lnLBPL lnLBPL+p P [day] p-value
(PS1 Only) (PS1 Only) (PS1 Only) (PS1 Only) (Ext.) (Ext.) (Ext.) (Ext.)
PSO J35.7068−4.2314 122.084 123.290 440.65+272.31−197.69 0.299 146.222 149.598 438.39+253.81−386.18 0.0341
PSO J35.8704−4.0263 135.932 139.463 933.95+37.15−42.84 0.0292 159.068 159.456 865.56+96.79−363.20 0.678
PSO J52.6172−27.6268 92.067 94.015 920.65+42.86−327.13 0.142 ... ... ... ...
PSO J129.4288+43.8234 149.500 151.079 39.67+286.67−23.32 0.206 142.371 144.485 315.03
+283.01
−276.98 0.120
PSO J130.9953+43.7685 144.857 146.891 676.61+125.09−94.90 0.130 167.262 171.668 719.63
+73.29
−66.70 0.0122
PSO J131.1273+44.8582 122.882 126.346 854.51+63.68−76.31 0.0313 139.890 142.284 951.58
+11.24
−28.75 0.0912
PSO J131.7789+45.0939 116.029 122.908 683.18+35.24−564.75 0.00102 144.822 147.025 791.39
+161.90
−138.09 0.110
PSO J148.8485+1.8124 118.235 121.033 843.44+124.29−115.70 0.0609 136.484 137.714 884.32
+82.72
−527.27 0.292
PSO J149.4989+2.7827 117.984 115.298 837.46+36.91−33.08 ... 133.597 134.509 970.23
+0
−187.79 0.401
PSO J149.2447+3.1393 103.452 107.834 937.46+25.36−54.63 0.0125 100.806 102.672 911.45
+53.11
−256.88 0.154
PSO J149.9400+1.5090 132.670 135.961 139.77+221.62 0.0372 150.163 150.793 410.55+260.84−58.38 −399.15 0.532
PSO J149.6873+1.7192 106.677 111.590 809.69+124.54−225.45 0.00735 138.831 145.153 882.80
+73.75
−36.24 0.00179
PSO J150.9191+3.3880 91.667 96.013 750.01+215.48−184.51 0.0129 108.393 109.759 714.21
+155.27
−444.73 0.255
161
Table 5.19: Re-analyses Using the Maximum Likelihood Method under the BPL Model (MD05-MD10)
PS1 Designation lnLBPL lnLBPL+p P [day] p-value lnLBPL lnLBPL+p P [day] p-value
(PS1 Only) (PS1 Only) (PS1 Only) (PS1 Only) (Ext.) (Ext.) (Ext.) (Ext.)
PSO J160.6037+56.9160 163.097 165.788 940.31+27.71−362.28 0.0678 159.657 163.796 887.59
+81.60
−138.39 0.0159
PSO J161.2980+57.4038 121.914 127.821 889.98+80.16−59.83 0.00272 102.539 105.230 749.69
+212.07
−247.93 0.0678
PSO J163.2331+58.8626 157.344 157.899 114.96+329.92−110.07 0.574 158.363 159.913 114.03
+344.93
−65.06 0.212
PSO J185.8689+46.9752 115.700 117.084 276.42+282.49−17.50 0.250 115.590 121.326 948.27
+16.31
−23.68 0.00322
PSO J213.9985+52.7527 149.063 151.434 252.65+591.97 +178.37−88.02 0.0933 151.416 153.182 691.74−261.62 0.171
PSO J214.9172+53.8166 120.828 121.999 89.51+283.66−46.33 0.310 119.487 122.671 859.68
+40.56
−289.43 0.0414
PSO J242.5040+55.4391 156.506 157.818 923.49+44.85−125.14 0.269 178.187 178.543 860.67
+102.24
−267.75 0.700
PSO J242.8039+54.0585 141.217 145.791 11.11+289.78−0.21 0.0103 176.463 178.068 693.01
+79.97
−280.02 0.200
PSO J243.5676+54.9741 152.160 152.454 957.06+7.42 0.745 184.377 186.607 127.27+288.24−262.57 −111.75 0.107
PSO J333.0298+0.9687 63.575 63.509 417.03+408.59 ... 109.700 110.376 62.03+297.64−131.40 −52.35 0.508
PSO J333.9832+1.0242 179.692 185.789 472.80+314.55−15.44 0.00224 246.109 249.901 443.50
+178.16
−401.83 0.0225
PSO J334.2028+1.4075 165.162 166.541 573.64+106.54−153.45 0.251 215.223 216.272 131.73
+132.76
−127.23 0.350
PSO J351.5679−1.6795 124.080 134.519 832.49+81.66 2.92×10−5−28.33 113.713 120.584 827.56+70.45−19.54 0.00103
2.5 0.5
This Work
Grahahm et al. (2015b)
Charisi et al. (2016) 0.4 This Work (significant after reanalysis)2.0
Charisi et al. (2016) (significant after reanalysis) Volonteri et al. (2009)
This Work
0.3 Haiman et al. (2009) (z = 2)1.5 This Work (significant after reanalysis)
Boroson & Lauer (2009)
1.0 0.2
0.5 0.1
0.0 0.0
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
z z
Figure 5.12: Left panel: the cumulative number of SMBHB candidates per 1, 000
quasars from this work (blue circles), C16 (green diamonds), and G15 (red squares).
Both this work and C16 re-analyzed their samples by adding extended data; the
rates calculated from the candidates that are ruled significant after re-analysis are
indicated with filled symbols. The cumulative rate inferred from C16 is inconsistent
with the rate measured in this work. We also compare the rates from this work,
C16, and G15 with the spectroscopic search by BL09 (orange star). Right panel:
We compare the cumulative number of SMBHB candidates per deg2 of sky area
from this work with the predicted rates by VMD09 (purple triangles) and HKM09
(asterisk). The cumulative rate measured from this work is marginally consistent
with the prediction by VMD09 or HKM09 (see the text for details).
162
−3 −1
N(<z) [10  quasar ]
−2
N(<z) [deg ]
−12
10
−14
10
−16
10
−18
10
−20
10 −8 −7
10 10
fGW (Hz)
Figure 5.13: We use the current upper limit on CW (Perera et al., 2018) as a proxy
for detection sensitivity (dashed line). The gravitational wave strain amplitudes of
the candidates (bars) are lower by a few orders of magnitude than the upper limit.
The best candidate (PSO J185) is indicated as a red bar. All chirp masses are
calculated for the range q = 0.1− 1.0.
Table 5.20: Comparing PS1 MDS and LSST Capabilities
PS1 MDS LSST
Single visit 5σ magnitude depth in g band [mag] 22.5 25.0
Expected photometric error at g = 17 mag [mag] 0.02 0.005
Sky coverage [deg2] 50 20,000
163
h
Figure 5.14: From an initial sample of 8, 996 quasars drawn from the quasar lumi-
nosity function, 1, 700 can be detected by LSST (dashed histograms). Assuming
they follow the variability amplitude-absolute magnitude relation in Heinis et al.
(2016b), 1, 199 can be detected as variable quasars (solid histograms). They span a
range of redshifts (0.3 < z < 3.1) and absolute magnitudes (−17 < M < −27).
164
165
Table 5.21: Overview of test cases used in this study.
zmin zmax logMch,min logMch,max log fGW,min log fGW,max number of realizations
LSST 1.6 1.6 6.1 9.2 -7.7 -5.9 ...
CRTS 1.0 1.0 7.1 9.8 -7.7 -5.7 ...
PS1 MDS 1.0 1.0 7.3 9.5 -7.3 -5.1 ...
20-yr IPTA 0.003 1.0 8.5 10.5 -8.7 -8.0 383
6-yr SKA1 0.003 1.0 8.5 10.0 -8.2 -7.5 395
6-yr SKA2 0.003 1.0 8.0 10.0 -8.2 -7.5 322
−8
10
20−yr IPTA
−10
10 CRTS (> 3 cycles)6−yr SKA1
PS1 MDS (> 3 cycles)
6−yr SKA2
−12
10
CRTS
−14
10
−16
10
S1 MD
S
P
−18
10
−20
10 LSST
−22
10 −9 −8 −7 −6 −5
10 10 10 10 10
fGW (Hz)
Figure 5.15: We calculate the GW strain amplitudes of SMBHBs that are most
likely to be detected as individual sources by the 20-yr IPTA (blue contours), 6-yr
SKA1 (green contours), and 6-yr SKA2 (red contours), adopting their z, Mch, and
fGW values from Rosado et al. (2015). The contours have been smoothed, and the
darker (lighter) shade of contour color represents 68% (99.7%) of the realizations.
The abrupt cut-offs reflect the ranges of parameters applied (Table 5.21). We also
estimate the range of h and fGW that PS1 MDS, CRTS, and LSST are sensitive
to (solid, dashed, and dotted box, respectively). Though the SMBHB candidates
from PS1 MDS or CRTS (Ncycle > 3 only) are unlikely to be detected as individual
CW sources by IPTA (grey boxes and black circles, respectively), those from LSST
could be detected by the SKA. The best candidate from this work, PSO J185, is
also included in figure (orange circle). For simplicity, we have adopted q = 1 for all
calculations.
166
h
Table 5.22: The SED of PSO J185
Catalog Filter/Band ν (Hz) νLν (erg s
−1)
FIRST 1.4 GHz 3.75×109 (2.50×1041)
AllWISE W1 2.40×1014 9.38×1044
AllWISE W2 1.75×1014 7.96×1044
AllWISE W3 6.93×1013 (1.68×1045)
AllWISE W4 3.64×1013 (4.29×1045)
SDSS u 2.27×1015 3.07×1045
SDSS g 1.69×1015 2.58×1045
SDSS r 1.29×1015 1.80×1045
SDSS i 1.05×1015 2.05×1045
SDSS z 8.81×1014 1.67×1045
GALEX NUV 3.55×1015 2.02×1045
XMM-Newton 1.5 keV 9.72×1017 (4.89×1046)
167
47
46
45
44 46.0
45.8
43
45.6
42 45.4
45.2
41
45.0
15.0 15.2 15.4 15.6 15.8 16.0
40
10 12 14 16 18
log(ν) (Hz)
Figure 5.16: We construct the SED of PSO J185 using multi-band archival data (or
upper limits; blue squares). We then compare with the mean SEDs of radio loud
and radio quiet quasars (solid and dashed curves, respectively) from Elvis et al.
(1994). The mean SED has been normalized to PSO J185 at λrest = 2000Å. PSO
J185 is consistent with the SED of a radio quiet quasar, while showing no evidence
for a spectral notch (dashed orange lines mark the range kT0− 15kT0, and the solid
orange line marks 4kT0). The inset shows the SED of PSO J185 (and the mean
quasar SEDs) in the frequency range kT0 − 15kT0.
168
−1
log(νL ) (erg s )
ν
Figure 5.17: The PS1-only light curve of the best candidate from PS1 MDS, PSO
J185.8689+46.9752 (hereafter PSO J185) (left), and its extended light curve, which
consists of new data from DCT/LMI (squares) and LCO/Spectral (diamond). Si-
nusoids of periods determined from the periodogram are imposed to guide the eye
(dashed lines).
169
20
C IV
15 PSO J185.8689+46.9752
z = 1.681
10
C III] Mg II
5
0
1500 2000 2500 3000 3500
rest wavelength [Å]
Figure 5.18: The SDSS spectrum of PSO J185 (also shown in Figure 5.27). Its broad
MgII line (a blowup is shown in Figure 5.27, last but one panel) shows no apparent
spectral features that deviate from those of a normal quasar.
Figure 5.19: The Lomb-Scargle periodogram of PSO J185 (PS1-only). The dashed
line marks the peak at 958 days, which is consistent with the best-fit period from
the maximum likelihood analysis of the extended light curve (arrow). The best-fit
period of ≈ 280 days from the PS1-only light curve is also seen as a possible peak
in the periodogram (at f ∼ 0.003 day−1).
170
−17 2
flux [10  erg/s/cm /Å]
Figure 5.20: The PS1 and extended light curves of the remaining candidates from
PS1 MDS. Sinusoids of periods determined from the periodogram are imposed to
guide the eye (dashed lines). Different sources of archival or new monitoring data
are represented by different symbols: GALEX – dots, SDSS/S82 – stars, PS1/MDS
– circles, DCT/LMI – squares, LCO/Spectral – diamonds.
171
Figure 5.21: Figure 5.20 continued.
172
Figure 5.22: Figure 5.21 continued.
173
Figure 5.23: Figure 5.22 continued.
174
Figure 5.24: Figure 5.23 continued.
175
Figure 5.25: Figure 5.24 continued.
176
Figure 5.26: Figure 5.25 continued.
177
40 20
C IV
30 PSO J160.6037+56.9160 15 PSO J242.8039+54.0585
z = 1.445 z = 0.961
Mg II
20 10
C III] Mg II
H β
10 5
0 0
1500 2000 2500 3000 3500 4000 2000 2500 3000 3500 4000 4500 5000
rest wavelength [Å] rest wavelength [Å]
20 20
C III]
15 PSO J129.4288+43.8234 15 PSO J148.8485+1.8123
z = 0.959 Lyα z = 2.378
Mg II C IV
10 10
Si IV + O IV] C III] Mg II
5 5
0 0
2000 2500 3000 3500 4000 4500 5000 1500 2000 2500 3000
rest wavelength [Å] rest wavelength [Å]
30 30
25 C IV 25 C III]
PSO J213.9985+52.7527 PSO J130.9952+43.7685
z = 1.867 z = 0.986
20 20
15 Si IV + O IV] C III] Mg II 15 Mg II
10 10
5 5
0 0
1500 2000 2500 3000 3500 2000 2500 3000 3500 4000 4500 5000
rest wavelength [Å] rest wavelength [Å]
Figure 5.27: We retrieved archival SDSS spectra from the SDSS Science Archive and
obtained spectroscopic observations with Gemini/GMOS or DCT/DeVeny. Promi-
nent emission lines, including the black hole mass estimators C IV and Mg II, are
indicated with red tick marks. The last two panels show example procedures of our
spectral fitting of the continuum and the broad MgII or CIV line.
178
−17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
−17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
30 20
25 C IV
PSO J149.9400+1.5090 15 PSO J185.8689+46.9752
C III] z = 1.106 z = 1.68120
15 Mg II 10
C III] Mg II
10
5
5
0 0
2000 2500 3000 3500 4000 1500 2000 2500 3000 3500
rest wavelength [Å] rest wavelength [Å]
20
14
C IV
12
PSO J214.9172+53.8166 15 PSO J131.1273+44.8582
z = 1.169 z = 2.011
10
8
C III] Mg II 10 C III] Mg II
6
4 5
2
0 0
2000 2500 3000 3500 4000 4500 1500 2000 2500 3000
rest wavelength [Å] rest wavelength [Å]
30 10
C III]
25 C IV
8
PSO J35.7068−4.2314 PSO J131.7789+45.0939
z = 1.564 Mg II z = 1.233
20
6
15 C III]
Mg II 4
10
2
5
0 0
1500 2000 2500 3000 3500 4000 2000 2500 3000 3500 4000 4500
rest wavelength [Å] rest wavelength [Å]
30 30
25 C IV 25
PSO J35.8704−4.0264 PSO J149.4989+2.7826
z = 1.916 Ly C IV z = 2.37620 20 α
15 C III] 15
Mg II Si IV + O IV] C III] Mg II
10 10
5 5
0 0
1500 2000 2500 3000 3500 1500 2000 2500 3000
rest wavelength [Å] rest wavelength [Å]
Figure 5.28: Figure 5.27 continued.
179
−17 2 −17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
−17 2 −17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
30 70
60 Lyα25
PSO J150.9191+3.3880 PSO J333.9833+1.0242
z = 0.719 50Mg II z = 2.23420 C IV
40
15 H β Si IV + O IV]
30 C III] Mg II
10
20
5 10
0 0
2500 3000 3500 4000 4500 5000 1500 2000 2500 3000
rest wavelength [Å] rest wavelength [Å]
10 10
C III]
8 8
PSO J161.2980+57.4038 PSO J149.6873+1.7192
z = 1.799 Mg II Mg II z = 1.350
6 6
4 4
2 2
0 0
2000 2500 3000 2500 3000 3500
rest wavelength [Å] rest wavelength [Å]
25 50
Lyman α
20 40
C III] PSO J243.5676+54.9741 PSO J163.2331+58.8626
z = 1.268 z = 2.165
15 Mg II 30 Si IV + O IV]
C IV C III]
10 20
5 10
0 0
2000 2500 3000 3500 1200 1400 1600 1800 2000 2200 2400
rest wavelength [Å] rest wavelength [Å]
60
14
50 C III]
12
PSO J242.5040+55.4391 PSO J351.5679−1.6795
z = 1.780 Mg II z = 1.15610 40
C III]
8
Mg II 30
6
20
4
10
2
0 0
2000 2500 3000 2000 2500 3000 3500
rest wavelength [Å] rest wavelength [Å]
Figure 5.29: Figure 5.28 continued.
180
−17 2 −17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
−17 2 −17 2 −17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
10
14
PSO J333.0298+0.9687 C III] PSO J334.2028+1.4075
12
C III] z = 1.284
8 z = 2.070±0.005
10 Mg II
Mg II Atmospheric
6 Atmospheric
8 A−band
A−band
6 4
4
2
2
0 0
2000 2500 3000 3500 1800 2000 2200 2400 2600 2800 3000 3200
rest wavelength [Å] rest wavelength [Å]
100 20
80 C III] PSO J149.2447+3.1392 15 C IV PSO J52.6172−27.6268
z = 1.859 z = 2.132
60 Mg II
C III]
10
Mg II
40
5
20
0 0
2000 2500 3000 1600 1800 2000 2200 2400 2600 2800 3000
rest wavelength [Å] rest wavelength [Å]
8
5
6
4
4
3
2
2
0
−2 1
0
2200 2400 2600 2800 3000 1450 1500 1550 1600 1650 1700
rest wavelength (Å) rest wavelength (Å)
Figure 5.30: Figure 5.29 continued.
181
−17 2
flux [10  erg/s/cm /Å] −17 2
flux (arbitrary units) flux [10  erg/s/cm /Å]
−17 2 −17 2
flux [10  erg/s/cm /Å] flux [10  erg/s/cm /Å]
flux (arbitrary units)
Chapter 6: Conclusions and Future Work
In this thesis, I have presented a systematic search for periodically varying
quasars, which are theoretically-predicted manifestations of SMBHBs at sub-pc sep-
arations. Using a custom-developed pipeline, I have mined the PS1 MDS dataset
and selected 26 candidates from ten MD fields. These candidates have at least 1.5
cycles of variation over the PS1 baseline, and their observed periods are between
∼ 400 days and ∼ 1000 days. Their PS1 magnitudes are between gP1 = 19 and
21, at least 1 mag deeper than previous systematic searches. By following up these
sources spectroscopically, I was able to measure their black hole masses and red-
shifts. The majority of the candidates have masses between 108M and 10
10M,
similar to previous work; however, our search is more sensitive at higher redshifts,
and the candidates have an almost flat distribution between 1 < z < 2. I then
followed up these candidates with extended-baseline imaging, in order to put their
periodicity to the test and break false positive signals due to the stochastic variabil-
ity of regular quasars. As evidence for a true periodic signal should strengthen over
a longer baseline, we are able to apply more rigorous down-selection criteria to the
sample of 26 candidates, as well as a periodic candidate reported in the literature,
the SMBHB candidate PG 1302−102. PG 1302−102 fails to demonstrate persistent
182
periodicity by our analysis, and at most six candidates from our PS1 MDS sam-
ple are statistically significant, assuming that normal stochastic quasar variability is
modeled by the Damped Radom Walk (DRW) process or a broken-power law (BPL)
PSD. We show that the down-selected sample infers an accumulative SMBHB de-
tection rate of < 0.6 per 1000 quasars, or < 0.12 per deg2 out to z = 2. These upper
limits are consistent with theoretical predictions, but the much-higher rates inferred
from previous work (Charisi et al., 2016; Graham et al., 2015b) are incompatible
with our upper limits.
Among the many questions that this thesis attempts to answer are the ex-
pectations for a true periodic signal (Chapter 4), the extended-baseline monitoring
of periodic candidates (Chapter 3 Section 3.4; Chapter 5 Section 5.5), and the
multi-wavelength (in particular, SED) properties of SMBHB candidates (Chapter 5
Section 5.7.4). However, a few questions remain and deserve further investigation:
1. What would be a convincing detection of periodicity?
2. Can we search for evidence for periodicity in the X-rays?
3. Can we find robust periodic candidates with extended-baseline imaging?
4. Can we search for complementary evidence for an SMBHB?
6.1 What Would be a Convincing Detection of Periodicity?
In Chapter 4, I have reanalyzed the putative periodicity of the SMBHB can-
didate PG 1302−102 and shown that it is inconsistent with the expectation that
183
the significance of a true periodic signal should increase over a longer baseline, even
though the new data have larger photometric uncertainties. While I have performed
the test for a simulated light curve generated from the DRW model plus a periodic
signal (“DRW+periodic”), it can be further improved and made more convincing in
the following ways:
• Generate a number of realizations of the DRW process, add periodic signals
of various amplitudes and periods, and repeat the analysis; how does their
significance increase with the length of the simulated data?
• Generate a number of “DRW-only” realizations using a Monte Carlo method,
and repeat the analysis; since a large fluctuation in the power spectrum can
mimic a periodic signal, does their significance decrease rapidly with data
length?
• Generate a number of realizations of BPL power spectra of various slopes;
compared with the DRW+periodic light curves, how does the signal degrade
with a different choice of the power-law slope?
In addition to applying the technique described in Chapter 4 to simulated
data, one can also use QPOs (a quasi-periodic oscillation is indicated by variability
power concentrated in a narrow range of frequencies, superimposed on a broad
continuum of red noise) reported in the literature as test cases. One of the highest-
significance X-ray QPO in an AGN is detected in the Seyfert galaxy RE J1034+396
(Gierliński et al., 2008), which shows 16 cycles of periodicity with a high quality
184
factor (Q = f/∆f , which indicates the coherence of the signal); recently, a candidate
optical QPO has been reported in a Seyfert galaxy in the Kepler field (Smith et al.,
2018b).
So far, the significance of an individual candidate is evaluated by assuming
a red noise model (i.e. DRW or BPL). Due to our poor knowledge of the red
noise continuum of a given candidate, in our analysis in Chapter 5 Section 5.6, a
candidate is considered significant only if it passes the test under both assumed
models. Thanks to the recent progress made on the analysis of high-quality Kepler
AGN light curves (Aranzana et al., 2018; Smith et al., 2018a), we now have a better
knowledge of the distribution of the (high-frequency) power-law slopes for a sample
of AGNs. Equipped with this knowledge, one can calibrate the false alarm rate
by simulating a population of variable AGNs or quasars that follow this power-law
slope distribution several times over; this statistical approach can then be applied
to a large sample of periodic candidates in order to determine their significance as
an ensemble (Section 6.2).
6.2 Can We Find Robust Periodic Candidates with Extended-baseline
Imaging?
I have demonstrated that extended-baseline monitoring is crucially important
in breaking the false positive signal due to stochastic variability and have shown
through simulated data that a true periodic signal can be identified with more
confidence on a longer baseline. Most of the candidates from the Graham et al.
185
(2015b) sample have extended data from LINEAR, but an extended-baseline analysis
has not been performed on those candidates (with the exception of PG 1302−102).
Charisi et al. (2016) have included data from CRTS and/or iPTF in their reanalysis
of PTF periodic candidates; however, as I have discussed in Chapter 4 Section
4.4, most of their candidates show a decrease in significance when extended data are
included, which may indicate that the detected periodicity does not persist. Despite
being shallower than PS1 MDS, the CRTS and PTF surveys are all-sky surveys and
provide a much larger sample of periodic candidates (∼ 200). We should perform a
reanalysis of the periodic candidates from Graham et al. (2015b) and Charisi et al.
(2016) and treat them with the same rigorous method as described in Chapter 4.
Combined with the statistical method proposed in Section 6.1, the two samples can
potentially yield robust periodic candidates for follow-up studies (e.g. Section 6.4).
Further, one can further extend the observational baseline of those candidates
with new data produced by the Zwicky Transient Facility (ZTF). The new survey
builds on the success of PTF and iPTF and is the third operating phase of the
Palomar family of telescopes. It will operate for three years (2018-2020) and survey
the northern sky with a three-day cadence. While the length of the ZTF survey
is comparatively short for the purpose of identifying robust periodic sources that
vary for many cycles, its magnitude depth is comparable to that of CRTS or PTF
(R ∼ 20.4) and is naturally suited to follow up those candidates on an even longer
baseline.
186
6.3 Can We Search for Evidence for Periodicity in the X-rays?
So far, systematic searches for periodic quasars have only been performed in
the optical time domain, thanks to the abundance, length, and cadence of light
curves produced by large optical time domain surveys. The recently-published
catalog of sources from the all-sky hard X-ray survey with the Swift Burst Alert
Telescope (BAT) can help change this landscape.
The 105-month catalog (Oh et al., 2018) contains 1632 hard X-ray sources at
14−195 keV from Swift BAT, about 900 of which are AGNs (not including blazars).
The catalog includes monthly-sampled light curves extracted from monthly mosaic
images of the sources, and it is worthwhile to apply the technique developed for this
thesis to search for periodic sources in the BAT catalog.
X-ray periodicity has been recently reported in an SMBHB candidate, MCG+11-
11-032, by Severgnini et al. (2018). The source was initially proposed as a dual
AGN identified through its double-peaked [OIII] emission lines in the SDSS spec-
trum (Wang et al., 2009), and follow-up observations have found two components
separated by 0.77 arcsec, or 0.55 kpc (z = 0.036) (Comerford et al., 2012). Sev-
ergnini et al. (2018) further analyzed the BAT light curve of the source and found
a tentative period of 25 months over a 123-month baseline, and the hardness ratio
of the source suggests that the flux variation is intrinsic and not caused by variable
absorption. Its X-ray Telescope (XRT) spectrum is best-fitted by an absorbed power
law plus a reflection component and two narrow emission lines at 6.16 keV and 6.56
keV, respectively (the rest-frame energy of the Fe Kα line is at 6.4 keV), which they
187
propose are produced at the inner edge of the circumbinary disk or by gas bound to
the individual black holes. Interestingly, under the SMBHB hypothesis, the orbital
velocity inferred from the ∆E of the two emission lines (∼ 0.06c) is consistent with
that inferred from the BAT light curve (Severgnini et al., 2018).
Other X-ray variability signatures of an SMBHB have also been explored in
recent simulations. Tang et al. (2018) followed an equal-mass 106M SMBHB from
60 rg to merger (which corresponds to ∼ 100 orbits, or weeks prior to merger) and
find clear periodicity at twice the binary orbital frequency throughout the process in
the soft X-ray band (at 2 keV), which is produced by gas flung outwards by the black
holes and periodically hitting the inner edge of the circumbinary disk. The hard
X-ray band at 10 keV, on the other hand, is dominated by the minidisk emission.
While variability in the hard X-ray band is more stochastic at the beginning of the
process, clear periodicity develops at the late stage of the inspiral and pulsates at
the same frequency as the soft X-ray band. This signature could be used to identify
the EM counterpart of a source for the Laser Interferometer Space Antenna (LISA),
the space-based GW detector which is sensitive to 0.1 mHz - 1 Hz frequencies.
The 105-month BAT catalog also provides spectra constructed from eight en-
ergy bands: 14 − 20 keV, 20 − 24 keV, 24 − 35 keV, 35 − 50 keV, 50 − 75 keV,
75− 100 keV, 100− 150 keV, and 150− 195 keV, in which one can look for excess
X-ray emission in a BAT periodic candidate and follow up with observations with
a higher energy resolution. Enhanced X-ray emission should be a unique signa-
ture of shock-heated minidisks, which would provide complementary evidence for
an SMBHB (Section 6.4).
188
6.4 Can We Search for Complementary Evidence for an SMBHB?
As several previous work have predicted (e.g. Farris et al. 2015; Roedig et al.
2014), the spectrum of an SMBHB system has contributions from three distinct
components — the minidisks, the accretion streams, and the circumbinary disk
— and should distinguish from a regular accreting SMBH. In the analytical work
by Roedig et al. (2014), the streams shock-heat the outer edge of the minidisks
and produce bright hard X-ray emission at ∼ 100 keV due to Compton cooling.
They also argue that a spectral “notch” should be another distinctive signature
of a binary, due to missing radiation from the cavity in the circumbinary disk;
the temperature of the notch lies between the lowest temperature in the minidisks
and the highest temperature in the circumbinary disk and depends on parameters
including black hole mass, binary separation, accretion rate, and radiative efficiency.
The 2D hydrodynamical simulation by Farris et al. (2015) finds similar enhancement
in the X-rays, but they expect that the notch is partially filled by emission from the
streams and unlikely to be noticeable. The recent 3D MHD simulation by d’Ascoli
et al. (2018) has also produced a spectrum of an accreting SMBHB from UV to the X-
rays, where the hard X-ray component has the most contribution from the minidisks.
I can search for this enhanced feature in the NuSTAR hard X-ray spectrum of an
SMBHB candidate down-selected from the PS1 MDS sample, or look for hints of
spectral hardening in the soft X-rays, which may be more observationally feasible.
These features are complementary to the variability signatures of an SMBHB and
would be strong supporting evidence for a binary.
189
Bibliography
Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12
Andrae, R., Kim, D.-W., & Bailer-Jones, C. A. L. 2013, AAP, 554, A137
Aranzana, E., Körding, E., Uttley, P., Scaringi, S., & Bloemen, S. 2018, MNRAS,
476, 2501
Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2014, ApJ, 794, 141
—. 2016, ApJ, 821, 13
Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47
Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665
Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ,
843, 14
Barrows, R. S., Lacy, C. H. S., Kennefick, D., Kennefick, J., & Seigar, M. S. 2011,
New Astronomy, 16, 122
Barth, A. J., & Stern, D. 2018, ArXiv e-prints, arXiv:1803.00691
Bauer, A., Baltay, C., Coppi, P., et al. 2009, ApJ, 696, 1241
Becker, R. H., White, R. L., Gregg, M. D., et al. 2001, ApJS, 135, 227
Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
Bertin, E. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351,
Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset,
D. Ponz, & S. Enrique, 112
Bertin, E., & Arnouts, S. 1996, AAPS, 117, 393
Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the
Pacific Conference Series, Vol. 281, Astronomical Data Analysis Software and
Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228
190
Bhatti, W. A., Richmond, M. W., Ford, H. C., & Petro, L. D. 2010, ApJS, 186, 233
Bianchi, L., Herald, J., Efremova, B., et al. 2011, APSS, 335, 161
Bond, J. R., Jaffe, A. H., & Knox, L. 1998, PRD, 57, 2117
Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53
Bowen, D. B., Campanelli, M., Krolik, J. H., Mewes, V., & Noble, S. C. 2017, ApJ,
838, 42
Bowen, D. B., Mewes, V., Campanelli, M., et al. 2018, ApJL, 853, L17
Britzen, S., Fendt, C., Witzel, G., et al. 2018, MNRAS, 478, 3199
Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2017, ApJ, 834, 111
Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145
Charisi, M., Haiman, Z., Schiminovich, D., & D’Orazio, D. J. 2018, MNRAS, 476,
4617
Chornock, R., Bloom, J. S., Cenko, S. B., et al. 2010, ApJL, 709, L39
Comerford, J. M., Gerke, B. F., Stern, D., et al. 2012, ApJ, 753, 42
Comerford, J. M., Griffith, R. L., Gerke, B. F., et al. 2009a, ApJL, 702, L82
Comerford, J. M., Pooley, D., Barrows, R. S., et al. 2015, ApJ, 806, 219
Comerford, J. M., Gerke, B. F., Newman, J. A., et al. 2009b, ApJ, 698, 956
Cutri, R. M., & et al. 2013, VizieR Online Data Catalog, 2328
d’Ascoli, S., Noble, S. C., Bowen, D. B., et al. 2018, ArXiv e-prints, arXiv:1806.05697
Deane, R. P., Paragi, Z., Jarvis, M. J., et al. 2014, Nature, 511, 57
Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
Detweiler, S. 1979, ApJ, 234, 1100
D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, MNRAS, 436, 2997
D’Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015, Nature, 525, 351
Dorn-Wallenstein, T., Levesque, E. M., & Ruan, J. J. 2017, ApJ, 850, 86
Dotti, M., Montuori, C., Decarli, R., et al. 2009, MNRAS, 398, L73
Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870
191
Edelson, R., Vaughan, S., Malkan, M., et al. 2014, ApJ, 795, 2
Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549
Eggers, D., Shaffer, D. B., & Weistrop, D. 2000, AJ, 119, 460
Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1
Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, ApJS, 201, 23
Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014, ApJ, 783, 134
—. 2015, MNRAS, 446, L36
Ferrarese, L., & Ford, H. 2005, SSR, 116, 523
Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9
Foord, A., Gültekin, K., Reynolds, M., et al. 2017, ApJ, 851, 106
Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220
Forster, K., Green, P. J., Aldcroft, T. L., et al. 2001, ApJS, 134, 35
Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300
Fried, J. W., & Schulz, H. 1983, AAP, 118, 166
Fu, H., Yan, L., Myers, A. D., et al. 2012, ApJ, 745, 67
Fu, H., Zhang, Z.-Y., Assef, R. J., et al. 2011, ApJL, 740, L44
Garcia, A., Sodré, L., Jablonski, F. J., & Terlevich, R. J. 1999, MNRAS, 309, 803
Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJL, 539, L13
Gezari, S., Martin, D. C., Forster, K., et al. 2013, ApJ, 766, 60
Gierliński, M., Middleton, M., Ward, M., & Done, C. 2008, Nature, 455, 369
Gillespie, D. T. 1996, American Journal of Physics, 64, 225
Gold, R., Paschalidis, V., Etienne, Z. B., Shapiro, S. L., & Pfeiffer, H. P. 2014,
PRD, 89, 064060
Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and
Computational Science, Vol. 5, No. 1, p. 65-80, 2010, 5, 65
Graham, A. W., Erwin, P., Caon, N., & Trujillo, I. 2001, ApJL, 563, L11
Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015a, Nature, 518, 74
—. 2015b, MNRAS, 453, 1562
192
Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198
Haiman, Z., Kocsis, B., & Menou, K. 2009, ApJ, 700, 1952
Heinis, S., Kumar, S., Gezari, S., et al. 2016a, ApJ, 821, 86
Heinis, S., Gezari, S., Kumar, S., et al. 2016b, ApJ, 826, 62
Hellings, R. W., & Downs, G. S. 1983, ApJL, 265, L39
Hobbs, G. 2013, Classical and Quantum Gravity, 30, 224007
Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, Classical and Quantum
Gravity, 27, 084013
Honeycutt, R. K. 1992, PASP, 104, 435
Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1
Horne, J. H., & Baliunas, S. L. 1986, ApJ, 302, 757
Ivezic, Z., Tyson, J. A., Acosta, E., et al. 2008, ArXiv e-prints, arXiv:0805.2366
Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, Advancing Astrophysics with
the Square Kilometre Array (AASKA14), 37
Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2004, ApJ, 606, 799
Jester, S., Schneider, D. P., Richards, G. T., et al. 2005, AJ, 130, 873
Jun, H. D., Stern, D., Graham, M. J., et al. 2015, ApJL, 814, L12
Kaiser, N., Burgett, W., Chambers, K., et al. 2010, in Society of Photo-Optical
Instrumentation Engineers (SPIE) Conference Series, Vol. 7733, Society of Photo-
Optical Instrumentation Engineers (SPIE) Conference Series, 77330E
Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
Kharb, P., Lal, D. V., & Merritt, D. 2017, Nature Astronomy, 1, 727
Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502
Kocsis, B., Haiman, Z., & Loeb, A. 2012, MNRAS, 427, 2680
Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJL, 582, L15
Kormendy, J., & Ho, L. C. 2013, ARAA, 51, 511
Kormendy, J., & Richstone, D. 1995, ARAA, 33, 581
Koss, M., Mushotzky, R., Treister, E., et al. 2011, ApJL, 735, L42
Kovacs, G. 1981, APSS, 78, 175
193
Kovačević, A. B., Pérez-Hernández, E., Popović, L. Č., et al. 2018, MNRAS, 475,
2051
Koz lowski, S. 2016, MNRAS, 459, 2787
Koz lowski, S., Kochanek, C. S., Udalski, A., et al. 2010, ApJ, 708, 927
Kramer, M., & Champion, D. J. 2013, Classical and Quantum Gravity, 30, 224009
Kun, E., Frey, S., Gabányi, K. É., et al. 2015, MNRAS, 454, 1290
Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, PASP, 121, 1395
Lehto, H. J., & Valtonen, M. J. 1996, ApJ, 460, 207
Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2015, MNRAS, 453, 2576
Liu, T., Gezari, S., & Miller, M. C. 2018, ApJL, 859, L12
Liu, T., Gezari, S., Heinis, S., et al. 2015, ApJL, 803, L16
Liu, T., Gezari, S., Burgett, W., et al. 2016, ApJ, 833, 6
Liu, X., Shen, Y., Strauss, M. A., & Greene, J. E. 2010, ApJ, 708, 427
Lobanov, A. P., & Roland, J. 2005, A&A, 431, 831
Lomb, N. R. 1976, APSS, 39, 447
Lu, J.-F., & Zhou, B.-Y. 2005, ApJL, 635, L17
MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83
MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, ApJ, 721, 1014
Magnier, E. 2006, in The Advanced Maui Optical and Space Surveillance Technolo-
gies Conference, E50
Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20
Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285
Maness, H. L., Taylor, G. B., Zavala, R. T., Peck, A. B., & Pollack, L. K. 2004,
ApJ, 602, 123
Martini, P. 2004, Coevolution of Black Holes and Galaxies, 169
Matthews, T. A., & Sandage, A. R. 1963, ApJ, 138, 30
McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184
McHardy, I. 2010, The Jet Paradigm: From Microquasars to Quasars, ed. T. Belloni
(Berlin, Heidelberg: Springer Berlin Heidelberg), 203–232
194
McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008
McLure, R. J., & Dunlop, J. S. 2004, MNRAS, 352, 1390
McLure, R. J., & Jarvis, M. J. 2002, MNRAS, 337, 109
Miller, L., Turner, T. J., Reeves, J. N., et al. 2010, MNRAS, 403, 196
Milosavljević, M., & Merritt, D. 2003, in American Institute of Physics Confer-
ence Series, Vol. 686, The Astrophysics of Gravitational Wave Sources, ed. J. M.
Centrella, 201–210
Moore, C. J., Taylor, S. R., & Gair, J. R. 2015, Classical and Quantum Gravity, 32,
055004
Morganson, E., Burgett, W. S., Chambers, K. C., et al. 2014, ApJ, 784, 92
Mushotzky, R. F., Edelson, R., Baumgartner, W., & Gandhi, P. 2011, ApJL, 743,
L12
Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51
Oh, K., Koss, M., Markwardt, C. B., et al. 2018, ApJS, 235, 4
Padmanabhan, N., Schlegel, D. J., Finkbeiner, D. P., et al. 2008, ApJ, 674, 1217
Peng, C. Y., Impey, C. D., Ho, L. C., Barton, E. J., & Rix, H.-W. 2006, ApJ, 640,
114
Perera, B. B. P., Stappers, B. W., Babak, S., et al. 2018, MNRAS, arXiv:1804.10571
Peters, P. C. 1964, Physical Review, 136, 1224
Peterson, B. M. 2001, in Advanced Lectures on the Starburst-AGN, ed. I. Aretxaga,
D. Kunth, & R. Mújica, 3
Pihajoki, P. 2016, MNRAS, 457, 1145
Pojmanski, G. 1997, ACTAA, 47, 467
Press, W. H. 1978, Comments on Astrophysics, 7, 103
Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751
Richards, G. T., Fan, X., Newberg, H. J., et al. 2002, AJ, 123, 2945
Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006a, ApJS, 166, 470
Richards, G. T., Strauss, M. A., Fan, X., et al. 2006b, AJ, 131, 2766
Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49
195
Roedig, C., Krolik, J. H., & Miller, M. C. 2014, ApJ, 785, 115
Rosado, P. A., Sesana, A., & Gair, J. 2015, MNRAS, 451, 2417
Ross, N. P., McGreer, I. D., White, M., et al. 2013, ApJ, 773, 14
Runnoe, J. C., Eracleous, M., Mathes, G., et al. 2015, ApJS, 221, 7
Runnoe, J. C., Eracleous, M., Pennell, A., et al. 2017, MNRAS, 468, 1683
Sazhin, M. V. 1978, Sov. Ast., 22, 36
Scargle, J. D. 1982, ApJ, 263, 835
Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
Schlafly, E. F., Finkbeiner, D. P., Jurić, M., et al. 2012, ApJ, 756, 158
Schmidt, K. B., Marshall, P. J., Rix, H.-W., et al. 2010, ApJ, 714, 1194
Sesana, A., Haiman, Z., Kocsis, B., & Kelley, L. Z. 2017, ArXiv e-prints,
arXiv:1703.10611
Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255
Sesar, B., Stuart, J. S., Ivezić, Ž., et al. 2011, AJ, 142, 190
Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, AJ, 134, 2236
Severgnini, P., Cicone, C., Della Ceca, R., et al. 2018, MNRAS, arXiv:1806.10150
Shannon, R. M., Ravi, V., Lentati, L. T., et al. 2015, Science, 349, 1522
Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48
Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T., & Schneider, D. P. 2008,
ApJ, 680, 169
Shen, Y., & Loeb, A. 2010, ApJ, 725, 249
Shi, J.-M., & Krolik, J. H. 2015, ApJ, 807, 131
Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118
Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ,
325, 628
Simm, T., Salvato, M., Saglia, R., et al. 2016, AAP, 585, A129
Smith, K. L., Mushotzky, R. F., Boyd, P. T., et al. 2018a, ArXiv e-prints,
arXiv:1803.06436
196
Smith, K. L., Mushotzky, R. F., Boyd, P. T., & Wagoner, R. V. 2018b, ApJL, 860,
L10
Smith, K. L., Shields, G. A., Bonning, E. W., et al. 2010, ApJ, 716, 866
Snodgrass, C., & Carry, B. 2013, The Messenger, 152, 14
Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
Sudou, H., Iguchi, S., Murata, Y., & Taniguchi, Y. 2003, Science, 300, 1263
Tang, Y., Haiman, Z., & MacFadyen, A. 2018, MNRAS, 476, 2249
Tang, Y., MacFadyen, A., & Haiman, Z. 2017, MNRAS, 469, 4258
The NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ, 813,
65
Timmer, J., & Koenig, M. 1995, AAP, 300, 707
Tonry, J. L., Stubbs, C. W., Kilic, M., et al. 2012a, ApJ, 745, 42
Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012b, ApJ, 750, 99
Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARAA, 35, 445
Valtonen, M. J., Lehto, H. J., Takalo, L. O., & Sillanpää, A. 2011, ApJ, 729, 33
Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851
Valtonen, M. J., Zola, S., Ciprini, S., et al. 2016, ApJL, 819, L37
Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549
Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692
Vasiliev, E., Antonini, F., & Merritt, D. 2015, ApJ, 810, 49
Vaughan, S., & Uttley, P. 2006, Advances in Space Research, 38, 1405
Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145
Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689
Vestergaard, M., & Wilkes, B. J. 2001, ApJS, 134, 1
Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559
Volonteri, M., Miller, J. M., & Dotti, M. 2009, ApJL, 703, L86
Wang, J.-M., Chen, Y.-M., Hu, C., et al. 2009, ApJL, 705, L76
White, R. L., Becker, R. H., Helfand, D. J., & Gregg, M. D. 1997, ApJ, 475, 479
197
White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341
Wright, E. L. 2006, PASP, 118, 1711
York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, AJ, 120, 1579
Zheng, Z.-Y., Butler, N. R., Shen, Y., et al. 2016, ApJ, 827, 56
Zhu, X.-J., Hobbs, G., Wen, L., et al. 2014, MNRAS, 444, 3709
Zoghbi, A., Reynolds, C., & Cackett, E. M. 2013, ApJ, 777, 24
198