The Astrophysics of Nanohertz Gravitational Waves

The Astronomy and Astrophysics Review (2019) 27:5
https://doi.org/10.1007/s00159-019-0115-7
REVIEW ARTICLE
The astrophysics of nanohertz gravitational waves
Sarah Burke-Spolaor, et al. [full author details at the end of the article]
Received: 8 November 2018 / Published online: 18 June 2019

© The Author(s) 2019
Abstract
Pulsar timing array (PTA) collaborations in North America, Australia, and Europe,
have been exploiting the exquisite timing precision of millisecond pulsars over decades
of observations to search for correlated timing deviations induced by gravitational
waves (GWs). PTAs are sensitive to the frequency band ranging just below 1 nanohertz
to a few tens of microhertz. The discovery space of this band is potentially rich with
populations of inspiraling supermassive black hole binaries, decaying cosmic string
networks, relic post-inflation GWs, and even non-GW imprints of axionic dark matter.
This article aims to provide an understanding of the exciting open science questions
in cosmology, galaxy evolution, and fundamental physics that will be addressed by
the detection and study of GWs through PTAs. The focus of the article is on providing
an understanding of the mechanisms by which PTAs can address specific questions
in these fields, and to outline some of the subtleties and difficulties in each case. The
material included is weighted most heavily toward the questions which we expect will
be answered in the near-term with PTAs; however, we have made efforts to include
most currently anticipated applications of nanohertz GWs.
Keywords Gravitational waves · Stars · Neutron · Galaxies · Evolution · Black hole

physics · Cosmology · Miscellaneous
Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Pulsar timing in brief . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Pulsar timing and timing residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 The pulsar term, the Earth term, and correlation analysis . . . . . . . . . . . . . . . . . . . 5
2.3 Types of gravitational-wave signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3 The population and evolution of supermassive black hole binaries . . . . . . . . . . . . . . . . . 10
3.1 GW emission from supermassive black hole binaries . . . . . . . . . . . . . . . . . . . . . 12
3.1.1 PTAs and the binary life cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.1.2 Continuous waves: binary inspiral . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.1.3 Memory: binary coalescence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.4 The stochastic binary gravitational-wave background . . . . . . . . . . . . . . . . . . 15
3.1.5 Gravitational-wave background anisotropy . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Supermassive binaries and their environments . . . . . . . . . . . . . . . . . . . . . . . . . 18
123
5 Page 2 of 78 S. Burke-Spolaor et al.
3.2.1 The mass function of supermassive black holes . . . . . . . . . . . . . . . . . . . . . 19

3.2.2 Dynamical friction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2.3 Stellar loss-cone scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.4 Viscous circumbinary disk interaction . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.5 Eccentricity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.6 Measuring the spin of supermassive black hole binaries . . . . . . . . . . . . . . . . . 25
4 Cosmic strings and cosmic superstrings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5 The nature of gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.1 How many gravitational wave polarization states exist? . . . . . . . . . . . . . . . . . . . . 32
5.2 Characterizing the dispersion relation of gravitational waves . . . . . . . . . . . . . . . . . 35
5.3 Examples of gravity theories that predict additional polarizations and modified dispersion
relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.3.1 Einstein Aether . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3.2 Massive gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3.3 f (R) Gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6 Relic gravitational waves and early-universe cosmology . . . . . . . . . . . . . . . . . . . . . . 39
6.1 Relic gravitational waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.2 Cosmological measurements and the cosmological constant . . . . . . . . . . . . . . . . . . 41
7 Additional science from “Hidden Planets” to dark matter . . . . . . . . . . . . . . . . . . . . . . 42
7.1 Stellar convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.2 Dark matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.2.1 Cold dark matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.2.2 Scalar-field dark matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
7.3 Primordial black holes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
7.4 Solar system ephemerides and wandering planets . . . . . . . . . . . . . . . . . . . . . . . 44
8 Gravitational-wave synergies: multi-band GW science . . . . . . . . . . . . . . . . . . . . . . . 46
9 Electromagnetic synergies: multi-messenger science . . . . . . . . . . . . . . . . . . . . . . . . 49
9.1 Identifying and tracking an SMBHB host . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
9.2 Direct searches for SMBHBs in electromagnetic bands . . . . . . . . . . . . . . . . . . . . 51
9.3 Topics in binary supermassive black hole multi-messenger science . . . . . . . . . . . . . . 52
9.3.1 Accretion dynamics: active nucleus geometry . . . . . . . . . . . . . . . . . . . . . . 52
9.3.2 Constraining black hole mass–host galaxy relations . . . . . . . . . . . . . . . . . . . 53
9.3.3 Active nuclei preceding and following SMBHB coalescence . . . . . . . . . . . . . . 53
9.3.4 SMBHB effects on galaxy central morphology and dynamics . . . . . . . . . . . . . . 54
9.3.5 Tidal disruption events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
9.3.6 Testing electromagnetic SMBHB emission models . . . . . . . . . . . . . . . . . . . 55
9.3.7 Constraining SMBHB populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
9.3.8 Tests of gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
10 Current outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
1 Introduction
The first direct detection of gravitational waves (GWs) in 2015 started a new era in
astrophysics, in which we can now use gravity itself as a unique messenger from
the cosmos (Abbott et al. 2016a, b, 2017a). The subsequent detection of a neutron
star merger with associated electromagnetic counterparts in 2017 marked the begin-
ning of gravitational waves’ contribution to the era of multi-messenger astronomy
(Abbott et al. 2017b). The ground-based laser interferometers that made those detec-
tions (LIGO/VIRGO Collaboration) are sensitive to sources that emit gravitational
radiation between about ten and a few thousand herz. In the coming decades, we will
open up additional bands in the GW spectrum, which will allow us to probe entirely
new astrophysical sources and physics.
123
The astrophysics of nanohertz gravitational waves Page 3 of 78 5
The detection of nanohertz GWs by Pulsar Timing Arrays (PTAs) is expected to be

the next major milestone in GW astrophysics. In the future, PTAs and ground-based
laser interferometry experiments will be complemented by space-based laser interfer-
ometers (Amaro-Seoane et al. 2017) and observations of primordial GWs, imprinted
in the polarization of the cosmic microwave background (e. g., BICEP2/Keck Collab-
oration et al. 2015), providing comprehensive access to the GW Universe. Current
PTA efforts are spearheaded by a number of groups worldwide, including the Euro-
pean Pulsar Timing Array (EPTA, Babak et al. 2016; Desvignes et al. 2016; Lentati
et al. 2015), the North American Nanohertz Observatory for Gravitational Waves
(NANOGrav, Arzoumanian et al. 2018b) and the Parkes Pulsar Timing Array (PPTA,
Lasky et al. 2016; Manchester et al. 2013; Shannon et al. 2015). The individual groups
are also the constituents of an international collaboration, known as the International
Pulsar Timing Array (IPTA, Hobbs et al. 2010; Verbiest et al. 2016).
This paper presents a comprehensive background in astrophysical theory that can
be addressed observationally by PTAs, and thus the science that will be extracted
from the detection of GWs at nanohertz frequencies. The immediate focus of PTAs
has been a stochastic GW background, hypothesized to result from the ensemble of
in-spiraling supermassive black hole binaries. However, the astrophysics resulting
from detection and study of GWs by PTAs is much richer, and some of it has been
developed alongside steady PTA sensitivity improvements over the past decade. We
limit this paper to describe the astrophysics that is related to GW detection in the PTA
band, and in Sect. 7 to gravitational effects on PTAs not due to the pulsars or their
companions. Throughout the present work, the “PTA band” refers to GW frequencies
of approximately 1–1000 nHz.
We do not aim to cover the rich (non-GW-related) astrophysics accessible by pulsar
timing. The prolific ancillary science from a PTA as a whole includes, but it not limited
to: neutron star population dynamics (Cordes and Chernoff 1998; Matthews et al.
2016), the formation histories of compact objects (Bassa et al. 2016, Fonseca et al.
2016, Kaplan et al. 2016), and the characterization of the ionized interstellar medium
(Jones et al. 2017; Keith et al. 2013; Lam et al. 2016; McKee et al. 2018), including
plasma lensing events (Coles et al. 2015; Lam et al. 2018), tests of general relativity
(Kramer et al. 2006; Taylor and Weisberg 1982; Zhu et al. 2015) and the physics of
nuclear matter (Demorest et al. 2010).
For readers seeking a brief summary, each section is led by an outline of the most
salient overview points from that section. The layout of the remainder of this paper is as
follows: Sect. 2 provides a backdrop of concepts in pulsar timing that are relevant to the
understanding of this review. In Sect. 3, we discuss PTA applications to supermassive
black hole binaries; in Sect. 4, we consider cosmic strings; in Sect. 5, we assess
whether nanohertz GWs can present unique tests of general relativity; in Sect. 6.1 we
consider topics in cosmology, and in Sect. 7 we consider other (potentially more exotic)
possibilities. In Sects. 8 and 9, we describe potential synergistic science in multi-band
GW studies (in particular with the Laser Interferometer Space Antenna, LISA), and
in multi-messenger studies (in particular with electromagnetic observations of binary
supermassive black holes), respectively. In Sect. 10, we summarize the current and
the expected near-future developments in this field.
123
2 Pulsar timing in brief
Here, we introduce the critical concepts for understanding how PTAs can
access their target science.
Timing residuals: Evidence of GWs can be seen by the influence they have on
the arrival time of pulsar signals at the Earth. The measured versus predicted
arrival time of pulses, as a function of time, is referred to as the timing residuals.
Pulsar versus Earth term: A propagating GW will pass both the pulsar and
Earth, affecting their local space–time at different times. Pulsar timing can
detect a GW’s passage through an individual pulsar (“pulsar term”), and can
detect a wave’s passage through the Earth (“Earth term”) as a signal correlated
between pulsars.
Correlation analysis: While we can detect the Earth or pulsar term in one
pulsar, a GW can only be confidently detected by observing the correlated influ-
ence of the GW on multiple pulsars, demonstrating a dominantly quadrupolar
signature.
GW signals:
• Continuous waves from orbiting binary black holes.
• GW bursts from single-encounter supermassive black hole (SMBH) pairs
and cosmic strings.
• Bursts with memory, singular, rapid and permanent step changes in space
time that can accompany SMBH binary (SMBHB) coalescence and cosmic
strings.
• GW background, the combined sum from all sources of GW emission.
2.1 Pulsar timing and timing residuals
We time a pulsar by building a “timing model”, which is a mathematical description

of anything we know about what will affect the arrival times of its pulses at the Earth
(for details on how precision pulse arrival times are measured, see e. g., Lorimer
and Kramer 2012). Effects we know about and model include (but are not limited to)
the time-dependent position of the Earth in the solar system, the natural slowing of a
pulsar’s period due to rotational energy loss, and any orbital motion of the pulsar, if it is
in a binary. The parameters of the timing model are iteratively refined to minimize the
root mean square (RMS) of the “timing residuals”, which are the difference between
the observed pulse arrival time and the arrival time expected based on the timing model.
As a GW moves between the Earth and a pulsar, it alters the local space–time, and
thus changes the effective path length light must travel. By this process, the pulses will
arrive slightly earlier or later than expected. GWs and any other processes influencing
pulse arrival times that are unaccounted for in the timing model will manifest as
structure in the pulsar’s timing residuals. Since a pulsar’s timing model is modified
over time to remove as much structure as possible from the timing residuals (forming
so-called “post-fit” timing residuals), some of the residual structure induced by a GW
will be effectively “absorbed” by the timing model.
123
30
60
20
40
Residual (ns)
Residual (ns)
20 10
0 0
-20 -10
-40
-20
-60
-30
0 500 1000 1500 2000 2500 3000 3500 0 500 1000 1500 2000 2500 3000 3500
Day since observing start Day since observing start
(a) Continuous wave (b) Background

40
20 30
20
Residual (ns)
Residual (ns)
10 10
0 0
-10
-10 -20
-30
-20
-40
0 500 1000 1500 2000 2500 3000 3500 0 500 1000 1500 2000 2500 3000 3500
Day since observing start Day since observing start
(c) Burst with Memory (d) Burst
Fig. 1 Each panel shows pulsar timing residuals for three pulsars (black triangles, red stars, blue circles)
simulated with weekly observing cadence and 1 ns of white noise in their arrival times. The pulsar-to-pulsar
variations demonstrate how the quadrupolar signature of GWs will manifest as correlated timing residuals
in distinct pulsars. Note that 1 ns is not a noise level yet achieved for any pulsar; however, here it allows us
to demonstrate each observable signal type with a high signal-to-noise ratio. Panels are: a continuous waves
from an equal mass 109 M SMBHB at redshift z = 0.01. The distortion from a perfect sinusoid is caused
by self-interference from the pulsar term (Sect. 2.2). In this case, the pulsar term has a lower frequency
because we see the effects on the pulsar from an earlier phase in the SMBHB’s inspiral evolution. This
interferes with the Earth term, which takes a direct path from the source to Earth and therefore is a view of
a more advanced stage of evolution. b A GW background with h c = 10−15 and α = −2/3. c A memory
event of h = 5 × 10−15 , whose wavefront passes the Earth on day 1500. d A burst source with an arbitrary
waveform
Simulated post-fit residuals influenced by a variety of GW signals, which PTAs

are poised to detect, are illustrated in Fig. 1. Residuals for three different pulsars are
shown to demonstrate how the GW signal can vary from pulsar to pulsar. As can be
easily seen, PTAs are sensitive to effects that last on timescales from ∼weeks, which
is the approximate cadence of pulsar observations, to decades, which is how long
PTA experiments have been running. We note that the scale of these GW effects is
realistic given the properties of SMBHBs, but the noise level is optimistically small
by a factor of 20 or more depending on the pulsar. Since realistic signals will not
have such a high signal-to-noise ratio, PTAs time dozens of pulsars to mitigate signal-
to-noise limitations in individual pulsars and search for correlations in their timing
residuals.
2.2 The pulsar term, the Earth term, and correlation analysis
A GW passing through the galaxy perturbs the local space–time at the Earth and at
the pulsar, but at different times. Pulsar timing can detect a GW’s passage through
123
an individual pulsar (pulsar term) and also through the Earth (Earth term). The Earth
term signal is correlated between different pulsars, while the pulsar term is not.1
Nanohertz GW sources of interest are thought to be tens to hundreds of megaparsecs
away and perhaps further (e. g., Rosado et al. 2015)—it is thus well justified to approx-
imate the GW as a plane wave. With this simplifying approximation, we can consider
the influence of a GW on the observed pulse arrival time as a Doppler shift between
the reference frame of the pulsar we observe and the solar system barycenter. The
Doppler-shifted pulsation frequency as measured by an observer at the quasi-inertial
solar system barycenter is given by f obs = (1 + z) f emit . The observed redshift varies
with time, depending on the time-varying influence of the GW on the local space–time
of the pulsar and the solar system barycenter. More specifically, at time t:
p̂i p̂ j
z(t) = [h i j (t) − h i j (t − tl )], (1)
2(1 + ˆ · p̂)
where h i j is the space–time perturbation (typically referred to as the “strain”), and

gives the metric perturbation describing the GW in transverse-traceless gauge. With
the solar system barycenter as a reference position, the parameter p̂ is a vector pointing
to the pulsar position, ˆ is a vector along the direction in which the GW propagates,
ˆ
tl = (l/c)(1 + · p̂), and l is the distance between the pulsar and the solar system
barycenter. The timing perturbation to pulse arrival times is the integral of the redshift
over time (Anholm et al. 2009; Detweiler 1979), which reduces to the difference
between the Earth term (evaluated at time t) and the pulsar term (evaluated at time
t − tl ).
Note the pulsar term, if observed, always depicts an earlier time in the evolution
of a GW signal. This is because the Earth term samples GWs arriving at the Earth
directly from the source, while the pulsar term is associated with a longer path length,
encompassing the GW’s trip to the pulsar from the source and then the traversal of
light from the pulsar to the Earth. This may be an important effect in studying SMBHB
evolution, as discussed in more detail in Sect. 3.2.6.
Figure 1 illustrates the simulated post-fit residuals for four types of GW signals. For
each signal, the residuals of three separate pulsars are shown. Figure 1a, b, representing
continuous GWs from an individual SMBHB and a stochastic GW background from
an ensemble of SMBHBs, respectively, correspond to long-duration signals, for which
both the Earth and pulsar terms are active simultaneously. In these cases, the pulsar
term interferes with the Earth term and lessens the extent to which the residuals
between different pulsars are correlated or anti-correlated. In Fig. 1a, we have modeled
the inspiral phase of an SMBHB, which over the course of its evolution changes its
frequency and phase. Because each pulsar has a different position relative to the Earth
and the GW source, the pulsar terms probe different stages of the SMBHB orbit and
the pulsar terms interfere with the Earth-term signal in slightly different ways. For
burst-like signals, Fig. 1c, d, the Earth term can be active, while the pulsar term is
1 Note, a more precise remark is that the Earth and pulsar terms are not correlated as long as two pulsars are
separated by many gravitational wavelengths, that is to say that f L 1, where f is the GW frequency and
L is distance to the pulsar. This assumption is called the short-wavelength approximation (e. g., Mingarelli
and Mingarelli 2018).
123
quiescent or vice versa. If the Earth term is active, but all pulsar terms are not, the
timing perturbation from the GW will be fully correlated or anti-correlated across all
pulsars in a PTA.
It is important to note that a GW can only be confidently detected by PTAs if
the correlated influence of the GW on multiple pulsars is observed (Taylor et al.
2017a; Tiburzi et al. 2016). Earth-based clock errors will influence all pulsar timing
residuals the same way (monopolar signature, e. g., Hobbs et al. 2012), and errors in our
solar system models will influence ecliptic pulsars more severely (dipolar signature,
e. g., Champion et al. 2010). GWs are expected to have a quadrupolar signature,
the directional correlations of which are expected to depend on the nature of gravity,
the polarization of the GWs, and also the nature of other noise sources affecting the
PTA (Taylor et al. 2017a; Tiburzi et al. 2016). Therefore, the relative positions of a
GW source and two pulsars will dictate how the residuals of those two pulsars are
correlated. How these correlations appear as a function of the angle between pulsars
depends on the nature of gravity and the polarization of a GW, and is commonly
shown for pulsar timing data as the “Hellings and Downs” curve (Hellings and Downs
1983). Section 5 discusses the correlation analysis and the Hellings and Downs curve
in more detail, and describes how various models of gravity will dictate the shape of
the correlations observed.
2.3 Types of gravitational-wave signal
Depending on the origin of the GW signal, the induced residuals might appear as
deterministic signals or a stochastic background. Here, we simply aim to set up a
reference point of what signal modes we expect to detect with PTAs. In the remainder
of the document, we further discuss what information can be extracted about the
universe, depending on the type of the detected GW signal.
Cyclic signals (continuous waves; Fig. 1a) can arise from objects in an actively
orbiting binary system. Bursts (Fig. 1d) represent rapid but temporally finite acceler-
ations, e.g., during the pericenter passage in a highly eccentric or unbound orbit of
two SMBHs (e. g., Finn and Lommen 2010). These classes can be detected by PTAs
as long as their characteristic timescale is between weeks and a few decades.
Bursts with memory (Fig. 1c) represent a rapid and permanent deformation in
spacetime. A burst with memory (BWM) event is generally expected to occur on
timescales less than 1 day (e. g., Favata 2010). Because the duration of memory’s
ramp-up time is relatively brief compared to PTA observing cadence, it is commonly
modeled as an instantaneous step function in strain. PTAs cannot typically detect the
memory event as it occurs, but they can see the resulting difference between the pre- and
post-event spacetimes; the BWM creates a sudden, long-term change in the apparent
period of a pulsar. This leaves a low-frequency ramp-like signature in the pulsar timing
residuals (e. g., Arzoumanian et al. 2015b; Madison et al. 2014; van Haasteren and
Levin 2010; Wang et al. 2015). In Fig. 1c, the memory event at day 1500 causes a
characteristic “ω” shape to be seen in the timing residuals, indicating a difference in
the spacetime before and after the event. The residual shape here is influenced by our
123
fit to the period and period derivative of each pulsar, which is required in the pulsar
timing model. See Madison et al. (2014) and Sect. 3.1.3 for more details on this effect.
Finally, all of nature’s deterministic signals can contribute en masse to a stochas-
tic GW background (Fig. 1b). The strain of the background is frequency dependent,
described by the characteristic strain spectrum, h c ( f ). This is calculated by integrat-
ing the squared GW strain signal over the entire emitting population (Phinney 2001).
For most GW sources in the nanohertz to microhertz frequency regime, the predicted
characteristic strain spectrum can be simplified as a single power law:
α
f
hc( f ) = A , (2)
f yr
or in terms of the energy density of GWs, GW ∝ f 2(1+α) . In this way, the background
can be characterized with an amplitude A, and a single spectral index α. The amplitude
A is commonly defined at a frequency f yr = 1 year−1 ∼ 32 nHz. As discussed
in future sections, the details of physical processes in galaxies, cosmic strings, and
inflation can potentially make the spectrum more complex than a single power law.
Nevertheless, as the PTA sensitivity to a GW background increases, they will be able to
detect the amplitude and spectral shape of the background. In Fig. 1b, the background
appears as a red noise process, because the index α is negative, leading to greater
variations at longer timescales/lower frequencies. For a good introductory review on
methods to detect the stochastic background, we refer readers to Romano and Cornish
(2017).
Figure 2 provides a bird’s-eye view of the PTA sensitivities required to success-
fully breach each area of science that we describe in the remainder of this document.
Regardless of the emission source, PTAs will grow in sensitivity by adding pulsars to
the array, by decreasing the average RMS residuals, and simply by timing pulsars for
a longer duration. Thus, in Fig. 2, we show how the requirements change as a function
of these parameters. As seen in the top (“now”) panel, PTAs have already breached the
expectations for the background of GWs from SMBHBs formed in galaxy mergers,
and are now setting increasingly stringent limits on galaxy/black hole co-evolution.
In the coming years to decades, we expect this to become first a detection of the back-
ground and then become an exploration of the physics of discrete SMBHB systems.
Deeper explorations of gravity, dark matter, and other effects should be soon to follow
thereafter.
123
1000
NANOGrav
Indirect evidence of
12.5 years
Gravitational Waves
es
av
nd sW
RMS Precision (ns)

ou ou
k gr n u
c i
Ba nt
B Co
M BH HB
S B
100 SM
Tests of Gravity
Primordial BHs
Dark Matter
10
1 10 100 1000
Number of Pulsars in Array
1000
Indirect evidence of
Gravitational Waves
und
ro es
c kg av
Ba u sW
RMS Precision (ns)
B uo
BH nt
in Estimated
SM Co
IPTA 2030
B HB
SM
100
Tests of Gravity
Primordial BHs
Dark Matter
10
1 10 100 1000
Number of Pulsars in Array
Fig. 2 Here, we outline the approximate number of pulsars and timing precision required to access various
science, based on current predictions for each signal. The upper and lower panels represent a 10- and 25-
year timing array, respectively. In the top plot, the black curve shows a representative PTA, reflecting the
upcoming NANOGrav 12.5-year data release. That data set contains approximately 70 pulsars; however
the timescale over which each pulsar has been timed ranges from ∼1 to 20 years. The lower plot shows
expectations for the future IPTA, assuming approximately 100 pulsars. Each curve shows pulsars that are
timed to a precision lower than or equal to the indicated RMS timing precision. The location and shape of
the SMBHB regions reflect the scaling relations of Siemens et al. (2013). These assume a detection signal-
to-noise ratio of at least five, and an SMBHB background of h c 1 × 10−15 , which is just below the most
recent limit placed independently by several PTAs on this background source of GWs. A longer-duration
PTA requires less precision and fewer pulsars for a detection because the signal-to-noise ratio scales with
total observing time
123
3 The population and evolution of supermassive black hole binaries
SMBHBs are expected to be the brightest sources of nanohertz GWs. They

form during major mergers between two galaxies, each of which contain their
own SMBH.
Stochastic background
• In the simplest model, an ensemble of SMBHBs should produce a
stochastic GW background with a characteristic strain h c = A f −2/3
(gw ( f ) ∝ f 2/3 ), with an amplitude O[10−15 ] (O[10−9 ]) at a GW
frequency of 1 year−1 .
• The GW background spectrum might be detected with a shallower slope
at frequencies 10 nHz, if SMBHB evolution is accelerated due to strong
interactions with their environments (stars and gas).
• Galaxy merger rates, SMBH-host co-evolution, dynamical relaxation
timescales, and whether SMBHBs might stall at wide separations can all
affect the amplitude scaling of this background.
• PTA constraints or measurement of the background’s amplitude or spectral
shape can give information on all of the above astrophysical uncertainties
for the ensemble SMBHB population.
• A detection of the SMBHB background would confirm the consensus view
that SMBHs reside in most or all massive galaxies.
• Constraints on background anisotropy may indicate local binary clustering
or large-scale cosmic features.
Discrete continuous-wave sources
• Massive or nearby systems may be individually resolvable from the back-
ground. Detections will provide the most direct probe of the early-inspiral
stage of an SMBHB merger, and can provide measurements of the binary’s
position, phase, and an entangled estimate of chirp mass and luminosity
distance (M/DL ). Chirp mass is defined in the section below.
• If detected, the pulsar term can permit M and DL to be disentangled.
Evolution of the waveform over PTA experimental durations is unlikely
for SMBHBs; however, this would also disentangle the M/DL term.
Bursts with memory
• PTA measurement of a burst with memory can provide the date, approx-
imate sky location, and reduced mass over co-moving distance of an
SMBHB coalescence.
It is now broadly accepted that SMBHs in a mass range around 106 –1010 M
reside at the centers of most galaxies (Kormendy and Richstone 1995; Magorrian
et al. 1998), with several scaling relations between the SMBH and galactic-bulge
properties (e.g. M• − σ , M• − Mbulge , Ferrarese and Merritt 2000; Gebhardt et al.
2000), indicating a potential co-evolution between the two. In the standard paradigm
123
Fig. 3 Binary SMBHs can form during a major merger. Pulsar timing arrays’ main targets are continuous-
wave binaries within ∼0.1 pc separation (second panel in the lower figure; Sect. 3.1.2), although we may
on rare occasion detect “GW memory” from a binary’s coalescence (Favata 2010, Sect. 3.1.3). Millions of
such binaries will contribute to a stochastic GW background, also detectable by PTAs (Sect. 3.1.4). A major
unknown in both binary evolution theory and GW prediction is the means by which a binary progresses
from ∼10 pc separations down to ∼0.1 pc, after which the binary can coalesce efficiently due to GWs (e. g.,
Begelman et al. 1980). If it cannot reach sub-parsec separations, a binary may “stall” indefinitely; such
occurrences en masse can cause a drastic reduction in the ensemble GWs from this population. Alternately,
if the binary interacts excessively with the environment within 0.1 pc orbital separations, the expected
strength and spectrum of the expected GWs will change. Image credits: Galaxies, Hubble/STSci; 4C37.11,
Rodriguez et al. (2006); Simulation visuals, C. Henze/NASA; Circumbinary accretion disk, C. Cuadra
of structure formation, galaxies and SMBHs grow through a continuous process of gas
and dark matter accretion, interspersed with major and minor mergers. Major galaxy
mergers form binary SMBHs, and these are currently the primary target for PTAs. In
this section, we lay out a detailed picture of what is not known about the SMBHB
population, how those unknowns influence GW emission from this population, and
what problems PTAs can solve in this area of study.
In Fig. 3, we summarize the life cycle of binary SMBHs. SMBHB formation begins
with a merger between two massive galaxies, each containing their own SMBH.
Through the processes of dynamical friction and mass segregation, the SMBHs will
sink to the center of the merger remnant through interactions with the galactic gas, stars,
and dark matter. Eventually, they will form a gravitationally bound SMBHB (Begel-
man et al. 1980). Through continued interaction with the environment, the binary orbit
will tighten, and GW emission will increasingly dominate their evolution.
Any detection of GWs in the nanohertz regime, either from the GW background or
from individual SMBHBs, will be by itself a great scientific accomplishment. Beyond
that first detection, however, there are a variety of scientific goals that can be attained
from detections of the various types of GW signals. The subsections below discuss
these in turn, first setting up GW emission from SMBHB systems and then detailing
the influence of environmental interactions. Each section describes a different aspect
of galaxy evolution that PTAs can access.
123
Throughout this document, we refer to SMBHB parameters using the following

conventions: SMBH masses m 1 and m 2 have a mass ratio q = m 2 /m 1 defined such
that 0 ≥ q ≥ 1. The total mass is M = m 1 + m 2 and the reduced mass is μ =
m 1 m 2 /(m 1 + m 2 ). Chirp mass is defined as M ≡ (m 1 m 2 )3/5 /(m 1 + m 2 )1/5 . The
binary inclination, eccentricity, and semi-major axis are given by the symbols i, e, and
a, respectively. The parameter D is the radial comoving distance to the binary system.
Other specific parameters will be defined in-line where relevant.
3.1 GW emission from supermassive black hole binaries
3.1.1 PTAs and the binary life cycle
As shown in Fig. 3, SMBHBs can emit discrete PTA-detectable GWs in two phases
of their life cycle. PTAs can detect continuous waves during SMBHBs’ active inspiral
phase, up to a few days before coalescence. PTAs can also detect the moment of
coalescence of an SMBHB by detecting its related burst with memory (covered in more
detail in Sect. 3.1.3). As noted earlier (Sect. 2.2), the “pulsar term” of a binary contains
information about an earlier stage of binary evolution. For a sufficiently strong GW
signal, the pulsar term can be measured in several pulsars, and thus multiple snapshots
of the evolutionary progression of the binary can be detected simultaneously (Corbin
and Cornish 2010; Ellis 2013; Mingarelli et al. 2012; Taylor et al. 2014).
Because binary inspiral is slower at larger separations, the number density is
much higher for discrete systems at the earlier stages of inspiral (that is, at low GW
frequency). At these stages, the binary may still be interacting closely with its environ-
ment. Here, we review deterministic and stochastic GW emission from binary SMBHs,
and in the next sub-sections we develop from this into how environments can influence
the GW signals—and how PTAs can uniquely explore these physical processes.
3.1.2 Continuous waves: binary inspiral
A PTA detection of the correlated signal from a continuous-wave SMBHB (Fig. 1a)
will produce constraints on a system’s orbital parameters (Arzoumanian et al. 2014;
Babak and Sesana 2012; Babak et al. 2016; Ellis 2013; Lee et al. 2011a; Taylor et al.
2016; Zhu et al. 2014), in much the same way as ground-based instruments can with
stellar-mass binaries. PTAs may have potentially poorer parameter estimation; this is
because PTAs will typically observe an early portion of the binary inspiral, and only
have a glimpse of this phase over the ∼1–2 decade observational time spans of PTAs.
During this time, we are unlikely to observe frequency evolution of the binary that
creates the information-rich “chirping” signal seen by ground-based laser interferom-
eters (Taylor et al. 2016), which can allow GW experiments to derive distances to a
binary and a detailed model of the system’s evolution. However, if the system is of
sufficiently high mass, has initially high eccentricity, or is detected in a late stage of
inspiral evolution, then chirping within the observational window may be detectable
(Lee et al. 2011b; Taylor et al. 2016).
123
As a binary evolves and accelerates in its orbit, it has a greater chance at decoupling
from the environment. Somewhere below separations of ∼1 pc, this may occur and the
binary can be considered as an isolated physical system. In this case, the dissipation
of orbital energy will depend only on the constituent SMBH masses, the orbital semi-
major axis, and the binary’s eccentricity. As explored in the sections below, the timing
of this decoupling has a distinct effect on the detectable GW signals from SMBHB
systems. Here, we lay out pure orbital evolution due to gravitational radiation.
The leading order equations for GW-driven orbital evolution are (Peters 1964; Peters
and Mathews 1963):

73 2 37 4
1+e + 96 e
da 64 m 1 m 2 (m 1 + m 2 ) 24
=− ,
dt 5 a 3 (1 − e2 )7/2

