The astrophysics of nanohertz gravitational waves

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.

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 beginning of gravitational waves’ contribution to the era of multi-messenger astronomy (Abbott et al. 2017b). The ground-based laser interferometers that made those detections (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.

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 interferometers (Amaro-Seoane et al. 2017) and observations of primordial GWs, imprinted in the polarization of the cosmic microwave background (e. g., BICEP2/Keck Collaboration et al. 2015), providing comprehensive access to the GW Universe. Current PTA efforts are spearheaded by a number of groups worldwide, including the European 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.

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

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 $\sim$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.

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 $10^9\,M_\odot$ 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_\mathrm{c}=10^{-15}$ and $\alpha =-2/3$. c A memory event of $h=5\times 10^{-15}$, whose wavefront passes the Earth on day 1500. d A burst source with an arbitrary waveform

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

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 approximate 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_{\mathrm{obs}}=(1+z)f_{\mathrm{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:

$$\begin{aligned} z(t)=\frac{{\hat{p}}^i{\hat{p}}^j}{2(1+{\hat{\Omega }}\cdot {\hat{p}})}[h_{ij}(t)-h_{ij}(t-t_l)], \end{aligned}$$

(1)

where $h_{ij}$ 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 ${\hat{p}}$ is a vector pointing to the pulsar position, ${\hat{\Omega }}$ is a vector along the direction in which the GW propagates, $t_l=(l/c)(1+{\hat{\Omega }}\cdot {\hat{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-t_l$).

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 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 accelerations, 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 “$\omega$” 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 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 stochastic GW background (Fig. 1b). The strain of the background is frequency dependent, described by the characteristic strain spectrum, $h_{\mathrm{c}}(f)$. This is calculated by integrating 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:

$$\begin{aligned} h_{\mathrm{c}} (f) = A\left( \frac{f}{f_{\mathrm {yr}}}\right) ^{\alpha } \ , \end{aligned}$$

(2)

or in terms of the energy density of GWs, $\Omega _{\mathrm{GW}} \propto f^{2(1+\alpha )}$. In this way, the background can be characterized with an amplitude A, and a single spectral index $\alpha$. The amplitude A is commonly defined at a frequency $f_{\mathrm{yr}}=1\, \hbox {year}^{-1} \sim 32\,\mathrm {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 $\alpha$ 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 successfully 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 background 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.

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 $\sim$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_\mathrm{c}\lesssim 1\times 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

3 The population and evolution of supermassive black hole binaries

It is now broadly accepted that SMBHs in a mass range around $10^6$–$10^{10}$ $\hbox {M}_\odot$ 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_\bullet -\sigma$, $M_\bullet -M_{\mathrm {bulge}}$, Ferrarese and Merritt 2000; Gebhardt et al. 2000), indicating a potential co-evolution between the two. In the standard paradigm 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 (Begelman 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.

Fig. 3

Binary SMBHs can form during a major merger. Pulsar timing arrays’ main targets are continuous-wave binaries within $\sim$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 $\sim$10 pc separations down to $\sim$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

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\ge q\ge 1$. The total mass is $M=m_1+m_2$ and the reduced mass is $\mu =m_1m_2/(m_1+m_2)$. Chirp mass is defined as ${\mathcal {M}}\equiv (m_1m_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 environment. 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 $\sim$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 interferometers (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).

As a binary evolves and accelerates in its orbit, it has a greater chance at decoupling from the environment. Somewhere below separations of $\sim$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):

$$\begin{aligned} \left\langle \frac{\mathrm{d}a}{\mathrm{d}t} \right\rangle= & {} -\dfrac{64}{5}\frac{m_1m_2(m_1+m_2)}{a^3}\dfrac{\left( 1+\dfrac{73}{24}e^2+\frac{37}{96}e^4\right) }{(1-e^2)^{7/2}}, \nonumber \\ \left\langle \frac{\mathrm{d}e}{\mathrm{d}t} \right\rangle= & {} -\dfrac{304}{15}\frac{m_1m_2(m_1+m_2)}{a^4}e\dfrac{\left( 1+\dfrac{121}{304}e^2\right) }{(1-e^2)^{5/2}}, \end{aligned}$$

(3)

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 eccentricity 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 ${\mathrm{d}}f/{\mathrm{d}}t\propto 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:

$$\begin{aligned} f_{\mathrm {peak}} \approx \frac{(1+e)^{1/2}}{(1-e)^{3/2}} \frac{f_{K,r}}{(1+z)}~, \end{aligned}$$

(4)

where

$$\begin{aligned} f_{K,r} = \frac{f(1+z)}{n}~. \end{aligned}$$

(5)

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

$$\begin{aligned} \frac{f_{K,r}(e)}{f_{K,r}(e_0)} = \left( \frac{\sigma (e_0)}{\sigma (e)}\right) ^{3/2}~, \end{aligned}$$

(6)

where $e_0$ is the eccentricity at some reference epoch, and,

$$\begin{aligned} \sigma (e) = \frac{e^{12/19}}{1-e^2}\left[ 1 + \frac{121}{304}e^2\right] ^{870/2299}~. \end{aligned}$$

(7)

Hence, a binary with (for example) $e_0=0.95$, when its orbital frequency is 1 nHz, will have $e\approx 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:

$$\begin{aligned} h(f_{K,r}) = \sqrt{\frac{32}{5}}\,\frac{(G{\mathcal {M}})^{5/3}}{c^4D}\,(2\pi f_{K,r})^{2/3}~. \end{aligned}$$

(8)

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 eccentricity. However, the strain amplitude is scaled by the degenerate parameters ${\mathcal {M}}^{5/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\lesssim 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 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

$$\begin{aligned} h \simeq h_+\simeq \frac{(1-\sqrt{8}/3)G\mu }{24 c^2 D}\mathrm{sin}^2(i)[17+\cos ^2(i)][1+{\mathcal {O}}(\mu ^2/M^2)] \simeq 0.02\frac{G\mu }{c^2D} \ , \end{aligned}$$

(9)

and the cross-hand polarization term $h_\times$ 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).

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

3.1.4 The stochastic binary gravitational-wave background

The superposition of GWs from the large population of SMBHBs predicted by hierarchical 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_{\mathrm{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 outline how many discrete continuous-wave sources can contribute to a GW background.

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

$$\begin{aligned} h_{\mathrm{c}}^2(f) =&\int _0^\infty \mathrm {d}z \int _0^\infty \mathrm {d}M_{1} \int _0^1 \mathrm {d}q \, {\frac{\mathrm {d}^4N}{\mathrm {d}z\,\mathrm {d}M_{1}\,\mathrm {d}q\,\mathrm {d}t_r}} \nonumber \\&\qquad \times \sum _{n=1}^\infty \left\{ {\frac{g[n,e(f_{K,r})]}{(n/2)^2}} {\frac{\mathrm {d}t_r}{\mathrm {d}\ln f_{K,r}}} {h^2(f_{K,r}) } \right\} , \end{aligned}$$

(10)

where the contributing factors are:

• (i) $\mathrm {d}^4N / \mathrm {d}z\mathrm {d}M_1\mathrm {d}q\mathrm {d}t_r$ is the comoving occupation function of binaries per redshift, primary mass, and mass–ratio interval, where $t_r$ 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 {$\cdots$} describes the distribution of GW strain over harmonics, 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 function g(n, e) is a distribution function whose form is given in Peters and Mathews (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) $\mathrm {d}t_r / \mathrm {d}\ln f_{K,r}$ describes the time each binary spends emitting in a particular logarithmic 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_{\mathrm{c}}(f)$ is easily deduced. In this case, $g(n=2,e=0)=1$ and $g(n \ne 2,e=0)=0$ such that $f=2f_{K,r}/(1+z)$, and the residence time $\mathrm {d}t_r / \mathrm {d}\ln f_{K,r} \propto f^{-8/3}$ is given by the quadrupole radiation formula (Peters 1964). Collecting terms in frequency gives $h_{\mathrm{c}}(f) \propto f^{-2/3}$, as per Eq. 2 (this single power law spectrum is shown as the black, dashed line in Fig. 4).²

The $h_{\mathrm{c}} \propto f^{-2/3}$ power law GW background spectrum assumes a continuous distribution 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 \gtrsim 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).³ 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 harmonics 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 frequencies 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_{\mathrm{c}}$ at low frequencies ($\lesssim 10^{-8}$ Hz) and increasing it at higher frequencies ($\gtrsim 10^{-8}$ Hz) (Enoki et al. 2007; Huerta et al. 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 creates 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 cluster 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 background 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 $\sim 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 $\sim 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_\bullet -M_{\mathrm {bulge}}$ 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).

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

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 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\propto {\mathcal {M}}^{5/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, the claims of a $\sim 2\times 10^{10}M_\odot$ 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 ($>10^9M_\odot$) end of the relation, demonstrating that $M_\bullet -M_{\mathrm{bulge}}$ 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 galactic bulge mass: $M_{\mathrm {BH}} = \beta M_{\mathrm {bulge}}$. 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 $\beta$ (the most widely used one is $\beta \approx 0.003$); that issue is analyzed in detail in the Section IIC of Rasskazov and Merritt (2017b). Given the observed mass and mass ratio distribution of merging galaxy pairs, $\beta$ 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):

$$\begin{aligned} h_{\mathrm{c}}(f)= & {} A\frac{(f/f_{\mathrm{yr}})^{-2/3}}{1+(f_{\mathrm{b}}/f)^{53/30}}, \end{aligned}$$

(11)

$$\begin{aligned} A= & {} 2.77\times 10^{-16}\, \left( \frac{\beta }{10^{-3}}\right) ^{0.83},\end{aligned}$$

(12)

$$\begin{aligned} f_{\mathrm{b}}= & {} 1.35\times 10^{-9}\,\mathrm{Hz}\, \left( \frac{\beta }{10^{-3}}\right) ^{-0.68}. \end{aligned}$$

(13)

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 Milosavljević 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 (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):

$$\begin{aligned} T_{\mathrm {DF}} \approx \frac{19}{\ln \Lambda }\ \left( \frac{R_e}{5\,\mathrm {kpc}}\right) ^2\left( \frac{\sigma }{200\,\mathrm {km\ s^{-1}}}\right) \left( \frac{10^8 M_\odot }{M_\bullet }\right) \ \mathrm {Gyr}, \end{aligned}$$

(14)

where $R_e$ and $\sigma$ are the galaxy’s effective radius and velocity dispersion, $M_\bullet$ is the SMBH’s mass and $\ln \Lambda$ is the Coulomb logarithm.⁴ However, the braking of the individual SMBHs in reality will be much shorter. In a genuine merger system, the $M_\bullet$ 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\,\mathrm {Gyr}$ (Dosopoulou and Antonini 2017; Kelley et al. 2017a; Yu 2002).⁵

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 ($\lesssim 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 $\sim$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 individual 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):

$$\begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&= -\frac{G\rho }{\sigma }H a^2 \ , \end{aligned}$$

(15)

$$\begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t}&= \frac{G\rho }{\sigma } H K a \ , \end{aligned}$$

(16)

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” (Milosavljević 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 numerical 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 eccentricity), 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 population⁶, we can re-arrange 15 to deduce the orbital frequency evolution, and hence the evolution of the emitted GW frequency, such that $\mathrm{d}f/\mathrm{d}t \propto f^{1/3}$. Inserting this into 10 and collecting terms in frequency gives $h_\mathrm{c}(f)\propto f$, which is markedly different from the fiducial $\propto 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 analytic 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 pressure (the so-called $\alpha$-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):

$$\begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t} = -\frac{2{\dot{m}}_1}{\mu }(aa_0)^{1/2}, \end{aligned}$$

(17)

where ${\dot{m}}_1$ is the mass accretion rate onto the primary BH, and $a_0$ 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

$$\begin{aligned} a_0 =&3\times 10^3 \left( \dfrac{\alpha }{10^{-2}}\right) ^{4/7} \left( \dfrac{m_2}{10^6 M_\odot }\right) ^{5/7} \nonumber \\&\times \left( \dfrac{m_1}{10^8 M_\odot }\right) ^{-11/7} \left( 100d\dfrac{{\dot{m}}_1}{{\dot{M}}_E}\right) ^{-3/7} r_g, \end{aligned}$$

(18)

where $\alpha$ is a disk viscosity parameter, ${\dot{M}}_E = 4\pi Gm_1 / c\kappa _{\mathrm {T}}$ is the Eddington accretion rate of the primary BH ($\kappa _{\mathrm {T}}$ is the Thompson opacity coefficient), and $r_g = 2Gm_1 / c^2$ 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 Teukolsky 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 prescriptions 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 eccentricity will experience a reduction in eccentricity, leading to the discovery of a limiting eccentricity for disk–binary interactions that falls in the interval $e_{\mathrm {crit}} \in [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).

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 surrounding disks, Kocsis and Sesana (2011) studied scaling relations for the strain spectrum from different disk–binary scenarios, varying from $h_\mathrm{c}(f)\propto f^{-1/6}$ for secondary-dominated type-II migration in a radiation-dominated $\alpha$-disk, to $h_\mathrm{c}(f)\propto 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 $\propto 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 standard $\alpha =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).

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)

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 (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 $10^9~M_\odot$ SMBHB system with an Earth term frequency of 100 nHz should produce 32.1 gravitational wave cycles, of which 31.7 are from the leading Newtonian order (or $\hbox {p}^{0}$N), 0.9 wave cycles are from $\hbox {p}^1$N order, and −0.7 are from $\hbox {p}^{1.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 $\sim 3000$ years ago (roughly equivalent to the typical light travel time between the Earth and the pulsar). Over this time, one expects 4305.1 wave cyles, of which 4267.8 are Newtonian, 77.3 come from $\hbox {p}^1$N order, $-45.8$ are from $\hbox {p}^{1.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\pi f L <1$ to not lose a single wave cycle, where f is the GW frequency and L is the distance to the pulsar. Therefore, it is in principal necessary to measure the pulsar distances to, e.g., $\sim 1.5$ pc for a GW signal at 1 nHz.

Many pulsar distances are poorly constrained, since most are estimated via the dispersion 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 measurement of $156.3 \pm 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\pm 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 Array⁷ 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

Topological defects are a generic prediction of grand unified theories. As the universe 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 attention of the theoretical cosmology community has focused on strings as the only viable candidate that could lead to potential observational signatures. The most simple symmetry breaking that leads to the formation of cosmic strings occurs in the Abelian Higgs model, where the symmetry group U(1) breaks

$$\begin{aligned} U(1) \rightarrow 1 \end{aligned}$$

at some temperature, or energy scale, $T_s$. They are characterized by their mass per unit length $\mu$, which in natural units is determined by the temperature at which the phase transition takes place, $\mu \sim T_s^2$. The tension of a string is normally given in terms of the dimensionless parameter $G\mu /c^2.$ 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 condensates, 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.

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

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 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 superstrings 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 theory, 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 \le 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 areas of the cosmic (super)string $G\mu /c^2$-p parameter space excluded by the present and potential future experiments, including LIGO and the three leading PTA experiments (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 scaling 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.

Fig. 9

Plot of regions of the cosmic (super)string $G\mu /c^2$-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)

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\mu /c^2 \lesssim 10^{-7}$, relegating them to at most a secondary role in structure formation. However, cosmic (super)strings are still potential candidates for the generation 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. Correspondingly, 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\mu /c^2 < 5.3(2) \times 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 $\Omega$ versus frequency for dimensionless string tensions $G\mu /c^2$ in the range $10^{-23}$-$10^{-9}$ for $p=1$. Overlaid are the current experimental constraints from PTAs⁸ 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.

Fig. 10

Plot of the gravitational wave spectrum in terms of the dimensionless parameter $\Omega$, 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\mu /c^2$ 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

General relativity has been an exceptionally successful physical theory and continues 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 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 observational 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_{\mu \nu } = \eta _{\mu \nu } + h_{\mu \nu }$, where $\eta _{\mu \nu }$ is the flat-space Minkowski metric and $|h_{\mu \nu }| \ll 1$. Since the metric is a symmetric $4 \times 4$ matrix, $h_{\mu \nu }$ 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 polarization 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_{\mu 0}=0$, which is called the synchronous 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):

$$\begin{aligned} e_{ij}^+= & {} \left( \begin{array}{ccc}1 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad -1 &{} \quad 0 \\ 0 &{} \quad 0 &{} 0\end{array}\right) , \ e_{ij}^{\times } = \left( \begin{array}{ccc}0 &{} 1 &{} \quad 0 \\ 1 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0\end{array}\right) ,\end{aligned}$$

(19)

$$\begin{aligned} e_{ij}^b= & {} \left( \begin{array}{ccc}1 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 1 &{} \quad 0 \\ 0 &{} \quad 0 &{} 0\end{array}\right) , \ e_{ij}^\ell = \left( \begin{array}{ccc}0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 1\end{array}\right) , \end{aligned}$$

(20)

$$\begin{aligned} e^x_{ij}= & {} \left( \begin{array}{ccc}0 &{} \quad 0 &{} \quad 1 \\ 0 &{} \quad 0 &{} \quad 0 \\ 1 &{} \quad 0 &{} \quad 0\end{array}\right) , \ e^y_{ij} = \left( \begin{array}{ccc}0 &{} \quad 0 &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad 1 \\ 0 &{} \quad 1 &{} \quad 0\end{array}\right) , \end{aligned}$$

(21)

where we have normalized these polarization tensors to satisfy $e^P_{ij} e_{P'}^{ij} = \delta _{P,P'}$. Three of these polarization tensors ($+$/$\times$, and b) are transverse; the other scalar mode and the two vector components are longitudinal, as shown in Fig. 11.

Fig. 11

The six possible GW polarizations in synchronous gauge where $h_{\mu 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)

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 $\theta$. 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.

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

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

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, $\theta$. The shape of this angular correlation function of pulse residuals, $C(\theta )$, 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 $\theta$ (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, 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 $\times$ 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 $\sim$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 binary system. Given that the scalar longitudinal mode produces the largest autocorrelation (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 parameterization in geometric units (Blas et al. 2016):

$$\begin{aligned} \omega ^2 = m_{\mathrm{g}}^2 + c_{g}^2 k^2, \end{aligned}$$

(22)

where $m_{\mathrm{g}}$ is the mass of the graviton and $c_{g}$ 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 suppression from which $m_{\mathrm{g}}$ can be constrained using a variety of observations. The most constraining limit comes from weak lensing observations which yield $m_{\mathrm{g}}< 6.9 \times 10^{-32}\ \mathrm{eV}$ and can be translated into a frequency $f_{g} = m_{\mathrm{g}}c^2/h > 1.7\times 10^{-17}\ \mathrm{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_{\mathrm{g}}< 1.27 \times 10^{-30}\ \mathrm{eV}$ (Gupta and Desai 2018). A massive graviton changes the functional form of the gravitational potential from a shape of $\sim 1/r$ to $\sim e^{-\mathrm{km}_{\mathrm{g}}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_{\mathrm{g}}< 4.4 \times 10^{-22}\ \mathrm{eV}$ or $f_{g} > 10^{-7}\ \mathrm{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 $c_{g} = 1$, constraints from GW170104 provide the best dynamical constraint to the graviton mass and yield $m_{\mathrm{g}}< 7.7 \times 10^{-23}\ \mathrm{eV}$ or $f_{g}>1.8 \times 10^{-8}\ \mathrm{Hz}$ (Abbott et al. 2017d).

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_\mathrm{g} = 3 \times 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 $\sim$10 pulsars with RMS accuracy of 100 ns, PTA measurements are not yet competitive with other experiments for probing graviton mass

Pulsars can also be sensitive to the graviton mass and, given that the best dynamical 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 graviton 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_\mathrm{g} = 3 \times 10^{-22}$ eV with 90% confidence level. A 10-year observation with the same RMS timing accuracy and cadence but with 300 pulsars would probe graviton masses down to $m_\mathrm{g} = 3 \times 10^{-23}$ eV.

In addition to using the correlated residuals between pulsars in response to a stochastic 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_{\mathrm{g}}\ll 10^{-22}\ \mathrm{eV}$ so that $f_{g} \gg \mathrm{year}^{-1}$, in the PTA band the gravitational dispersion relation becomes $\omega ^2 \simeq c_g^2 k^2$. If $c_g<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-c_g<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 analysis of GW150914 allowed the first direct constraint on $c_g$ giving $c_g < 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\times {10}^{-15}<c_g-1<7\times {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 (Jacobson and Mattingly 2001), massive gravity (de Rham et al. 2011), and f(R) gravity (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 gravity (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, $\{{\mathcal {C}}_1,{\mathcal {C}}_2,{\mathcal {C}}_3,{\mathcal {C}}_4\}$ (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_{\mu 0}=0$ and taking the ${\mathcal {C}}_i \rightarrow 0$ limit, the theory predicts: plus/cross tensor polarizations traveling at the speed of light; two longitudinal vector modes traveling at a speed possibly different from the speed of light; and a linear combination of the two scalar modes

$$\begin{aligned} h_{ij}^s=v_0\left[ \frac{3}{c_s^2}e_{ij}^\ell -({\mathcal {C}}_1+{\mathcal {C}}_4)e_{ij}^b\right] , \end{aligned}$$

(23)

where $v_0$ is the linear perturbation to the time component of the Aether field, $c_s$ is the speed of the scalar GWs, and $e_{ij}$ 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 discontinuity, 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 freedom. This is known as the Boulware–Deser ghost (Boulware and Deser 1972), proof 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 $c_s$ with a polarization

$$\begin{aligned} h^{s}_{ij} = \frac{\pi }{6M_{\mathrm{pl}}}\left[ e^b_{ij} + \left( 1-\frac{1}{c_s^2}\right) e^\ell _{ij}\right] , \end{aligned}$$

(24)

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+\mu /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

$$\begin{aligned} m^2 = \frac{1}{3} \left[ \frac{1}{f^{\prime \prime }(R_0)} - \frac{R_0}{f^{\prime }(R_0)} \right] , \end{aligned}$$

(25)

where $R_0$ 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)

$$\begin{aligned} h^s_{ij}= -\frac{1}{2}f^{\prime }(R_0)\left( e^b_{ij}+\frac{m^2}{\omega ^2} e^\ell _{ij}\right) , \end{aligned}$$

(26)

and follow the dispersion relation $m^2 = \omega ^2 - k^2$.

6 Relic gravitational waves and early-universe cosmology

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.

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 $\Omega _{\mathrm{gw}}(f)<2.3\times 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 interferometers 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 $\Omega _{\mathrm{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_\odot$, 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).

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)

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 background. 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 \sim 10^{10}$. 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 surpass 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 > 10^{17}$) and the GWs produced during that epoch. Such a phase transition could produce GWs with $f \sim 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 $\Lambda$ (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 constraints on the properties of dark energy. Bernabeu et al. (2011) and Espriu and Puigdomènech (2013) describe an approach in which measurements of pulsar timing residuals due to GWs emitted in the nearby universe ($z \ll 1$) could place a complementary constraint on $\Lambda$. 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.

7 Additional science from “Hidden Planets” to dark matter

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 $\Lambda$-Cold-Dark-Matter ($\Lambda$CDM) model for the evolution of the universe predicts that the mass function of dark matter halos may extend to very small masses ($\ll M_\odot$), 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 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_\odot$. 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_\odot$.

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 $\Lambda$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,

$$\begin{aligned} \lambda _{dB}\approx 600~\mathrm{pc}\left( \frac{10^{-23}~\mathrm{eV}}{m}\right) \left( \frac{10^{-3}}{v}\right) , \end{aligned}$$

(27)

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

$$\begin{aligned} f\approx 5\cdot 10^{-9}~\mathrm{Hz}\left( \frac{m}{10^{-23}~\mathrm{eV}}\right) \end{aligned}$$

(28)

and an amplitude

$$\begin{aligned} \Psi (\mathbf{x})=\pi \frac{G\rho (\mathbf{x})}{m^2}, \end{aligned}$$

(29)

where $\rho (\mathbf{x})$ is the density of dark matter particles at position $\mathbf{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 over 12 years to compute rigorous Bayesian and frequentist constraints on the local density of ultralight scalar-field dark matter. They found that $\rho _{\mathrm {SF}}\le 6$ Gev $\hbox {cm}^{-3}$ at $95\%$ confidence for ultralight bosons with $m\le 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 frequency of the bump that can be produced in the GW background spectrum depends on the primordial BH mass, and the 10–100 $M_\odot$ 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_\odot$). 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):

$$\begin{aligned} \Delta t_{\mathrm {SSB}} = \frac{{\mathbf {R}}_{\mathrm {SSB}}\cdot \hat{{\mathbf {n}}}}{c}, \end{aligned}$$

(30)

where ${\mathbf {R}}_{\mathrm {SSB}}$ is the vector connecting the Earth to the solar system barycenter, $\hat{{\mathbf {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 ${\mathbf {R}}_{\mathrm {SSB}}$ 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 statistics 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 perturbations 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 $\lesssim 10^{-4} M_\odot$, 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).

8 Gravitational-wave synergies: multi-band GW science

The LIGO/Virgo GW detections (Abbott et al. 2016a, b, 2017a) in the high-frequency GW regime ($\gtrsim$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 ($\sim$ mHz) expected to be launched in the 2030s (Amaro-Seoane et al. 2012).

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)

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 ${\mathcal {M}} \sim 10^9 \, M_\odot$ (Mingarelli et al. 2017). In contrast, LISA is sensitive to the inspiral, merger, and ringdown of SMBHBs with masses from $10^4 - 10^7~M_\odot$ 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\times 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 the probability of a sequential detection being extremely low ($4.7\times 10^{-4}$– $3.3\times 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 gravitational 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 circumstances, 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, 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, several 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 conceptually 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) overlapping 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 $\mu$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 detection (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, astrometric 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 \gtrsim 3\times 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).

9 Electromagnetic synergies: multi-messenger science

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 $\sim 7\sigma$) will be up to a few times 1000 $\hbox {deg}^2$ (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 ($\sim 10^{10}\,M_\odot$). 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 determine 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 photometry, or the morphology of their host. We summarize the methods (and efforts) to electromagnetically identify SMBHBs in Sect. 9.2.

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 ($\theta ,\phi$) 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 $\times$ 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)

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 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 science, 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). Additionally, 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 information 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.

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 identification 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 systematic 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 SMBHBs. 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 ($\sim$150) have emerged from systematic 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 compared 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 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 ${\mathcal {M}}_c$ (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 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_\bullet$ 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 $M_c$. As we write this, PTAs are sensitive to SMBHBs of $M_c>10^9\,M_\odot$ out to approximately 200 Mpc, or a redshift of $z\simeq 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 $M_c$ 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\gtrsim 0.1$, implying an error when inferring M of $\sim$0.5 dex. Even in the more conservative limit of $q\gtrsim 0.01$, the error in the total mass inferred from PTA measurement of ${\mathcal {M}}$ would be only $\sim 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) demonstrated 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 predictions 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:

$$\begin{aligned} t_{\mathrm{X-ray}} \simeq 7(1+z)\bigg (\frac{M}{10^6\,M_\odot }\bigg )^{1.32}\,\mathrm{yr} \end{aligned}$$

(31)

(Milosavljević and Phinney 2005). For the high-mass systems that PTAs are expected to discover ($M\gtrsim 10^8\,M_\odot$), 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 evolution 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. However, given the broad coverage of current and upcoming synoptic sky surveys, like the Catalina Sky Survey (CRTS, Drake et al. 2009), the Panoramic Survey Telescope 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 $\sim$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; Milosavljević 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.

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 coordinated 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 astrophysical application. Most famously, a ${\mathcal {M}}\simeq 10^{10}\,M_\odot$ SMBHB was hypothesized to exist in galaxy 3C66B due to an elliptical motion detected through long-baseline interferometry; 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 ${\mathcal {M}}>0.5\times 10^{8}\,M_\odot$ 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 $\sim$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 $\sim 10^4$ 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 covariant 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. However, the predictive power of these discoveries currently suffer from small number of statistics. Future large-scale systematic surveys, such as LSST, or the Square Kilometer 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 propagate at the same speed. This agreement can be tested if a resolvable continuous wave 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 ($\lambda _{g}=h/(m_\mathrm{g} \mathrm{c})$) of such a system (Will 1998):

$$\begin{aligned} \lambda _{g} > 3 \times 10^{12} \mathrm{km} \left( \frac{D}{200 \mathrm{Mpc}} \frac{100 \mathrm{Hz}}{f} \right) ^{\frac{1}{2}}\left( \frac{1}{f \Delta t} \right) ^{\frac{1}{2}} \end{aligned}$$

(32)

and choosing nominal values for the parameters of an SMBHB, of $f=1\times 10^{-9}$ Hz, $D=1$ Gpc and an observing time of $\Delta t=5$ years, one obtains a lower limit of $\lambda _{g} \sim 10^{19}$ km. This limit is three orders of magnitude less precise than the current Particle Data Group limit of $10^{22}$ 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_\mathrm{g}$, if one assumes a Yukawa-type potential ($V_g\sim e^{-\mathrm{km}_\mathrm{g}}/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.

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

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 is no longer applicable if the mass of the lens is too small compared to the wavelength of the radiation; less than $10^5 M_\odot (f/\mathrm{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 \sim 10$ nHz, there is a reasonable probability of obtaining a few lenses, for which the time difference between the GW and electromagnetic 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 American 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 Timing 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 precision, 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 timing more than 70 pulsars. In addition to new pulsar discoveries, improvements to radio observing systems will continue to improve mitigation of intrinsic pulse-to-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 improvements 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, simply 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 maintain unique access to exploring the inspiral of high-mass SMBHB systems, and will uniquely probe the physical properties of strings.