121 2
1+ e
de 304 m 1 m 2 (m 1 + m 2 ) 304
=− e , (3)
dt 15 a4 (1 − e2 )5/2
where the derivatives of the orbital separation and eccentricity are averaged over an
orbital period. One should note from 3 that GW emission always causes the eccen-
tricity to decrease, i.e., the binary will become more circular as it inspirals toward
coalescence. For purely circular systems, the GW emission frequency will be twice
that of the orbital frequency, and will evolve as df /dt ∝ f 11/3 .
An important concept here is that of residence times; because the binary’s inspiral
evolves more rapidly fast at smaller separations, it spends less time residing in a high-
frequency state once its inspiral is GW dominated (and accordingly, it spends less time
residing in a state of small separation). Thus, we would naturally expect fewer binary
systems emitting at high GW frequencies, and many more binary systems emitting at
low GW frequencies. As you will see in the next section, this becomes a critical point
in assessing the shape of the GW background.
For a population of eccentric binaries, the GW emission will be distributed across
a spectrum of harmonics of each binary’s orbital frequency. At higher eccentricities,
the frequency of peak emitted GW power shifts to higher and higher harmonics. This
peak will coincide approximately with the pericenter frequency (Kocsis et al. 2012),
such that:
(1 + e)1/2 f K ,r
f peak ≈ , (4)
(1 − e)3/2 (1 + z)
where
f (1 + z)
f K ,r = . (5)
n
Here, f is the (Keplerian, observed, Earth-reference-frame) GW frequency, and z is
source redshift. The parameter n describes the harmonic of the orbital frequency at
which the GWs are emitted; for a circular system, n = 2. Eccentric systems emit at
the orbital frequency itself as well as at higher harmonics (e. g., Wen 2003).
Binary eccentricity and the Keplerian orbital frequency co-evolve in the following
mass-independent way (Peters 1964; Peters and Mathews 1963; Taylor et al. 2016);
123
3/2
f K ,r (e) σ (e0 )
= , (6)
f K ,r (e0 ) σ (e)
where e0 is the eccentricity at some reference epoch, and,

e12/19 121 2 870/2299
σ (e) = 1+ e . (7)
1 − e2 304
Hence, a binary with (for example) e0 = 0.95, when its orbital frequency is 1 nHz,
will have e ≈ 0.3 by the time its orbital frequency has evolved to 100 nHz.
Finally, the strain at which GWs are emitted is often quoted as the “RMS strain”
averaged over orbital orientations:

32 (GM)5/3
h( f K ,r ) = (2π f K ,r )2/3 . (8)
5 c4 D
This equation and those above highlight the fact that continuous-wave detection by
PTAs will enable a measurement of the system’s orbital frequency and eccentric-
ity. However, the strain amplitude is scaled by the degenerate parameters M5/3 /D;
therefore, chirp mass and source distance cannot be directly measured unless there
is orbital frequency evolution observed over the course of the PTA observations, or
unless the host galaxy of the continuous-wave source is identified (Sect. 9). Some
loose constraints on the mass and distance of the continuous wave might be inferred
simply based on the fact that statistically, we expect the first few continuous-wave
detections to be of the heaviest, relatively equal-mass systems, at low to moderate
redshifts (z 1), as demonstrated originally in Sesana et al. (2009).
It is worth noting here that, until now in this section, we have ignored the pulsar
term (Sect. 2.2). Because the Earth term is correlated between different pulsars, it
will always be discovered at a higher S/N than the pulsar term. If the pulsar term
can be measured in multiple pulsars, however, we can map multiple phases of the
binary’s inspiral history. We can exploit this information through a technique known
as temporal aperture synthesis to disentangle chirp mass and distance, as well as
improve the precision of other parameters (Corbin and Cornish 2010; Ellis 2013; Ellis
et al. 2012; Taylor et al. 2014; Zhu et al. 2016). Likewise, if we have many pulsars
and excellent timing precision, then we can potentially place constraints on BH spin
terms in the waveform (Mingarelli et al. 2012). This is discussed further in Sect. 3.2.6.
3.1.3 Memory: binary coalescence
SMBHBs are one of the two leading sources that we expect to produce detectable
GW memory events (the other potential source, as noted later in this paper, is cosmic
string loops). In the case of SMBHBs, the inspiral and even the coalescence produce
the oscillatory waveform that we see as continuous waves. However, the stress tensor
after the binary’s coalescence will differ from its mean before the coalescence; this is
apparent in the “BURST” panel in Fig. 3, where the waveform is offset from zero after
the SMBHB coalescence’s ring-down period. That offset, which grows over the entire
123
past history of the binary’s evolution, but most precipitously during the coalescence, is
the non-oscillatory term we call memory (Braginskii and Thorne 1987; Christodoulou
1991; Favata 2010; Thorne 1992; Zel’dovich and Polnarev 1974). All SMBHBs are
expected to produce a GW memory signal. SMBHBs may produce both linear and non-
linear memory, where the former is related to the SMBHB motion in the final moments
of coalescence, whereas the non-linear signal is produced by the GWs themselves (see
the discussion in Favata 2011). The memory strain from a coalescing binary depends
on the binary parameters and the black hole spin. For a circular binary, the order of
the strain can reasonably be approximated as
√
(1 − 8/3)Gμ 2 Gμ
h h+ sin (i)[17 + cos2 (i)][1 + O(μ2 /M 2 )] 0.02 , (9)
24c2 D c2 D
and the cross-hand polarization term h × goes to zero (Madison et al. 2014).
Note that due to its non-oscillatory nature, this strain deformation is permanent.
As previously noted, this leads to a sudden observed change in the observed period
of pulsars. After this event, timing residuals will begin to deviate from zero with a
linear upward or downward trend. We would observe that ramp as such, if we knew
the intrinsic spin period and spin-down rate of pulsars (instead, we measure these
from timing data). After fitting the timing residuals for the pulsar’s period and period
derivative, one finds the signature shown in Fig. 1c, with the sharp peak of the curve
representing the moment of coalescence. Based on the time and the amplitude of this
signature in the residuals, PTAs can measure the date of coalescence and obtain a
covariant measurement of the SMBHB’s reduced mass and co-moving distance. If the
signature is detected in the Earth term (i.e., correlated between multiple pulsars), a
position of the memory event can be loosely inferred, likely to a few thousand square
degrees depending on the number of pulsars in the array and how well they are timed.
Predictive simulations have found that PTAs are highly unlikely to detect GW
memory from SMBHB mergers due to the extreme rarity of bright events (Cordes
and Jenet 2012; Ravi et al. 2015; Seto 2009). Nonetheless, Cutler et al. (2014) find
memory to be increasingly important for investigating phenomena at high redshift and
discuss its potential for discovering unexpected phenomena. Techniques to search for
memory in PTAs have been developed (e.g., Madison et al. 2014; Pshirkov et al. 2010;
van Haasteren and Levin 2010) and applied to place limits with PTAs (Arzoumanian
et al. 2015b; Wang et al. 2015).
3.1.4 The stochastic binary gravitational-wave background
The superposition of GWs from the large population of SMBHBs predicted by hierar-
chical galaxy formation (Volonteri et al. 2003a) will produce a stochastic background.
The GW background has greater power at lower frequencies (Fig. 1b); we typically
visualize this as a plot of characteristic strain, h c ( f ). Figure 4 demonstrates the effect
on the GW spectrum of variations in assumptions about the SMBHB population. This
connection between PTA characterization of the background and SMBHB evolution
and galaxy dynamical evolution is the focus of the following subsections. Here, we
123
Frequency [yr−1 ]
10−2 10−1 100 101
Fig. 4 The GW spectrum at nanohertz frequencies from supermassive black hole binaries. We adapted the
data from the SMBH binary populations and evolutionary models of Kelley et al. (2017a) and Kelley et al.
(2017b), highlighting the effects of variations in particular binary model parameters on the resulting GW
spectrum. The dashed black line is the spectrum using only the population mass distribution and assuming
GW-driven evolution, and the gray-shaded region represents the uncertainty in the overall distribution of
SMBHB in the universe. The cyan (orange) line is the GW background from a particular realization of an
SMBHB population using a high (low) eccentricity model. The time sampling corresponds to a PTA with
duration of 20 years and a cadence of 0.05 year. The NANOGrav 11 year detection sensitivity and GW
background upper limits (Arzoumanian et al. 2018a) are illustrated with a gray dotted line and triangle,
respectively, while the EPTA (Lentati et al. 2015) and PPTA (Shannon et al. 2015) upper limits are denoted
by a square and circle, respectively. We note that the PPTA limit appears to be most constraining; however,
it is known to be sensitive to the choice of planetary ephemeris; this effect has been accounted for in
subsequent analysis of other PTA data and results in less constraining limits (Arzoumanian et al. 2018a).
Note: the shaded regions are schematic
outline how many discrete continuous-wave sources can contribute to a GW back-

ground.
The calculation reveals the cosmological and astrophysical factors that can influence
the spectral amplitude and shape of the GW background (Rajagopal and Romani 1995;
Sesana et al. 2004; Wyithe and Loeb 2003). In particular, we typically calculate the
characteristic strain spectrum from an astrophysical population of inspiraling compact
binaries (e.g., that shown in Fig. 4) as
∞ ∞ 1 d4 N
h 2c ( f ) = 0 dz 0 dM1 0 dq dz dM1 dq dtr
∞ g[n,e( f K ,r )] dtr
× n=1 (n/2)2 d ln f K ,r h 2 ( f K ,r ) , (10)
where the contributing factors are:

(i) d4 N /dzdM1 dqdtr is the comoving occupation function of binaries per redshift,
primary mass, and mass–ratio interval, where tr measures time in the binary’s
rest frame. (Uncertainties in this quantity are illustrated as the gray shaded region
in Fig. 4.)
(ii) The expression within {· · · } describes the distribution of GW strain over harmon-
ics, n, of the binary orbital frequency when the system is eccentric. As previously
noted, a circular system will emit at twice the orbital frequency, while eccentric
systems emit at the orbital frequency itself as well as higher harmonics. The func-
tion g(n, e) is a distribution function whose form is given in Peters and Mathews
123
(1963). Thus, non-zero eccentricity in a binary redistributes the power of the GW

spectrum, as shown by green-shaded regions in Fig. 4.
(iii) dtr /d ln f K ,r describes the time each binary spends emitting in a particular log-
arithmic frequency interval, the “residence time”, where f K ,r is the Keplerian
orbital frequency in the binary’s rest frame. This is largely controlled by the
impact of the direct SMBHB environment, as shown by red- and blue-shaded
regions in Fig. 4. The effects of environment are explored in much greater detail
in Sect. 3.2 below.
(iv) h( f K ,r ) is the orientation-averaged GW strain amplitude of a single binary as
given by Eq. 8. Note that in Fig. 4, sharp peaks in the orange and cyan lines indicate
contributions from single sources that may be loud enough to be “resolved” from
the background.
In the simple case of a population of circular binaries whose orbital evolution is driven
entirely by the emission of GWs, the form of h c ( f ) is easily deduced. In this case,
g(n = 2, e = 0) = 1 and g(n = 2, e = 0) = 0 such that f = 2 f K ,r /(1 + z), and the
residence time dtr /d ln f K ,r ∝ f −8/3 is given by the quadrupole radiation formula
(Peters 1964). Collecting terms in frequency gives h c ( f ) ∝ f −2/3 , as per Eq. 2 (this
single power law spectrum is shown as the black, dashed line in Fig. 4). 2
The h c ∝ f −2/3 power law GW background spectrum assumes a continuous distri-
bution of circular SMBHBs evolving purely due to GW emission over an infinite range
of frequencies. As the residence time decreases with frequency, the probability that a
binary still exists (i. e., has not coalesced) also decreases at higher frequencies—that
is, there are far fewer binary systems with a high-frequency orbit. At f 10 nHz, the
Poisson noise in the number of binaries contributing significantly to a given frequency
bin becomes important, and realistic GW spectra become ‘jagged’ (e.g., orange and
cyan lines in Fig 4). At even higher frequencies, the probability of a given frequency bin
containing any binaries approaches zero, and the spectrum steepens sharply relative to
the power law estimate in response (e.g., Sesana et al. 2008).3 At the same time, with
fewer sources contributing substantially to the GW spectrum at higher frequencies, the
chance of finding a discrete system that outshines the combined background becomes
larger; in this case, we say that the binary can be “resolved” from the background as
a continuous-wave source.
Binaries with non-zero eccentricity emit GW radiation over a spectrum of harmon-
ics of the orbital frequency, rather than just the second harmonic as in the circular
case. For large eccentricities, this can significantly shift energy from lower frequen-
cies to higher ones, and change the numbers of binaries contributing energy in a given
observed frequency bin. This can thus substantially change the shape of the spectrum
(i.e., green shaded region in Fig 4), decreasing h c at low frequencies ( 10−8 Hz)
and increasing it at higher frequencies ( 10−8 Hz) (Enoki et al. 2007; Huerta et al.
2 It is important to note that most PTA information on GW backgrounds is derived from the lowest few
accessible frequencies in the data sets (the most sensitive point is typically at frequencies ∼ 2/T ). However,
published upper limits often quote a standardized constraint on Aα within a fiducial single-power law model,
for instance in the SMBH binary case by projecting the limit to a frequency of f = 1 year−1 using α = 2/3.
If there are spectral turnovers as described herein, such projections are non-physical, even if practical to
use for comparison of two PTAs.
3 The power law average effectively includes the contribution from fractional binary systems.
123
2015; Kelley et al. 2017b; Rasskazov and Merritt 2017b; Ravi et al. 2014; Taylor et al.
2017b).
Here again, it is worth explicitly tying these ideas back to the effect of binary
residence times on the GW spectrum; while the strain of an individual SMBHB rises
with frequency (Eq. 8), the number of binary systems contributing to the background
falls with frequency, leading to the generally downward-sloped GW spectrum at high
frequencies. The turnover seen at low frequencies for the case of eccentric binaries
and strong environmental influence (green, blue, and red curves) can likewise be
interpreted in part as due to these effects decreasing the residence time of the binaries
at those frequencies: the systems are pushed to evolve much faster through that phase
than a circular, purely GW-driven binary would, therefore fewer systems contribute.
The “turnover frequency” that is seen at low GW frequencies for an eccentric
and/or environmentally influenced population, as well as the shape of the spectrum
before and after that turnover in the spectrum, is rich with information about nuclear
environments, binary eccentricities, and the influence of gas on binary evolution, as
will be explored more fully in Sect. 3.2.
3.1.5 Gravitational-wave background anisotropy
The incoherent superposition of GWs from the cosmic merger history of SMBHBs cre-
ates a GW background, as we have discussed. However, some of these SMBHBs may
be nearby, but not quite resolvable as continuous waves (Sect. 3.1.2). This can induce
departures from isotropy in the GW background (e.g., Mingarelli et al. 2017; Roebber
and Holder 2017; Simon et al. 2014). Moreover, it may be possible for a galaxy clus-
ter to host more than one inspiralling SMBHB, and thus this sky region may present
excess stochastic GW power. Indeed, many physical processes can induce GW back-
ground anisotropy, which can be characterized and detected using methods developed
in numerous works (Conneely et al. 2019; Cornish and van Haasteren 2014; Gair et al.
2014; Mingarelli and Sidery 2014; Mingarelli et al. 2013; Taylor and Gair 2013).
Importantly, GW background anisotropy will influence the shape of the observed
Hellings and Downs curve (Sects. 2.2, 5), leading to different correlation functions
between pulsar residuals than would be observed for an isotropic background (Gair
et al. 2014; Mingarelli et al. 2013). The current limit on GW background anisotropy
is ∼ 40% of the isotropic component (Taylor et al. 2015). Indeed, the expected level
of anisotropy due to Poisson noise in the GW background is expected to be ∼ 20% of
the monopole signal (Mingarelli et al. 2013; Taylor and Gair 2013), in agreement with
current astrophysical predictions based on nearby galaxies (Mingarelli et al. 2017).
These estimates, however, assume a specific M• − Mbulge relation (McConnell and
Ma 2013) for the prediction of anisotropy levels (Fig. 5).
The detection of the isotropic GW background may follow after a GW background
detection (Taylor et al. 2016).
3.2 Supermassive binaries and their environments
We now discuss in detail factors which influence both the number and waveforms
of continuous-wave sources, and the amplitude and shape of the characteristic strain
123
1.0
Angular Power Spectrum
Median Anisotropy
0.8
C /C0
0.6
0.4
0.2
0 10 20 30 40 50
Fig. 5 A model of anisotropy induced by local nearby SMBHBs. Left: a view of the GW intensity induced
by the superposition of many SMBHBs in a random realization of the local universe from Mingarelli et al.
(2017). Here, there are 111 GW sources in the PTA band (at a frequency of 1 nHz for the sake of this image),
with the color scale representing the characteristic strain as a function of sky position. The level of this
anisotropy is determined by taking the angular power spectrum of the background and normalizing it to the
isotropic component, which we have done on the right. Right: angular power spectrum the same GW sky.
Though there are fluctuations, the median anisotropy as a fraction of the monopole is 0.19
spectrum for the gravitational wave background from SMBHBs. PTAs will ultimately
measure at least both continuous waves and the amplitude and spectral shape of the
GW background at various frequencies. Thus, these measurements have the potential
to explore the factors discussed below.
The influence of several of the factors below, such as binary inspiral rate and their
average eccentricity in early evolutionary phases, has covariant effects on the expected
GW signals (Fig 4). Thus, the measurements of PTAs for some of the effects discussed
below can be enhanced by constraints on these factors from other areas of astrophysics,
e.g., through electromagnetic observation of the SMBHB population and through
numerical simulations (Sect. 9).
This section is laid out as follows. Section 3.2.1 discusses how PTAs can directly
measure the SMBH mass function via the influence of this parameter on the GW
background. Except for the SMBH mass, the strain of the GW background at various
frequencies depends on the residence time of the binaries, which in turn depends
on the physical mechanisms that drive SMBHBs to coalescence. As illustrated in
Fig. 3, binaries may be influenced by several external effects, particularly in the phase
immediately preceding continuous-wave GW emission in the PTA band. These effects
are discussed in Sects. 3.2.2, 3.2.3 and 3.2.4. Finally, Sects. 3.2.5 and 3.2.6 review
how PTAs might access information about the eccentricity of binary systems and the
spin of individual SMBHs in a binary.
3.2.1 The mass function of supermassive black holes
The GW amplitude depends strongly on the masses of SMBHB components: h ∝

M5/3 (Eq. 8). Because of that, errors in the assumed SMBH mass distribution may
lead to significant errors in GW background level estimates. Unfortunately, there is
a tendency for different SMBH mass-estimation techniques (stellar kinematics, gas
kinematics, reverberation mapping, AGN emission lines) to systematically disagree
with each other, with stellar kinematics usually giving the highest values. For example,
123
the claims of a ∼ 2 × 1010 M SMBH in NGC 1277 (van den Bosch et al. 2012) were
subsequently found to be too large by factors of 3–5 (Emsellem 2013). Such a bias
is unsurprising: beyond the Local Group, only a few galaxies show evidence for a
central increase in the RMS stellar velocities (see Fig. 2.5 in Merritt 2013) expected in
the presence of an SMBH. That implies the SMBH sphere of influence is unresolved
and we can only measure an upper limit on its mass. Other methods have their own
biases, e.g., different AGN emission lines give different SMBH mass estimates (Shen
et al. 2008). These biases were highlighted in the “bias-corrected” SMBH mass–host
galaxy relations of Shankar et al. (2016); PTA limits on the SMBH background have
also supported the possibility of biased SMBH masses at the upper (> 109 M ) end of
the relation, demonstrating that M• − Mbulge relations must be constrained to below
certain values, otherwise we should have already detected the GW background (Simon
and Burke-Spolaor 2016).
Unlike galaxy masses, only a handful of SMBH masses are directly measured.
Therefore, when constructing an SMBHB population we have to rely on various
SMBH–galaxy scaling relations, such as the relation between SMBH mass and galac-
tic bulge mass: MBH = β Mbulge . While these types of relations have been thoroughly
studied, there may be systematic biases in the SMBH populations they measure (e.g.,
Shen et al. 2008). Unsurprisingly, different studies give different values of β (the most
widely used one is β ≈ 0.003); that issue is analyzed in detail in the Section IIC
of Rasskazov and Merritt (2017b). Given the observed mass and mass ratio distribu-
tion of merging galaxy pairs, β is the only free parameter defining the SMBHB mass
distribution. Rasskazov and Merritt (2017b) came up with the following analytical
approximation for the GW background strain spectrum (assuming zero eccentricity
and triaxial galaxies):
( f / f yr )−2/3
hc( f ) = A , (11)
1 + ( f b / f )53/30
0.83
−16 β
A = 2.77 × 10 , (12)
10−3
−0.68
β
f b = 1.35 × 10−9 Hz . (13)
10−3
As one can see, lower SMBH masses not only lower the GW background amplitude, but
also increase the turnover frequency, since lighter SMBHBs enter the GW-dominated
regime at higher orbital frequencies.
3.2.2 Dynamical friction
When their host galaxies merge together, the resident SMBHs sink to the center of
the resulting galactic remnant through dynamical friction (Antonini and Merritt 2012;
Merritt and Milosavljević 2005). This is the consequence of many weak and long-range
gravitational scattering events within the surrounding stellar, gas, and dark matter
distributions, creating a drag that causes the SMBHs to decelerate and transfer energy
to the ambient media (Chandrasekhar 1943). The most simple treatment of this phase
123
(resulting in a gross overestimate of the dynamical friction timescale) considers a point

mass (the black hole) travelling in a single isothermal sphere (the galaxy). In this case,
the inspiral timescale is in the order of 10 Gyr (Binney and Tremaine 1987):
2 108 M
19 Re σ
TDF ≈ Gyr, (14)
ln 5 kpc 200 km s−1 M•
where Re and σ are the galaxy’s effective radius and velocity dispersion, M• is the
SMBH’s mass and ln is the Coulomb logarithm.4 However, the braking of the
individual SMBHs in reality will be much shorter. In a genuine merger system, the M•
in the denominator cannot be modeled with just the SMBHs, as they will initially be
surrounded by their constituent galaxies, and later by nuclear stars. After the galaxies
begin to interact, each SMBH is carried by its parent galaxy through the early stages
of merger as the galaxies are stripped and mixed into one. The in-spiral timescale
is dominated by the lower-mass black hole (in the case of an unequal mass–ratio
merger), which along with a dense core of stars and gas within the SMBH influence
radius will be left to inspiral on its own. Extending the above equation to include a
more realistic model of tidal stripping, the resultant timescale can often be shorter
than 1 Gyr (Dosopoulou and Antonini 2017; Kelley et al. 2017a; Yu 2002).5
By extracting energy from the SMBHs on the kiloparsec separation scale, dynamical
friction is a critical initial step toward binary hardening and coalescence. For systems
with extreme mass ratios ( 10−2 ) or very low total masses, dynamical friction may
not be effective at forming a bound binary from the two SMBH within a Hubble time.
In this case, the pair might become “stalled” at larger separations, with one of the two
SMBH left to wander the galaxy at ∼kpc separations (Dvorkin and Barausse 2017;
Kelley et al. 2017a; Kulier et al. 2015; McWilliams et al. 2014; Yu 2002). It is possible
that a non-negligible fraction of galaxies may have such wandering SMBH, some of
which may be observable as offset-AGN, discussed in Sect. 9.
3.2.3 Stellar loss-cone scattering
At parsec separations, dynamical friction becomes an inefficient means of further

binary hardening. At this stage, the dominant hardening mechanism results from indi-
vidual three-body scattering events between stars in the galactic core and the SMBH
binary (Begelman et al. 1980). Stars slingshot off the binary which can extract orbital
energy from the system (Mikkola and Valtonen 1992; Quinlan 1996). The ejection of
stars by the binary leads to hardening of the semi-major axis, and eccentricity evolution
is usually parametrized as (Quinlan 1996):
4 This is related to the impact factor of galactic material and the binary. For circular binaries to moderate
eccentricities, the value should be around 2 ln 5, where ln 5 applies to the case of a circular
binary (Gualandris and Merritt 2008).
5 On the other hand, some simulations have shown that in galaxies that are rapidly stripped or natively
lack a dense stellar core, these inspiral timescales may become indefinite (Tremmel et al. 2018). Note
that pessimistic assumptions about these and other evolution uncertainties still give predictions resulting in
detectable GW signals (e. g., Taylor et al. 2016).
123
da
= − Gρ
σ Ha ,
2 (15)
dt
de
= Gρ
σ HKa , (16)
dt
where H is a dimensionless hardening rate, and K is a dimensionless eccentricity

growth rate. Both of these parameters can be computed from numerical scattering
experiments (e.g., Sesana et al. 2006).
However, only stars in centrophilic orbits with very low angular momentum have
trajectories which bring them deep enough into the galactic center to interact with the
binary. The region of stellar-orbit phase space that is occupied by these types of stars
is known as the “loss cone” (LC; Frank and Rees 1976). Stars which extract energy
from the binary in a scattering event tend to be ejected from the core, depleting the
LC. In general, the steady-state scattering rate of stars is expected to be relatively
low as stars are resupplied to the LC at larger radii where relaxation from star–star
scatterings is slow. Like with dynamical friction at larger scales, binaries can also
stall here, at parsec scales, due to inefficacy of the LC, which is typically known as
the “final parsec problem” (Milosavljević and Merritt 2003b, a). Generally, binaries
that do not reach sub-parsec separations will be unable to merge via GW emission
within a Hubble time (Dvorkin and Barausse 2017; McWilliams et al. 2014). Some
simulations suggest, however, that even in the case of a depleted, steady-state LC, the
most massive SMBHBs, which dominate in GW energy production and also tend to
carry the largest stellar masses, may still be able to reach the GW-dominated regime
and eventually coalesce (Kelley et al. 2017a).
Various mechanisms have been explored to see whether the LC can be efficiently
refilled or populated to ensure continuous hardening of the binary down to milliparsec
separations. In general, any form of bulge morphological triaxiality will ensure a
continually refilled LC that can mitigate the final parsec problem (Khan et al. 2013;
Vasiliev et al. 2014, 2015; Vasiliev and Merritt 2013). Isolated galaxies often exhibit
triaxiality, and given that the SMBH binaries of interest are the result of galactic
mergers, triaxiality and general asymmetries can be expected as a natural post-merger
by-product. Also, post-merger galaxies often harbor large, dense molecular clouds
that can be channeled into the galactic center, acting as a perturber for the stellar
distribution that will refill the LC (Young and Scoville 1991), or even directly hardening
the binary (Goicovic et al. 2017). Finally, since binary coalescence times can be of the
order of Gyrs, while galaxies can undergo numerous merger events over cosmic time
(e.g., Rodriguez-Gomez et al. 2015), subsequent mergers can lead to the formation
of hierarchical SMBH systems (Amaro-Seoane et al. 2010; Bonetti et al. 2018; Ryu
et al. 2018). In this scenario, a third SMBH can not only stir the stellar distribution
for LC refilling, but may also act as a perturber through the Kozai–Lidov mechanism
(Kozai 1962; Lidov 1962) wherein orbital inclination can be exchanged for eccentricity
(Amaro-Seoane et al. 2010; Blaes et al. 2002; Makino and Ebisuzaki 1994). Not only
could a third SMBH more effectively refill the LC, but it could also increase the SMBH
binary eccentricity which speeds up GW inspiral (see Sect. 3.1.2).
Even in the absence of a third perturbing SMBH, binary eccentricity can be
enhanced through stellar LC scattering. This has been observed in many numeri-
123
cal scattering experiments (Quinlan 1996; Sesana et al. 2006), where the general trend
appears to be that equal-mass binaries (most relevant for PTAs) with very low initial
eccentricity will maintain this or become slightly more eccentric. For binaries with
moderate to large eccentricity (or simply with extreme mass ratios at any initial eccen-
tricity), the eccentricity can grow significantly such that high values are maintained
even into the PTA band (Roedig and Sesana 2012; Sesana 2010).
The rotation of the stellar distribution (when the stars have nonzero total angular
momentum) can impact the evolution of the binary’s orbital elements. In particular,
a stellar distribution that is co-rotating with the binary will tend to circularize its
orbit. But if the stellar distribution is counter-rotating with respect to the binary, then
interactions with stars in individual scattering events are more efficient at extracting
angular momentum from the binary, enhancing the eccentricity to potentially quite
high values (e > 0.9) (Mirza et al. 2017; Rasskazov and Merritt 2017a; Sesana et al.
2011). However, the binary’s orbital inclination (with respect to the stellar rotation axis)
also changes: it is always decreases, so that in the end the initially counter-rotating
binaries tend to become co-rotating with the stellar environment (Gualandris et al.
2012; Rasskazov and Merritt 2017a). In most cases, this joint evolution of eccentricity
and inclination leads to the eccentricities at PTA orbital frequencies being much lower
than we would expect from a non-rotating stellar environment model (Rasskazov and
Merritt 2017b).
Interaction of a binary with its surrounding galactic stellar distribution will lead to
attenuation of the characteristic strain spectrum of GWs at low frequencies (i.e., blue
shaded region in Fig. 4). This can be separated into two distinct effects: (1) the direct
coupling leads to an accelerated binary evolution, such that the amount of time spent
by each binary at low frequencies is reduced; (2) extraction of angular momentum
by stellar slingshots can excite eccentricity, which leads to faster GW-driven inspiral
and (again) lower residence time at low frequencies. Assuming an isothermal density
profile for the stellar population6 , we can re-arrange 15 to deduce the orbital frequency
evolution, and hence the evolution of the emitted GW frequency, such that df /dt ∝
f 1/3 . Inserting this into 10 and collecting terms in frequency gives h c ( f ) ∝ f , which
is markedly different from the fiducial ∝ f −2/3 circular GW-driven behavior. The
excitation of binary eccentricity by interaction with stars will further attenuate the
strain spectrum at low frequencies, leading to an even sharper turnover (e.g., Taylor
et al. 2017b).
3.2.4 Viscous circumbinary disk interaction
At centiparsec to milliparsec separations, viscous dissipation of angular momentum

to a gaseous circumbinary disk may play an important role in hardening the binary
(Begelman et al. 1980; Kocsis and Sesana 2011). This influence will depend on the
details of the dissipative physics of the disk; however, the simple case of a binary
exerting torques on a coplanar prograde disk has a self-consistent non-stationary ana-
lytic solution (Ivanov et al. 1999). These studies have assumed a geometrically thin
optically thick disk whose viscosity is proportional to the sum of thermal and radiation
6 stellar density ρ(r ) ∝ r −1 , and constant velocity dispersion σ (r ) = σ .

0
123
pressure (the so-called α-disk, Shakura and Sunyaev 1973). The binary torque will
dominate over the viscous torque in the disk, leading to the formation of a cavity in
the gas distribution and the accumulation of material at the outer edge of this cavity
(i.e., Type II migration). The excitation of a spiral density wave in the disk torques
the binary and leads to hardening through the following semi-major axis evolution
(Haiman et al. 2009; Ivanov et al. 1999):
da 2ṁ 1
=− (aa0 )1/2 , (17)
dt μ
where ṁ 1 is the mass accretion rate onto the primary BH, and a0 is the semi-major
axis, at which the disk mass enclosed is equal to the mass of the secondary BH, given
by Ivanov et al. (1999) as
5/7
α 4/7 m2
a0 = 3 × 10 3
10−2 106 M
−11/7
m1 ṁ 1 −3/7
× 100d rg , (18)
108 M Ṁ E
where α is a disk viscosity parameter, Ṁ E = 4π Gm 1 /cκT is the Eddington accretion

rate of the primary BH (κT is the Thompson opacity coefficient), and r g = 2Gm 1 /c2
is the Schwarzchild radius of the primary BH.
The disk–binary dynamics may be much more complicated, for example, the
disk may be composed of several physically distinct regions (Shapiro and Teukol-
sky 1986), discriminated by the dominant pressure (thermal or radiation) and opacity
contributions (Thompson or free-free). Additionally, high-density disks (equivalently:
high-accretion rates) may provide rapid hardening, but may also be unstable due to self-
gravity. Furthermore, while the system will initially pass through “disk-dominated”
Type II migration (where the secondary BH is carried by the disk like a cork floating
in a water drain), it will generally transition to “planet-dominated” migration (where
the secondary BH is dynamically dominant over the disk) which can be significantly
slower.
Eccentricity growth may be significant during this disk-coupled phase (Armitage
and Natarajan 2005; Cuadra et al. 2009), although there are no generalized prescrip-
tions of the form of 16. The growth of eccentricity is driven by outer Lindblad resonant
interaction of the binary with gas in the disk at large distances (Goldreich and Sari
2003). However, Roedig et al. (2011) found that binaries with high initial eccentric-
ity will experience a reduction in eccentricity, leading to the discovery of a limiting
eccentricity for disk–binary interactions that falls in the interval ecrit ∈ [0.6, 0.8].
The emerging picture is that in a low-eccentricity orbit, the density wake excited by
the secondary BH will lag behind it at apoapsis, causing deceleration and increasing
eccentricity. Whereas in a high-eccentricity orbit, the density wake instead advances
ahead of the secondary BH, causing a net acceleration and reduction in eccentricity.
All previously mentioned studies considered prograde disks, but if a retrograde disk
forms around the binary then the eccentricity can grow rapidly, leading to significantly
diminished GW-inspiral time (Schnittman and Krolik 2015).
123
Coupling of a viscous circumbinary disk with a SMBH binary, like in the stellar
LC scattering scenario, will lead to attenuation of the characteristic strain spectrum
of GWs through both direct coupling and excitation of eccentricity (e.g., red shaded
region in Fig. 4). Focusing on direct coupling of a circular binary to its surround-
ing disks, Kocsis and Sesana (2011) studied scaling relations for the strain spectrum
from different disk–binary scenarios, varying from h c ( f ) ∝ f −1/6 for secondary-
dominated type-II migration in a radiation-dominated α-disk, to h c ( f ) ∝ f 1/2 for
the Ivanov et al. (1999)-model in Eq. 17. Across all models, the characteristic strain
spectrum can be flattened or even increasing due to disk coupling (very different from
the fiducial ∝ f −2/3 circular GW-driven behavior), and spectral attenuation is further
enhanced through disk excitation of binary eccentricity. When the disk–binary models
are applied to populations of SMBHB, the overall GW background spectra tend to
more closely resemble the canonical −2/3 power law, because each disk regime only
applies to a fairly narrow frequency range for a given binary mass (Kelley et al. 2017a;
Kocsis and Sesana 2011).
3.2.5 Eccentricity
The influence of an initial eccentricity on SMBH evolution, without the influence

of external driving factors as in the previous subsection, is shown in Eqs. 6 and 7
There have been several studies of the influence of non-zero binary eccentricity on the
characteristic strain spectrum of nanohertz GWs (Enoki et al. 2007; Huerta et al. 2015;
Kelley et al. 2017b; Rasskazov and Merritt 2017b; Ravi et al. 2014; Taylor et al. 2017b).
The exact shape and amplitude of the spectrum will depend on the detailed interplay
of direct environmental couplings with eccentricity, and, in the case of the latter, the
distribution of eccentricities at binary formation (Ravi et al. 2014; Ryu et al. 2018).
Broadly speaking, (1) eccentricity increases the GW luminosity of the binary, meaning
it evolves faster and thus spends less time emitting in each frequency resolution bin; (2)
eccentricity distributes the strain preferentially toward higher harmonics of the orbital
frequency. These effects lead to the strain being diminished at lower observed GW
frequencies, but also somewhat enhanced at higher frequencies—i.e., the spectrum
can exhibit a turnover to a positive slope at low frequencies, but then a small “bump”
enhancement at the turnover transition.
In simulated SMBHB populations that include eccentricity, non-zero eccentricities
tend to reduce the mean occupation number at lower frequencies, thus making the
stochastic background appear to have a flatter (or inverted) spectrum than the stan-
dard α = 2/3. However, these eccentricities also work to redistribute the power to
higher frequencies, where in a circular population the background would be otherwise
dominated by small numbers of binaries (Kelley et al. 2017b).
3.2.6 Measuring the spin of supermassive black hole binaries
When GWs transit our galaxy, they perturb pulsar signals both at the Earth and at the
pulsar, i.e., they affect both the “Earth term” and the “pulsar term”; as a reminder, the
pulsar term shows a GW signal that is delayed in time by an amount proportional to
the light travel time between the Earth and the pulsar (Sect. 2.2). That is, if a source
123
2 kpc
Earth term
3 kpc
Galaxy merger and
1 kpc subsequent supermassive
black hole merger
pulsar terms Gravitational waves
Fig. 6 Gravitational waves spanning thousands of years in a binary’s evolutionary cycle can be detected
from a continuous GW source by using the pulsar term. As an example, we have drawn a few pulsars with
line-of-sight path length differences to the Earth. These relative time delays between the pulsar terms can
be used to probe the evolution and the dynamics of an SMBHB systems over these many thousands of
years. Right: a major galaxy merger leads to the creation of an SMBHB, emitting nanohertz GWs. Left: the
pulsar term from each pulsar probes a different part of the SMBHB evolution, since they are all at different
distances from the Earth. The blue sinusoid is a cartoon of the GW waveform and shows that the pulsar
terms can be coherently concatenated to probe the binary’s evolution, allowing one to measure, e.g., the
spin of the SMBHB (Mingarelli et al. 2012)
(such as a SMBHB) is evolving, the pulsar term encodes information about earlier
phases in the SMBH evolution. We can use this information to our advantage: when
a continuous GW signal is detected, one can look for the perturbation caused by the
same source in the pulsars, but thousands of years ago. These pulsar terms can be
used to map the evolution of an SMBHB system over many thousands of years: each
pulsar term is a snapshot of the binary during a different point in the history of its
evolution (Fig. 6; Mingarelli et al. 2012), and the phase evolution of the SMBHB can
thus be measured. This is important, since SMBH spins affect the phase evolution of
the binary, thus, constraining the phase evolution allows one to constrain the SMBH
spins (Mingarelli et al. 2012).
One estimates the number of expected gravitational wave cycles observed at the
Earth via the post-Newtonian expansion (Blanchet 2006), which is a function of the
SMBH mass and spin. For example, over a 10-year observation, an equal-mass 109 M
SMBHB system with an Earth term frequency of 100 nHz should produce 32.1 gravi-
tational wave cycles, of which 31.7 are from the leading Newtonian order (or p0 N), 0.9
wave cycles are from p1 N order, and −0.7 are from p1.5 N order. This last term is from
spin–orbit coupling and depends on the SMBH spins. Accessing the pulsar term when
it arrives at the Earth gives information about the SMBHB system over ∼ 3000 years
ago (roughly equivalent to the typical light travel time between the Earth and the pul-
sar). Over this time, one expects 4305.1 wave cyles, of which 4267.8 are Newtonian,
77.3 come from p1 N order, −45.8 are from p1.5 N order, etc... One can therefore see
at a glance that spinning binaries evolve more quickly, which in turn affects the phase
evolution of the waveform. This signal is imprinted in the pulsar terms of the pulsars
in the array, and is therefore only accessible via PTA observations of the pulsar terms.
However, to do pulsar term phase matching, we require that 2π f L < 1 to not lose
a single wave cycle, where f is the GW frequency and L is the distance to the pulsar.
123
Therefore, it is in principal necessary to measure the pulsar distances to, e.g., ∼ 1.5 pc
for a GW signal at 1 nHz.
Many pulsar distances are poorly constrained, since most are estimated via the dis-
persion measure of the pulse (e. g., Cordes and Lazio 2002; Yao et al. 2017). However,
a number of nearby, well-timed pulsars have accurate position measurements based on
parallax measured by very long baseline interferometry (e. g., Deller et al. 2019). For
instance, the well-timed millisecond pular PSR J0437–4715 has a distance measure-
ment of 156.3 ± 1.3 pc (Deller et al. 2008), and is thus suitable for this measurement.
Pulsar timing can also be used to obtain a parallax measurement to the pulsar, as in
Lyne and Graham-Smith (2012), but with relatively large errors. Optical surveys such
as Gaia (Gaia Collaboration et al. 2018) can be used to measure parallaxes to some
pulsars’ white dwarf companions (Jennings et al. 2018), though again with limited
accuracy due to the low brightness of the white dwarfs. The independent distance
measurements to the pulsar’s binary companion can also be used in combination with
the pulsar distance measurements to improve this estimate (Mingarelli et al. 2018). The
most precise way to constrain pulasr distances is through measuring the binary’s orbital
period derivative—a dynamical distance measurement (Shklovskii 1970)—which in
the case of PSR J0437–4715 furthers its distance constraints to 156.79 ± 0.25 pc
(Reardon et al. 2016).
Thus, while current pulsar distances are generally not suitable for an extensive
study of this effect, future instruments such as the Square Kilometre Array (e. g., Smits
et al. 2011) or Next-generation Very Large Array7 and NASA’s WFIRST telescope
(The WFIRST Astrometry Working Group et al. 2017) hold tremendous promise for
enabling precise pulsar (and binary companion) distance measurements, which will
also enable more tests of fundamental physics.
4 Cosmic strings and cosmic superstrings
Cosmic strings and superstrings produce GWs in the nanohertz band. Cosmic
strings are topological defects that can form during phase transitions in the
early universe, while cosmic superstrings are fundamental strings stretched
to cosmological scales due to the expansion of the universe. In a cosmolog-
ical setting, and for the most simple superstring models, both cosmic string
and superstring networks evolve in the same way. Cosmic (super)strings can
exchange partners when they meet and produce loops when they self-intersect.
These loops then oscillate and lose energy to GWs generating a stochastic back-
ground, with a power law spectrum having a spectral index close to that of
SMBHBs. Detection of such a background, or GWs from an individual cosmic
(super)string loop, would therefore provide a unique window into high-energy
physics. Currently, pulsar timing experiments produce the most constraining
bounds on the energy scale and other model parameters of cosmic strings and
superstrings.
7 http://ngvla.nrao.edu/page/scibook
123
Topological defects are a generic prediction of grand unified theories. As the uni-
verse expands and cools, the symmetries of the grand unified theory are broken down
into the Standard Model in one or more stages called phase transitions. At each of
these phase transitions, topological defects generically form, with the type of defect
depending on what symmetry is being broken. Cosmic strings are a one-dimensional
(or line-like) type of topological defect that can form in the early universe during one
(or more) of these phase transitions. The other common types of topological defects
are monopoles and domain walls. Both of these are ruled out, however, because they
lead to cosmological disasters (e.g., the monopole problem), and much of the atten-
tion of the theoretical cosmology community has focused on strings as the only viable
candidate that could lead to potential observational signatures. The most simple sym-
metry breaking that leads to the formation of cosmic strings occurs in the Abelian
Higgs model, where the symmetry group U (1) breaks
U (1) → 1
at some temperature, or energy scale, Ts . They are characterized by their mass per unit
length μ, which in natural units is determined by the temperature at which the phase
transition takes place, μ ∼ Ts2 . The tension of a string is normally given in terms
of the dimensionless parameter Gμ/c2 . which is just the ratio of the string energy
scale to the Planck scale squared. It is worth pointing out that analogous processes
abound in condensed matter systems such as superfluid helium, Bose–Einstein con-
densates, superconductors, and liquid crystals, which can also lead to the production
of topological defects.
Phase transitions in the early universe can therefore lead to the formation of a
network of cosmic strings. Analytic work and numerical simulations show that the
network quickly evolves toward an attractor called the scaling solution. In this regime,
the statistical properties of the system—such as the correlation lengths of long strings
and the size of loops—scale with the cosmic time, and the energy density of the
string network becomes a small constant fraction of the radiation or matter density.
This attractor solution is possible due to reconnections: when two strings meet they
exchange partners (“intercommute”), and when a string self-intersects it chops off a
loop (see Fig. 7). The loops produced by the network oscillate, generate gravitational
waves, and shrink, gradually decaying away. This process removes energy from the
string network, converting it to gravitational waves, and providing the very signal we
seek to detect. The way the scaling solution works is as follows: if the density of
strings in the network becomes large, then strings will meet more often and reconnect,
producing extra loops which then decay gravitationally, removing the surplus string
from the network. If, on the other hand, the density of strings becomes too low, strings
will not meet often enough to produce loops, and their density will start to grow. In
this way, the cosmic string network finds a stable equilibrium density and a Hubble
volume of the universe with a string network statistically always looking like that
shown in Fig. 8, stretched by the cosmic time.
Superstrings refer to the fundamental strings of string theory that in some models
can be stretched to cosmological scales by the expansion of the universe. Fortunately,
much of what we have learned about the evolution of cosmic string networks can be
123
Fig. 7 A depiction of the production of a cosmic string loop from a self-intersecting string. The loop
oscillates and produces gravitational waves slowly decaying away. This process allows the string network
to reach the scaling regime
Fig. 8 Simulation of a cosmic

string network in the matter era.
Long strings are shown in blue
and loops are color coded to
show their size, with red being
the shortest to yellow being the
largest. Image credit K.D. Olum
applied to the evolution of cosmic superstrings. Aside from the possibility of forming
more than one type of string, the most significant difference is that cosmic super-
strings do not always reconnect when they meet. This occurs for two reasons: (i)
because these string theory models require more than three dimensions, and though
strings may appear to meet in three dimensions they can miss each other in the extra
dimensions, and (ii) because superstrings, being the fundamental strings of string the-
ory, interact with a probability proportional to the string coupling squared. The net
effect is to lower the reconnection probability from p = 1 for cosmic strings to p ≤ 1
for cosmic superstrings. A network of cosmic superstrings still enters the scaling
regime, albeit at a density higher by a factor inversely proportional to the reconnection
probability (because strings need to interact all that more often to produce loops that
gravitationally radiate away). The smaller reconnection probability of superstrings
therefore actually increases the chances of finding observational signatures because
the equilibrium string density of the scaling solution is higher. Figure 9 shows the
123
Fig. 9 Plot of regions of the

cosmic (super)string Gμ/c2 - p
plane excluded by present
experiments LIGO and PTAs.
PTAs currently place the
strongest constraints on cosmic
strings. LISA has the potential to
improve these limits (or provide
a detection). The excluded areas
are to the right of each curve.
The figure is from
Blanco-Pillado et al. (2018)
areas of the cosmic (super)string Gμ/c2 - p parameter space excluded by the present
and potential future experiments, including LIGO and the three leading PTA experi-
ments (PPTA, NANOGrav, and EPTA) as of 2017, as well as the planned LISA space
mission. As the reconnection probability decreases, the density of strings in the scal-
ing regime increases, and thus experiments are sensitive to lower and lower string
tensions. Following our previous discussion, the limits at p = 1 represent the limits
specifically for cosmic strings.
In the 1990s, cosmic microwave background data ruled out cosmic (super)strings
as the primary source of density perturbations, placing constraints on the string tension
at the level of Gμ/c2 10−7 , relegating them to at most a secondary role in structure
formation. However, cosmic (super)strings are still potential candidates for the gen-
eration of a host of interesting astrophysical phenomena: including ultrahigh-energy
cosmic rays, fast radio bursts, gamma ray bursts, and of course gravitational radiation.
Clearly, any positive detection of an observational signature of cosmic (super)strings
would have profound consequences for theoretical physics.
PTAs are currently the most sensitive experiment for the detection of cosmic
(super)strings and will remain so for more than a decade and a half. Correspond-
ingly, the most constraining upper limits on the energy scale of cosmic (super)strings
come from PTA analyses. As of the writing of this paper, the most constraining upper
limit published by a PTA collaboration (for p = 1) is Gμ/c2 < 5.3(2) × 10−11 from
the the NANOGrav collaboration (Arzoumanian et al. 2018a). Later, Blanco-Pillado
et al used results from all PTAs and recalculated upper limits on the string tension;
Fig. 10 shows the stochastic background spectrum produced by cosmic strings in terms
of the dimensionless density parameter versus frequency for dimensionless string
tensions Gμ/c2 in the range 10−23 -10−9 for p = 1. Overlaid are the current experi-
mental constraints from PTAs1 and ground-based GW detectors, and future constraints
from spaced-based detectors. PTA sensitivity will not be superseded until the LISA
mission which is scheduled for launch in 2034.
1 While in Blanco-Pillado et al. (2018) the PPTA data gave the most constraining upper limit, the PPTA
result used by Blanco-Pillado et al. did not account for uncertainties in the ephemerides which would likely
weaken the given constraint (Arzoumanian et al. 2018a).
123
Fig. 10 Plot of the gravitational

wave spectrum in terms of the
dimensionless parameter , as a
function of frequency in hertz.
The figure shows cosmic
(super)string spectra for p = 1
for values of the (dimensionless)
string tension Gμ/c2 in the
range of 10−23 –10−9 , as well as
the spectrum produced by
supermassive binary black holes
(SMBBH), along with the
current and future experimental
constraints. The figure is from
Blanco-Pillado et al. (2018)
5 The nature of gravity
Understanding gravity—and testing general relativity as a theory of gravity—

is a leading pursuit of most GW detectors, and PTAs are no exception to this.
Polarization modes: Upon a high-significance detection of gravitational radia-
tion, PTAs will be able to assess its polarization; this is done by correlating the
residuals of pairs of pulsars and constructing an angular correlation function
based on the angle between the pair. This can be measured regardless of the
type of radiation (continuous, memory, background, etc.). By characterizing
the shape of the correlation function, the polarization mode of the signal can
be measured. Various classes of gravity models, including General Relativity,
will be supported or ruled out. PTAs will explore targets other than ground-
and space-based interferometry, permitting competitive and independent con-
straints on theories of gravity.
Graviton mass and gravity group velocity: Variations in the graviton mass
will cause slight variations to the correlation function which may be detected
if PTAs acquire sufficient measurement precision. Graviton mass can also be
constrained or measured by the temporal offset of any co-detected electro-
magnetic counterpart to a single GW source. If the graviton mass is minimal
(m g 10−22 eV), which appears to be based on recent LIGO measure-
ments, PTAs will instead produce precise limits on the group velocity of GWs.
For context, we finish this section by summarizing three examples of theories
that PTAs can test, which make several predictions that differ from general
relativity.
General relativity has been an exceptionally successful physical theory and contin-
ues to stand up to over 100 years of tests of its accuracy (Will 2014). Pulsar timing of
a pulsar–neutron star binary provided the first indirect proof of gravitational radiation
(Taylor and Weisberg 1982), while subsequent studies of the pulsar–pulsar binary have
led to tests of GR to exquisite precision (e. g., Kramer and Stairs 2008). This includes
the discovery of what had been the remaining unobserved prediction of Einstein for
123
general relativity: the detection of GWs from two black holes by LIGO. Over the
years, there have been many other theories of gravity put forward. Some of these are
based on physically esthetic motivations, for instance a change in symmetry or adding
a generalization, such as the scalar field in Brans–Dicke theory (Brans and Dicke
1961). Others, especially in recent years, strive to explain an observed phenomenon
that is unexplained by current physical theories, such as dark matter or dark energy
(Akrami et al. 2013). Whatever the nature of the theory, it must pass all of the observa-
tional tests that ave made general relativity such a successful theory. Myriad new tests
of gravity have become possible with GW observations (e. g., Abbott et al. 2017g;
Chatziioannou et al. 2012; Eardley et al. 1973b; Yunes and Pretorius 2009; Yunes and
Siemens 2013).
5.1 How many gravitational wave polarization states exist?
General relativity predicts the existence of GWs which travel at the speed of light, are
transverse, and have two polarizations, the standard plus and cross. Other theories of
gravity may predict the existence of GWs with different properties.
For any theory in which a metric encodes the dynamics of spacetime, there are
six distinct GW polarizations possible (Eardley et al. 1973b). A linear GW is a small
deviation from flat spacetime with a metric given by gμν = ημν + h μν , where ημν
is the flat-space Minkowski metric and |h μν | 1. Since the metric is a symmetric
4 × 4 matrix, h μν has ten independent components: using our freedom to choose the
coordinate system (i.e., gauge freedom), we can remove four of these components,
leaving only six. These components can be grouped into the way each transforms
under rotations, giving us two scalar components, two vector components, and two
tensor components (Eardley et al. 1973a, b). The tensor components are those most
commonly pursued and are commonly referred to as the plus and cross polarizations.
There has recently been considerable interest in using observations of gravitational
waves to look for evidence of these non-Einsteinian polarizations. For example, Isi et al.
(2015) determined that the signal to noise of match-filtered searches for gravitational
waves with aLIGO can be greatly impacted if they consist of non-Einsteinian polar-
ization modes. Analysis of the detection of GW170814 (the first GW event detected
by the two LIGO observatories and Virgo) favored tensor polarization modes over a
pure vector or scalar modes by a Bayes factor of more than 200 and 1000, respectively
(Abbott et al. 2017e; Isi and Weinstein 2017).
We will work in a coordinate system in which h μ0 = 0, which is called the syn-
chronous gauge. If we consider a plane wave traveling in the z-direction, then a generic
GW is described by six polarization tensors (Eardley et al. 1973a, b):
⎛ ⎞ ⎛ ⎞
1 0 0 0 1 0
ei+j = ⎝ 0 −1 0 ⎠ , ei×j = ⎝ 1 0 0⎠, (19)
0 0 0 0 0 0
⎛ ⎞ ⎛ ⎞
1 0 0 0 0 0
eibj = ⎝ 0 1 0 ⎠ , ei j = ⎝ 0 0 0⎠, (20)
0 0 0 0 0 1
123
(a) (c) (e)
(b) (d) (f)
Fig. 11 The six possible GW polarizations in synchronous gauge where h μ0 = 0. The solid and dotted
line in each case represents the maxima and minima in the strain variations induced by an oscillatory GW.
General relativity predicts only plus and cross modes; however, some theories outlined in Sect. 5.3 give rise
to the other modes. Reproduced from Nishizawa et al. (2009)
⎛ ⎞ ⎛ ⎞
0 0 1 0 0 0
eixj = ⎝ 0 0 ⎠ , ei j = ⎝ 0 1⎠,
y
0 0 (21)
1 0 0 0 1 0
ij
where we have normalized these polarization tensors to satisfy eiPj e P = δ P,P . Three
of these polarization tensors (+/×, and b) are transverse; the other scalar mode and
the two vector components are longitudinal, as shown in Fig. 11.
PTAs are sensitive to the polarization of GWs of any sufficiently bright source,
including single sources, the stochastic background, etc. (Chatziioannou et al. 2012).
For illustration, let us imagine a single, high-intensity plane GW traveling along the
positive z-axis as we observe a PTA with pulsar pairs scattered across the sky, with any
two pulsars separated by some arbitrary sky angle θ . For our single wave with various
polarization modes induced, we show the predicted modulation of the observed pulse
phase at various sky positions in Fig. 12. Here, we can see the quadrupolar response to
the usual plus and cross polarization, the dipolar response to the vector polarizations,
and the monopolar response to the scalar polarizations.
For each pulsar pair, we can measure the correlated response of their timing residuals
and plot this response as a function of a pairs’ angular separation, θ . The shape of this
angular correlation function of pulse residuals, C(θ ), will be different depending on
the type of polarization in the GW being observed. In Fig. 13, we show the correlation
of residuals for a collection of pulsar pairs separated by an angle θ (i.e., the Hellings
and Downs curve, first formulated by Hellings and Downs 1983). We can see that each
of the different polarization states generates a distinct correlation structure which can,
123
Fig. 12 The induced variations in the observed pulse phase as a function of pulsar sky location for a long-
wavelength GW traveling along the positive z-axis. The scalar/vector/tensor and longitudinal/transverse
nature of each polarization mode is apparent
0.5
C(θ )
0.0
0 50 100 150
θ
Fig. 13 The angular correlation functions (Hellings and Downs curves; Sect. 5.1) as a function of the sky
angle between two pulsars. We show the six possible polarizations: plus/cross (solid blue), vector x/y
(dotted red), scalar breathing (dot-dashed yellow), scalar longitudinal (dashed green). Note that there are
only four curves because the plus/cross and vector x/y modes have been polarization averaged, respectively.
See Maggiore (2007) for more details. Lee et al. (2008) found that a PTA can discriminate between the
tensor and non-tensor modes with a PTA with bi-weekly observations of 500 pulsars for 5 years and an
RMS timing accuracy of 100 ns
given sufficient PTA sensitivity, be measured (Gair et al. 2015). As discussed in Lee
et al. (2008), biweekly observations for 5 years with RMS timing accuracy of 100 ns can
detect non-Einsteinian polarization (i.e., other than the + and × modes) with 60 pulsars
for the longitudinal scalar mode and the vector modes; 40 pulsars for the transverse
scalar mode. To discriminate non-Einsteinian modes from Einsteinian modes, the
PTA needs to monitor 40 pulsars for the transverse scalar mode, 100 pulsars for the
longitudinal scalar mode, and 500 pulsars for the vector modes. As such, currently the
detection of such modes is believed to be at least in the intermediate future; currently,
up to ∼10 pulsars have the required residual RMS levels. However, this will change
rapidly with the timing programs beginning on next-generation radio facilities.
As an example, Cornish et al. (2018) established the first PTA upper limits on non-
Einsteinian polarizations. Currently, these limits are derived from the auto-correlation
of the residuals of each pulsar with itself. This analysis assumed a gravitational wave
background produced by a large collection of unresolved binaries and took into account
the fact that dipole and quadrupole radiation will be emitted at different rates from a
123
binary system. Given that the scalar longitudinal mode produces the largest autocor-
relation (see Fig. 13), it has the most stringent PTA upper limit.
5.2 Characterizing the dispersion relation of gravitational waves
Theories which predict these novel polarization states generically also predict non-
standard GW dispersion relations which can result in a number of effects measurable
by PTAs. We can model a range of modified dispersion relations using the parameter-
ization in geometric units (Blas et al. 2016):
ω2 = m 2g + cg2 k 2 , (22)
where m g is the mass of the graviton and cg gives a group velocity different from the
speed of light. Theories of massive gravity (de Rham et al. 2011) allow for a non-zero
graviton mass and some quantum gravity theories predict a modified group velocity
and dispersive effects (Blas and Lim 2015; Liberati 2013).
A number of experiments have already placed limits on the graviton mass. In the
presence of a graviton mass, the Newtonian potential has a length-dependent sup-
pression from which m g can be constrained using a variety of observations. The
most constraining limit comes from weak lensing observations which yield m g <
6.9×10−32 eV and can be translated into a frequency f g = m g c2 /h > 1.7×10−17 Hz
(Choudhury et al. 2004). Recent independent analyses have focused on using galaxy
clusters and the Sunyaev–Zeldovich effect to limit the mass to m g < 1.27 × 10−30 eV
(Gupta and Desai 2018). A massive graviton changes the functional form of the grav-
itational potential from a shape of ∼ 1/r to ∼ e−kmg r /r ; this is the so-called Yukawa
potential. This change in the potential leads to a change in density fluctuations of
matter, which in turn affects the power spectrum of fluctuations measured in weak
gravitational lensing surveys. It should be noted that these constraints are subject to
model-dependent uncertainties in the exact distribution of dark matter, and so should
be used with caution. Other, less model-dependent, constraints from the dynamics of
objects in the solar system yield a limit m g < 4.4 × 10−22 eV or f g > 10−7 Hz
(Talmadge et al. 1988; Will 1998).
With a non-zero graviton mass, GWs acquire a frequency-dependent phase velocity
leading to an additional phase term in GW signals that would be detected by LIGO (Will
1998). Again, fixing cg = 1, constraints from GW170104 provide the best dynamical
constraint to the graviton mass and yield m g < 7.7 × 10−23 eV or f g > 1.8 × 10−8 Hz
(Abbott et al. 2017d).
Pulsars can also be sensitive to the graviton mass and, given that the best dynam-
ical limit to the graviton mass lies in the middle of the PTA frequency band, it is
unsurprising that PTAs may complement the limits set by LIGO. A massive gravi-
ton would be detectable through slight variations in the Hellings and Downs curve,
and an example of this is shown in Fig. 14. As discussed in Lee et al. (2010b), for
a PTA with bi-weekly observations of 60 pulsars for 5 years with an RMS timing
accuracy of 100 ns, massless gravitons can be distinguished from gravitons heavier
than m g = 3 × 10−22 eV with 90% confidence level. A 10-year observation with the
123
Fig. 14 The Hellings and Downs curves for the standard polarization modes (plus/cross), but with different
values of the graviton mass. Reproduced from Lee et al. (2010a). A PTA with bi-weekly observations
of 300 pulsars for 10 years with an RMS timing accuracy of 100 ns can probe graviton masses as low
as m g = 3 × 10−23 eV. This would be competitive with the current LIGO limit and thus could possibly
provide independent confirmation of a LIGO measurement. However, given that we currently only have ∼10
pulsars with RMS accuracy of 100 ns, PTA measurements are not yet competitive with other experiments
for probing graviton mass
same RMS timing accuracy and cadence but with 300 pulsars would probe graviton
masses down to m g = 3 × 10−23 eV.
In addition to using the correlated residuals between pulsars in response to a stochas-
tic GW background, if a PTA detects a single SMBH binary merger which has a
time-tagged electromagnetic counterpart, a comparison between the arrival time of
the two signals allows for a complementary probe of the graviton mass. We discuss
this possibility further in Sect. 9.3.8.
Pulsar timing can also place a limit on the group velocity of GWs. Assuming that the
graviton mass m g 10−22 eV so that f g year −1 , in the PTA band the gravitational
dispersion relation becomes ω2 cg2 k 2 . If cg < 1, the pulse train from a pulsar travels
faster than the group velocity of GWs, thus “surfing” on the wave. These interactions
can lead to phase changes in the waves that significantly amplify the observed timing
residuals. Assuming a reasonable model for the SMBH binary background and taking
the timing residuals for PSR B1937+21, Baskaran et al. (2008) find that 1−cg < 10−2 .
This limit will improve significantly as PTAs become more sensitive.
In other GW bands, the GW group velocity can be estimated through measurements
of the propagation time between the two LIGO detectors. Multi-messenger detections
at Hz-frequencies can also provide constraints (Lombriser and Taylor 2016); an anal-
ysis of GW150914 allowed the first direct constraint on cg giving cg < 1.7 (Blas
et al. 2016). An analysis of the neutron star binary merger GW170817 (Abbott et al.
2017f) with its electromagnetic counterpart placed the most restrictive constraint on
the relative propagation speed of gravitational waves and electromagnetic waves of
−3 × 10−15 < cg − 1 < 7 × 10−16 (Abbott et al. 2017c).
5.3 Examples of gravity theories that predict additional polarizations and

modified dispersion relations
There are several gravity theories which predict additional polarizations and modified
GW dispersion relations. Here, we will briefly describe three of them: Einstein Aether
123
(Jacobson and Mattingly 2001), massive gravity (de Rham et al. 2011), and f (R) grav-
ity (Carroll et al. 2004). We note that these three theories have been chosen for having
a complete and consistent theoretical description and that each face some challenges
when confronted with observations. See Isi and Stein (2018) for a detailed discussion
on how to probe alternative gravity theories with gravitational wave measurements.
We also note that the binary neutron star merger observed by LIGO/VIRGO (Abbott
et al. 2017f) places strong constraints on the properties of alternative theories of grav-
ity (Baker et al. 2017; Boran et al. 2018; Lombriser and Taylor 2016; Sakstein and
Jain 2017).
5.3.1 Einstein Aether
Einstein Aether posits the existence of a Lorentz symmetry-violating, dynamical,

time-like vector field. This is in addition to the gravitational metric, a two-tensor field.
This arrangement preserves three-dimensional rotational symmetry while generating
a preferred rest frame.
The dynamical time-like vector field introduces a preferred frame arising from local
physics leading to a generally covariant theory. Einstein Aether is an example of the
more general vector-tensor theories (Nordtvedt and Will 1972; Will and Nordtvedt
1972) and can be thought of as a simplified model of spontaneous Lorentz symmetry
breaking which may occur in string theory (Kostelecky and Samuel 1989).
Considering a theory that includes all covariant couplings which is quadratic
in all derivatives, the theory can be specified by the value of four new constants,
{C1 , C2 , C3 , C4 } (Jacobson and Mattingly 2004). In the limit where these constants
vanish, this theory is indistinguishable from general relativity. As shown in Jacobson
and Mattingly (2004), there are five GW polarizations in this theory. Transforming to
a gauge where h μ0 = 0 and taking the Ci → 0 limit, the theory predicts: plus/cross
tensor polarizations traveling at the speed of light; two longitudinal vector modes trav-
eling at a speed possibly different from the speed of light; and a linear combination of
the two scalar modes
3
h is j = v0 2 ei j − (C1 + C4 )eibj , (23)
cs
where v0 is the linear perturbation to the time component of the Aether field, cs is the
speed of the scalar GWs, and ei j is as given in Eqs. 19–21.
5.3.2 Massive gravity
Building a well-behaved model of gravity in which the graviton has a non-zero mass
has been an issue for decades; however, recent attempts at understanding the nature of
dark energy have reinvigorated the conversation. In such models, a massive graviton
might introduce a scale that could explain the observed long-range behavior of the
gravitational field. Since the early 1970s it has been known that adding a mass in the
absence of non-linear interactions gives rise to an observationally ruled out discon-
tinuity, the dDVZ discontinuity (van Dam and Veltman 1970; Zakharov 1970). The
addition of non-linear interactions can cure this, but at the cost of a degree of free-
dom. This is known as the Boulware–Deser ghost (Boulware and Deser 1972), proof
123
that only five propagating degrees of freedom can exist in any massive, interacting
extension of general relativity.
Only very recently was there a solution to this issue; the imposition of a symmetry,
the Gallileon symmetry (de Rham et al. 2011; Hassan and Rosen 2012b), restricts the
form of the non-linear interactions and projects out the Boulware–Deser ghost. Later,
it was shown that bimetric theories can be derived from these new massive gravity
theories, and so the theories have been connected (Hassan and Rosen 2012a).
Just as in the case of Einstein–Aether theory, the stable propagating mode is a
superposition of the 6 degrees of freedom identified in Sect. 5.1. Transforming the
metric far away from the source (where the non-linear terms in the field equations
vanish) into synchronous gauge, massive gravity generates the standard plus/cross
tensor polarizations which travel at the speed of light and a scalar wave which travels
at a speed cs with a polarization

π 1
h is j = e + 1 − 2 ei j ,
b
(24)
6Mpl i j cs
which is a linear combination of the breathing and longitudinal modes.
5.3.3 f (R) Gravity
Another commonly considered alternative theory of gravity imagines a gravitational

Lagrangian which is a function of the Ricci curvature scalar, f (R) where f (R) = R
gives the Lagrangian for GR. The original motivation for f (R) gravity was to provide
an alternative gravity theory that could account for the observed accelerated expansion
of the universe (Carroll et al. 2004). In its original formulation, the gravitational
Lagrangian was modified to include a term inversely proportional to the Ricci scalar:
R+μ/R. As the universe expands, the homogeneous value of the Ricci scalar decreases
until, at a late enough time, this new term becomes dynamically important. Although
it was shown that this specific modification was at odds with measurements of the
deflection of light around the Sun (Chiba 2003; Erickcek et al. 2006), there are specific
forms of the function f (R) which pass the solar system tests and lead to a late period
of accelerated expansion (Chiba et al. 2007; Hu and Sawicki 2007).
f (R) gravity is an example of a scalar–tensor theory (see, e.g., Will 1993): by
performing a conformal transformation, the theory can be described as a scalar–tensor
theory with a scalar mass

1 1 R0
m = 2
− , (25)
3 f (R0 ) f (R0 )
where R0 is the Ricci scalar for the background over which the GWs propagate.
An analysis of the polarization modes shows that this theory predicts the standard
plus/cross tensor polarizations propagating at the speed of light and, again, a linear
combination of the breathing and longitudinal scalar modes (Eqs. 20 and 21; Liang
et al. 2017)
123

1 m2
h is j = − f (R0 ) eibj + 2 ei j , (26)
2 ω
and follow the dispersion relation m 2 = ω2 − k 2 .
6 Relic gravitational waves and early-universe cosmology
Leading theories for the evolution of the early universe invoke an inflationary
epoch to explain the isotropy and other broadly observed properties of the
universe.
If this epoch occurred, it would have produced GWs by rapidly expand-

ing quantum fluctuations present in the pre-inflation epoch. The nature of
inflationary expansion will be encoded in the strength and spectrum of
the produced broad-band GW background. For standard inflation models,
the inflationary background will likely be fainter than that of SMBHBs;
however, some inflationary modes may have a spectrum that rises into the
PTA/LISA/LIGO bands. PTAs have already begun to place stringent limits
on those models. PTAs will most sharply probe the spectrum (scalar spectral
index), while cosmic microwave background (CMB) experiments will probe
the relative strength of these “relic” GWs and the more standard density waves
(the latter have already been detected and characterized by CMB probes).
If primordial GWs can be detected, PTAs may also be sensitive to phase transi-
tions in the early universe, allowing (among other possibilities) an independent
constraint on the existence of the cosmological constant.
6.1 Relic gravitational waves
GW signals can arise from the early universe if a period of rapid inflation occurred
(Fig. 15). Quantum fluctuations from early in the universe would have been amplified
by inflation, and, like strings, would produce a broad-band signal detectable by multiple
instruments. These “relic” GWs are a long-standing target for current and future CMB
experiments, looking for the tensor modes induced by these waves ( Kamionkowski and
Kovetz 2016, and references therein). CMB experiments will constrain only the long-
wavelength (low-frequency) portion of the spectrum of inflationary GW signals, and
can most effectively constrain the scalar-to-tensor ratio. However, higher-frequency
experiments, like PTAs and space- and ground-based laser interferometry, are able
to more effectively constrain the spectrum of these early fluctuations by looking
at the scalar spectral index (the spectral index of the detected background, which
reflects the spectrum of the fluctuation scales and the mode of their amplification).
Various models of inflation predict differing values for the observed scalar spectral
index.
123
Fig. 15 The quantum

fluctuations in the very early
universe, if rapidly amplified by
inflation, could cause a highly
broad-band GW background
signal in addition to affecting the
cosmic microwave background.
While cosmic microwave
background experiments can
constrain the scalar-to-tensor
ratio, higher-frequency
experiments (as PTAs, LISA,
LIGO/VIRGO shown here) will
most sharply probe the spectrum
of these early fluctuations.
Figure reproduced from (Battye
and Shellard 1996)
Only weak—and often highly model-dependent—constraints exist on the shape

of the inflationary GW spectrum. Following the recognition that primordial black
holes could contribute to the LIGO detections, there has been a renewed interest in
examining what kind of constraints could be placed on the spectrum of inflationary
perturbations at shorter wavelengths, where PTAs can contribute broad-band limits.
In fact, the most stringent limit on relic GWs, before LIGO’s first observation, had
been set at gw ( f ) < 2.3 × 10−10 by the combination of PPTA, LIGO, and CMB
limits (Lasky et al. 2016). This marked the first time the most optimistic predictions
of various models were impinged upon by observations.
Given the rapid progress across the GW spectrum (from ground-based interferom-
eters to PTAs to CMB experiments), our assessment is that this area is likely to remain
observationally driven. That is, while models can always be constructed to evade the
observational limits, the extent to which a blue inflationary GW spectrum (i.e., one
where the energy density rises with GW frequency) remains viable will be determined
by the observations across the GW spectrum.
While we can currently place constraints on inflationary GWs, we state with some
confidence that PTAs will not be the primary driver for the study of inflationary GWs.
This is because it is highly unlikely that the inflationary epoch signal will dominate
over the much brighter GW background from inspiraling SMBHBs. For instance, the
most widely accepted “slow-roll inflation” model has an expected gw of 10−15 , while
the most pessimistic backgrounds from SMBHBs are expected to be around two orders
of magnitude above that level (Lasky et al. 2016).
Certain axion inflation models also predict that a large bump in the GW background
spectrum may be produced by the production of stellar-mass primordial BHs in the
early universe (e. g., Namba et al. 2016). If the primordial BHs are produced within
a narrow mass range with a peak at a mass of in the order of a few to 100 M , this
GW background bump may peak in the pulsar timing band, significantly enhancing the
inflationary signal in a limited gravitational bandwidth that is defined by the primordial
BH mass distribution (García-Bellido et al. 2017).
123
6.2 Cosmological measurements and the cosmological constant
As the early universe cooled, successive constituents decoupled and began to evolve
(more or less) independently. The most notable such example is the decoupling of
photons, which led to the formation of the CMB. Prior to photon decoupling, it is
expected that neutrinos underwent the same process, leading to a relic neutrino back-
ground. Lattanzi et al. (2010) and Benini et al. (2011) note that cosmic neutrino
decoupling should have happened a few seconds after the Big Bang, at a redshift
z ∼ 1010 . A GW entering the horizon at this time currently has frequency of order
1 nHz (e.g., see Fig. 15), suggesting that the neutrino decoupling time can be probed
by nanohertz GWs. Specifically, they predict that the effective viscosity from the
neutrinos will suppress such GWs. They also conclude that a multi-decadal data
set would be required to detect this effect. With current PTAs beginning to sur-
pass the decadal observational duration, it seems that this effect still remains out
of reach.
The early universe may have also experienced phase transitions as it cooled, likely
before the epoch of neutrino decoupling. Caprini et al. (2010) consider a quantum
chromodynamics phase transition (which they take to occur at z > 1017 ) and the
GWs produced during that epoch. Such a phase transition could produce GWs with
f ∼ 1 Hz. Their assessment was that the NANOGrav sensitivity, at the time, was
insufficient to detect GWs from this phase transition, but predicted that PTAs might
be able to detect these GWs on the timescale of the year 2020. However, their prediction
was based on an idealized PTA (with higher precision timing residuals, and on a smaller
number of pulsars than currently timed). For a realistic estimate, this calculation would
need to be reformulated to encompass the actual sensitivity of current PTAs, i.e., taking
into account the true number of pulsars being timed by the various PTAs, and realistic
timing residuals.
The standard cosmological model includes a “dark energy” component (Planck
Collaboration et al. 2016; and references therein), which may be equivalent to
the “cosmological constant.” The expectation is that dark energy becomes relevant
on large scales, and significant international efforts are devoted to placing con-
straints on the properties of dark energy. Bernabeu et al. (2011) and Espriu and
Puigdomènech (2013) describe an approach in which measurements of pulsar tim-
ing residuals due to GWs emitted in the nearby universe (z 1) could place a
complementary constraint on . One (acknowledged) caution about their results is
that their estimates are based on a somewhat idealized PTA with a relatively short
observation duration (3 year). Given that current PTAs have durations in excess of
a decade, and substantially more pulsars, it is likely that different (and probably
more strigent) constraints could be obtained by repeating their analysis for current
PTAs.
123
7 Additional science from “Hidden Planets” to dark matter
There are numerous scientific endeavors that may be attainable with PTAs,
including those for which: (1) our theoretical understanding is still under active
development or is speculative in nature, or (2) the pursuit/detection of the
science is likely much further in the future than the topics described in previous
sections. This includes GW studies of stellar convection and dark matter, as
well as the effect of primordial black holes to produce a GW background or
direct gravitational perturbations on Earth or the pulsar, or direct perturbations
on Earth caused by our own solar system bodies.
The discussion here is structured in terms of specific sources and predictions.
We also recommend reading Cutler et al. (2014), which considers the blind
discovery space of PTAs. That work showed that, whatever the source, GW
memory at high redshift is a potential discovery area.
7.1 Stellar convection
Although the common assumption for nanohertz GWs is that they are produced by
extragalactic or cosmological sources, Bennett and Melatos (2014) evaluate the mass
quadrupole produced by turbulence within a convective star, such as the Sun. They
show that the ensemble of stars should produce an irreducible background. While their
primary focus is LISA, they note that the amplitude of the GW signal is amplified in the
near-field zone, which is certainly the case for the Earth term for frequencies relevant
to PTAs. A simple extrapolation of their results into the nanohertz GW frequency
band suggests that the Sun may be a contributor to the PTA noise budget, and thus
potentially detectable. Such an extrapolation relies upon assumptions regarding the
nature of turbulence within a star, and that the actual magnitude of any GW signal
may be substantially lower; thus, with increasingly stringent limits, PTAs may place
constraints on solar and stellar convection.
7.2 Dark matter
Several authors have considered various classes of dark matter candidates that may
produce observable signatures in PTA data.
7.2.1 Cold dark matter
The concordance -Cold-Dark-Matter (CDM) model for the evolution of the uni-
verse predicts that the mass function of dark matter halos may extend to very small
masses ( M ), with the mass function cutoff related to the nature of dark matter.
Thus, detection of small-scale dark matter clumps can both provide support for the
cold dark matter scenario and constrain the mass of dark matter particles. There are two
possible influences on pulsar timing (Baghram et al. 2011; Dror et al. 2019; Ishiyama
et al. 2010; Kashiyama and Oguri 2018; Siegel et al. 2007): (i) Shapiro delay of the
radio pulses propagating through a distribution of dark-matter sub-structure (i.e., the
123
integrated Sachs–Wolfe effect), and (ii) Doppler delay due to an acceleration of the
Earth or pulsar caused by the nearby passage of a dark matter clump.
The Shapiro effect is line-of-sight dependent, leading to delays in pulsars that are
related only by the statistical properties of dark matter sub-structure. This is expected
to be challenging to disentangle from intrinsic pulsar noise processes, but could give
access to clumps in the mass range 10−4 − 10−3 M . In the Doppler delay, a clump
passing by a pulsar produces a timing-residual influence only in that pulsar, whereas
a clump passing by Earth produces a dipolar-correlated timing delay in all pulsars in
the array (similar to unmodeled solar system ephemeris systematics; Baghram et al.
2011). The Doppler delay effect is expected to dominate over the Shapiro delay from
sub-structure (Dror et al. 2019), potentially constraining clumps in the mass range
10−9 − 10−8 M .
7.2.2 Scalar-field dark matter
Scalar-field dark matter models (sometimes referred to as “fuzzy dark matter“) address
some of the issues that CDM has in reproducing the observed number density of
sub-galactic scale structures in the universe (Hu et al. 2000; Hui et al. 2017). Structure
formation is suppressed below the de Broglie wavelength,

10−23 eV 10−3
λd B ≈ 600 pc , (27)
m v
where m is the mass of the dark matter particles and v is their characteristic velocity
in units of c. Such a population of particles will produce an oscillating pressure that,
though averaging to zero, will cause sinusoidal variations in the local gravitational
potential with a frequency
m
f ≈ 5 · 10−9 Hz (28)
10−23 eV
and an amplitude
Gρ(x)
(x) = π , (29)
m2
where ρ(x) is the density of dark matter particles at position x. While it is not a GW
effect, pulsar signals propagating through such a time-variable gravitational potential
will have their pulsation frequencies sinusoidally modulated. The perturbation to the
observed times of arrival could be as large as several hundred nanoseconds, which
would be a signal that is accessible to a number of currently timed pulsars (Khmelnitsky
and Rubakov 2014).
Porayko and Postnov (2014) searched for this signature of scalar-field dark matter
and placed the first observational constraints on the mass of such particles using the
NANOGrav 5-Year Data Release (Demorest et al. 2013). Their upper limits on the
local oscillating gravitational potential were an order of magnitude above the current
predictions. In follow-up work, Porayko et al. (2018) used 26 PPTA pulsars observed
123
over 12 years to compute rigorous Bayesian and frequentist constraints on the local
density of ultralight scalar-field dark matter. They found that ρSF ≤ 6 Gev cm−3
at 95% confidence for ultralight bosons with m ≤ 10−23 eV. This improves upon
previous constraints by a factor of 2–5.
The ability to measure this oscillating gravitational potential improves as more
pulsars are added to a PTA. The current IPTA contains more than twice the number
of pulsars than were used by Porayko et al. (2018), and more continue to be added
to the constituent PTA programs on an annual basis. The prospects for near-future
constraints from the IPTA on ultralight scalar-field dark matter are thus very exciting.
7.3 Primordial black holes
The first detection of GWs from merging stellar-mass BHs by LIGO/Virgo ( Abbott
et al. 2016c), as well as the subsequent mergers detected, raised the question of the
formation channel that produced these BHs. While a full discussion of this topic is
beyond the scope of this paper, it has been proposed that these BHs are primordial,
produced during the inflationary period in the early universe. As previously discussed
in the last paragraph of Sect. 6.1, some of these inflation models predict solar mass
and larger primordial BHs that may be produced as a by-product of inflation. Inflation
itself, in addition to the production of primordial BHs, contributes to a potential GW
background signal (e. g., Cheng et al. 2017; García-Bellido et al. 2016). The peak fre-
quency of the bump that can be produced in the GW background spectrum depends on
the primordial BH mass, and the 10–100 M range precisely maps into the nanohertz
GW band (García-Bellido et al. 2017). Thus, PTA searches of the GW background
also serve as probes of primordial BH mass distributions in the critical range that
LIGO detects.
Much smaller primordial BHs have also been proposed as relevant to PTAs, by
Seto and Cooray (2007) and Kashiyama and Seto (2012), due to the perturbations
they can cause when passing close to the Earth. In contrast to the BHs responsible for
the LIGO events, the primordial BHs considered by these authors have much smaller
masses (< 10−10 M ). Nonetheless, from the perspective of an observationally driven
constraint on primordial BHs, these constraints remain of potential interest. These
authors show that a primordial BH may produce perturbations of order 20 ns over a
duration of the order 15 years. With the various PTA data sets now exceeding a decade,
this signal may become feasible to detect in the next decade, if other contributions
to pulsar timing residuals at similar levels (tens of nanoseconds) can be sufficiently
controlled or modeled.
7.4 Solar system ephemerides and wandering planets
A fundamental term in pulsar timing is an astrometric term, designed to transfer the

pulsar TOAs into the frame of the solar system barycenter, which is assumed to be
(quasi-)inertial. The time delay between the Earth and barycentric reference frames is
given by (Lorimer and Kramer 2012):
123
RSSB · n̂
tSSB = , (30)
c
where RSSB is the vector connecting the Earth to the solar system barycenter, n̂ is the
unit vector in the direction of the pulsar, and c is the speed of light. The barycenter
position is the center of mass of the solar system, involving weighted contributions
from all planets and important dynamical objects. Any uncertainties in the position or
masses of solar system bodies will create a corresponding uncertainty in the barycenter
position.
In Eq. 30, it is clear that uncertainties in the knowledge of the position of the
solar system barycenter will affect the accuracy to which pulsars can be timed. As the
position of the solar system barycenter RSSB is determined from the orbits and masses
of the solar system planets (and minor bodies), there is a long history of assessing
whether, and to what extent, precision pulsar timing can be used to constrain properties
of the solar system (e.g., Li et al. 2016; Mulholland 1971). These efforts generally
find that the position of the solar system barycenter is not known to better than about
100 m, a value that is consistent with expectations from the spacecraft data used to
construct the solar system ephemeris (W. M. Folkner 2017, private communication).
There have also been efforts using PPTA and IPTA data to measure the masses of
the solar system planets by using precision pulsar timing (Champion et al. 2010 and
Caballero et al. 2018, respectively).
Recent efforts to constrain the nanohertz stochastic GW background with the
NANOGrav 11-year data set uncovered differences in upper limits and detection statis-
tics caused by systematic errors in published solar system ephemerides (Arzoumanian
et al. 2018a). Upper limit variations were small, but Bayes factors for a long timescale
red noise process shared by all pulsars (a potential first signature of GWs) varied by
over an order of magnitude, from compelling to insignificant evidence. To mitigate
these systematic errors, Arzoumanian et al. (2018a) added gas-giant mass perturba-
tions and Jupiter orbital-element perturbations to the global PTA signal and noise
model. This led to consistent constraints and detection statistics, regardless of the
assumed baseline ephemeris version. Work is ongoing to expand this model into a full
Bayesian solar system ephemeris.
There is also a long history of assessing whether precision pulsar timing could
reveal currently unknown members of the solar system or a distant companion star to
the Sun (e.g., Guo et al. 2018; Harrison 1977; Zakamska and Tremaine 2005). Though
some of the early efforts may have been affected by severe selection effects (Henrichs
and Staller 1978), the existence of distant members of the solar system is a question
that has recently attracted considerably more attention, given the suggestion of a
previously unknown “Planet Nine” (Batygin and Brown 2016). While we are unaware
of any large-scale effort to constrain the existence of this body with the modern PTAs,
Guo et al. (2018) estimate that PTAs could constrain the mass to be 10−4 M ,
and Zakamska and Tremaine (2005) show that measuring the acceleration of the solar
system’s barycenter with precision pulsar timing may be highly competitive with other
methods for constraining the existence of distant companions (> 300 au).
123
8 Gravitational-wave synergies: multi-band GW science
PTAs have some target science in common with other GW detectors at higher
frequencies (LIGO, LISA) and lower frequencies (CMB experiments). In some
cases, this allows for synergistic or even direct multi-band GW studies.
PTAs and higher-frequency millihertz experiments like LISA uniquely share
SMBHBs as a direct common target. PTAs and LISA will make synergistic
measurements of galaxy and SMBHB evolution. However, LISA will probe
lower-mass and higher-redshift systems, while PTAs will preferentially detect
the largest systems in the nearest few Gpc. These experiments will work
together to fully characterize black hole growth over cosmic time. For a small
subset of discrete systems, both experiments may co-detect the selfsame sys-
tem at different phases of evolution.
Cosmic strings and relic GW signals are also intrinsically broad-band GW
emitters, and where relevant we have summarized the links between exper-
iments sensitive to those phenomena in Sects. 4 and 6.1, respectively.
Complementary constraints from multiple wave bands on the nature of gravity
are described in Sect. 5.
The LIGO/Virgo GW detections (Abbott et al. 2016a, b, 2017a) in the high-

frequency GW regime (Hz) have ushered us into the era of GW astrophysics. In
addition to PTAs, there are many other instruments and techniques currently in the
design, planning and prototyping phases for future GW observatories. For instance,
LISA is a planned space-based GW detector, sensitive to intermediate frequencies (∼
mHz) expected to be launched in the 2030s (Amaro-Seoane et al. 2012).
SMBHBs are among the primary targets of both PTAs and LISA. However, the
two experiments probe different stages of SMBHB evolution and they are sensitive
to SMBHBs in different mass ranges. PTAs are most sensitive to the early inspiral
(orbital periods of years or longer) of nearby sources with M ∼ 109 M (Mingarelli
et al. 2017). In contrast, LISA is sensitive to the inspiral, merger, and ringdown of
SMBHBs with masses from 104 − 107 M at a wide range of redshifts (Amaro-
Seoane et al. 2012). The two populations of SMBHBs probed by PTAs and LISA
are linked via the growth and evolution of SMBH across cosmic time, as shown in
Fig. 16. Given that the same fundamental physics is required to produce and evolve
both populations of BH binaries, there exists a strong link between the methodology
of evolutionary models used to study them, and the insights that observations of their
GW signatures will provide. In particular, PTA observations of the GW background,
and measurement of its spectral index, will provide valuable constraints on the mass
function and eccentricity distribution of SMBH (e.g., Kelley et al. 2017b) which will
better constrain the detection rates for LISA.
Under the right circumstances, an individual source could be observed by PTAs
and LISA at different stages of its evolution (Pitkin et al. 2008; Spallicci 2013). For
a PTA with high-frequency cutoff of 4 × 10−7 Hz, an SMBHB will transition from
being observable by PTAs to being observable by LISA in 1–50 years. This transition
time can be reduced by increasing the pulsar observing cadence, which improves the
sensitivity of PTAs at high frequencies. However, astrophysical rate estimates suggest
123
Fig. 16 LIGO, LISA, and PTAs have complementary coverage to study the full range of black hole
masses at various stages of the Universe. Here we show the approximate signal-to-noise ratio for the
complementary wavebands of these three instruments as they are currently (darker shading/black contours)
and in the early- to mid-2030’s era (lighter shading). This plot focuses only on individual (rather than
stochastic) black hole detections. All curves assume instrument-limited sensitivity, without an astrophysical
background. Individual inspiral/coalescence events at high redshift will be detectable by LISA, while
systems in the extended inspiral phase at higher masses and lower redshift are detectable by PTAs as
continuous gravitational waves. The source classes of LISA and PTAs are particularly linked through the
evolution of MBHs across cosmic time. Understanding the growth of MBHs will require the contributions
of both PTA and LISA data. Figure produced by Andrew Kaiser and Sean McWilliams (WVU); a more
rigorous version will be published in Kaiser & McWilliams (in prep)
the probability of a sequential detection being extremely low (4.7 × 10−4 – 3.3 × 10−6
per year to merger per year of survey) due to the small number of individual sources
observable by both detectors (Spallicci 2013).
It is also possible to use ringdown observations made by LISA as triggers to search
PTA data for past continuous waves, or for memory-inducing SMBHB coalescence
events (Sect. 3.1.3). LISA can observe the ringdown of higher mass sources, even
when the inspiral and merger happen outside of the LISA band, meaning there is better
overlap for direct observations of these sources with both PTAs and LISA. Parameter
constraints from observing the ringdown can be used to improve the search for the
inspiral, extrapolating a SMBHB model back in time to predict the expected gravita-
tional waveform throughout the previous years of pulsar observations. Currently, the
planned launch date for LISA is 2034, at which point PTAs will have accumulated
over 30 years of data that can be used for such a search.
Galactic sources of GWs that LISA will be able to study may, under certain cir-
cumstances, be possible to investigate with pulsar timing. Globular clusters (GCs)
likely host GW sources detectable by LISA (Kremer et al. 2018). GCs are also known
to host large populations of pulsars (Freire et al. 2017; Ransom et al. 2005). Pulsars
in a GC will be within a few parsecs of GW sources in that cluster and could act as
sensitive probes of those GWs (Jenet et al. 2005; Madison et al. 2017). Pulsars in GCs
have some limitations in sensitivity due to accelerations from intra-cluster dynamics,
which tend to cause low-frequency structure in the timing residuals of those pulsars,
123
possibly masking or mimicking low-frequency GWs. Still, several GC pulsars are

currently timed as part of the PPTA and EPTA, and could serve as probes of Galactic
LISA targets.
Alternatively, a pulsar behind a GC, significantly more distant than the cluster,
would be extremely sensitive to GWs from sources near its line of sight. However,
such an arrangement is unlikely. In the absence of multi-messenger information, sev-
eral such pulsars would also be necessary to distinguish GWs from unmodeled effects
in an individual pulsar. Additionally, for pulsars within or behind a GC, the con-
ceptually straightforward division of the timing signal into an Earth term and pulsar
term does not apply because the GWs cannot be approximated as plane waves. This
fact will require PTAs to adopt currently undeveloped detection techniques. Finally,
several pulsar timing GW limits at high GW frequencies (10−6 –10−3 Hz) over-
lapping with the LISA band have been published (Dolch et al. 2016; Perera et al.
2018; Yardley et al. 2010; Yi et al. 2014). This is possible because of high-cadence
observations of selected pulsars. These limits are many orders of magnitude above
the Doppler tracking GW limits from the Cassini spacecraft and the expected GW
sensitivity of LISA (Armstrong et al. 2003). However, both Cassini and LISA are
entirely solar system bound; one possible advantage of pulsar timing in the μHz
to mHz band is its sensitivity to a galactic GW source close to a pulsar line of
sight. Studying galactic sources is an area of investigation that needs to be further
developed by the PTA community, but it could yield exciting overlap with LISA
science.
Recently, there has also been work on the possibility of astrometric GW detec-
tion (Book and Flanagan 2011; Schutz 2010), for example using an instrument like
GAIA that can accurately measure the positions of over a billion stars in the Milky
Way (Gaia Collaboration et al. 2016). This technique uses correlated deviations in
the observed positions of astrophysical objects (most likely stars, or possibly quasars)
to measure passing GWs. PTA use a similar procedure based on time delays which
are determined by integrating a GW induced redshift over the path of each photon.
This integration introduces a frequency dependence in which the characteristic strain
sensitivity of the array is linearly related to the GW frequency, and thus is worse at
high frequencies. Astrometric detection, on the other hand, performs an effectively
instantaneous measurement of each photon’s source direction. Within the Nyquist
band, determined by the observational cadence and duration of the observatory, astro-
metric detection has a sensitivity relatively independent of frequency. This detection
method is still in its early days of study; however, astrometric GW detection may
be able to complement PTAs at the high-frequency end of the PTA sensitivity range
( f 3 × 10−8 Hz; Moore et al. 2017) and could potentially be used in a coherent
multi-messenger PTA+Gaia analysis to validate the detection of nanohertz GWs, or
even constrain beyond-GR GW polarizations (Mihaylov et al. 2018; O’Beirne and
Cornish 2018).
123
9 Electromagnetic synergies: multi-messenger science
Almost all of the science cases described in the previous sections can be signif-
icantly augmented with electromagnetic studies. However, for the PTA GW
band, the science that would benefit most clearly from a concerted multi-
messenger and multi-wavelength effort is the identification of the inspiral
and/or coalescence from discrete SMBHBs. Here, we outline a selection of
studies related to SMBHs that can be uniquely achieved with a host galaxy
identification and the direct measurement and tracking of an SMBHB elec-
tromagnetic counterpart (note, the latter implies the host galaxy has been
identified).
We also summarize a few practical aspects of multi-messenger detection in
the nanohertz waveband, which may include both new and archival data.
9.1 Identifying and tracking an SMBHB host
To identify the host of a discrete (continuous wave/memory/burst) GW source, we

need to be able to localize the object on the sky. The size of a PTA’s localization error
depends on the signal-to-noise ratio of the GW detection, and the relative position of
pulsars with respect to the GW source (e. g., Ellis 2013). However, typical localizations
for at least the near future (detection significance of ∼ 7σ ) will be up to a few times
1000 deg2 (Fig. 17). Although the chirp mass of and distance to a SMBH binary are
covariant in a PTA detection (Eq. 8), we can estimate an upper limit on the distance
by assuming a maximum chirp mass (∼ 1010 M ). Still, the large number of resulting
galaxies in this error region will require us to turn to electromagnetic information to
identify the most likely host.
There are two main goals of multi-messenger astronomy in the nanohertz regime:
(1) We aim to simply identify the likely galactic host of the GW signal, and deter-
mine its distance. This information will allow us to determine the GW source’s chirp
mass, by breaking the chirp-mass/luminosity-distance degeneracy, mentioned above.
(2) We aim to perform true “multi-messenger” monitoring of a target, by detecting
electromagnetic signals from the vicinity of a SMBH binary, during the inspiral or
coalescence. This will be possible, if the SMBHs are illuminated as AGN, or if the
binary leaves a distinct observable signature on the stellar or gas dynamics, the pho-
tometry, or the morphology of their host. We summarize the methods (and efforts) to
electromagnetically identify SMBHBs in Sect. 9.2.
SMBHBs are formed in galaxy mergers, which for much of their duration are
ostentatious events; they exhibit large-scale asymmetries, tidal tails, sudden bursts of
star-formation (e.g., ULIRG stage), and potentially an abundance of tidal disruption
events (Burke-Spolaor 2013). However, since, as discussed extensively in Sect. 3, the
efficiency of SMBHB formation and inspiral is highly uncertain, it is also unclear
whether these easily observable features would be present by the time a binary SMBH
enters the GW-dominated phase. If the final parsec problem is resolved efficiently, it
123
Fig. 17 Representative localization and parameter measurements for a continuous-wave detection. These
images show the marginalized 2-D posterior probability density functions in the sky coordinates (θ, φ)
and the log of the chirp mass and distance for injected signal-to-noise ratios of 7, 14, and 20 shown from
top to bottom. The × symbol indicates the injected parameters and the solid, dashed and dot-dashed lines
represent the 1, 2, and 3 sigma credible regions, respectively. Figure from Ellis (2013)
is possible that these large-scale indicators of a galaxy merger persist, until the binary
reaches the inspiral phase. Therefore, we can easily survey for this type of signature
in the GW error region, and potentially narrow down the list of candidate hosts.
However, if the archival data prove inadequate, new dedicated observing campaigns
will also be required to systematically catalog and identify outstanding features of
galaxies in the error region of PTAs.
Before we move on to discuss direct SMBHB detection and multi-messenger sci-
ence, it is worth noting an important practical point about electromagnetic SMBHB
emission. The evolution of SMBHBs within the PTA band is slow, and the inpsiral
phase may last for long periods of time (of order years to millennia). Addition-
ally, SMBHBs at these stages may have long-enduring electromagnetic signatures.
Therefore, search efforts for multi-messenger signals in the nanohertz regime have a
tremendous advantage (e.g., compared to LIGOs electromagnetic counterparts); for
the majority of objects, tracking of a GW target will not necessarily be time critical.
Electromagnetic data collected years before the GW detection might reveal key infor-
mation for a binary SMBH, its host, and its evolution. Therefore, the vast availability
of archival images, and spectra, along with time series from the ongoing time-domain
surveys, as well as data from the emergent highly sensitive instruments will allow
great potential for discovery. They may identify any electromagnetic signatures on the
large scale (i.e., unique identifiers in the properties of the host galaxy), or they may
allow the tracking of any multi-messenger counterparts that mark the GW-emitting
binary itself.
123
9.2 Direct searches for SMBHBs in electromagnetic bands
Given the significance of SMBHBs in galaxy evolution and the prospects for GW
astrophysics, over the past few decades, several groups have searched for compact
sub-pc binaries using electromagnetic data across the electromagnetic spectrum. In
principle, direct imaging of material near the SMBHB would be the most direct identi-
fication method. However, only radio interferometry can provide sufficient resolution
to probe sub-parsec scales (Burke-Spolaor 2011; Kharb et al. 2017); although this
provides a convenient observation method to test some binary SMBHB candidates, it
can only reach those up to moderate redshifts at the earlier stages of binary evolution.
Therefore, many recent searches have relied on the effects of the binary on its
environment (e.g., Doppler-shifted broad emission lines, periodic variability, distorted
morphology of radio jets, etc; see Burke-Spolaor 2013 for a discussion of multi-
wavelength counterparts to SMBHBs). However, since these signatures are indirect,
other physical processes can also explain the observed features.
For instance, in the standard picture of AGN, the broad emission lines arise in
regions close to the central SMBH, whereas the narrow emission lines are associated
with large-scale features of the host galaxy. In the presence of a compact binary, the
orbital motion of the SMBHB will be imprinted in the broad emission lines, producing
noticeable shifts compared to the narrow lines. Alternatively, if both SMBHs have their
own broad emission line region, the broad lines should be double peaked, reflecting
the motion of the binary. Several candidates have been identified, especially from sys-
tematic searches in the spectroscopic database of SDSS (Eracleous et al. 2012; Ju et al.
2013). However, similar features can be produced by inflows/outflows, the geometry
of the broad emission line region, as well as by a special type of AGN, the so-called
double-peaked emitters. Therefore, only long-term monitoring can prove/disprove the
binary nature of these candidates, if coherent changes in the spectra are observed, as
expected due to the orbit of a binary.
Periodic variability in AGN and quasars can also signify the presence of SMB-
HBs. The binary is expected to periodically perturb the circumbinary disk and several
hydrodynamical simulations have found that SMBHBs can produce bright quasar-like
luminosities, with the mass accretion rate onto the BHs periodically modulated at the
orbital period of the binary. Several candidates (∼150) have emerged from system-
atic searches in analysis of time-resolved photometry (Charisi et al. 2016; Graham
et al. 2015b, a; Liu et al. 2015; Zheng et al. 2016). However, the identified periods
are relatively long, and only a few cycles have been observed. This, combined with
the stochastic underlying variability of quasars, can lead to false detections (Vaughan
et al. 2016). Again long-term monitoring can test whether the detected periodicity
is persistent and thus attributed to a SMBHB. For instance, follow-up observations
showed that the periodicity of quasar PSO J334.2028+01.4075 is not persistent (Liu
et al. 2016), whereas for quasar PG1302-102, the periodic model is preferred com-
pared to pure stochastic noise (Liu et al. 2018). Additionally, several independent
signatures have been proposed, which, if observed along with the periodicity, will
increase our confidence in the binary nature of the candidates. These include multi-
wavelength observations of Doppler-boosted emission (Charisi et al. 2018; D’Orazio
123
and Haiman 2017; D’Orazio et al. 2015b), periodicity with a characteristic frequency
pattern (Charisi et al. 2015; D’Orazio et al. 2015a), and periodic self-lensing flares
(D’Orazio and Di Stefano 2018).
Last but not least, the orbit of the binary can affect the morphology of radio jets. If
one of the BHs is emitting a radio jet, the orbital motion of the base of the jet will result
in a jet with helical structure (Kaastra and Roos 1992). Furthermore, the presence of
a secondary BH can lead to jet precession, producing a jet with conical or helical
morphology (e. g., Britzen et al. 2017). As before, these signatures are not unique to
SMBHBs, as instabilities in the jet or the misalagnment of the accretion disk with the
jet axis can also produce similar features. Independent lines of evidence are required
to test whether these galaxies host SMBHBs. We note that the above signatures, albeit
indirect, will be significant in searches for the host galaxy, once PTAs detect a GW
source.
9.3 Topics in binary supermassive black hole multi-messenger science
9.3.1 Accretion dynamics: active nucleus geometry
In gas-rich mergers, copious amounts of gas will be funneled into the central regions of
the post-merger galaxy. Once the two SMBHs become gravitationally bound, the gas
will settle in a circumbinary accretion disk. The disk can dissipate the SMBHB energy
and angular momentum, driving the binary to the GW regime. Since the existence of
ambient gas can catalyze the binary evolution, this has been suggested as one solution
to the final parsec problem (Sect. 3). Aditionally, torques from a circumbinary disk
could influence the inspiral evolution of the binary significantly enough for this effect
to be visible in the gravitational waveform. Theoretical work on circumbinary disks
covers a variety of topics, including the timescales of SMBHB evolution, whether
the disks will be retrograde or prograde with respect to the binary orbit (Schnittman
and Krolik 2015), whether such disks will support accretion onto one SMBH or both
SMBHs (Cuadra et al. 2009), the effect of the gaseous disk on the binary eccentricity,
etc. (see also Sections 3.1.2 and 3.2.4).
More importantly, the presence of gas is crucial in providing bright electromagnetic
counterparts, which will allow us to pinpoint the tentative SMBHBs. This is an area of
intense investigation, with several remaining open questions regarding the expected
emission signatures. For instance:
• Will the binary have one or two radio jets?
• Will the circumbinary disk itself support a stable or unstable jet?
• Will a secondary black hole clear a gap or a cavity (Hayasaki et al. 2007), and if
so, can we use this as an observable to infer its mass ratio?
The GW signal will provide an independent measurement of the chirp mass Mc
(provided that the distance is constrained electromagnetically from the redshift of the
source). It will also allow for the measurement of the separation of the two black holes
and time-tagged tracking of the orbital motion of the binary. These, in combination
with any observed electromagnetic signatures, will allow us to explore in detail the
interplay between the disk and the binary orbit, the accretion rate onto the individual
123
BHs and the overall disk geometry, as well as the processes involved in jet generation,
thus providing unprecedented insight into the physics of active galactic nuclei.
9.3.2 Constraining black hole mass–host galaxy relations
The mass of the central SMBH M• is strongly correlated with several properties of
the host galaxy, e.g., bulge mass (Magorrian et al. 1998), galactic velocity dispersion
(Ferrarese and Merritt 2000; Gebhardt et al. 2000), and Sérsic index (Graham and
Driver 2007), among other properties. This suggests that SMBHs co-evolve with their
host galaxies. These relations have been critical throughout astronomy because of
their ability to predict SMBH masses where there is no direct measurement. However,
the reliability of the SMBH masses derived from these relations has recently come
into question, and it has been suggested that SMBH masses might be biased upward
(Shankar et al. 2016). As mentioned above, PTAs can only measure a degenerate chirp
mass and distance to the source. However, with the identification of the galactic host,
this degeneracy can be broken, permitting a measurement of Mc . As we write this,
PTAs are sensitive to SMBHBs of Mc > 109 M out to approximately 200 Mpc,
or a redshift of z 0.05 (Babak et al. 2016; NANOGrav Collaboration 2018). As
the PTA horizon grows and more distant sources are detected, PTAs will have some
limited capability to probe the redshift dependence of these relations, particularly at
the highest-mass end.
Note that, while Mc does not reflect the total black hole mass directly, numerous
simulations have demonstrated that all SMBHBs detected by PTAs will likely have
mass ratios q 0.1, implying an error when inferring M of ∼0.5 dex. Even in the
more conservative limit of q 0.01, the error in the total mass inferred from PTA
measurement of M would be only ∼ 1 dex.
As discussed here, PTAs can probe the SMBH mass–host relations through the
detection of continuous waves, however Simon and Burke-Spolaor (2016) demon-
strated that they can also probe these relations by looking at the amplitude of the GW
background. Because SMBH masses are currently linked to GW background predic-
tions through modeling of host galaxy populations, that study was able to use limits on
the GW background to place constraints on the high-mass slope, scatter, and intercept
of the relationship between SMBH mass and host bulge mass.
9.3.3 Active nuclei preceding and following SMBHB coalescence
In the case of a GW memory event, identification of the host will allow us to “wait
and watch” for a post-coalescence ignition of an active galactic nucleus. If there
is a circumbinary disk, it will closely track the shrinking orbit of the SMBHB. If
circumbinary disks are highly viscous, at some point it is possible that the binary and
circumbinary disk will “decouple” if the GW-induced inspiral becomes so rapid that
the viscosity of the disk no longer allows it to track the binary. After the coalescence,
there may be a delay in any nuclear activation while the disk settles; this delay depends
both on the disk viscosity, which governs the infall rate of the disk, and the SMBHB
masses, which govern the rapidity of the final break-away inspiral. The delay between
coalescence and subsequent turn-on of a luminous X-ray source can be modeled as:
123
1.32
M
tX−ray 7(1 + z) yr (31)
106 M
(Milosavljević and Phinney 2005). For the high-mass systems that PTAs are expected
to discover (M 108 M ), this implies timescales that are impractically long, in the
order of hundreds to thousands of years. However, Tanaka and Menou (2010) proposed
a much more rapid temporal evolution of the circumbinary disk, finding that there may
be a rapid brightening to super-Eddington luminosities in the few years following a
merger. These signatures would be detectable in X-rays, with potential signatures in
optical and infrared bands. Subsequent ignition of a radio jet might also arise once the
accretion disk is settled.
Another promising prospect is to search through archival data to uncover the evo-
lution of the system leading up to and during the merger in electromagnetic bands.
In practice, the extent of such studies will be highly dependent on the position of
the source on the sky, and the prevalence of historical data in that direction. How-
ever, given the broad coverage of current and upcoming synoptic sky surveys, like
the Catalina Sky Survey (CRTS, Drake et al. 2009), the Panoramic Survey Tele-
scope and Rapid Response System (Pan-STARRS, Kaiser et al. 2010), the Zwicky
Transient Facility (ZTF, Bellm 2014), the Large Synoptic Survey Telescope (LSST,
LSST Science Collaboration et al. 2017), and the Extended Roentgen Survey with
an Imaging Telescope Array (eROSITA, Merloni et al. 2012), we expect PTAs to
be detecting SMBHBs in an era rich with time-domain data over broad areas of the
sky.
9.3.4 SMBHB effects on galaxy central morphology and dynamics
Under certain circumstances (e. g., dry/low gas fraction mergers), the stellar scattering
described in Sect. 3.2.3 can work to flatten the density profile of a galaxy away from
a central stellar cusp, leading to a “core” in a post-merger system. Such cores are
commonly observed in the most luminous galaxies, and are marked by a flat stellar
photometric profile of the central ∼kpc of a galaxy (Faber et al. 1997), as opposed
to a much steeper profile external to the core. Simulations have shown direct links
between the core phenomenon and SMBHB stellar scattering (Ebisuzaki et al. 1991;
Makino 1997; Milosavljević and Merritt 2001), and some work has indicated these
cores may even survive subsequent mergers (Volonteri et al. 2003b). However, direct
observational proof of these links does not yet exist.
SMBHBs discovered by PTAs will have already undergone the stellar scattering
phase. With the identification of multiple galactic hosts of continuous-wave sources,
the presence of cores in those hosts could indicate direct links between the dynamical
actions of SMBHBs and the dynamics and morphology of the larger-scale core.
Taking a more general view the central dynamics of galaxies containing a
continuous-wave emitter, compared to those without a binary or at earlier stages of
orbital inspiral, may serve to reveal unanticipated stellar or gas dynamical interactions
during a merger.
123
9.3.5 Tidal disruption events
PTAs will not directly detect tidal disruption events, however a number of studies have
postulated that tidal disruption events may occur at much higher rates in SMBHB
systems. This is largely caused by the chaotic three-body interaction of the SMBHB
with stars in the stellar core, which can cause an increase of a few to several thousand
over the rate of disruption events around an isolated black hole (Chen et al. 2011;
Darbha et al. 2018). It has also been postulated that tidal disruption events might
encode electromagnetic information about the binary’s orbit in their time-resolved
signatures (Liu et al. 2014).
LSST is expected to detect thousands of tidal disruption events in the coming decade
(Wegg and Nate Bode 2011). This may highlight individual galaxies with excess event
rates. Such galaxies could signify the host of a potential continuous-wave GW source,
which would assist in host identification for any future continuous-wave detections by
PTAs.
Similarly, if we can otherwise identify the host of a GW, this detection can be
cross-matched with transient catalogs to test for excess events, or to allow direct
multi-messenger studies of any tidal disruption event caused by the continuous-wave
emitter. Any binary modeling extracted from a tidal disruption event will provide PTAs
greater sensitivity to the continuous waves from that event.
Extreme-mass-ratio inspiral events (commonly termed EMRIs; in this case, the
inspiral of the object before it is disrupted) may also be detected by LISA. The coor-
dinated detection of GWs from a binary SMBH by PTAs, an EMRI signature from
the low-mass object, and the electromagnetic observation of a tidal disruption event
could provide a fascinating set of studies of distinct processes in a common object.
9.3.6 Testing electromagnetic SMBHB emission models
Several binary SMBHB candidates have emerged, but it has proven difficult to find
a smoking-gun signature that can conclusively confirm the binary nature of these
sources. PTAs, on the other hand, can provide a straight-forward way to test the
binary hypothesis for such candidate systems identified in electromagnetic bands.
Given the current PTA capabilities, this test can be applied for candidates that are
sufficiently nearby, massive, and at the right orbital phase to be emitting GWs in
nanohertz frequencies, i.e., they need to be at sub-parsec orbital separations.
Even without a direct detection of nanohertz GWs, the PTA upper limits already
provide astrophysically relevant constraints on the SMBHB population. For instance,
recently, Sesana et al. (2018) and Holgado et al. (2018) examined the GW background
from the population of SMBHBs inferred from the samples of periodic quasars and
blazars, respectively. They found that, even if the individual candidates are below the
sensitivity of PTAs, the inferred population is in moderate tension with the current
limits on the GW background, which suggests that the samples may be contaminated
by false detections.
Additionally, in some cases, electromagnetic constraints on a binary system may
raise the PTA sensitivity for the same source by constraining the searchable parameter
space. There are already several examples of this type of multi-messenger astrophys-
123
ical application. Most famously, a M 1010 M SMBHB was hypothesized to exist

in galaxy 3C66B due to an elliptical motion detected through long-baseline interfer-
ometry; it was suggested that this implied a binary-modulated jet precession (Sudou
et al. 2003). However, Jenet et al. (2004) used only 7 years of data on PSR 1855+09
to place limits on the allowed chirp mass and eccentricity of this system, effectively
ruling out the purported binary parameters. More general applications of ruling out
binaries in specific systems involve targeted campaigns to place limits on binaries in
nearby, massive systems (Lommen and Backer 2001; Schutz and Ma 2016).
These types of studies are constantly improving, as the PTA horizon grows with
time, and with the addition of more pulsars. Based on the non-detection of continuous
waves in PTA data thus far, we can already en masse rule out nanohertz-frequency
binaries with M > 0.5 × 108 M in the Virgo cluster (NANOGrav Collaboration
2018).
9.3.7 Constraining SMBHB populations
The previous sections discussed electromagnetic signatures of SMBHB systems that

are primarily already in the microhertz to nanohertz gravitational waveband. However,
we note that there is also a direct benefit to systematically search for “GW precursor”
binaries—systems easily resolved as two AGN but well within the ∼few kpc virial
radius of a typical merger. In part because they likely have a higher residence time
at very large separations (Eq. 14), and in part because they are directly resolvable,
these systems may be more readily accessible to large electromagnetic surveys. These
can provide direct statistics on the systems that will make up our binary population
(e. g., Wen et al. 2009). These systems also complement GW-based studies of massive
merger environments and of the efficiency of SMBHB evolution (Sect. 3.2). So far,
there have been dozens of such discoveries; a few examples include the resolved dual
AGN in NGC6240 (Komossa et al. 2003), the spatially resolved broad line regions in
the galaxies of Comerford et al. (2013), and the record-holding binary at a separation of
7 pc in 0402+379 (Rodriguez et al. 2006). Recently, astrometric tracking of 0402+379
led Bansal et al. (2017) to suggest an orbital period for this binary of ∼ 104 year. For a
summary on the detections of dual AGN, we refer the reader to Rubinur et al. (2018)
Statistics on the population of dual AGN can directly constrain some of the covari-
ant parameters contributing to the GW background that are discussed in Sect. 3.2 and
summarized in Fig. 4. This includes eccentricity distributions, inferred inspiral rates,
and observed evidence of ongoing interaction with the galactic environment. How-
ever, the predictive power of these discoveries currently suffer from small number of
statistics. Future large-scale systematic surveys, such as LSST, or the Square Kilo-
meter Array, coupled with detailed observations using instruments with specialized
capabilities, such as high resolution, rapid spectrum acquisition, and/or time-resolved
imaging, will significantly enhance our understanding of the binary population.
9.3.8 Tests of gravity
As mentioned in Sect. 5.2, in general relativity electromagnetic and GW signals prop-

agate at the same speed. This agreement can be tested if a resolvable continuous wave
123
Fig. 18 Limits on the mass of

the graviton, given a
characteristic noise strain
amplitude for a pulsar timing
array with a continuous wave
signal of given chirp mass
(Hazboun et al. 2013). The lines
and coloring represent contours
of the graviton mass limit. A
fiducial frequency of 10−9 Hz is
used for this illustration
source also happens to be an eclipsing electromagnetic binary (Cutler et al. 2003; Will
1998). Such a system allows for a differential test of the propagation speed of GWs
and light, and the mass of the graviton. Using the order-of-magnitude limit on the
graviton Compton wavelength (λg = h/(m g c)) of such a system (Will 1998):
1 1
D 100Hz 2 1 2
λg > 3 × 10 km 12
(32)
200Mpc f f t
and choosing nominal values for the parameters of an SMBHB, of f = 1 × 10−9

Hz, D = 1 Gpc and an observing time of t = 5 years, one obtains a lower limit of
λg ∼ 1019 km. This limit is three orders of magnitude less precise than the current
Particle Data Group limit of 1022 km, suggesting that today’s PTAs might not be able
to produce competitive limits without more sophisticated GW data analysis schemes
(Choudhury et al. 2004).
The approaches in Will (1998), Cutler et al. (2003), and Hazboun et al. (2013)
are framed in the context of a null test, in which the assumed observation is the
simultaneous arrival, within the margin of error, of some fiducial feature of both
the electromagnetic and GW signal. While the error in arrival time is that of both
the electromagnetic and GW signals, in practice the electromagnetic arrival time is
much more precise, and is often ignored in terms of these calculations. The error,
related to the signal-to-noise ratio in general, and to the strain sensitivity for a specific
detector, is then used to set a limit on the difference in propagation speed between
the electromagnetic and GW signals. This velocity limit can then be used to calculate
a lower limit on the Compton wavelength and an upper limit on the mass of the
graviton m g , if one assumes a Yukawa-type potential (Vg ∼ e−kmg /r ) for the graviton.
This type of potential is a phenomenological model which acts to restrict a interaction’s
ability to act at a distance. Figure 18 shows a contour plot of the limits on m g for a
generic PTA at f gw = 10−9 Hz.
Takahashi (2017) considers an alternate approach to testing the concordance of
electromagnetic and GW propagation using gravitational lensing. Standard analyses
assume that gravitational lensing can be treated with geometric optics. This assumption
123
is no longer applicable if the mass of the lens is too small compared to the wavelength
of the radiation; less than 105 M ( f /Hz)−1 .
Both of these tests of gravitational/electromagnetic wave propagation concordance
would require the detection of an individual SMBH binary with a high enough signal-
to-noise ratio in GWs and some electromagnetic measure of the orbiting SMBHs
(e.g., optical or X-ray light curve). Takahashi (2017) shows that, for signal-to-noise
ratios of 10–100 and GW frequencies f ∼ 10 nHz, there is a reasonable probability of
obtaining a few lenses, for which the time difference between the GW and electromag-
netic signals would be of order 1–10 days, if of order 100 individual SMBH binaries
could be detected as (continuous wave) GW sources. (Longer time delays are required
for detection at lower signal-to-noise ratio.) At higher (lower) GW frequencies, shorter
(longer) time delays result in similar predictions for comparable signal-to-noise ratios.
10 Current outlook
The focus of this paper has been on the possible sources that could generate nanohertz
GW signals. We now consider briefly the observational status of pulsar timing arrays,
both current and near future, as a means of providing some measure of the likelihood
of detecting these signals. As indicated in the first two sections above and referenced
throughout, the sensitivity of pulsar timing arrays depends on length of data set, RMS
timing precision of pulsars, and the total number of contributing pulsars.
The three large national and regional PTAs are the European Pulsar Timing Array
(EPTA, Babak et al. 2016; Desvignes et al. 2016; Lentati et al. 2015), the North Ameri-
can Nanohertz Observatory for Gravitational Waves (NANOGrav, Arzoumanian et al.
2015a, 2016), and the Parkes Pulsar Timing Array (PPTA, Manchester et al. 2013).
These efforts, in addition to other emerging timing arrays in other nations (China,
India), are currently working to combine their data into an International Pulsar Tim-
ing Array (IPTA; e. g., Hobbs et al. 2010). In the future, it is expected that new facilities
such as MeerKAT and the Square Kilometre Array will increase the sensitivity of PTAs
manifold through extensive pulsar searching and timing programs (e. g., Janssen et al.
2015).
As an approximate illustration of the performance of these PTAs, recent published
work has reported the timing of 20–40 millisecond pulsars with sub-microsecond pre-
cision, while some PTAs have begun to time > 70 pulsars to this precision. Currently,
the total number of pulsars timed worldwide does not yet exceed 100, because there
is overlap between the PTAs. Some of this overlap is intentional, for the purposes
of calibration between the different telescopes, in addition to providing the ability to
independently verify any reported GW detection.
Millisecond pulsars continue to be discovered (e.g., Cromartie et al. 2016; Pleunis
et al. 2017), and, while not all are suitable for inclusion into a PTA, many are, so these
PTAs can and are still expanding. As a specific example, the NANOGrav 9 Year Data
Release contained 37 pulsars (Arzoumanian et al. 2015a), its 11 Year Data Release
contained 45 pulsars (Arzoumanian et al. 2018b), and that collaboration is now tim-
ing more than 70 pulsars. In addition to new pulsar discoveries, improvements to
radio observing systems will continue to improve mitigation of intrinsic pulse-to-
123
pulse variations (“jitter”; Liu et al. 2012), which dominates the noise budget of many
PTA pulsars. We also have an ever-increasing understanding of other contributions to
the PTA noise budget, including long-term timing noise (“red noise”); some of this
is now commonly mitigated by tracking radio-frequency-dependent variations in the
interstellar medium, although the origin of long-term noise in some pulsars remains
yet unclear (e. g., Keith et al. 2013). These improvements can provide sudden improve-
ments in PTA sensitivity, as the data are re-processed using novel noise suppression
techniques (e. g., Jones et al. 2017). Given these ongoing improvements, we consider
it plausible that together, worldwide PTAs will include approximately 100 pulsars at
typical precisions of approximately 500 ns by 2020. With the use of MeerKAT and
SKA, these numbers will increase even more rapidly. Of course, with all PTAs, sim-
ply keeping continuous timing programs for ongoing pulsars to obtain a longer time
baseline will also afford increased sensitivity.
Many of the possibilities discussed in later sections are somewhat more speculative
(Sects. 6.1, 5, 7). In many cases, these are likely to be “observationally driven,” with
models that are updated to evade new limits as they are placed. Most of these models are
capable of producing data over many orders of magnitude across many parameters,
thus updated models will continually be available, however, they will likely cover
smaller and smaller allowed parameter space.
Current PTAs are already placing astrophysical limits on the existence of SMBH
binaries and their interactions with their environments (Sect. 3) and the existence of
cosmic strings (Sect. 4). Based on projections of SMBHB populations, we expect a
GW background detection to be imminent from that population (Vigeland and Siemens
2016), with the detection of continuous-wave SMBHBs soon to follow (Mingarelli
et al. 2017). The detection of these systems will be occurring concurrently with a
number of planned synoptic sky surveys, enabling a much greater scope for multi-
messenger studies (Sect. 9) in addition to potential multi-band GW science with LISA
when it flies in the mid-2030s (Sect. 8). For the forseeable future, PTAs will main-
tain unique access to exploring the inspiral of high-mass SMBHB systems, and will
uniquely probe the physical properties of strings.
Acknowledgements We thank Andrew Kaiser and Sean Williams for providing Fig. 16; for questions
regarding this Figure please contact [email protected]. NANOGrav is supported by NSF Physics
Frontier Center award #1430284. SBS is supported by NSF award #1458952 and by the CIFAR Azrieli
Global Scholars program. Portions of this research were carried out at the Jet Propulsion Laboratory,
California Institute of Technology, under a contract with the National Aeronautics and Space Administration.
DRM was, until recently, a Jansky Fellow of the National Radio Astronomy Observatory, a facility of the NSF
operated under cooperative agreement by Associated Universities, Inc. The Flatiron Institute is supported
by the Simons Foundation. AR received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 Programme for Research and Innovation ERC-2014-STG under grant
agreement No. 638435 (GalNUC).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna-
tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons license, and indicate if changes were made.
123
References
Abbott BP, Abbott R, Abbott TD, Abernathy MR, Acernese F, Ackley K, Adams C, Adams T, Addesso
P, Adhikari RX et al (2016a) GW151226: observation of gravitational waves from a 22-solar-mass
binary black hole coalescence. Phys Rev Lett 116(24):241103. https://doi.org/10.1103/PhysRevLett.
116.241103. arXiv:1606.04855
Abbott BP, Abbott R, Abbott TD, Abernathy MR, Acernese F, Ackley K, Adams C, Adams T, Addesso P,
Adhikari RX et al (2016b) Observation of gravitational waves from a binary black hole merger. Phys
Rev Lett 116(6):061102. https://doi.org/10.1103/PhysRevLett.116.061102 arXiv:1602.03837
Abbott BP, Abbott R, Abbott TD, Acernese F, Ackley K, Adams C, Adams T, Addesso P, Adhikari RX,
Adya VB et al (2017a) GW170104: observation of a 50-solar-mass binary black hole coalescence
at Redshift 0.2. Phys Rev Lett 118(22):221101. https://doi.org/10.1103/PhysRevLett.118.221101.
arXiv:1706.01812
Abbott BP, Abbott R, Abbott TD, Acernese F, Ackley K, Adams C, Adams T, Addesso P, Adhikari RX,
Adya VB et al (2017b) Multi-messenger observations of a binary neutron star merger. ApJ 848:L12.
https://doi.org/10.3847/2041-8213/aa91c9 arXiv:1710.05833
Abbott BP et al (2017c) Gravitational waves and gamma-rays from a binary neutron star merger:
GW170817 and GRB 170817A. Astrophys J 848(2):L13. https://doi.org/10.3847/2041-8213/aa920c
arXiv:1710.05834
Abbott BP et al (2017d) GW170104: observation of a 50-solar-mass binary black hole coalescence
at redshift 0.2. Phys Rev Lett 118(22):221,101. https://doi.org/10.1103/PhysRevLett.118.221101.
arXiv:1706.01812
Abbott BP et al (2017e) GW170814: a three-detector observation of gravitational waves from a binary black
hole coalescence. Phys Rev Lett 119(14):141,101. https://doi.org/10.1103/PhysRevLett.119.141101.
arXiv:1709.09660
Abbott BP et al (2017f) GW170817: observation of gravitational waves from a binary neutron star inspiral.
Phys Rev Lett 119(16):161,101. https://doi.org/10.1103/PhysRevLett.119.161101. arXiv:1710.05832
Abbott BP et al (2017g) Multi-messenger observations of a binary neutron star merger. Astrophys J
848(2):L12. https://doi.org/10.3847/2041-8213/aa91c9 arXiv:1710.05833
Akrami Y, Koivisto TS, Sandstad M (2013) Accelerated expansion from ghost-free bigravity: a statis-
tical analysis with improved generality. JHEP 03:099. https://doi.org/10.1007/JHEP03(2013)099.
arXiv:1209.0457
Amaro-Seoane P, Sesana A, Hoffman L, Benacquista M, Eichhorn C, Makino J, Spurzem R (2010) Triplets
of supermassive black holes: astrophysics, gravitational waves and detection. MNRAS 402(4):2308–
2320
Amaro-Seoane P, Aoudia S, Babak S, Binétruy P, Berti E, Bohé A, Caprini C, Colpi M, Cornish NJ,
Danzmann K, Dufaux JF, Gair J, Jennrich O, Jetzer P, Klein A, Lang RN, Lobo A, Littenberg T,
McWilliams ST, Nelemans G, Petiteau A, Porter EK, Schutz BF, Sesana A, Stebbins R, Sumner
T, Vallisneri M, Vitale S, Volonteri M, Ward H (2012) Low-frequency gravitational-wave science
with eLISA/NGO. Class Quantum Grav 29(12):124016. https://doi.org/10.1088/0264-9381/29/12/
124016. arXiv:1202.0839
Amaro-Seoane P, Audley H, Babak S, Baker J, Barausse E, Bender P, Berti E, Binetruy P, Born M, Bortoluzzi
D, Camp J, Caprini C, Cardoso V, Colpi M, Conklin J, Cornish N, Cutler C, Danzmann K, Dolesi R,
Ferraioli L, Ferroni V, Fitzsimons E, Gair J, Gesa Bote L, Giardini D, Gibert F, Grimani C, Halloin
H, Heinzel G, Hertog T, Hewitson M, Holley-Bockelmann K, Hollington D, Hueller M, Inchauspe H,
Jetzer P, Karnesis N, Killow C, Klein A, Klipstein B, Korsakova N, Larson SL, Livas J, Lloro I, Man
N, Mance D, Martino J, Mateos I, McKenzie K, McWilliams ST, Miller C, Mueller G, Nardini G,
Nelemans G, Nofrarias M, Petiteau A, Pivato P, Plagnol E, Porter E, Reiche J, Robertson D, Robertson
N, Rossi E, Russano G, Schutz B, Sesana A, Shoemaker D, Slutsky J, Sopuerta CF, Sumner T, Tamanini
N, Thorpe I, Troebs M, Vallisneri M, Vecchio A, Vetrugno D, Vitale S, Volonteri M, Wanner G, Ward
H, Wass P, Weber W, Ziemer J, Zweifel P (2017) Laser interferometer space antenna. ArXiv e-prints
arXiv:1702.00786
Anholm M, Ballmer S, Creighton JDE, Price LR, Siemens X (2009) Optimal strategies for gravitational
wave stochastic background searches in pulsar timing data. Phys Rev D 79(8):084030. https://doi.org/
10.1103/PhysRevD.79.084030. arXiv:0809.0701
Antonini F, Merritt D (2012) Dynamical friction around supermassive black holes. ApJ 745:83. https://doi.
org/10.1088/0004-637X/745/1/83. arXiv:1108.1163
123
Armitage PJ, Natarajan P (2005) Eccentricity of supermassive black hole binaries coalescing from gas-rich
mergers. ApJ 634(2):921–927
Armstrong JW, Iess L, Tortora P, Bertotti B (2003) Stochastic Gravitational wave background: upper limits
in the 10−6 to 10−3 Hz band. ApJ 599:806–813. https://doi.org/10.1086/379505
Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin SJ, Chatterjee S, Cordes JM, Demorest PB,
Deng X, Dolch T, Ellis JA, Ferdman RD, Garver-Daniels N, Jenet F, Jones G, Kaspi VM, Koop M,
Lam MT, Lazio TJW, Lommen AN, Lorimer DR, Luo J, Lynch RS, Madison DR, McLaughlin MA,
McWilliams ST, Nice DJ, Palliyaguru N, Pennucci TT, Ransom SM, Sesana A, Siemens X, Stairs IH,
Stinebring DR, Stovall K, Swiggum J, Vallisneri M, van Haasteren R, Wang Y, Zhu WW, NANOGrav
Collaboration (2014) Gravitational waves from individual supermassive black hole binaries in circular
orbits: limits from the north american nanohertz observatory for gravitational waves. ApJ 794:141.
https://doi.org/10.1088/0004-637X/794/2/141. arXiv:1404.1267
Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin S, Chatterjee S, Christy B, Cordes JM, Cornish
N, Crowter K, Demorest PB, Dolch T, Ellis JA, Ferdman RD, Fonseca E, Garver-Daniels N, Gonzalez
ME, Jenet FA, Jones G, Jones ML, Kaspi VM, Koop M, Lam MT, Lazio TJW, Levin L, Lommen AN,
Lorimer DR, Luo J, Lynch RS, Madison D, McLaughlin MA, McWilliams ST, Nice DJ, Palliyaguru N,
Pennucci TT, Ransom SM, Siemens X, Stairs IH, Stinebring DR, Stovall K, Swiggum JK, Vallisneri
M, van Haasteren R, Wang Y, Zhu W (2015a) The NANOGrav nine-year data set: observations, arrival
time measurements, and analysis of 37 millisecond pulsars. ApJ 813:65. https://doi.org/10.1088/0004-
637X/813/1/65. arXiv:1505.07540
Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin SJ, Chatterjee S, Christy B, Cordes JM, Cornish
NJ, Demorest PB, Deng X, Dolch T, Ellis JA, Ferdman RD, Fonseca E, Garver-Daniels N, Jenet F,
Jones G, Kaspi VM, Koop M, Lam MT, Lazio TJW, Levin L, Lommen AN, Lorimer DR, Luo J, Lynch
RS, Madison DR, McLaughlin MA, McWilliams ST, Nice DJ, Palliyaguru N, Pennucci TT, Ransom
SM, Siemens X, Stairs IH, Stinebring DR, Stovall K, Swiggum J, Vallisneri M, van Haasteren R, Wang
Y, Zhu WW, NANOGrav Collaboration (2015b) NANOGrav constraints on gravitational wave bursts
with memory. ApJ 810:150. https://doi.org/10.1088/0004-637X/810/2/150 arXiv:1501.05343
Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin SJ, Chatterjee S, Christy B, Cordes JM, Cornish
NJ, Crowter K, Demorest PB, Deng X, Dolch T, Ellis JA, Ferdman RD, Fonseca E, Garver-Daniels
N, Gonzalez ME, Jenet F, Jones G, Jones ML, Kaspi VM, Koop M, Lam MT, Lazio TJW, Levin
L, Lommen AN, Lorimer DR, Luo J, Lynch RS, Madison DR, McLaughlin MA, McWilliams ST,
Mingarelli CMF, Nice DJ, Palliyaguru N, Pennucci TT, Ransom SM, Sampson L, Sanidas SA, Sesana
A, Siemens X, Simon J, Stairs IH, Stinebring DR, Stovall K, Swiggum J, Taylor SR, Vallisneri M, van
Haasteren R, Wang Y, Zhu WW, NANOGrav Collaboration (2016) The NANOGrav nine-year data
set: limits on the isotropic stochastic gravitational wave background. ApJ 821:13. https://doi.org/10.
3847/0004-637X/821/1/13, arXiv:1508.03024
Arzoumanian Z, Baker PT, Brazier A, Burke-Spolaor S, Chamberlin SJ, Chatterjee S, Christy B, Cordes
JM, Cornish NJ, Crawford F, Thankful Cromartie H, Crowter K, DeCesar M, Demorest PB, Dolch T,
Ellis JA, Ferdman RD, Ferrara E, Folkner WM, Fonseca E, Garver-Daniels N, Gentile PA, Haas R,
Hazboun JS, Huerta EA, Islo K, Jones G, Jones ML, Kaplan DL, Kaspi VM, Lam MT, Lazio TJW,
Levin L, Lommen AN, Lorimer DR, Luo J, Lynch RS, Madison DR, McLaughlin MA, McWilliams ST,
Mingarelli CMF, Ng C, Nice DJ, Park RS, Pennucci TT, Pol NS, Ransom SM, Ray PS, Rasskazov A,
Siemens X, Simon J, Spiewak R, Stairs IH, Stinebring DR, Stovall K, Swiggum J, Taylor SR, Vallisneri
M, van Haasteren R, Vigeland S, Zhu WW, The NANOGrav Collaboration (2018a) The NANOGrav 11
year data set: pulsar-timing constraints on the stochastic gravitational-wave background. ApJ 859:47.
https://doi.org/10.3847/1538-4357/aabd3b, arXiv:1801.02617
Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin S, Chatterjee S, Christy B, Cordes JM, Cornish
NJ, Crawford F, Thankful Cromartie H, Crowter K, DeCesar ME, Demorest PB, Dolch T, Ellis JA,
Ferdman RD, Ferrara EC, Fonseca E, Garver-Daniels N, Gentile PA, Halmrast D, Huerta EA, Jenet
FA, Jessup C, Jones G, Jones ML, Kaplan DL, Lam MT, Lazio TJW, Levin L, Lommen A, Lorimer
DR, Luo J, Lynch RS, Madison D, Matthews AM, McLaughlin MA, McWilliams ST, Mingarelli
C, Ng C, Nice DJ, Pennucci TT, Ransom SM, Ray PS, Siemens X, Simon J, Spiewak R, Stairs IH,
Stinebring DR, Stovall K, Swiggum JK, Taylor SR, Vallisneri M, van Haasteren R, Vigeland SJ, Zhu
W, NANOGrav Collaboration (2018b) The NANOGrav 11-year data set: high-precision timing of 45
millisecond pulsars. ApJS 235:37. https://doi.org/10.3847/1538-4365/aab5b0, arXiv:1801.01837
Babak S, Sesana A (2012) Resolving multiple supermassive black hole binaries with pulsar timing arrays.
Phys Rev D 85(4):044034. https://doi.org/10.1103/PhysRevD.85.044034. arXiv:1112.1075
123
Babak S, Petiteau A, Sesana A, Brem P, Rosado PA, Taylor SR, Lassus A, Hessels JWT, Bassa CG, Burgay
M, Caballero RN, Champion DJ, Cognard I, Desvignes G, Gair JR, Guillemot L, Janssen GH, Karup-
pusamy R, Kramer M, Lazarus P, Lee KJ, Lentati L, Liu K, Mingarelli CMF, Osłowski S, Perrodin
D, Possenti A, Purver MB, Sanidas S, Smits R, Stappers B, Theureau G, Tiburzi C, van Haasteren
R, Vecchio A, Verbiest JPW (2016) European pulsar timing array limits on continuous gravitational
waves from individual supermassive black hole binaries. MNRAS 455:1665–1679. https://doi.org/10.
1093/mnras/stv2092 arXiv:1509.02165
Baghram S, Afshordi N, Zurek KM (2011) Prospects for detecting dark matter halo substructure with pulsar
timing. Phys Rev D 84:043511. https://doi.org/10.1103/PhysRevD.84.043511. arXiv:1101.5487
Baker T, Bellini E, Ferreira PG, Lagos M, Noller J, Sawicki I (2017) Strong constraints on cosmological
gravity from GW170817 and GRB 170817A. Phys Rev Lett 119(25):251, 301. https://doi.org/10.
1103/PhysRevLett.119.251301. arXiv:1710.06394
Bansal K, Taylor GB, Peck AB, Zavala RT, Romani RW (2017) Constraining the orbit of the
supermassive black hole binary 0402+379. ApJ 843:14. https://doi.org/10.3847/1538-4357/aa74e1.
arXiv:1705.08556
Baskaran D, Polnarev AG, Pshirkov MS, Postnov KA (2008) Limits on the speed of gravitational
waves from pulsar timing. Phys Rev D 78(044):018. https://doi.org/10.1103/PhysRevD.78.044018.
arXiv:0805.3103
Bassa CG, Janssen GH, Stappers BW, Tauris TM, Wevers T, Jonker PG, Lentati L, Verbiest JPW, Desvignes
G, Graikou E, Guillemot L, Freire PCC, Lazarus P, Caballero RN, Champion DJ, Cognard I, Jessner
A, Jordan C, Karuppusamy R, Kramer M, Lazaridis K, Lee KJ, Liu K, Lyne AG, McKee J, Osłowski S,
Perrodin D, Sanidas S, Shaifullah G, Smits R, Theureau G, Tiburzi C, Zhu WW (2016) A millisecond
pulsar in an extremely wide binary system. MNRAS 460:2207–2222. https://doi.org/10.1093/mnras/
stw1134. arXiv:1604.00129
Battye RA, Shellard EPS (1996) Primordial gravitational waves: a probe of the early universe. ArXiv e-prints
arXiv:astro-ph/9604059
Batygin K, Brown ME (2016) Evidence for a distant giant planet in the solar system. AJ 151:22. https://
doi.org/10.3847/0004-6256/151/2/22. arXiv:1601.05438
Begelman MC, Blandford RD, Rees MJ (1980) Massive black hole binaries in active galactic nuclei. Nature
287:307–309. https://doi.org/10.1038/287307a0
Bellm E (2014) The Zwicky transient facility. In: Wozniak PR, Graham MJ, Mahabal AA, Seaman R
(eds) Proceedings of the third hot-wiring transient universe workshop, SLAC, Stanford, eConf, vol
C131113.1, pp 27–33, arXiv:1410.8185
Benini R, Lattanzi M, Montani G (2011) Signatures of the neutrino thermal history in the spectrum of
primordial gravitational waves. Gen Relat Gravity 43:945–958. https://doi.org/10.1007/s10714-010-
0994-4. arXiv:1009.6110
Bennett MF, Melatos A (2014) Stochastic microhertz gravitational radiation from stellar convection. ApJ
792:55. https://doi.org/10.1088/0004-637X/792/1/55. arXiv:1407.2281
Bernabeu J, Espriu D, Puigdomènech D (2011) Gravitational waves in the presence of a cosmological
constant. Phys Rev D 84(6):063523. https://doi.org/10.1103/PhysRevD.84.063523. arXiv:1106.4511
BICEP2/Keck Collaboration, Planck Collaboration, Ade PAR, Aghanim N, Ahmed Z, Aikin RW, Alexan-
der KD, Arnaud M, Aumont J, Baccigalupi C, et al (2015) Joint analysis of BICEP2/keck array
and planck data. Phys Rev Lett 114(10):101301. https://doi.org/10.1103/PhysRevLett.114.101301.
arXiv:1502.00612
Binney J, Tremaine S (1987) Galactic dynamics. Princeton University Press, Princeton, NJ
Blaes O, Lee MH, Socrates A (2002) The kozai mechanism and the evolution of binary supermassive black
holes. ApJ 578(2):775–786
Blanchet L (2006) Gravitational radiation from post-Newtonian sources and inspiralling compact binaries.
Living Rev Relat 9:4. https://doi.org/10.12942/lrr-2006-4
Blanco-Pillado JJ, Olum KD, Siemens X (2018) New limits on cosmic strings from gravitational wave obser-
vation. Phys Lett B 778:392–396. https://doi.org/10.1016/j.physletb.2018.01.050. arXiv:1709.02434
Blas D, Lim E (2015) Phenomenology of theories of gravity without Lorentz invariance: the preferred frame
case. Int J Mod Phys D 23(1443):009. https://doi.org/10.1142/S0218271814430093. arXiv:1412.4828
las D, Ivanov MM, Sawicki I, Sibiryakov S (2016) On constraining the speed of gravitational waves following
GW150914. JETP Lett 103(10):624–626, https://doi.org/10.1134/S0021364016100040, https://doi.
org/10.7868/S0370274X16100039, [Pisma Zh Eksp Teor Fiz 103:708 (2016)], arXiv:1602.04188
123
Bonetti M, Haardt F, Sesana A, Barausse E (2018) Post-Newtonian evolution of massive black hole triplets
in galactic nuclei - II. Survey of the parameter space. MNRAS 477:3910–3926. https://doi.org/10.
1093/mnras/sty896. arXiv:1709.06088
Book LG, Flanagan ÉÉ (2011) Astrometric effects of a stochastic gravitational wave background. Phys Rev
D 83(2):024024. https://doi.org/10.1103/PhysRevD.83.024024. arXiv:1009.4192
Boran S, Desai S, Kahya EO, Woodard RP (2018) GW170817 falsifies dark matter emulators. Phys Rev D
97:041501. https://doi.org/10.1103/PhysRevD.97.041501 arXiv:1710.06168
Boulware DG, Deser S (1972) Can gravitation have a finite range? Phys Rev D 6:3368–3382. https://doi.
org/10.1103/PhysRevD.6.3368
Braginskii VB, Thorne KS (1987) Gravitational-wave bursts with memory and experimental prospects.
Nature 327:123–125. https://doi.org/10.1038/327123a0
Brans C, Dicke RH (1961) Mach’s principle and a relativistic theory of gravitation. Phys Rev 124:925–935.
https://doi.org/10.1103/PhysRev.124.925
Britzen S, Qian SJ, Steffen W, Kun E, Karouzos M, Gergely L, Schmidt J, Aller M, Aller H, Krause M,
Fendt C, Böttcher M, Witzel A, Eckart A, Moser L (2017) A swirling jet in the quasar 1308+326.
A&A 602:A29. https://doi.org/10.1051/0004-6361/201629999
Burke-Spolaor S (2011) A radio census of binary supermassive black holes. MNRAS 410:2113–2122.
https://doi.org/10.1111/j.1365-2966.2010.17586.x. arXiv:1008.4382
Burke-Spolaor S (2013) Multi-messenger approaches to binary supermassive black holes in the
‘continuous-wave’ regime. Class Quantum Grav 30(22):224013. https://doi.org/10.1088/0264-9381/
30/22/224013. arXiv:1308.4408
Caballero RN, Guo YJ, Lee KJ, Lazarus P, Champion DJ, Desvignes G, Kramer M, Plant K, Arzoumanian Z,
Bailes M, Bassa CG, Bhat NDR, Brazier A, Burgay M, Burke-Spolaor S, Chamberlin SJ, Chatterjee S,
Cognard I, Cordes JM, Dai S, Demorest P, Dolch T, Ferdman RD, Fonseca E, Gair JR, Garver-Daniels
N, Gentile P, Gonzalez ME, Graikou E, Guillemot L, Hobbs G, Janssen GH, Karuppusamy R, Keith
MJ, Kerr M, Lam MT, Lasky PD, Lazio TJW, Levin L, Liu K, Lommen AN, Lorimer DR, Lynch
RS, Madison DR, Manchester RN, McKee JW, McLaughlin MA, McWilliams ST, Mingarelli CMF,
Nice DJ, Osłowski S, Palliyaguru NT, Pennucci TT, Perera BBP, Perrodin D, Possenti A, Ransom
SM, Reardon DJ, Sanidas SA, Sesana A, Shaifullah G, Shannon RM, Siemens X, Simon J, Spiewak
R, Stairs I, Stappers B, Stinebring DR, Stovall K, Swiggum JK, Taylor SR, Theureau G, Tiburzi
C, Toomey L, van Haasteren R, van Straten W, Verbiest JPW, Wang JB, Zhu XJ, Zhu WW (2018)
Studying the solar system with the international pulsar timing array. MNRAS 481:5501–5516. https://
doi.org/10.1093/mnras/sty2632. arXiv:1809.10744
Caprini C, Durrer R, Siemens X (2010) Detection of gravitational waves from the QCD phase transition
with pulsar timing arrays. Phys Rev D 82(6):063511. https://doi.org/10.1103/PhysRevD.82.063511.
arXiv:1007.1218
Carroll SM, Duvvuri V, Trodden M, Turner MS (2004) Is cosmic speed-up due to new gravitational physics?
Phys Rev D 70(043):528. https://doi.org/10.1103/PhysRevD.70.043528. arXiv:astro-ph/0306438
Champion DJ, Hobbs GB, Manchester RN, Edwards RT, Backer DC, Bailes M, Bhat NDR, Burke-Spolaor
S, Coles W, Demorest PB, Ferdman RD, Folkner WM, Hotan AW, Kramer M, Lommen AN, Nice DJ,
Purver MB, Sarkissian JM, Stairs IH, van Straten W, Verbiest JPW, Yardley DRB (2010) Measuring
the mass of solar system planets using pulsar timing. ApJ 720:L201–L205. https://doi.org/10.1088/
2041-8205/720/2/L201. arXiv:1008.3607
Chandrasekhar S (1943) Dynamical friction. I. General considerations: the coefficient of dynamical friction.
ApJ 97:255. https://doi.org/10.1086/144517
Charisi M, Bartos I, Haiman Z, Price-Whelan AM, Márka S (2015) Multiple periods in the variability of
the supermassive black hole binary candidate quasar PG1302–102 MNRAS 454:L21–L25, https://
doi.org/10.1093/mnrasl/slv111, arXiv:1502.03113
Charisi M, Bartos I, Haiman Z, Price-Whelan AM, Graham MJ, Bellm EC, Laher RR, Márka S (2016) A
population of short-period variable quasars from PTF as supermassive black hole binary candidates.
MNRAS 463:2145–2171. https://doi.org/10.1093/mnras/stw1838. arXiv:1604.01020
Charisi M, Haiman Z, Schiminovich D, D’Orazio DJ (2018) Testing the relativistic Doppler boost hypothesis
for supermassive black hole binary candidates. MNRAS 476:4617–4628. https://doi.org/10.1093/
mnras/sty516. arXiv:1801.06189
Chatziioannou K, Yunes N, Cornish N (2012) Model-independent test of general relativity: an extended post-
Einsteinian framework with complete polarization content. Phys Rev D 86(022):004. https://doi.org/
123
10.1103/PhysRevD.86.022004, https://doi.org/10.1103/PhysRevD.95.129901, [Erratum: Phys. Rev.

D95, no.12,129901(2017)] arXiv:1204.2585
Chen X, Sesana A, Madau P, Liu FK (2011) Tidal stellar disruptions by massive black hole pairs. II.
Decaying binaries. ApJ 729:13. https://doi.org/10.1088/0004-637X/729/1/13. arXiv:1012.4466
Cheng SL, Lee W (2017) Ng KW (2017) Production of high stellar-mass primordial black holes in trapped
inflation. JHEP 02:8. https://doi.org/10.1007/JHEP02(2017)008. arXiv:1606.00206
Chiba T (2003) 1/R gravity and scalar-tensor gravity. Phys Lett B 575:1–3. https://doi.org/10.1016/j.
physletb.2003.09.033. arXiv:astro-ph/0307338
Chiba T, Smith TL, Erickcek AL (2007) Solar System constraints to general f (R) gravity. Phys Rev D
75(124):014. https://doi.org/10.1103/PhysRevD.75.124014. arXiv:astro-ph/0611867
Choudhury SR, Joshi GC, Mahajan S, McKellar BHJ (2004) Probing large distance higher dimensional
gravity from lensing data. Astropart Phys 21:559–563. https://doi.org/10.1016/j.astropartphys.2004.
04.001. arXiv:hep-ph/0204161
Christodoulou D (1991) Nonlinear nature of gravitation and gravitational-wave experiments. Phys Rev Lett
67:1486–1489. https://doi.org/10.1103/PhysRevLett.67.1486
Coles WA, Kerr M, Shannon RM, Hobbs GB, Manchester RN, You XP, Bailes M, Bhat NDR, Burke-Spolaor
S, Dai S, Keith MJ, Levin Y, Osłowski S, Ravi V, Reardon D, Toomey L, van Straten W, Wang JB,
Wen L, Zhu XJ (2015) Pulsar observations of extreme scattering events. ApJ 808:113. https://doi.org/
10.1088/0004-637X/808/2/113. arXiv:1506.07948
Comerford JM, Schluns K, Greene JE, Cool RJ (2013) Dual supermassive black hole candidates in
the AGN and galaxy evolution survey. ApJ 777:64. https://doi.org/10.1088/0004-637X/777/1/64.
arXiv:1309.2284
Conneely C, Jaffe AH, Mingarelli CMF (2019) On the amplitude and stokes parameters of a stochastic
gravitational-wave background. MNRAS https://doi.org/10.1093/mnras/stz1022, arXiv:1808.05920
Corbin V, Cornish NJ (2010) Pulsar timing array observations of massive black hole binaries. ArXiv e-prints
arXiv:1008.1782
Cordes JM, Chernoff DF (1998) Neutron star population dynamics. II. Three-dimensional space velocities
of young pulsars. ApJ 505:315–338. https://doi.org/10.1086/306138. arXiv:astro-ph/9707308
Cordes JM, Jenet FA (2012) Detecting gravitational wave memory with pulsar timing. ApJ 752:54. https://
doi.org/10.1088/0004-637X/752/1/54
Cordes JM, Lazio TJW (2002) NE2001.I. A new model for the galactic distribution of free electrons and
its fluctuations. ArXiv e-prints arXiv:astro-ph/0207156
Cornish NJ, van Haasteren R (2014) Mapping the nano-Hertz gravitational wave sky. ArXiv e-prints
arXiv:1406.4511
Cornish NJ, O’Beirne L, Taylor SR, Yunes N (2018) Constraining alternative theories of gravity using pul-
sar timing arrays. Phys Rev Lett 120(18):181, 101. https://doi.org/10.1103/PhysRevLett.120.181101.
arXiv:1712.07132
Cromartie HT, Camilo F, Kerr M, Deneva JS, Ransom SM, Ray PS, Ferrara EC, Michelson PF, Wood KS
(2016) Six new millisecond pulsars from arecibo searches of fermi gamma-ray sources. ApJ 819:34.
Cuadra J, Armitage PJ, Alexander RD, Begelman MC, Alexander RD (2009) Massive black hole binary
mergers within subparsec scale gas discs. MNRAS 393(4):1423–1432
Cutler C, Hiscock WA, Larson SL (2003) LISA, binary stars, and the mass of the graviton. Phys Rev D
67(024):015. https://doi.org/10.1103/PhysRevD.67.024015. arXiv:gr-qc/0209101
Cutler C, Burke-Spolaor S, Vallisneri M, Lazio J, Majid W (2014) The gravitational-wave discovery space
of pulsar timing arrays. Phys Rev D 89(4):042003. https://doi.org/10.1103/PhysRevD.89.042003.
arXiv:1309.2581
van Dam H, Veltman MJG (1970) Massive and massless Yang–Mills and gravitational fields. Nucl Phys
B22:397–411. https://doi.org/10.1016/0550-3213(70)90416-5
Darbha S, Coughlin ER, Kasen D, Quataert E (2018) Gravitational interactions of stars with supermassive
black hole binaries - I. Tidal disruption events. MNRAS 477:4009–4034. https://doi.org/10.1093/
mnras/sty822. arXiv:1802.07850
Deller AT, Verbiest JPW, Tingay SJ, Bailes M (2008) Extremely high precision VLBI astrometry of PSR
J0437–4715 and implications for theories of gravity. ApJ 685:L67. https://doi.org/10.1086/592401.
arXiv:0808.1594
123
Deller AT, Goss WM, Brisken WF, Chatterjee S, Cordes JM, Janssen GH, Kovalev YY, Lazio TJW, Petrov
L, Stappers BW, Lyne A (2019) Microarcsecond VLBI pulsar astrometry with PSRπ II. Parallax
distances for 57 pulsars. ApJ 875:100. https://doi.org/10.3847/1538-4357/ab11c7. arXiv:1808.09046
Demorest PB, Pennucci T, Ransom SM, Roberts MSE, Hessels JWT (2010) A two-solar-mass neu-
tron star measured using Shapiro delay. Nature 467:1081–1083. https://doi.org/10.1038/nature09466.
arXiv:1010.5788
Demorest PB, Ferdman RD, Gonzalez ME, Nice D, Ransom S, Stairs IH, Arzoumanian Z, Brazier A,
Burke-Spolaor S, Chamberlin SJ, Cordes JM, Ellis J, Finn LS, Freire P, Giampanis S, Jenet F, Kaspi
VM, Lazio J, Lommen AN, McLaughlin M, Palliyaguru N, Perrodin D, Shannon RM, Siemens X,
Stinebring D, Swiggum J, Zhu WW (2013) Limits on the stochastic gravitational wave background
from the North American nanohertz observatory for gravitational waves. ApJ 762:94. https://doi.org/
10.1088/0004-637X/762/2/94. arXiv:1201.6641
Desvignes G, Caballero RN, Lentati L, Verbiest JPW, Champion DJ, Stappers BW, Janssen GH, Lazarus
P, Osłowski S, Babak S, Bassa CG, Brem P, Burgay M, Cognard I, Gair JR, Graikou E, Guillemot
L, Hessels JWT, Jessner A, Jordan C, Karuppusamy R, Kramer M, Lassus A, Lazaridis K, Lee KJ,
Liu K, Lyne AG, McKee J, Mingarelli CMF, Perrodin D, Petiteau A, Possenti A, Purver MB, Rosado
PA, Sanidas S, Sesana A, Shaifullah G, Smits R, Taylor SR, Theureau G, Tiburzi C, van Haasteren R,
Vecchio A (2016) High-precision timing of 42 millisecond pulsars with the European Pulsar Timing
Array. MNRAS 458:3341–3380. https://doi.org/10.1093/mnras/stw483. arXiv:1602.08511
Detweiler S (1979) Pulsar timing measurements and the search for gravitational waves. ApJ 234:1100–1104.
https://doi.org/10.1086/157593
Dolch T, NANOGrav Collaboration, Ellis JA, Chatterjee S, Cordes JM, Lam MT, Bassa C, Bhattacharyya
B, Champion DJ, Cognard I, Crowter K, Demorest PB, Hessels JWT, Janssen G, Jenet FA, Jones G,
Jordan C, Karuppusamy R, Keith M, Kondratiev VI, Kramer M, Lazarus P, Lazio TJW, Lorimer DR,
Madison DR, McLaughlin MA, Palliyaguru N, Perrodin D, Ransom SM, Roy J, Shannon RM, Smits
R, Stairs IH, Stappers BW, Stinebring DR, Stovall K, Verbiest JPW, Zhu WW (2016) Single-source
gravitational wave limits from the J1713+0747 24-hr global campaign. J Phys Conf Ser 716:012014.
https://doi.org/10.1088/1742-6596/716/1/012014. arXiv:1509.05446
D’Orazio DJ, Di Stefano R (2018) Periodic self-lensing from accreting massive black hole binaries. MNRAS
474:2975–2986. https://doi.org/10.1093/mnras/stx2936. arXiv:1707.02335
D’Orazio DJ, Haiman Z (2017) Lighthouse in the dust: infrared echoes of periodic emission from
massive black hole binaries. MNRAS 470:1198–1217. https://doi.org/10.1093/mnras/stx1269.
arXiv:1702.01219
D’Orazio DJ, Haiman Z, Duffell P, Farris BD, MacFadyen AI (2015a) A reduced orbital period for the
supermassive black hole binary candidate in the quasar PG 1302–102? MNRAS 452:2540–2545.
https://doi.org/10.1093/mnras/stv1457. arXiv:1502.03112
D’Orazio DJ, Haiman Z, Schiminovich D (2015b) Relativistic boost as the cause of periodicity in
a massive black-hole binary candidate. Nature 525:351–353. https://doi.org/10.1038/nature15262.
arXiv:1509.04301
Dosopoulou F, Antonini F (2017) Dynamical friction and the evolution of supermassive black hole
binaries: the final hundred-parsec problem. ApJ 840:31. https://doi.org/10.3847/1538-4357/aa6b58.
arXiv:1611.06573
Drake AJ, Djorgovski SG, Mahabal A, Beshore E, Larson S, Graham MJ, Williams R, Christensen E,
Catelan M, Boattini A, Gibbs A, Hill R, Kowalski R (2009) First results from the catalina real-time
transient survey. ApJ 696:870–884. https://doi.org/10.1088/0004-637X/696/1/870. arXiv:0809.1394
Dror JA, Ramani H, Trickle T, Zurek KM (2019) Pulsar timing probes of primordial black holes and
subhalos. ArXiv e-prints arXiv:1901.04490
Dvorkin I, Barausse E (2017) The nightmare scenario: measuring the stochastic gravitational wave back-
ground from stalling massive black hole binaries with pulsar timing arrays. MNRAS 470(4):4547–4556
Eardley DM, Lee DL, Lightman AP (1973a) Gravitational-wave observations as a tool for testing relativistic
gravity. Phys Rev D 8:3308–3321. https://doi.org/10.1103/PhysRevD.8.3308
Eardley DM, Lee DL, Lightman AP, Wagoner RV, Will CM (1973b) Gravitational-wave observations as a
tool for testing relativistic gravity. Phys Rev Lett 30:884–886. https://doi.org/10.1103/PhysRevLett.
30.884
Ebisuzaki T, Makino J, Okumura SK (1991) Merging of two galaxies with central black holes. Nature
354:212–214. https://doi.org/10.1038/354212a0
123
Ellis JA (2013) A Bayesian analysis pipeline for continuous GW sources in the PTA band. Class Quantum
Grav 30(22):224004. https://doi.org/10.1088/0264-9381/30/22/224004. arXiv:1305.0835
Ellis JA, Siemens X, Creighton JDE (2012) Optimal strategies for continuous gravitational wave detection in
pulsar timing arrays. ApJ 756:175. https://doi.org/10.1088/0004-637X/756/2/175. arXiv:1204.4218
Emsellem E (2013) Is the black hole in NGC 1277 really overmassive? MNRAS 433:1862–1870. https://
doi.org/10.1093/mnras/stt840. arXiv:1305.3630
Enoki M, Nagashima M, Nagashima M (2007) The effect of orbital eccentricity on gravitational wave
background radiation from supermassive black hole binaries. Prog Theor Phys 117(2):241–256
Eracleous M, Boroson TA, Halpern JP, Liu J (2012) A large systematic search for close supermassive
binary and rapidly recoiling black holes. ApJS 201:23. https://doi.org/10.1088/0067-0049/201/2/23.
arXiv:1106.2952
Erickcek AL, Smith TL, Kamionkowski M (2006) Solar system tests do rule out 1/R gravity. Phys Rev D
74(121):501. https://doi.org/10.1103/PhysRevD.74.121501. arXiv:astro-ph/0610483
Espriu D, Puigdomènech D (2013) Local measurement of using pulsar timing arrays. ApJ 764:163.
Faber SM, Tremaine S, Ajhar EA, Byun YI, Dressler A, Gebhardt K, Grillmair C, Kormendy J, Lauer TR,
Richstone D (1997) The centers of early-type galaxies with HST. IV. Central parameter relations. AJ
114:1771. https://doi.org/10.1086/118606. arXiv:astro-ph/9610055
Favata M (2010) The gravitational-wave memory effect. Class Quantum Gravity 27(8):084036. https://doi.
org/10.1088/0264-9381/27/8/084036. arXiv:1003.3486
Favata M (2011) The gravitational-wave memory from eccentric binaries. Phys Rev D 84(12):124013.
https://doi.org/10.1103/PhysRevD.84.124013. arXiv:1108.3121
Ferrarese L, Merritt D (2000) A fundamental relation between supermassive black holes and their host
galaxies. ApJ 539:L9–L12. https://doi.org/10.1086/312838. arXiv:astro-ph/0006053
Finn LS, Lommen AN (2010) Detection, localization, and characterization of gravitational wave
bursts in a pulsar timing array. ApJ 718:1400–1415. https://doi.org/10.1088/0004-637X/718/2/1400.
arXiv:1004.3499
Fonseca E, Pennucci TT, Ellis JA, Stairs IH, Nice DJ, Ransom SM, Demorest PB, Arzoumanian Z, Crowter
K, Dolch T, Ferdman RD, Gonzalez ME, Jones G, Jones ML, Lam MT, Levin L, McLaughlin MA,
Stovall K, Swiggum JK, Zhu W (2016) The NANOGrav nine-year data set: mass and geometric
measurements of binary millisecond pulsars. ApJ 832:167. https://doi.org/10.3847/0004-637X/832/
2/167. arXiv:1603.00545
Frank J, Rees MJ (1976) Effects of massive central black holes on dense stellar systems. MNRAS
176(3):633–647
Freire PCC, Ridolfi A, Kramer M, Jordan C, Manchester RN, Torne P, Sarkissian J, Heinke CO, D’Amico
N, Camilo F, Lorimer DR, Lyne AG (2017) Long-term observations of the pulsars in 47 Tucanae
- II. Proper motions, accelerations and jerks. MNRAS 471:857–876. https://doi.org/10.1093/mnras/
stx1533. arXiv:1706.04908
Collaboration Gaia, Prusti T, de Bruijne JHJ, Brown AGA, Vallenari A, Babusiaux C, Bailer-Jones CAL,
Bastian U, Biermann M, Evans DW et al (2016) The Gaia mission. A&A 595:A1. https://doi.org/10.
1051/0004-6361/201629272. arXiv:1609.04153
Collaboration Gaia, Brown AGA, Vallenari A, Prusti T, de Bruijne JHJ, Babusiaux C, Bailer-Jones CAL,
Biermann M, Evans DW, Eyer L et al (2018) Gaia data release 2. Summary of the contents and survey
properties. A&A 616:A1. https://doi.org/10.1051/0004-6361/201833051. arXiv:1804.09365
Gair J, Romano JD, Taylor S, Mingarelli CMF (2014) Mapping gravitational-wave backgrounds using
methods from CMB analysis: application to pulsar timing arrays. Phys Rev D 90(8):082001. https://
doi.org/10.1103/PhysRevD.90.082001. arXiv:1406.4664
Gair JR, Romano JD, Taylor SR (2015) Mapping gravitational-wave backgrounds of arbitrary polarisa-
tion using pulsar timing arrays. Phys Rev D 92(10):102,003. https://doi.org/10.1103/PhysRevD.92.
102003. arXiv:1506.08668
García-Bellido J, Peloso M, C Unal (2016) Gravitational waves at interferometer scales and primordial
black holes in axion inflation. JCAP 2016(12):031. https://doi.org/10.1088/1475-7516/2016/12/031.
arXiv:1610.03763
García-Bellido J, Peloso M (2017) Unal C (2017) Gravitational wave signatures of inflationary models
from primordial black hole dark matter. JCAP 09:013. https://doi.org/10.1088/1475-7516/2017/09/
013. arXiv:1707.02441
123
Gebhardt K, Bender R, Bower G, Dressler A, Faber SM, Filippenko AV, Green R, Grillmair C, Ho LC,
Kormendy J, Lauer TR, Magorrian J, Pinkney J, Richstone D, Tremaine S (2000) A relationship
between nuclear black hole mass and galaxy velocity dispersion. ApJ 539:L13–L16. https://doi.org/
10.1086/312840. arXiv:astro-ph/0006289
Goicovic FG, Sesana A, Cuadra J (2017) Infalling clouds on to supermassive black hole binaries—II. Binary
evolution and the final parsec problem. MNRAS
Goldreich P, Sari R (2003) Eccentricity evolution for planets in gaseous disks. ApJ 585(2):1024–1037
Graham AW, Driver SP (2007) A log-quadratic relation for predicting supermassive black hole masses from
the host Bulge Sérsic index. ApJ 655:77–87. https://doi.org/10.1086/509758. arXiv:astro-ph/0607378
Graham MJ, Djorgovski SG, Stern D, Drake AJ, Mahabal AA, Donalek C, Glikman E, Larson S, Christensen
E (2015a) A systematic search for close supermassive black hole binaries in the catalina real-time
transient survey. MNRAS 453:1562–1576. https://doi.org/10.1093/mnras/stv1726. arXiv:1507.07603
Graham MJ, Djorgovski SG, Stern D, Glikman E, Drake AJ, Mahabal AA, Donalek C, Larson S, Christensen
E (2015b) A possible close supermassive black-hole binary in a quasar with optical periodicity. Nature
518:74–76. https://doi.org/10.1038/nature14143. arXiv:1501.01375
Gualandris A, Merritt D (2008) Ejection of supermassive black holes from galaxy cores. ApJ 678:780–797.
https://doi.org/10.1086/586877. arXiv:0708.0771
Gualandris A, Dotti M, Sesana A (2012) Massive black hole binary plane reorientation in rotating stellar sys-
tems. MNRAS 420:L38–L42. https://doi.org/10.1111/j.1745-3933.2011.01188.x. arXiv:1109.3707
Guo YJ, Lee KJ, Caballero RN (2018) A dynamical approach in exploring the unknown mass in the Solar
system using pulsar timing arrays. MNRAS 475:3644–3653. https://doi.org/10.1093/mnras/stx3326.
arXiv:1802.05452
Gupta S, Desai S (2018) Limit on graviton mass using stacked galaxy cluster catalogs from SPT-SZ,
Planck-SZ and SDSS-redMaPPer. Ann Phys 399:85–92. https://doi.org/10.1016/j.aop.2018.09.017.
arXiv:1810.00198
Haiman Z, Kocsis B, Menou K (2009) The population of viscosity- and gravitational wave-driven super-
massive black hole binaries among luminous active galactic nuclei. ApJ 700(2):1952–1969
Harrison ER (1977) Has the Sun a companion Star. Nature 270:324–326. https://doi.org/10.1038/270324a0
Hassan SF, Rosen RA (2012a) Bimetric gravity from ghost-free massive gravity. JHEP 02:126. https://doi.
org/10.1007/JHEP02(2012)126. arXiv:1109.3515
Hassan SF, Rosen RA (2012b) Resolving the ghost problem in non-linear massive gravity. Phys Rev Lett
108(041):101. https://doi.org/10.1103/PhysRevLett.108.041101. arXiv:1106.3344
Hayasaki K, Mineshige S, Sudou H (2007) Binary black hole accretion flows in merged galactic nuclei.
PASJ 59:427–441. https://doi.org/10.1093/pasj/59.2.427. arXiv:astro-ph/0609144
Hazboun JS, Marcano MP, Larson SL (2013) Limiting alternative theories of gravity using gravitational
wave observations across the spectrum. ArXiv e-prints arXiv:1311.3153
Hellings RW, Downs GS (1983) Upper limits on the isotropic gravitational radiation background from
pulsar timing analysis. ApJ 265:L39–L42. https://doi.org/10.1086/183954
Henrichs HF, Staller RFA (1978) Has the sun really got a companion star. Nature 273:132–134. https://doi.
org/10.1038/273132a0
Hobbs G, Archibald A, Arzoumanian Z, Backer D, Bailes M, Bhat NDR, Burgay M, Burke-Spolaor S,
Champion D, Cognard I, Coles W, Cordes J, Demorest P, Desvignes G, Ferdman RD, Finn L, Freire
P, Gonzalez M, Hessels J, Hotan A, Janssen G, Jenet F, Jessner A, Jordan C, Kaspi V, Kramer M,
Kondratiev V, Lazio J, Lazaridis K, Lee KJ, Levin Y, Lommen A, Lorimer D, Lynch R, Lyne A,
Manchester R, McLaughlin M, Nice D, Oslowski S, Pilia M, Possenti A, Purver M, Ransom S,
Reynolds J, Sanidas S, Sarkissian J, Sesana A, Shannon R, Siemens X, Stairs I, Stappers B, Stinebring
D, Theureau G, van Haasteren R, van Straten W, Verbiest JPW, Yardley DRB, You XP (2010) The
international pulsar timing array project: using pulsars as a gravitational wave detector. Class Quantum
Grav 27(8):084013. https://doi.org/10.1088/0264-9381/27/8/084013. arXiv:0911.5206
Hobbs G, Coles W, Manchester RN, Keith MJ, Shannon RM, Chen D, Bailes M, Bhat NDR, Burke-Spolaor
S, Champion D, Chaudhary A, Hotan A, Khoo J, Kocz J, Levin Y, Oslowski S, Preisig B, Ravi V,
Reynolds JE, Sarkissian J, van Straten W, Verbiest JPW, Yardley D, You XP (2012) Development of a
pulsar-based time-scale. MNRAS 427:2780–2787. https://doi.org/10.1111/j.1365-2966.2012.21946.
x. arXiv:1208.3560
Holgado AM, Sesana A, Sandrinelli A, Covino S, Treves A, Liu X, Ricker P (2018) Pulsar timing constraints
on the Fermi massive black hole binary blazar population. MNRAS 481:L74–L78. https://doi.org/10.
1093/mnrasl/sly158. arXiv:1806.11111
123
Hu W, Sawicki I (2007) Models of f (R) cosmic acceleration that evade solar-system tests. Phys Rev D
76(064):004. https://doi.org/10.1103/PhysRevD.76.064004. arXiv:0705.1158
Hu W, Barkana R, Gruzinov A (2000) Fuzzy cold dark matter: the wave properties of ultralight particles.
Phys Rev Lett 85:1158–1161. https://doi.org/10.1103/PhysRevLett.85.1158. arXiv:astro-ph/0003365
Huerta E, McWilliams ST, Gair JR, Taylor S (2015) Detection of eccentric supermassive black hole binaries
with pulsar timing arrays: signal-to-noise ratio calculations. Phys Rev D 92(6):063010. https://doi.
org/10.1103/PhysRevD.92.063010. arXiv:1504.00928
Hui L, Ostriker JP, Tremaine S, Witten E (2017) Ultralight scalars as cosmological dark matter. Phys Rev
D 95:043541. https://doi.org/10.1103/PhysRevD.95.043541. arXiv:1610.08297
Ishiyama T, Makino J, Ebisuzaki T (2010) Gamma-ray signal from Earth-mass dark matter microhalos.
ApJ 723:L195–L200. https://doi.org/10.1088/2041-8205/723/2/L195. arXiv:1006.3392
Isi M, Stein LC (2018) Measuring stochastic gravitational-wave energy beyond general relativity. Phys Rev
Isi M, Weinstein AJ (2017) Probing gravitational wave polarizations with signals from compact binary
coalescences. ArXiv e-prints arXiv:1710.03794
Isi M, Weinstein AJ, Mead C, Pitkin M (2015) Detecting beyond-Einstein polarizations of continu-
ous gravitational waves. Phys Rev D 91(8):082002. https://doi.org/10.1103/PhysRevD.91.082002.
arXiv:1502.00333
Ivanov PB, Papaloizou JCB, Polnarev AG (1999) The evolution of a supermassive binary caused by an
accretion disc. MNRAS 307(1):79–90
Jacobson T, Mattingly D (2001) Gravity with a dynamical preferred frame. Phys Rev D 64(024):028. https://
doi.org/10.1103/PhysRevD.64.024028. arXiv:gr-qc/0007031
Jacobson T, Mattingly D (2004) Einstein-Aether waves. Phys Rev D 70(024):003. https://doi.org/10.1103/
PhysRevD.70.024003. arXiv:gr-qc/0402005
Janssen G, Hobbs G, McLaughlin M, Bassa C, Deller A, Kramer M, Lee K, Mingarelli C, Rosado P,
Sanidas S, Sesana A, Shao L, Stairs I, Stappers B, Verbiest JPW (2015) Gravitational wave astronomy
with the SKA. In: Advancing astrophysics with the Square Kilometre Array (AASKA14), PoS, p 37,
arXiv:1501.00127
Jenet FA, Lommen A, Larson SL, Wen L (2004) Constraining the properties of supermassive black hole
systems using pulsar timing: application to 3C 66B. ApJ 606:799. https://doi.org/10.1086/383020.
Jenet FA, Creighton T, Lommen A (2005) Pulsar timing and the detection of black hole binary systems in
globular clusters. ApJ 627:L125–L128. https://doi.org/10.1086/431949. arXiv:astro-ph/0505585
Jennings RJ, Kaplan DL, Chatterjee S, Cordes JM, Deller AT (2018) Binary pulsar distances and velocities
from gaia data release 2. ApJ 864:26. https://doi.org/10.3847/1538-4357/aad084. arXiv:1806.06076
Jones ML, McLaughlin MA, Lam MT, Cordes JM, Levin L, Chatterjee S, Arzoumanian Z, Crowter K,
Demorest PB, Dolch T, Ellis JA, Ferdman RD, Fonseca E, Gonzalez ME, Jones G, Lazio TJW, Nice
DJ, Pennucci TT, Ransom SM, Stinebring DR, Stairs IH, Stovall K, Swiggum JK, Zhu WW (2017)
The NANOGrav nine-year data set: measurement and analysis of variations in dispersion measures.
ApJ 841:125. https://doi.org/10.3847/1538-4357/aa73df. arXiv:1612.03187
Ju W, Greene JE, Rafikov RR, Bickerton SJ, Badenes C (2013) Search for supermassive black hole binaries
in the sloan digital sky survey spectroscopic sample. ApJ 777:44. https://doi.org/10.1088/0004-637X/
777/1/44. arXiv:1306.4987
Kaastra JS, Roos N (1992) Massive binary black-holes and wiggling jets. A&A 254:96
Kaiser N, Burgett W, Chambers K, Denneau L, Heasley J, Jedicke R, Magnier E, Morgan J, Onaka P, Tonry
J (2010) The Pan-STARRS wide-field optical/NIR imaging survey. In: Stepp LM, Gilmozzi R, Hall HJ
(eds) Ground-based and Airborne Telescopes III, Bellingham, WA, Proc. SPIE, vol 7733, p 77330E,
https://doi.org/10.1117/12.859188
Kamionkowski M, Kovetz ED (2016) The quest for B modes from inflationary gravitational waves. ARA&A
54:227–269. https://doi.org/10.1146/annurev-astro-081915-023433. arXiv:1510.06042
Kaplan DL, Kupfer T, Nice DJ, Irrgang A, Heber U, Arzoumanian Z, Beklen E, Crowter K, DeCesar ME,
Demorest PB, Dolch T, Ellis JA, Ferdman RD, Ferrara EC, Fonseca E, Gentile PA, Jones G, Jones ML,
Kreuzer S, Lam MT, Levin L, Lorimer DR, Lynch RS, McLaughlin MA, Miller AA, Ng C, Pennucci
TT, Prince TA, Ransom SM, Ray PS, Spiewak R, Stairs IH, Stovall K, Swiggum J, Zhu W (2016)
PSR J1024–0719: a millisecond pulsar in an unusual long-period orbit. ApJ 826:86. https://doi.org/
10.3847/0004-637X/826/1/86. arXiv:1604.00131
123
Kashiyama K, Oguri M (2018) Detectability of small-scale dark matter clumps with pulsar timing arrays.
ArXiv e-prints arXiv:1801.07847
Kashiyama K, Seto N (2012) Enhanced exploration for primordial black holes using pulsar timing arrays.
MNRAS 426:1369–1373. https://doi.org/10.1111/j.1365-2966.2012.21935.x. arXiv:1208.4101
Keith MJ, Coles W, Shannon RM, Hobbs GB, Manchester RN, Bailes M, Bhat NDR, Burke-Spolaor S,
Champion DJ, Chaudhary A, Hotan AW, Khoo J, Kocz J, Osłowski S, Ravi V, Reynolds JE, Sarkissian
J, van Straten W, Yardley DRB (2013) Measurement and correction of variations in interstellar disper-
sion in high-precision pulsar timing. MNRAS 429:2161–2174. https://doi.org/10.1093/mnras/sts486.
arXiv:1211.5887
Kelley LZ, Blecha L, Hernquist L (2017a) Massive black hole binary mergers in dynamical galactic envi-
ronments. MNRAS 464:3131–3157. https://doi.org/10.1093/mnras/stw2452. arXiv:1606.01900
Kelley LZ, Blecha L, Hernquist L, Sesana A, Taylor S (2017b) The gravitational wave background from
massive black hole binaries in illustris: spectral features and time to detection with pulsar timing
arrays. arXivorg p arXiv:1702.02180
Khan FM, Holley-Bockelmann K, Berczik P, Just A (2013) Supermassive black hole binary evolution in
axisymmetric galaxies: the final parsec problem is not a problem. ApJ 773:100. https://doi.org/10.
1088/0004-637X/773/2/100. arXiv:1302.1871
Kharb P, Lal DV, Merritt D (2017) A candidate sub-parsec binary black hole in the Seyfert galaxy NGC
7674. Nat Astron 1:727–733. https://doi.org/10.1038/s41550-017-0256-4. arXiv:1709.06258
Khmelnitsky A, Rubakov V (2014) Pulsar timing signal from ultralight scalar dark matter. J Cosmol
Astropart Phys 2:019. https://doi.org/10.1088/1475-7516/2014/02/019. arXiv:1309.5888
Kocsis B, Sesana A (2011) Gas-driven massive black hole binaries: signatures in the nHz gravitational wave
background. MNRAS 411(3):1467–1479
Kocsis B, Ray A, Zwart SP (2012) Mapping the galactic center with gravitational wave measurements using
pulsar timing. ApJ
Komossa S, Burwitz V, Hasinger G, Predehl P, Kaastra JS, Ikebe Y (2003) Discovery of a binary active
galactic nucleus in the ultraluminous infrared galaxy NGC 6240 using chandra. ApJ 582:L15–L19.
https://doi.org/10.1086/346145. arXiv:astro-ph/0212099
Kormendy J, Richstone D (1995) Inward bound-the search for supermassive black holes in galactic nuclei.
ARA&A 33:581. https://doi.org/10.1146/annurev.aa.33.090195.003053
Kostelecky VA, Samuel S (1989) Spontaneous breaking of Lorentz symmetry in string theory. Phys Rev D
39:683. https://doi.org/10.1103/PhysRevD.39.683
Kozai Y (1962) Secular perturbations of asteroids with high inclination and eccentricity. AJ 67:591. https://
doi.org/10.1086/108790
Kramer M, Stairs IH (2008) The double pulsar. ARA&A 46:541–572. https://doi.org/10.1146/annurev.
astro.46.060407.145247
Kramer M, Stairs IH, Manchester RN, McLaughlin MA, Lyne AG, Ferdman RD, Burgay M, Lorimer DR,
Possenti A, D’Amico N, Sarkissian JM, Hobbs GB, Reynolds JE, Freire PCC, Camilo F (2006) Tests of
general relativity from timing the double pulsar. Science 314:97–102. https://doi.org/10.1126/science.
1132305. arXiv:astro-ph/0609417
Kremer K, Chatterjee S, Breivik K, Rodriguez CL, Larson SL, Rasio FA (2018) LISA sources in Milky Way
globular clusters. Phys Rev Lett 120(19):191103. https://doi.org/10.1103/PhysRevLett.120.191103.
arXiv:1802.05661
Kulier A, Ostriker JP, Natarajan P, Lackner CN, Cen R (2015) Understanding black hole mass assembly
via accretion and mergers at late times in cosmological simulations. ApJ 799(2):178
Lam MT, Cordes JM, Chatterjee S, Jones ML, McLaughlin MA, Armstrong JW (2016) Systematic and
stochastic variations in pulsar dispersion measures. ApJ 821:66. https://doi.org/10.3847/0004-637X/
821/1/66. arXiv:1512.02203
Lam MT, Ellis JA, Grillo G, Jones ML, Hazboun JS, Brook PR, Turner JE, Chatterjee S, Cordes JM, Lazio
TJW, DeCesar ME, Arzoumanian Z, Blumer H, Cromartie HT, Demorest PB, Dolch T, Ferdman RD,
Ferrara EC, Fonseca E, Garver-Daniels N, Gentile PA, Gupta V, Lorimer DR, Lynch RS, Madison
DR, McLaughlin MA, Ng C, Nice DJ, Pennucci TT, Ransom SM, Spiewak R, Stairs IH, Stinebring
DR, Stovall K, Swiggum JK, Vigeland SJ, Zhu WW (2018) A second chromatic timing event of
interstellar origin toward PSR J1713+0747. ApJ 861:132. https://doi.org/10.3847/1538-4357/aac770.
arXiv:1712.03651
Lasky PD, Mingarelli CMF, Smith TL, Giblin JT, Thrane E, Reardon DJ, Caldwell R, Bailes M, Bhat NDR,
Burke-Spolaor S, Dai S, Dempsey J, Hobbs G, Kerr M, Levin Y, Manchester RN, Osłowski S, Ravi
123
V, Rosado PA, Shannon RM, Spiewak R, van Straten W, Toomey L, Wang J, Wen L, You X, Zhu
X (2016) Gravitational-wave cosmology across 29 decades in frequency. Phys Rev X 6(1):011035.
https://doi.org/10.1103/PhysRevX.6.011035. arXiv:1511.05994
Lattanzi M, Benini R, Montani G (2010) A possible signature of cosmic neutrino decoupling in the nHz
region of the spectrum of primordial gravitational waves. Class Quantum Gravity 27(19):194008.
Lee K, Jenet FA, Price RH, Wex N, Kramer M (2010a) Detecting massive gravitons using pul-
sar timing arrays. Astrophys J 722:1589–1597. https://doi.org/10.1088/0004-637X/722/2/1589.
arXiv:1008.2561
Lee K, Jenet FA, Price RH, Wex N, Kramer M (2010b) Detecting massive gravitons using pulsar timing
arrays. ApJ 722:1589–1597. https://doi.org/10.1088/0004-637X/722/2/1589. arXiv:1008.2561
Lee KJ, Jenet FA, Price RH (2008) Pulsar timing as a probe of non-Einsteinian polarizations of gravitational
waves. ApJ 685:1304–1319. https://doi.org/10.1086/591080
Lee KJ, Wex N, Kramer M, Stappers BW, Bassa CG, Janssen GH, Karuppusamy R, Smits R (2011a)
Gravitational wave astronomy of single sources with a pulsar timing array. MNRAS 414:3251–3264.
Lee KJ, Wex N, Kramer M, Stappers BW, Bassa CG, Janssen GH, Karuppusamy R, Smits R (2011b)
Gravitational wave astronomy of single sources with a pulsar timing array. MNRAS 414:3251–3264.
Lentati L, Taylor SR, Mingarelli CMF, Sesana A, Sanidas SA, Vecchio A, Caballero RN, Lee KJ, van
Haasteren R, Babak S, Bassa CG, Brem P, Burgay M, Champion DJ, Cognard I, Desvignes G, Gair
JR, Guillemot L, Hessels JWT, Janssen GH, Karuppusamy R, Kramer M, Lassus A, Lazarus P, Liu
K, Osłowski S, Perrodin D, Petiteau A, Possenti A, Purver MB, Rosado PA, Smits R, Stappers B,
Theureau G, Tiburzi C, Verbiest JPW (2015) European pulsar timing array limits on an isotropic
stochastic gravitational-wave background. MNRAS 453:2576–2598. https://doi.org/10.1093/mnras/
stv1538. arXiv:1504.03692
Li L, Guo L, Wang GL (2016) Detecting the errors in solar system ephemeris by pulsar timing. Res Astron
Astrophys 16:58. https://doi.org/10.1088/1674-4527/16/4/058
Liang D, Gong Y, Hou S, Liu Y (2017) Polarizations of gravitational waves in f (r ) gravity. Phys Rev D
95(104):034. https://doi.org/10.1103/PhysRevD.95.104034
Liberati S (2013) Tests of Lorentz invariance: a 2013 update. Class Quantum Grav 30(133):001. https://
doi.org/10.1088/0264-9381/30/13/133001. arXiv:1304.5795
Lidov ML (1962) The evolution of orbits of artificial satellites of planets under the action of gravita-
tional perturbations of external bodies. Planet Space Sci 9:719–759. https://doi.org/10.1016/0032-
0633(62)90129-0
Liu FK, Li S, Komossa S (2014) A milliparsec supermassive black hole binary candidate in the
galaxy SDSS J120136.02+300305.5. ApJ 786:103. https://doi.org/10.1088/0004-637X/786/2/103.
arXiv:1404.4933
Liu K, Keane EF, Lee KJ, Kramer M, Cordes JM, Purver MB (2012) Profile-shape stability and phase-jitter
analyses of millisecond pulsars. MNRAS 420:361–368. https://doi.org/10.1111/j.1365-2966.2011.
20041.x. arXiv:1110.4759
Liu T, Gezari S, Heinis S, Magnier EA, Burgett WS, Chambers K, Flewelling H, Huber M, Hodapp KW,
Kaiser N, Kudritzki RP, Tonry JL, Wainscoat RJ, Waters C (2015) 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. ApJ 803:L16. https://doi.org/10.1088/2041-8205/803/
2/L16. arXiv:1503.02083
Liu T, Gezari S, Burgett W, Chambers K, Draper P, Hodapp K, Huber M, Kudritzki RP, Magnier E, Metcalfe
N, Tonry J, Wainscoat R, Waters C (2016) A systematic search for periodically varying quasars in
pan-STARRS1: an extended baseline test in medium deep survey field MD09. ApJ 833:6. https://doi.
org/10.3847/0004-637X/833/1/6. arXiv:1609.09503
Liu T, Gezari S, Miller MC (2018) Did ASAS-SN kill the supermassive black hole binary candidate
PG1302-102? ApJ 859:L12. https://doi.org/10.3847/2041-8213/aac2ed. arXiv:1803.05448
Lombriser L (2016) Taylor A (2016) Breaking a dark degeneracy with gravitational waves. JCAP 03:031.
Lommen AN, Backer DC (2001) Using pulsars to detect massive black hole binaries via gravitational
radiation: Sagittarius A* and nearby galaxies. ApJ 562:297–302. https://doi.org/10.1086/323491.
123
Lorimer DR, Kramer M (2012) Handbook of pulsar astronomy. Cambridge University Press, Cambridge
LSST Science Collaboration, Marshall P, Anguita T, Bianco FB, Bellm EC, Brandt N, Clarkson W, Connolly
A, Gawiser E, Ivezic Z, Jones L, Lochner M, Lund MB, Mahabal A, Nidever D, Olsen K, Ridgway
S, Rhodes J, Shemmer O, Trilling D, Vivas K, Walkowicz L, Willman B, Yoachim P, Anderson S,
Antilogus P, Angus R, Arcavi I, Awan H, Biswas R, Bell KJ, Bennett D, Britt C, Buzasi D, Casetti-
Dinescu DI, Chomiuk L, Claver C, Cook K, Davenport J, Debattista V, Digel S, Doctor Z, Firth RE,
Foley R, Fong Wf, Galbany L, Giampapa M, Gizis JE, Graham ML, Grillmair C, Gris P, Haiman Z,
Hartigan P, Hawley S, Hlozek R, Jha SW, Johns-Krull C, Kanbur S, Kalogera V, Kashyap V, Kasliwal
V, Kessler R, Kim A, Kurczynski P, Lahav O, Liu MC, Malz A, Margutti R, Matheson T, McEwen JD,
McGehee P, Meibom S, Meyers J, Monet D, Neilsen E, Newman J, O’Dowd M, Peiris HV, Penny MT,
Peters C, Poleski R, Ponder K, Richards G, Rho J, Rubin D, Schmidt S, Schuhmann RL, Shporer A,
Slater C, Smith N, Soares-Santos M, Stassun K, Strader J, Strauss M, Street R, Stubbs C, Sullivan M,
Szkody P, Trimble V, Tyson T, de Val-Borro M, Valenti S, Wagoner R, Wood-Vasey WM, Zauderer BA
(2017) Science-driven optimization of the LSST observing strategy. ArXiv e-prints arXiv:1708.04058
Lyne A, Graham-Smith F (2012) Pulsar astronomy. Cambridge University Press, Cambridge
Madison DR, Cordes JM, Chatterjee S (2014) Assessing pulsar timing array sensitivity to gravitational wave
bursts with memory. ApJ 788:141. https://doi.org/10.1088/0004-637X/788/2/141. arXiv:1404.5682
Madison DR, Chernoff DF, Cordes JM (2017) Pulsar timing perturbations from Galactic gravitational
wave bursts with memory. Phys Rev D 96(12):123016. https://doi.org/10.1103/PhysRevD.96.123016.
arXiv:1710.04974
Maggiore M (2007) Gravitational waves. Vol. 1: Theory and experiments. Oxford Master Series in Physics,
Oxford University Press, Oxford
Magorrian J, Tremaine S, Richstone D, Bender R, Bower G, Dressler A, Faber SM, Gebhardt K, Green R,
Grillmair C, Kormendy J, Lauer T (1998) The demography of massive dark objects in galaxy centers.
AJ 115:2285–2305. https://doi.org/10.1086/300353. arXiv:astro-ph/9708072
Makino J (1997) Merging of galaxies with central black holes. II. Evolution of the black hole binary and
the structure of the core. ApJ 478:58–65. https://doi.org/10.1086/303773. arXiv:astro-ph/9608161
Makino J, Ebisuzaki T (1994) Triple black holes in the cores of galaxies. Astrophys J 436:607–610
Manchester RN, Hobbs G, Bailes M, Coles WA, van Straten W, Keith MJ, Shannon RM, Bhat NDR, Brown
A, Burke-Spolaor SG, Champion DJ, Chaudhary A, Edwards RT, Hampson G, Hotan AW, Jameson
A, Jenet FA, Kesteven MJ, Khoo J, Kocz J, Maciesiak K, Oslowski S, Ravi V, Reynolds JR, Sarkissian
JM, Verbiest JPW, Wen ZL, Wilson WE, Yardley D, Yan WM, You XP (2013) The parkes pulsar
timing array project. PASA 30:e017. https://doi.org/10.1017/pasa.2012.017. arXiv:1210.6130
Matthews AM, Nice DJ, Fonseca E, Arzoumanian Z, Crowter K, Demorest PB, Dolch T, Ellis JA, Ferdman
RD, Gonzalez ME, Jones G, Jones ML, Lam MT, Levin L, McLaughlin MA, Pennucci TT, Ransom
SM, Stairs IH, Stovall K, Swiggum JK, Zhu W (2016) The NANOGrav nine-year data set: astrometric
measurements of 37 millisecond pulsars. ApJ 818:92. https://doi.org/10.3847/0004-637X/818/1/92.
arXiv:1509.08982
McConnell NJ, Ma CP (2013) Revisiting the scaling relations of black hole masses and host galaxy prop-
erties. ApJ 764:184. https://doi.org/10.1088/0004-637X/764/2/184. arXiv:1211.2816
McKee JW, Lyne AG, Stappers BW, Bassa CG, Jordan CA (2018) Temporal variations in scattering and
dispersion measure in the Crab Pulsar and their effect on timing precision. MNRAS 479:4216–4224.
https://doi.org/10.1093/mnras/sty1727. arXiv:1806.10486
McWilliams ST, Ostriker JP, Pretorius F (2014) Gravitational waves and stalled satellites from massive
galaxy mergers at z ≤ 1. ApJ 789(2):156
Merloni A, Predehl P, Becker W, Böhringer H, Boller T, Brunner H, Brusa M, Dennerl K, Freyberg M,
Friedrich P, Georgakakis A, Haberl F, Hasinger G, Meidinger N, Mohr J, Nandra K, Rau A, Reiprich
TH, Robrade J, Salvato M, Santangelo A, Sasaki M, Schwope A, Wilms J, German eROSITA Consor-
tium t (2012) eROSITA science book: mapping the structure of the energetic universe. ArXiv e-prints
arXiv:1209.3114
Merritt D (2013) Dynamics and evolution of galactic nuclei. Princeton University Press, Princeton
Merritt D, Milosavljević M (2005) Massive black hole binary evolution. Living Rev Relativ 8:8.
10.12942/lrr-2005-8, arXiv:astro-ph/0410364
Mihaylov DP, Moore CJ, Gair JR, Lasenby A, Gilmore G (2018) Astrometric effects of gravitational wave
backgrounds with non-Einsteinian polarizations. ArXiv e-prints arXiv:1804.00660
Mikkola S, Valtonen MJ (1992) Evolution of binaries in the field of light particles and the problem of two
black holes. MNRAS 259(1):115–120
123
Milosavljević M, Merritt D (2001) Formation of galactic nuclei. ApJ 563:34–62. https://doi.org/10.1086/

323830. arXiv:astro-ph/0103350
Milosavljević M, Merritt D (2003a) Long-term evolution of massive black hole binaries. ApJ 596(2):860–
878
Milosavljević M, Merritt D (2003b) The final parsec problem. In: Centrella JM, Barnes S (eds) The astro-
physics of gravitational wave sources. American Institute of Physics, Melville, NY, AIP Conf. Proc.,
vol 686, pp 201–210, https://doi.org/10.1063/1.1629432
Milosavljević M, Phinney ES (2005) The afterglow of massive black hole coalescence. ApJ 622:L93–L96.
https://doi.org/10.1086/429618. arXiv:astro-ph/0410343
Mingarelli CMF, Mingarelli AB (2018) Proving the short-wavelength approximation in pulsar timing array
gravitational-wave background searches. J Phys Commun 2(10):105,002
Mingarelli CMF, Sidery T (2014) Effect of small interpulsar distances in stochastic gravitational wave
background searches with pulsar timing arrays. Phys Rev D 90(6):062011. https://doi.org/10.1103/
PhysRevD.90.062011. arXiv:1408.6840
Mingarelli CMF, Grover K, Sidery T, Smith RJE, Vecchio A (2012) Observing the dynamics of supermassive
black hole binaries with pulsar timing arrays. Phys Rev Lett 109(8):081104. https://doi.org/10.1103/
PhysRevLett.109.081104. arXiv:1207.5645
Mingarelli CMF, Sidery T, Mandel I, Vecchio A (2013) Characterizing gravitational wave stochastic back-
ground anisotropy with pulsar timing arrays. Phys Rev D 88(6):062005. https://doi.org/10.1103/
PhysRevD.88.062005. arXiv:1306.5394
Mingarelli CMF, Lazio TJW, Sesana A, Greene JE, Ellis JA, Ma CP, Croft S, Burke-Spolaor S, Taylor SR
(2017) The local nanohertz gravitational-wave landscape from supermassive black hole binaries. Nat
Astron 1:886–892. https://doi.org/10.1038/s41550-017-0299-6. arXiv:1708.03491
Mingarelli CMF, Anderson L, Bedell M, Spergel DN (2018) Improving binary millisecond pulsar distances
with gaia. ArXiv e-prints arXiv:1812.06262
Mirza MA, Tahir A, Khan FM, Holley-Bockelmann H, Baig AM, Berczik P, Chishtie F (2017) Galaxy
rotation and supermassive black hole binary evolution. MNRAS 470:940–947. https://doi.org/10.
1093/mnras/stx1248. arXiv:1704.03490
Moore CJ, Mihaylov D, Lasenby A, Gilmore G (2017) An astrometric search method for individually
resolvable gravitational wave sources with Gaia. ArXiv e-prints arXiv:1707.06239
Mulholland JD (1971) The system of planetary masses as error sources in pulsar timings. ApJ 165:105.
https://doi.org/10.1086/150879
Namba R, Peloso M, Shiraishi M, Sorbo L, Unal C (2016) Scale-dependent gravitational waves from a rolling
axion. JCAP 2016(01):041. https://doi.org/10.1088/1475-7516/2016/01/041. arXiv:1509.07521
NANOGrav Collaboration (2018) The NANOGrav 11-year data set: limits on gravitational waves from
individual supermassive black hole binaries. submitted
Nishizawa A, Taruya A, Hayama K, Kawamura S, Sakagami MA (2009) Probing nontensorial polarizations
of stochastic gravitational-wave backgrounds with ground-based laser interferometers. Phys Rev D
79:082002. https://doi.org/10.1103/PhysRevD.79.082002. arXiv:0903.0528
Nordtvedt K Jr, Will CM (1972) Conservation laws and preferred frames in relativistic gravity. II. Experi-
mental evidence to rule out preferred-frame theories of gravity. ApJ 177:775. https://doi.org/10.1086/
151755
O’Beirne L, Cornish NJ (2018) Constraining the polarization content of gravitational waves with astrometry.
ArXiv e-prints arXiv:1804.03146
Perera BBP, Stappers BW, Babak S, Keith MJ, Antoniadis J, Bassa CG, Caballero RN, Champion DJ,
Cognard I, Desvignes G, Graikou E, Guillemot L, Janssen GH, Karuppusamy R, Kramer M, Lazarus
P, Lentati L, Liu K, Lyne AG, McKee JW, Osłowski S, Perrodin D, Sanidas SA, Sesana A, Shaifullah G,
Theureau G, Verbiest JPW, Taylor SR (2018) Improving timing sensitivity in the microhertz frequency
regime: limits from PSR J1713+0747 on gravitational waves produced by supermassive black hole
binaries. MNRAS 478:218–227 arXiv:1804.10571
Peters PC (1964) Gravitational radiation and the motion of two point masses. Phys Rev 136:1224–1232.
https://doi.org/10.1103/PhysRev.136.B1224
Peters PC, Mathews J (1963) Gravitational radiation from point masses in a Keplerian Orbit. Phys Rev
131:435–440. https://doi.org/10.1103/PhysRev.131.435
Phinney ES (2001) A practical theorem on gravitational wave backgrounds. ArXiv e-prints
123
Pitkin M, Clark J, Hendry MA, Heng IS, Messenger C, Toher J, Woan G (2008) Is there potential com-
plementarity between LISA and pulsar timing? J Phys Conf Ser 122:012004. https://doi.org/10.1088/
1742-6596/122/1/012004. arXiv:0802.2460
Collaboration Planck, Adam R, Ade PAR, Aghanim N, Akrami Y, Alves MIR, Argüeso F, Arnaud M, Arroja
F, Ashdown M et al (2016) Planck 2015 results. I. Overview of products and scientific results. A&A
594:A1. https://doi.org/10.1051/0004-6361/201527101. arXiv:1502.01582
Pleunis Z, Bassa CG, Hessels JWT, Kondratiev VI, Camilo F, Cognard I, Grießmeier JM, Stappers
BW, van Amesfoort AS, Sanidas S (2017) A millisecond pulsar discovery in a survey of uniden-
tified Fermi γ -Ray sources with LOFAR. ApJ 846:L19. https://doi.org/10.3847/2041-8213/aa83ff.
arXiv:1709.01452
Porayko NK, Postnov KA (2014) Constraints on ultralight scalar dark matter from pulsar timing. Phys Rev
Porayko NK, Zhu X, Levin Y, Hui L, Hobbs G, Grudskaya A, Postnov K, Bailes M, Bhat NDR, Coles W, Dai
S, Dempsey J, Keith MJ, Kerr M, Kramer M, Lasky PD, Manchester RN, Osłowski S, Parthasarathy
A, Ravi V, Reardon DJ, Rosado PA, Russell CJ, Shannon RM, Spiewak R, van Straten W, Toomey L,
Wang J, Wen L, You X, Collaboration PPTA (2018) Parkes pulsar timing array constraints on ultra-
light scalar-field dark matter. Phys Rev D 98:102002. https://doi.org/10.1103/PhysRevD.98.102002.
arXiv:1810.03227
Pshirkov MS, Baskaran D, Postnov KA (2010) Observing gravitational wave bursts in pulsar timing measure-
ments. MNRAS 402:417–423. https://doi.org/10.1111/j.1365-2966.2009.15887.x. arXiv:0909.0742
Quinlan GD (1996) The dynamical evolution of massive black hole binaries I. Hardening in a fixed stellar
background. New Astron 1(1):35–56
Rajagopal M, Romani RW (1995) Ultra-low-frequency gravitational radiation from massive black hole
binaries. ApJ 446:543. https://doi.org/10.1086/175813. arXiv:astro-ph/9412038
Ransom SM, Hessels JWT, Stairs IH, Freire PCC, Camilo F, Kaspi VM, Kaplan DL (2005) Twenty-one
millisecond pulsars in terzan 5 using the green bank telescope. Science 307:892–896. https://doi.org/
10.1126/science.1108632. arXiv:astro-ph/0501230
Rasskazov A, Merritt D (2017a) Evolution of binary supermassive black holes in rotating nuclei. ApJ
837(2):135
Rasskazov A, Merritt D (2017b) Evolution of massive black hole binaries in rotating stellar nuclei: Implica-
tions for gravitational wave detection. Phys Rev D 95(8):084032. https://doi.org/10.1103/PhysRevD.
95.084032. arXiv:1606.07484
Ravi V, Wyithe JSB, Shannon RM, Hobbs GB, Manchester RN (2014) Binary supermassive black hole
environments diminish the gravitational wave signal in the pulsar timing band. MNRAS 442(1):56–68
Ravi V, Wyithe JSB, Shannon RM, Hobbs G (2015) Prospects for gravitational-wave detection and super-
massive black hole astrophysics with pulsar timing arrays. MNRAS 447:2772–2783. https://doi.org/
10.1093/mnras/stu2659. arXiv:1406.5297
Reardon DJ, Hobbs G, Coles W, Levin Y, Keith MJ, Bailes M, Bhat NDR, Burke-Spolaor S, Dai S, Kerr
M, Lasky PD, Manchester RN, Osłowski S, Ravi V, Shannon RM, van Straten W, Toomey L, Wang J,
Wen L, You XP, Zhu XJ (2016) Timing analysis for 20 millisecond pulsars in the parkes pulsar timing
array. MNRAS 455:1751–1769. https://doi.org/10.1093/mnras/stv2395. arXiv:1510.04434
de Rham C, Gabadadze G, Tolley AJ (2011) Resummation of massive gravity. Phys Rev Lett 106(231):101.
https://doi.org/10.1103/PhysRevLett.106.231101. arXiv:1011.1232
Rodriguez C, Taylor GB, Zavala RT, Peck AB, Pollack LK, Romani RW (2006) A compact supermassive
binary black hole system. ApJ 646:49–60. https://doi.org/10.1086/504825. arXiv:astro-ph/0604042
Rodriguez-Gomez V, Genel S, Vogelsberger M, Sijacki D, Pillepich A, Sales LV, Torrey P, Snyder G, Nelson
D, Springel V, Ma CP, Hernquist L (2015) The merger rate of galaxies in the illustris simulation: a
comparison with observations and semi-empirical models. MNRAS 1:17–64 arXiv:1502.01339
Roebber E, Holder G (2017) Harmonic space analysis of pulsar timing array redshift maps. ApJ 835:21.
Roedig C, Sesana A (2012) Origin and implications of high eccentricities in massive black hole binaries at
sub-pc scales. J Phys Conf Ser 363(1):012,035
Roedig C, Dotti M, Sesana A, Cuadra J, Colpi M, Colpi M (2011) Limiting eccentricity of subparsec massive
black hole binaries surrounded by self-gravitating gas discs. MNRAS 415(4):3033–3041
Romano JD, Cornish NJ (2017) Detection methods for stochastic gravitational-wave backgrounds: a unified
treatment. Living Rev Relat 20:2. https://doi.org/10.1007/s41114-017-0004-1. arXiv:1608.06889
123
Rosado PA, Sesana A, Gair J (2015) Expected properties of the first gravitational wave signal
detected with pulsar timing arrays. MNRAS 451:2417–2433. https://doi.org/10.1093/mnras/stv1098.
arXiv:1503.04803
Rubinur K, Das M, Kharb P (2018) Searching for dual active galactic nuclei. J Astrophys Astron 39:8.
https://doi.org/10.1007/s12036-018-9512-y. arXiv:1801.03451
Ryu T, Perna R, Haiman Z, Ostriker JP, Stone NC (2018) Interactions between multiple supermassive black
holes in galactic nuclei: a solution to the final parsec problem. MNRAS 473:3410–3433. https://doi.
org/10.1093/mnras/stx2524. arXiv:1709.06501
Sakstein J, Jain B (2017) Implications of the neutron star merger GW170817 for cosmological scalar-
tensor theories. Phys Rev Lett 119(25):251,303. https://doi.org/10.1103/PhysRevLett.119.251303.
arXiv:1710.05893
Schnittman JD, Krolik JH (2015) Evjisk. arXivorg arXiv:1504.00311v1
Schutz BF (2010) Astrometric and timing effects of gravitational waves. In: Klioner SA, Seidelmann PK,
Soffel MH (eds) Relativity in fundamental astronomy: dynamics, reference frames, and data analysis.
Cambridge University Press, Cambridge, IAU Symposium, vol 261, pp 234–239, https://doi.org/10.
1017/S1743921309990457
Schutz K, Ma CP (2016) Constraints on individual supermassive black hole binaries from pulsar timing
array limits on continuous gravitational waves. MNRAS 459:1737–1744. https://doi.org/10.1093/
mnras/stw768. arXiv:1510.08472
Sesana A (2010) Self consistent model for the evolution of eccentric massive black hole binaries in stellar
environments: implications for Gravitational Wave Observations. ApJ 719(1):851–864
Sesana A, Haardt F, Madau P, Volonteri M (2004) Low-frequency gravitational radiation from coalescing
massive black hole binaries in hierarchical cosmologies. ApJ 611:623–632. https://doi.org/10.1086/
422185. arXiv:astro-ph/0401543
Sesana A, Haardt F, Madau P (2006) Interaction of massive black hole binaries with their stellar environment.
I. Ejection of hypervelocity stars. ApJ 651(1):392–400
Sesana A, Vecchio A, Colacino CN (2008) The stochastic gravitational-wave background from mas-
sive black hole binary systems: implications for observations with Pulsar Timing Arrays. MNRAS
390(1):192–209
Sesana A, Vecchio A, Volonteri M (2009) Gravitational waves from resolvable massive black hole binary
systems and observations with Pulsar Timing Arrays. MNRAS 394:2255–2265. https://doi.org/10.
1111/j.1365-2966.2009.14499.x. arXiv:0809.3412
Sesana A, Gualandris A, Dotti M (2011) Massive black hole binary eccentricity in rotating stellar systems.
MNRAS: Lett 415(1):L35–L39
Sesana A, Haiman Z, Kocsis B, Kelley LZ (2018) Testing the binary hypothesis: pulsar timing constraints on
supermassive black hole binary candidates. ApJ 856:42. https://doi.org/10.3847/1538-4357/aaad0f
Seto N (2009) Search for memory and inspiral gravitational waves from supermassive binary black holes
with pulsar timing arrays. MNRAS 400:L38–L42. https://doi.org/10.1111/j.1745-3933.2009.00758.
x. arXiv:0909.1379
Seto N, Cooray A (2007) Searching for primordial black hole dark matter with pulsar timing arrays. ApJ
659:L33–L36. https://doi.org/10.1086/516570. arXiv:astro-ph/0702586
Shakura NI, Sunyaev RA (1973) Black holes in binary systems. Observational appearance. A&A 24:337–
355
Shankar F, Bernardi M, Sheth RK, Ferrarese L, Graham AW, Savorgnan G, Allevato V, Marconi A, Läsker
R, Lapi A (2016) Selection bias in dynamically measured supermassive black hole samples: its con-
sequences and the quest for the most fundamental relation. MNRAS 460:3119–3142. https://doi.org/
10.1093/mnras/stw678. arXiv:1603.01276
Shannon RM, Ravi V, Lentati LT, Lasky PD, Hobbs G, Kerr M, Manchester RN, Coles WA, Levin Y, Bailes
M, Bhat NDR, Burke-Spolaor S, Dai S, Keith MJ, Osłowski S, Reardon DJ, van Straten W, Toomey
L, Wang JB, Wen L, Wyithe JSB, Zhu XJ (2015) Gravitational waves from binary supermassive
black holes missing in pulsar observations. Science 349:1522–1525. https://doi.org/10.1126/science.
aab1910. arXiv:1509.07320
Shapiro SL, Teukolsky SA (1986) Black holes, white dwarfs and neutron stars: the physics of compact
objects. Wiley-VCH, Weinheim
Shen Y, Greene JE, Strauss MA, Richards GT, Schneider DP (2008) Biases in virial black hole masses: an
SDSS perspective. ApJ 680:169–190. https://doi.org/10.1086/587475. arXiv:0709.3098
Shklovskii IS (1970) Possible causes of the secular increase in pulsar periods. Soviet Ast. 13:562
123
Siegel ER, Hertzberg MP, Fry JN (2007) Probing dark matter substructure with pulsar timing. MNRAS
382:879–885. https://doi.org/10.1111/j.1365-2966.2007.12435.x. arXiv:astro-ph/0702546
Siemens X, Ellis J, Jenet F, Romano JD (2013) The stochastic background: scaling laws and time to detection
for pulsar timing arrays. Class Quantum Grav 30(22):224,015
Simon J, Burke-Spolaor S (2016) Constraints on black hole/host galaxy co-evolution and binary
stalling using pulsar timing arrays. ApJ 826:11. https://doi.org/10.3847/0004-637X/826/1/11.
arXiv:1603.06577
Simon J, Polin A, Lommen A, Stappers B, Finn LS, Jenet FA, Christy B (2014) Gravitational wave hotspots:
ranking potential locations of single-source gravitational wave emission. ApJ 784:60. https://doi.org/
10.1088/0004-637X/784/1/60. arXiv:1402.1140
Smits R, Tingay SJ, Wex N, Kramer M, Stappers B (2011) Prospects for accurate distance measurements
of pulsars with the Square Kilometre Array: enabling fundamental physics. A&A 528:A108. https://
doi.org/10.1051/0004-6361/201016141. arXiv:1101.5971
Spallicci ADAM (2013) On the complementarity of pulsar timing and space laser interferometry for the
individual detection of supermassive black hole binaries. ApJ 764:187. https://doi.org/10.1088/0004-
637X/764/2/187. arXiv:1107.5984
Sudou H, Iguchi S, Murata Y, Taniguchi Y (2003) Orbital motion in the radio galaxy 3C 66B: evidence
for a supermassive black hole binary. Science 300:1263. https://doi.org/10.1126/science.1082817.
Takahashi R (2017) Arrival time differences between gravitational waves and electromagnetic signals due to
gravitational lensing. ApJ 835:103. https://doi.org/10.3847/1538-4357/835/1/103. arXiv:1606.00458
Talmadge C, Berthias JP, Hellings RW, Standish EM (1988) Model independent constraints on possible mod-
ifications of Newtonian gravity. Phys Rev Lett 61:1159–1162. https://doi.org/10.1103/PhysRevLett.
61.1159
Tanaka T, Menou K (2010) Time-dependent models for the afterglows of massive black hole mergers. ApJ
714:404–422. https://doi.org/10.1088/0004-637X/714/1/404. arXiv:0912.2054
Taylor JH, Weisberg JM (1982) A new test of general relativity: gravitational radiation and the binary pulsar
PSR 1913+16. ApJ 253:908–920. https://doi.org/10.1086/159690
Taylor S, Ellis J, Gair J (2014) Accelerated Bayesian model-selection and parameter-estimation in contin-
uous gravitational-wave searches with pulsar-timing arrays. Phys Rev D 90(10):104028. https://doi.
org/10.1103/PhysRevD.90.104028. arXiv:1406.5224
Taylor S, Huerta E, Gair JR, McWilliams ST (2016) Detecting eccentric supermassive black hole binaries
with pulsar timing arrays: resolvable source strategies. ApJ 817(1):70
Taylor SR, Gair JR (2013) Searching for anisotropic gravitational-wave backgrounds using pulsar timing
arrays. Phys Rev D 88(8):084001. https://doi.org/10.1103/PhysRevD.88.084001. arXiv:1306.5395
Taylor SR, Mingarelli CMF, Gair JR, Sesana A, Theureau G, Babak S, Bassa CG, Brem P, Burgay M,
Caballero RN, Champion DJ, Cognard I, Desvignes G, Guillemot L, Hessels JWT, Janssen GH,
Karuppusamy R, Kramer M, Lassus A, Lazarus P, Lentati L, Liu K, Osłowski S, Perrodin D, Petiteau
A, Possenti A, Purver MB, Rosado PA, Sanidas SA, Smits R, Stappers B, Tiburzi C, van Haasteren R,
Vecchio A, Verbiest JPW, Collaboration EPTA (2015) Limits on anisotropy in the nanohertz stochastic
gravitational wave background. Phys Rev Lett 115(4):041101. https://doi.org/10.1103/PhysRevLett.
115.041101. arXiv:1506.08817
Taylor SR, Vallisneri M, Ellis JA, Mingarelli CMF, Lazio TJW, van Haasteren R (2016) Are we there yet?
Time to detection of nanohertz gravitational waves based on pulsar-timing array limits. ApJ 819:L6.
https://doi.org/10.3847/2041-8205/819/1/L6. arXiv:1511.05564
Taylor SR, Lentati L, Babak S, Brem P, Gair JR, Sesana A, Vecchio A (2017a) All correlations must die:
assessing the significance of a stochastic gravitational-wave background in pulsar timing arrays. Phys
Rev D 95(4):042002. https://doi.org/10.1103/PhysRevD.95.042002. arXiv:1606.09180
Taylor SR, Simon J, Sampson L (2017b) Constraints on the dynamical environments of supermassive
black-hole binaries using pulsar-timing arrays. Phys Rev Lett 118(18):181102. https://doi.org/10.
1103/PhysRevLett.118.181102. arXiv:1612.02817
The WFIRST Astrometry Working Group, Sanderson RE, Bellini A, Casertano S, Lu JR, Melchior P,
Bennett D, Shao M, Rhodes J, Malhotra S, Gaudi S, Fall M, Nelan E, Guhathakurta P, Anderson J, Ho S,
Libralato M (2017) Astrometry with the WFIRST wide-field imager. ArXiv e-prints arXiv:1712.05420
Thorne KS (1992) Gravitational-wave bursts with memory: the Christodoulou effect. Phys. Rev. D 45:520–
524. https://doi.org/10.1103/PhysRevD.45.520
123
Tiburzi C, Hobbs G, Kerr M, Coles WA, Dai S, Manchester RN, Possenti A, Shannon RM, You XP (2016)
A study of spatial correlations in pulsar timing array data. MNRAS 455:4339–4350. https://doi.org/
10.1093/mnras/stv2143. arXiv:1510.02363
Tremmel M, Governato F, Volonteri M, Quinn TR, Pontzen A (2018) Dancing to CHANGA: a self-consistent
prediction for close SMBH pair formation time-scales following galaxy mergers. MNRAS 475:4967–
4977. https://doi.org/10.1093/mnras/sty139. arXiv:1708.07126
van den Bosch RCE, Gebhardt K, Gültekin K, van de Ven G, van der Wel A, Walsh JL (2012) An over-
massive black hole in the compact lenticular galaxy NGC 1277. Nature 491:729–731. https://doi.org/
10.1038/nature11592. arXiv:1211.6429
van Haasteren R, Levin Y (2010) Gravitational-wave memory and pulsar timing arrays. MNRAS 401:2372–
2378. https://doi.org/10.1111/j.1365-2966.2009.15885.x. arXiv:0909.0954
Vasiliev E, Merritt D (2013) The loss-cone problem in axisymmetric nuclei. ApJ 774(1):87
Vasiliev E, Antonini F, Merritt D (2014) The final-parsec problem in nonspherical galaxies revisited. ApJ
785(2):163
Vasiliev E, Antonini F, Merritt D (2015) The final-parsec problem in the collisionless limit. ApJ 810(1):49
Vaughan S, Uttley P, Markowitz AG, Huppenkothen D, Middleton MJ, Alston WN, Scargle JD, Farr WM
(2016) False periodicities in quasar time-domain surveys. MNRAS 461:3145–3152. https://doi.org/
10.1093/mnras/stw1412. arXiv:1606.02620
Verbiest JPW, Lentati L, Hobbs G, van Haasteren R, Demorest PB, Janssen GH, Wang JB, Desvignes G,
Caballero RN, Keith MJ, Champion DJ, Arzoumanian Z, Babak S, Bassa CG, Bhat NDR, Brazier A,
Brem P, Burgay M, Burke-Spolaor S, Chamberlin SJ, Chatterjee S, Christy B, Cognard I, Cordes JM,
Dai S, Dolch T, Ellis JA, Ferdman RD, Fonseca E, Gair JR, Garver-Daniels NE, Gentile P, Gonzalez
ME, Graikou E, Guillemot L, Hessels JWT, Jones G, Karuppusamy R, Kerr M, Kramer M, Lam MT,
Lasky PD, Lassus A, Lazarus P, Lazio TJW, Lee KJ, Levin L, Liu K, Lynch RS, Lyne AG, Mckee J,
McLaughlin MA, McWilliams ST, Madison DR, Manchester RN, Mingarelli CMF, Nice DJ, Osłowski
S, Palliyaguru NT, Pennucci TT, Perera BBP, Perrodin D, Possenti A, Petiteau A, Ransom SM, Reardon
D, Rosado PA, Sanidas SA, Sesana A, Shaifullah G, Shannon RM, Siemens X, Simon J, Smits R,
Spiewak R, Stairs IH, Stappers BW, Stinebring DR, Stovall K, Swiggum JK, Taylor SR, Theureau G,
Tiburzi C, Toomey L, Vallisneri M, van Straten W, Vecchio A, Wang Y, Wen L, You XP, Zhu WW, Zhu
XJ (2016) The international pulsar timing array: first data release. MNRAS 458:1267–1288. https://
doi.org/10.1093/mnras/stw347. arXiv:1602.03640
Vigeland SJ, Siemens X (2016) Supermassive black hole binary environments: effects on the scaling laws
and time to detection for the stochastic background. Phys Rev D 94(12):123003. https://doi.org/10.
1103/PhysRevD.94.123003. arXiv:1609.03656
Volonteri M, Haardt F, Madau P (2003a) The assembly and merging history of supermassive black
holes in hierarchical models of galaxy formation. ApJ 582:559–573. https://doi.org/10.1086/344675.
Volonteri M, Madau P, Haardt F (2003b) The formation of galaxy stellar cores by the hierar-
chical merging of supermassive black holes. ApJ 593:661–666. https://doi.org/10.1086/376722.
Wang JB, Hobbs G, Coles W, Shannon RM, Zhu XJ, Madison DR, Kerr M, Ravi V, Keith MJ, Manchester
RN, Levin Y, Bailes M, Bhat NDR, Burke-Spolaor S, Dai S, Osłowski S, van Straten W, Toomey L,
Wang N, Wen L (2015) Searching for gravitational wave memory bursts with the parkes pulsar timing
array. MNRAS 446:1657–1671. https://doi.org/10.1093/mnras/stu2137. arXiv:1410.3323
Wegg C, Nate Bode J (2011) Multiple tidal disruptions as an indicator of binary supermassive black hole
systems. ApJ 738:L8. https://doi.org/10.1088/2041-8205/738/1/L8. arXiv:1011.5874
Wen L (2003) On the eccentricity distribution of coalescing black hole binaries driven by the Kozai mecha-
nism in globular clusters. ApJ 598:419–430. https://doi.org/10.1086/378794. arXiv:astro-ph/0211492
Wen ZL, Liu FS, Han JL (2009) Mergers of luminous early-type galaxies in the local universe and
gravitational wave background. ApJ 692:511–521. https://doi.org/10.1088/0004-637X/692/1/511.
arXiv:0810.5200
Will CM (1993) Theory and experiment in gravitational physics. Cambridge University Press, Cambridge
Will CM (1998) Bounding the mass of the graviton using gravitational wave observations of inspi-
ralling compact binaries. Phys Rev D 57:2061–2068. https://doi.org/10.1103/PhysRevD.57.2061.
arXiv:gr-qc/9709011
Will CM (2014) The confrontation between general relativity and experiment. Living Rev Rel 17:14. https://
doi.org/10.12942/lrr-2014-4. arXiv:1403.7377
123
Will CM, Nordtvedt K Jr (1972) Conservation laws and preferred frames in relativistic gravity. I. Preferred-
frame theories and an extended PPN formalism. ApJ 177:757. https://doi.org/10.1086/151754
Wyithe JSB, Loeb A (2003) Low-frequency gravitational waves from massive black hole binaries: pre-
dictions for LISA and pulsar timing arrays. ApJ 590:691–706. https://doi.org/10.1086/375187.
Yao JM, Manchester RN, Wang N (2017) A new electron-density model for estimation of pulsar and FRB
distances. ApJ 835:29. https://doi.org/10.3847/1538-4357/835/1/29. arXiv:1610.09448
Yardley DRB, Hobbs GB, Jenet FA, Verbiest JPW, Wen ZL, Manchester RN, Coles WA, van Straten W,
Bailes M, Bhat NDR, Burke-Spolaor S, Champion DJ, Hotan AW, Sarkissian JM (2010) The sensitivity
of the Parkes Pulsar Timing array to individual sources of gravitational waves. MNRAS 407:669–680.
Yi S, Stappers BW, Sanidas SA, Bassa CG, Janssen GH, Lyne AG, Kramer M, Zhang SN (2014) Limits
on the strength of individual gravitational wave sources using high-cadence observations of PSR
B1937+21. MNRAS 445:1245–1252. https://doi.org/10.1093/mnras/stu1826. arXiv:1409.2296
Young JS, Scoville NZ (1991) Molecular gas in galaxies. ARA&A 29:581–625. https://doi.org/10.1146/
annurev.aa.29.090191.003053
Yu Q (2002) Evolution of massive binary black holes. MNRAS 331(4):935–958
Yunes N, Pretorius F (2009) Fundamental theoretical bias in gravitational wave astrophysics and the param-
eterized post-Einsteinian framework. Phys Rev D 80(122):003. https://doi.org/10.1103/PhysRevD.
80.122003. arXiv:0909.3328
Yunes N, Siemens X (2013) Gravitational-wave tests of general relativity with ground-based detectors and
pulsar timing-arrays. Living Rev Rel 16:9. https://doi.org/10.12942/lrr-2013-9. arXiv:1304.3473
Zakamska NL, Tremaine S (2005) Constraints on the acceleration of the solar system from high-precision
timing. AJ 130:1939–1950. https://doi.org/10.1086/444476. arXiv:astro-ph/0506548
Zakharov VI (1970) Linearized gravitation theory and the graviton mass. JETP Lett 12:312 [Pisma Zh Eksp
Teor Fiz 12:447 (1970)]
Zel’dovich YB, Polnarev AG (1974) Radiation of gravitational waves by a cluster of superdense stars.
Soviet Ast 18:17
Zheng ZY, Butler NR, Shen Y, Jiang L, Wang JX, Chen X, Cuadra J (2016) SDSS J0159+0105: a radio-quiet
quasar with a centi-parsec supermassive black hole binary candidate. ApJ 827:56. https://doi.org/10.
3847/0004-637X/827/1/56. arXiv:1512.08730
Zhu WW, Stairs IH, Demorest PB, Nice DJ, Ellis JA, Ransom SM, Arzoumanian Z, Crowter K, Dolch T,
Ferdman RD, Fonseca E, Gonzalez ME, Jones G, Jones ML, Lam MT, Levin L, McLaughlin MA,
Pennucci T, Stovall K, Swiggum J (2015) Testing theories of gravitation using 21-year timing of pulsar
binary J1713+0747. ApJ 809:41. https://doi.org/10.1088/0004-637X/809/1/41. arXiv:1504.00662
Zhu XJ, Hobbs G, Wen L, Coles WA, Wang JB, Shannon RM, Manchester RN, Bailes M, Bhat NDR,
Burke-Spolaor S, Dai S, Keith MJ, Kerr M, Levin Y, Madison DR, Osłowski S, Ravi V, Toomey L, van
Straten W (2014) An all-sky search for continuous gravitational waves in the Parkes Pulsar Timing
Array data set. MNRAS 444:3709–3720. https://doi.org/10.1093/mnras/stu1717. arXiv:1408.5129
Zhu XJ, Wen L, Xiong J, Xu Y, Wang Y, Mohanty SD, Hobbs G, Manchester RN (2016) Detection and
localization of continuous gravitational waves with pulsar timing arrays: the role of pulsar terms.
MNRAS 461:1317–1327. https://doi.org/10.1093/mnras/stw1446. arXiv:1606.04539
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps
and institutional affiliations.
123
Aﬃliations
Sarah Burke-Spolaor1,2,3 · Stephen R. Taylor4,5 · Maria Charisi4 ·

Timothy Dolch6 · Jeﬀrey S. Hazboun7 · A. Miguel Holgado8 ·
Luke Zoltan Kelley9,10 · T. Joseph W. Lazio11 · Dustin R. Madison1,12 ·
Natasha McMann13,14 · Chiara M. F. Mingarelli15 · Alexander Rasskazov16 ·
Xavier Siemens17 · Joseph J. Simon11 · Tristan L. Smith18
B Sarah Burke-Spolaor
[email protected]
B Stephen R. Taylor
1 Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown,
WV 26506, USA
2 Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA
3 Canadian Institute for Advanced Research, CIFAR Azrieli Global Scholar, MaRS Centre West
Tower, 661 University Ave. Suite 505, Toronto, ON M5G 1M1, Canada
4 TAPIR, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
5 Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
6 Department of Physics, Hillsdale College, 33 E. College St., Hillsdale, MI 49242, USA
7 Physical Sciences Division, University of Washington Bothell, 18115 Campus Way NE, Bothell,
WA 98011-8246, USA
8 Department of Astronomy and National Center for Supercomputing Applications, University of
Illinois at Urbana-Champaign, 1002 W. Green St., Urbana, IL 61801, USA
9 Harvard–Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
10 Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA
11 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena,
CA 91109, USA
12 The National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903,
USA
13 Department of Life and Physical Sciences, Fisk University, 1000 17th Ave N., Nashville, TN
37208, USA
14 Department of Physics and Astronomy, Vanderbilt University, PMB 401807, Nashville, TN
37206, USA
15 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010,
USA
16 Eötvös University, Institute of Physics, Pázmány P. s. 1/A, Budapest 1117, Hungary
17 Department of Physics, Center for Gravitation, Cosmology and Astrophysics, University of
Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
18 Department of Physics and Astronomy, Swarthmore College, Swarthmore, PA 19081, USA
123

The Astrophysics of Nanohertz Gravitational Waves

Uploaded by

Copyright:

Available Formats

The Astrophysics of Nanohertz Gravitational Waves

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

The Astrophysics of Nanohertz Gravitational Waves

Uploaded by

Copyright:

Available Formats

The Astronomy and Astrophysics Review (2019) 27:5

The astrophysics of nanohertz gravitational waves

Received: 8 November 2018 / Published online: 18 June 2019

Keywords Gravitational waves · Stars · Neutron · Galaxies · Evolution · Black hole

3.2.1 The mass function of supermassive black holes . . . . . . . . . . . . . . . . . . . . . 19

The detection of nanohertz GWs by Pulsar Timing Arrays (PTAs) is expected to be

2 Pulsar timing in brief

2.1 Pulsar timing and timing residuals

We time a pulsar by building a “timing model”, which is a mathematical description

(a) Continuous wave (b) Background

Simulated post-fit residuals influenced by a variety of GW signals, which PTAs

where h i j is the space–time perturbation (typically referred to as the “strain”), and

2.3 Types of gravitational-wave signal

RMS Precision (ns)

3 The population and evolution of supermassive black hole binaries

SMBHBs are expected to be the brightest sources of nanohertz GWs. They

Throughout this document, we refer to SMBHB parameters using the following

3.1 GW emission from supermassive black hole binaries

3.1.1 PTAs and the binary life cycle

3.1.2 Continuous waves: binary inspiral

where e0 is the eccentricity at some reference epoch, and,

3.1.3 Memory: binary coalescence

3.1.4 The stochastic binary gravitational-wave background

outline how many discrete continuous-wave sources can contribute to a GW back-

where the contributing factors are:

(1963). Thus, non-zero eccentricity in a binary redistributes the power of the GW

3.1.5 Gravitational-wave background anisotropy

3.2 Supermassive binaries and their environments

3.2.1 The mass function of supermassive black holes

The GW amplitude depends strongly on the masses of SMBHB components: h ∝

3.2.2 Dynamical friction

(resulting in a gross overestimate of the dynamical friction timescale) considers a point

3.2.3 Stellar loss-cone scattering

At parsec separations, dynamical friction becomes an inefficient means of further

where H is a dimensionless hardening rate, and K is a dimensionless eccentricity

3.2.4 Viscous circumbinary disk interaction

At centiparsec to milliparsec separations, viscous dissipation of angular momentum

6 stellar density ρ(r ) ∝ r −1 , and constant velocity dispersion σ (r ) = σ .

where α is a disk viscosity parameter, Ṁ E = 4π Gm 1 /cκT is the Eddington accretion

The influence of an initial eccentricity on SMBH evolution, without the influence

3.2.6 Measuring the spin of supermassive black hole binaries

pulsar terms Gravitational waves

4 Cosmic strings and cosmic superstrings

Fig. 8 Simulation of a cosmic

Fig. 9 Plot of regions of the

Fig. 10 Plot of the gravitational

5 The nature of gravity

Understanding gravity—and testing general relativity as a theory of gravity—

5.1 How many gravitational wave polarization states exist?

(a) (c) (e)

(b) (d) (f)

5.2 Characterizing the dispersion relation of gravitational waves

5.3 Examples of gravity theories that predict additional polarizations and

5.3.1 Einstein Aether

Einstein Aether posits the existence of a Lorentz symmetry-violating, dynamical,

5.3.2 Massive gravity

which is a linear combination of the breathing and longitudinal modes.

5.3.3 f (R) Gravity

Another commonly considered alternative theory of gravity imagines a gravitational