The Local Nanohertz Gravitational-Wave Landscape From Supermassive Black Hole Binaries

# The Local Nanohertz Gravitational-Wave Landscape From Supermassive Black Hole Binaries

Chiara M. F. Mingarelli11affiliationmark: 22affiliationmark: 33affiliationmark: $\dagger$ $\dagger$affiliationmark: , T. Joseph W. Lazio33affiliationmark: , Alberto Sesana44affiliationmark: , Jenny E. Greene55affiliationmark: , Justin A. Ellis33affiliationmark: , Chung-Pei Ma66affiliationmark: , Steve Croft66affiliationmark: , Sarah Burke-Spolaor77affiliationmark: 88affiliationmark: , Stephen R. Taylor33affiliationmark: Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, D-53121 Bonn, Germany TAPIR, MC 350-17, California Institute of Technology, Pasadena, California 91125, USA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA School of Physics and Astronomy and Institute of Gravitational Wave Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA Department of Astronomy, University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA Department of Physics and Astronomy, West Virginia University, Morgantown, West Virginia 26506, USA Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505
###### Abstract

Supermassive black hole binaries (SMBHBs) in the 10 million to 10 billion range form in galaxy mergers, and live in galactic nuclei with large and poorly constrained concentrations of gas and stars. There are currently no observations of merging SMBHBs— it is in fact possible that they stall at their final parsec of separation and never merge. While LIGO has detected high frequency GWs, SMBHBs emit GWs in the nanohertz to millihertz band. This is inaccessible to ground-based interferometers, but possible with Pulsar Timing Arrays (PTAs). Using data from local galaxies in the 2 Micron All-Sky Survey, together with galaxy merger rates from Illustris, we find that there are on average sources emitting GWs in the PTA band, and binaries which will never merge. Local unresolved SMBHBs can contribute to GW background anisotropy at a level of , and if the GW background can be successfully isolated, GWs from at least one local SMBHB can be detected in 10 years.

$\dagger$$\dagger$affiliationtext: mingarelli@mpifr-bonn.mpg.de

Supermassive black holes (SMBHs) are widely held to exist at the heart of massive galaxies (Ferrarese & Ford, 2005). Galaxy mergers should form SMBH binary (SMBHB) systems, which eventually emit gravitational waves (GWs) and merge (Begelman et al., 1980). Galaxy mergers are a fundamental part of hierarchical assembly scenarios, forming the backbone of current structure formation models. Thus, the detection of GWs from merging SMBHs would have fundamental and far-reaching importance in cosmology, galaxy evolution, and fundamental physics, not accessible by any other means.

Pulsar Timing Arrays (PTAs) can detect nanohertz GWs by monitoring radio pulses between millisecond pulsars, which are highly stable clocks. GWs change the proper distance between the pulsars and the Earth, thus inducing a delay or advance of the pulse arrival times. The difference between the expected and actual arrival times of the pulses – the timing residuals – carries information about the GWs that is extracted by cross-correlating the pulsar residuals (Estabrook & Wahlquist, 1975; Sazhin, 1978; Detweiler, 1979; Hellings & Downs, 1983). Current PTAs include European PTA (Lentati et al., 2015) (EPTA), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Arzoumanian et al., 2016) the Parkes PTA (Shannon et al., 2015), and the International PTA (IPTA) (Verbiest et al., 2016), the latter being the union of the former three.

Here we introduce a bottom-up approach to constructing both realistic GW skies and future IPTA projections: we use IPTA pulsars with their real noise properties, and galaxies from the 2 Micron All Sky Survey (2MASS) (Skrutskie et al., 2006), together with galaxy merger rates from Illustris (Rodriguez-Gomez et al., 2015; Genel et al., 2014), to form multiple probabilistic realizations of the local GW Universe. In each realization, we search for SMBHB systems which emit continuous GWs (CGWs) in the PTA band, and also compute their contribution to the nanohertz GW background (GWB) and its anisotropy (Mingarelli et al., 2013; Schutz & Ma, 2016). We report on the physical properties of the most frequently selected SMBHBs and their host galaxies, and estimate their time to detection.

## Galaxy Selection

SMBHB merger timescales can be of the order  yrs after the galaxy merger, and therefore morphological merger signatures can be difficult to identify. Our focus here is on massive early-type galaxies, as these are likely to have formed from major mergers and would therefore host SMBHBs with approximate mass ratios of .

To approximate a mass selection, we select in the band using the 2MASS (Skrutskie et al., 2006) Extended Source Catalog (XSC) (Jarrett et al., 2000), following the procedure outlined in detail in Ref (Ma et al., 2014), but to a distance of 225 Mpc and over the full sky. We do not excise the galactic plane, but there are fewer reliable sources there. To approximate a mass selection, we select in the band using the 2MASS (Skrutskie et al., 2006) Extended Source Catalog (XSC) (Jarrett et al., 2000). We convert from the 2MASS K-band luminosity to the stellar mass via , appropriate for early-type galaxies (Cappellari, 2013), and apply a K-band cut to select galaxies with stellar mass of , since these are likely to host SMBHBs (Rosado & Sesana, 2014). At distances  Mpc, 2MASS itself begins to become incomplete at . Although our sample is not formally volume-complete within our chosen mass and distance cuts, it includes the majority of massive, nearby galaxies, which are expected to dominate the signal for PTAs.

This process creates a galaxy catalog with 5,110 early-type galaxies. In Ref (Ma et al., 2014; Schutz & Ma, 2016), we find that 33 of these galaxies contain dynamically measured SMBHs, Figure 1(a). Moreover, we manually add a further 9 galaxies from 2MASS: NGC 4889, NGC 4486a, NGC 1277, NGC 1332, NGC 3115, NGC 1550, NGC 1600, NGC 7436 and A1836 BCG, which did not make the luminosity cut, but also host dynamically measured SMBHs.

The SMBHB total mass, , is estimated by taking the stellar mass to be the bulge mass, and applying the empirical scaling relation from Ref (McConnell & Ma, 2013), with scatter , described in Methods. Dynamically measured SMBH masses fixed.

## Probability of a Galaxy Hosting a SMBHB in the PTA Band

We attack this problem by estimating two quantities: the probability that a SMBHB is in the PTA band, and the probability that a given galaxy has a SMBHB. The former is estimated via quantities which are derived from 2MASS, e.g. each galaxy’s stellar mass, , and through , the SMBHB mass M. The latter is obtained from the cumulative galaxy-galaxy merger rates from the Illustris cosmological simulation project (Rodriguez-Gomez et al., 2015; Genel et al., 2014).

A realization of a the local nanohertz GW sky is created as follows. The first quantity is the probability that a SMBHB is in the PTA band, , where is the time to coalescence of the binary and is the chirp mass.

In each realization, all galaxies in the sample are assigned a SMBHB with total mass , and mass ratio drawn from a log-normal distribution in . The time to coalescence is computed from  nHz, which is the beginning of the PTA band. is the effective lifetime of the binary: this is the sum of the dynamical friction (Binney & Tremaine, 2008) () and stellar hardening (Sesana & Khan, 2015) () timescales (details in Methods).

The second quantity is the probability that a galaxy hosts SMBHB. This is computed via the Illustris cumulative galaxy-galaxy merger rate, taken at the beginning of the binary evolution. This redshift is calculated at with Planck cosmological parameters (Planck Collaboration et al., 2016). The merger rate is then multiplied by time elapsed since then, . Thus, probability of galaxy hosting a SMBHB in the PTA band,

 pi=tc,iTz∫10.25dμ∗dNdt(M∗,z,μ∗)Tz. (1)

Calculating is essentially rewinding the SMBHB evolution: starting in the GW emission phase, we compute how long the binary spent in a stellar hardening phase, and then in a dynamical friction phase, for binary separations out to the effective radius of the galaxy. We can only rewind 12.5 Gyrs , as this is the maximum from Ref (Rodriguez-Gomez et al., 2015). If  Gyr (equivalently ), binary has likely stalled, and we set .

The total number of SMBHBs emitting GWs in the PTA band, for each Monte Carlo realization, is the sum of all these probabilities: . This number varies from realization to realization. We draw galaxies from the galaxy catalog according to the probability distribution, .

For each of these selected galaxies, we compute the inclination and polarization-averaged strain,

 (2)

assuming circular binary orbits, and is the GW frequency. Each of the probabilistically selected galaxies hosts a SMBHB in the PTA band with  nHz, where the evolution is assumed to be dominated by GW emission. We assign each SMBHB a GW frequency (Peters, 1964), :

 f=π−1M−5/8c[2565(tc−t)]−3/8, (3)

by drawing from a uniform distribution in , which is the time to coalescence of a SMBHB with from 1 nHz and 100 nHz, respectively111Formally we set and sample in . Galaxy distances are estimated via techniques outlined in the Methods section. We use the approximation only in Equation (2), since GW sources are all  Mpc.

## Projections for the IPTA

Continuous nanohertz GW upper limits come in two flavors: 1) as a function of GW source sky location, GW frequency, and pulsar sky location (Figure 1(b)), and 2) a sky-averaged strain upper limit as a function of GW frequency, which averages over all pulsar sky locations, Figure 1(d). Since upper limit curves, and not detection curves, were reported by Ref (Babak et al., 2016), we use the upper limit as a proxy for detection to underline the importance of the pulsar and SMBHB sky location. Indeed, the PTA response to CGWs is maximal when the GW source lies very close to the pulsar (Detweiler, 1979; Mingarelli & Sidery, 2014).

For IPTA projections, we are more rigorous, constructing an IPTA-like array to estimate the time to detection of CGW sources in 2MASS. We start with the 47 IPTA pulsars (Verbiest et al., 2016), and from 2016 onward we add four pulsars per year from sky locations accessible with IPTA telescopes, using the median white noise uncertainty of 300 ns.

Detection probability curves are computed using the statistic (Ellis et al., 2012) and False Alarm Probabilities (Arzoumanian et al., 2014) (FAPs) of , , and . For Gaussian distributions, these 1-FAPs are , and respectively.

## Gravitational Wave Backgrounds

The search for the GWB is ongoing (The NANOGrav Collaboration et al., 2015; Desvignes et al., 2016; Verbiest et al., 2016; Arzoumanian et al., 2016; Lentati et al., 2015; Shannon et al., 2015) with detection expected in 5 to 7 years (Taylor et al., 2016). While most PTAs publish limits on isotropic GWBs, searching for and characterizing GWB anisotropy is gaining momentum (Taylor et al., 2015). This is done by decomposing the sky on a basis of spherical harmonics (Mingarelli et al., 2013; Taylor & Gair, 2013; Gair et al., 2014), , and with pixel-bases (Cornish & van Haasteren, 2014).

In an effort to understand the contribution of these CGW sources to the GWB, we transform a GW sky realization to a GWB. This is done by casting each individual binary’s strain contribution to the closest pixel in a HEALPix sky map, and computing the characteristic strain , where is the strain of source , is its GW frequency, and is the inverse observation time (Babak et al., 2016).

The total power on the sky is , and is decomposed as , where is the direction of GW propagation. Anisotropy in described terms of , normalized to the isotropic component, . We calculate the angular power spectrum of a GW sky and take the monopole value as the contribution to the overall isotropic GWB.

## Results

Galaxies hosting SMBHBs We compute the probability of each galaxy in the catalog containing a SMBHB emitting GWs in the PTA band,  nHz. We carry out multiple realizations of the galaxy catalog, sampling over BH mass, mass ratio, and time to coalescence. We find that on average, there are galaxies hosting SMBHBs in the PTA band and stalled SMBHBs, despite the inclusion of a stellar hardening phase (Quinlan, 1996) to overcome the final parsec problem (Supplementary Figure 4).

Over multiple realizations, fewer than of GW skies hosted a currently detectable SMBHB. We also note that for 25 year datasets, binaries with are the likeliest to be detected, with more massive binaries being disfavored, Supplementary Figure 5(d).

Six pulsars dominate the EPTA’s sensitivity to CGWs, leading to an upper-limit map of the sky that is approximately dipolar, with a factor of 4 more in strain sensitivity behind the best pulsars, Figure 1. While we use the EPTA as an example, results are similar for NANOGrav, PPTA and the IPTA.

Anisotropy and contributions to the GW background In Figure 2(a), we show an example realization of the local nanohertz GW sky with a loud SMBHB, chosen at random. In this realization, NGC 4472 hosted a PTA-detectable SMBHB. Of course, NGC 4472 was only one of 87 galaxies in this realization, but it was the only one which contained a binary that was sufficiently loud to be detected. Assuming an isotropic GWB with an amplitude of a few times and a 25 year dataset, we find that such a single CGW source contributes less than 1% to the overall strain budget. This is to be expected, since a CGW is the ultimate anisotropy, and therefore contributes very little to an isotropic GWB. The strong CGW source does, however, dominate the angular power spectrum of the sky, Figure 2(d), where we measure GWB anisotropy. When the strong CGW source is removed, Figure 2(c), we find anisotropy from the superposition of GWs from undetected CGW sources at the level of , dominating the angular power spectrum of the GWB up to , Figure 2(d).

Time to Detection Evidence for nanohertz GWs will increase slowly and continuously. We estimate the time to detection of nearby GW sources with the IPTA, requiring a 95% detection probability under different false alarm probabilities (FAP): , , (, and for Gaussian distributions). This is done for 15, 20, and 25-year datasets, noting that current datasets are 10-15 years. Hence, 25-year predictions are for 10 years from now. Results for detections with FAP are in Figure 3, and summarized in Table 1.

We find that strong red noise in pulsar residuals greatly diminishes the chance of detecting CGW sources in the next 10 years. However, if the pulsar noise is white or if CGW sources can be extracted from an unresolvable GWB, then there is a 50% chance of detecting a local CGW source with a FAP in the next decade. Even more encouragingly, when one uses the all-sky GW strain sensitivity map for detections, a detection with FAP is possible in 10 years, Table 1.

## Discussion

Using the sky location and noise properties of IPTA pulsars, massive galaxies in 2MASS and galaxy merger rates from Illustris, we estimate when and where PTAs are likely to detect CGW sources in the local Universe. Over multiple realizations of the local GW sky, we find that of GW sources would have been detected with current PTA data (Desvignes et al., 2016), Figure 1(c), supporting the conclusions that current non-detection is unsurprising, Ref (Babak et al., 2016) .

In making IPTA predictions, we did not include new telescopes which will come on line in the next 10 years, such as MeerKAT (Booth et al., 2009) and possibly SKA Phase 1 (Janssen et al., 2015). Both telescopes will greatly increase PTA sensitivity in the Southern hemisphere.

Importantly, future IPTA detections depend on how successfully the GWB can be subtracted. The red noise in the pulsars is meant to emulate an unresolved GWB with , consistent with Ref (Sesana et al., 2016). While some pulsars exhibit intrinsic red noise, the highest precision timers – which contribute most to CGW sensitivity – are broadly consistent with being white noise dominated.

An overview of the detected CGW parameters is given in Figure 5. Interestingly, massive galaxies such as M87 have a lower probability of being selected, since this depends on . Therefore, binaries in e.g. M104 are more likely to host SMBHBs in the PTA band, Figure 5.

We performed a brief literature search on the likeliest galaxies to host SMBHBs (Supplementary Figures 5(a) and 5(b)) to assess whether they showed signs of merger or a current candidate binary. We find that NGC 3115 is the only object currently under investigation as a candidate binary or recoiling black hole (Wrobel & Nyland, 2012). While many of the other candidates in the top SMBH red-noise candidate list show signs of recent or ongoing merging activity, many galaxies in this mass range are involved in merging, and a more complete comparison between the merger properties of this sample with the general population is beyond the scope of this work.

## Methods

Galaxy Selection from 2MASS

We select our initial sample from the Two Mass Redshift Survey (2MRS (Huchra et al., 2012)). We first make a very broad selection of objects with and radial velocity distance  Mpc. We then cross-correlate with the Crook group catalog (Crook et al., 2007) to correct all radial velocities for galaxies in our sample to the radial velocity of the group. In practice, this has a negligible impact on the majority of galaxies in the sample, which live at the centers of their groups. Finally, with group distances and magnitudes in hand, our sample has 5,110 massive early-type galaxies with mag and Mpc.

To make our sample cleaner, and enable the clean conversion of stellar mass to an inferred black hole mass, we select only early-type galaxies (elliptical or S0), using morphologies from Hyperleda (Paturel et al., 2003). We visually inspected of the galaxies, and found the sample to be very clean.

The galaxy catalog is augmented by adding nearby galaxies which host SMBHs with dynamical mass measurements (Schutz & Ma, 2016; Ma et al., 2014). 33 such galaxies are added, 9 of which were not previously in the 2MASS sample (due to e.g. low luminosity). For these galaxies, we use the measured SMBH mass instead of inferring it from an empirical scaling relation.

Supermassive Black Hole Masses We estimate the SMBHB mass by converting the 2MASS K-band luminosity to the stellar mass, take this stellar mass to be the bulge mass, and further apply the scaling relation from Ref (McConnell & Ma, 2013) to obtain the total BH mass. When computing the total SMBHB mass this way, we incorporate an intrinsic scatter in the relation, , as follows: each logarithmic realization of is a random draw from a normal distribution with mean from , and . The answer is exponentiated. Since the exponent is drawn from a normal distribution, it follows that the total BH mass follows a log-normal distribution over many realizations, thus favoring smaller masses. We therefore consider the binary mass to be conservative.

IPTA predictions We add four pulsars per year to the existing 47 IPTA pulsars, chosen from sky locations in the field of view of the current IPTA telescopes. The white noise level is modeled as a combination of radiometer noise and jitter noise. The radiometer noise is estimated as the harmonic mean of the measured error bars (for each backend and observing frequency) to avoid overestimation due to low signal-to-noise ratio time of arrivals (TOAs), which would have large error bars. Jitter noise is obtained for each observing frequency (lcc+16). We then compute an infinite frequency TOA uncertainty from the low and high frequency noise estimates in order to simulate dispersion measure fitting. From 2016 onward we add four pulsars per year using the median white noise uncertainty of the existing pulsars in the array, typically around 300 ns. Further, we assume a new wide-band timing backend installed in 2018 at Arecibo and the Green Bank Telescope, which reduces the white noise rms by a factor of .

This backend upgrade is the dominant factor in the improved white noise detection curves, Figure 3(b), due to the fact that the signal-to-noise, , of a CGW detection is roughly , assuming the N pulsars have identical intrinsic properties, is the length of the dataset, is the cadence of the observation and is the white noise rms (Arzoumanian et al., 2014). For red noise with spectral index , at . There are of course other factors which motivate a large and expanding PTA, including the geometric term from the antenna beam pattern , where is the direction of propagation of the GW, is the direction to the pulsar (Detweiler, 1979; Arzoumanian et al., 2014; Mingarelli & Sidery, 2014), and is the GW polarization. Therefore, when (when the direction to the source, , is aligned with the pulsar) the response is maximal222When this alignment is exact there is no GW signal due to surfing effects but when the alignment is almost perfect the response is maximal, as observed in Figure 1(b).

Projections are made for 15, 20 and 25 year datasets with various false alarm probabilities. One can covert from the FAP, , to multiples of the standard deviation via erf. For example, a FAP of , assuming a Gaussian distribution.

Currently, a full Bayesian analysis is computationally intractable when performing these kinds of detection sensitivity analyses. Since we are performing these analyses as a function of CGW frequency and PTA configuration, we must compute the detection statistic (Bayes factor for Bayesian analysis, FAP for frequentist analysis) millions of times. In comparison to the frequentist FAP statistic, which takes fractions of a second to compute, the Bayes factor computation requires several hours to compute. Thus, for the number of simulations required in this work, a full Bayesian analysis for all simulations is intractable at this time, and we have instead used a frequentist proxy assuming only white noise in order to emulate the possible resolving capability of the Bayesian analysis.

Note that for pulsars with strong red noise, 100 times fewer sources are detected with white-noise dominated pulsars, Figure 3 and Table 1. The presence of unresolved red noise alters the position of the minimum frequency and achievable strain, shifting the lowest GW frequency accessible by PTAs to higher frequencies. Hence, the sources must be much closer (we find that the majority of these are within 20 Mpc) in order to be detected. We also find that, on average, the galaxies in our catalog underwent a major merger at , Figure 4(c).

Generating GW sky maps GW sky maps are created by interpolating a set of 128 original data points at 87 GW frequencies. We interpolate between the points using a bivariate spline approximation over a rectangular mesh on a sphere, and project the resulting sky onto a healpix map. This is done for 87 GW frequencies, with the resulting in 87 healpix maps saved as FITS files which are freely available. It is possible that more sources could be detected if future data points are extended to a greater range in declination. Interpolation errors close to the poles required us to make a hard cut at declinations of  degrees, eliminating potential galaxies as CGW sources.

Galaxies Hosting SMBHBs in the PTA band

A preliminary study of local potential nanohertz CGW sources was carried out by Ref (Simon et al., 2014). The authors assembled a complete galaxy catalog out to  Mpc, including galaxies with SMHBs with . However, in this study it was assumed that there was an equal probability for all galaxies to host a SMBHB, and that this was an equal mass binary. Moreover, the authors of Ref (Rosado & Sesana, 2014) employed a top-down approach to predict GW skies by creating a simulated galaxy catalog, and matching SMBHB merger rates from the Millennium Simulation (Springel et al., 2005) to the Sloan Digital Sky Survey (Abazajian et al., 2009)(SDSS), in an effort to identify the characteristics of potential SMBHB host galaxies. However, SDSS has a limited sky coverage, and does not allow a full-sky investigation of the loudest potential GW sources in the nearby Universe. The work presented here marks the use a galaxy survey to identify local massive galaxies, to assign each galaxy a probability that it hosts a SMBHB based on Illustris galaxy-galaxy merger rates, and to estimate the time to detection of these sources with PTAs.

Here we give an example of how probabilities are assigned to galaxies in the catalog derived from 2MASS, in our bottom-up approach. Consider, for example, galaxy NGC 4594333NGC 4594 has a dynamically measured SMBH mass, however, for illustration purposes we show how its mass would be assigned via the empirical relation, which has . Through the empirical scaling relation (Cappellari, 2013), the stellar mass is , and via with scatter , a SMBHB with total mass . The mass ratio of the binary, , is drawn from a log-normal distribution in [0.25, 1.0], from which we randomly draw . The chirp mass is therefore .

The dynamical friction timescale (Binney & Tremaine, 2008) is computed assuming the Coulomb logarithm is :

 tdf=264 Myr(a2 kpc)2(vc250 km/s)(108 M⊙M2), (4)

where is the mass of the secondary BH, with (Zahid et al., 2016), and is the galaxy’s effective radius, , Eq. 4 Ref (Dabringhausen et al., 2008), Figure 7.

We scale the stellar mass of the galaxy, , by a factor of 0.7, Ref (De Lucia & Blaizot, 2007; McBride et al., 2009) (see Supplementary Figure 4), to estimate the mass of the descendent galaxy when we begin the SMBHB evolution, in the dynamical friction phase. This scaled is also used to estimate the velocity dispersion . The parameters drawn are therefore: ,  km/s,  kpc, which when input in Eq. (4) yield a dynamical friction timescale of  Gyr.

The stellar hardening timescale (Quinlan, 1996) is computed as in Ref (Sesana & Khan, 2015), Eqs (6) and (7):

 a∗,gw=(64σinfM1M2M5Hρinf)1/5 ; tsh=σinfHρinfa∗,gw, (5)

where is computed via (Zahid et al., 2016) as above, and is the density profile evaluated at the influence radius, with (corresponding to a Hernquist profile (Dehnen, 1993)), more details given in Ref (Sesana & Khan, 2015). Here we assume that the binary eccentricity and the hardening constant , and find that the hardening timescale is  Gyr.

The sum of the dynamical friction and hardening timescales is  Gyr, or using Planck (Planck Collaboration et al., 2016) cosmological parameters. The cumulative galaxy-galaxy merger rate (Ref (Rodriguez-Gomez et al., 2015) Table 1) requires only and as inputs, since the dependence on is removed by integrating it between ; see Equation (1). We scale by 0.7, and find that is .

Finally, we compute the time this binary will spend in the PTA band. The time to coalescence of the binary, , is taken from the lower limit of the PTA band:  Hz, resulting in  Myr. By Eq. (1), the probability of NGC 4594 hosting a SMBHB in the PTA band in this particular realization is  mergers/Gyr) .

Computer Code and Data Availability A series of Jupyter Notebooks which reproduce our figures and all our results are freely available https://github.com/ChiaraMingarelli/nanohertz_GWs. Through this git repository, we also share the underlying detection curves and GW sensitivity sky maps.

Correspondence and requests for materials should be addressed to Chiara Mingarelli.

The authors thank S. Babak, J. Verbiest, D. Kaplan, E. Barr, K. Górski and E. Sheldon for useful discussions. The authors are grateful for open source scientific tools for Python (Jones et al., 2001–; The Astropy Collaboration et al., 2013; van der Walt et al., 2011; Barrett et al., 2005; Górski et al., 2005; Perez & Granger, 2007) and L. Singer’s open source plot.py (Singer et al., 2014). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by NASA and the NSF. C.M.F.M. was supported by a Marie Curie International Outgoing Fellowship within the European Union Seventh Framework Programme. S.R.T was partly supported by appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities and the Universities Space Research Association through a contract with NASA. AS is supported by a URF of the Royal Society. Part of these computations were performed on the Zwicky cluster at Caltech, which is supported by the Sherman Fairchild Foundation and NSF award PHY-0960291. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. This work has also been supported by NSF award 1458952. The NANOGrav project receives support from National Science Foundation Physics Frontier Center award number 1430284. Author contributions C.M.F.M. modeled the supermassive black hole evolution, developed and ran the Monte Carlo simulations used here to explore their evolution, analyzed the resulting data, produced all of the figures and table, and was the primary author of this paper. C.M.F.M, T.J.W.L and S.B.S. developed the concept of this work. A.S., C.P.M., S.C. and T.J.W.L. advised on supermassive black hole astrophysics and helped to interpret the results. J.E.G. and S.C. assembled and inspected the galaxy catalog. J.A.E. developed the time to detection methods for the IPTA. S.R.T. helped to develop the methods used to turn continuous gravitational-wave sources into a gravitational wave background, and explore its angular power spectrum.

## References

• Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
• Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, ArXiv e-prints, arXiv:1702.00786
• Arzoumanian et al. (2014) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2014, ApJ, 794, 141
• Arzoumanian et al. (2016) —. 2016, ApJ, 821, 13
• Babak et al. (2016) Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665
• Barrett et al. (2005) Barrett, P., Hunter, J., Miller, J. T., Hsu, J.-C., & Greenfield, P. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 91
• Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
• Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
• Booth et al. (2009) Booth, R. S., de Blok, W. J. G., Jonas, J. L., & Fanaroff, B. 2009, ArXiv e-prints, arXiv:0910.2935
• Cappellari (2013) Cappellari, M. 2013, The Astrophysical Journal Letters, 778, L2
• Cornish & van Haasteren (2014) Cornish, N. J., & van Haasteren, R. 2014, ArXiv e-prints, arXiv:1406.4511
• Crook et al. (2007) Crook, A. C., Huchra, J. P., Martimbeau, N., et al. 2007, ApJ, 655, 790
• Dabringhausen et al. (2008) Dabringhausen, J., Hilker, M., & Kroupa, P. 2008, MNRAS, 386, 864
• De Lucia & Blaizot (2007) De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2
• Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
• Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
• Detweiler (1979) Detweiler, S. 1979, 234, 1100
• Ellis et al. (2012) Ellis, J. A., Siemens, X., & Creighton, J. D. E. 2012, 756, 175
• Estabrook & Wahlquist (1975) Estabrook, F. B., & Wahlquist, H. D. 1975, General Relativity and Gravitation, 6, 439
• Ferrarese & Ford (2005) Ferrarese, L., & Ford, H. 2005, Space Science Reviews, 116, 523
• Gair et al. (2014) Gair, J., Romano, J. D., Taylor, S., & Mingarelli, C. M. F. 2014, Phys. Rev. D, 90, 082001
• Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175
• Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
• Hellings & Downs (1983) Hellings, R. W., & Downs, G. S. 1983, 265, L39
• Huchra et al. (2012) Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26
• Janssen et al. (2015) Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 37
• Jarrett et al. (2000) Jarrett, T. H., Chester, T., Cutri, R., et al. 2000, AJ, 119, 2498
• Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python, [Online; accessed 2016-08-24]
• Lentati et al. (2015) Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2015, ArXiv e-prints, arXiv:1504.03692
• Ma et al. (2014) Ma, C.-P., Greene, J. E., McConnell, N., et al. 2014, ApJ, 795, 158
• McBride et al. (2009) McBride, J., Fakhouri, O., & Ma, C.-P. 2009, MNRAS, 398, 1858
• McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, 764, 184
• Mingarelli & Sidery (2014) Mingarelli, C. M. F., & Sidery, T. 2014, Phys. Rev. D, 90, 062011
• Mingarelli et al. (2013) Mingarelli, C. M. F., Sidery, T., Mandel, I., & Vecchio, A. 2013, 88, 062005
• Paturel et al. (2003) Paturel, G., Petit, C., Prugniel, P., et al. 2003, A&A, 412, 45
• Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science Engineering, 9, 21
• Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
• Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
• Quinlan (1996) Quinlan, G. D. 1996, New Astronomy, 1, 35
• Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49
• Rosado & Sesana (2014) Rosado, P. A., & Sesana, A. 2014, MNRAS, 439, 3986
• Sazhin (1978) Sazhin, M. V. 1978, Soviet Ast., 22, 36
• Schutz & Ma (2016) Schutz, K., & Ma, C.-P. 2016, MNRAS, 459, 1737
• Sesana & Khan (2015) Sesana, A., & Khan, F. M. 2015, ArXiv e-prints, arXiv:1505.02062
• Sesana et al. (2016) Sesana, A., Shankar, F., Bernardi, M., & Sheth, R. K. 2016, MNRAS, 463, L6
• Shannon et al. (2015) Shannon, R. M., Ravi, V., Lentati, L. T., et al. 2015, Science, 349, 1522
• Simon et al. (2014) Simon, J., Polin, A., Lommen, A., et al. 2014, ApJ, 784, 60
• Singer et al. (2014) Singer, L. P., Price, L. R., Farr, B., et al. 2014, ApJ, 795, 105
• Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
• Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
• Taylor & Gair (2013) Taylor, S. R., & Gair, J. R. 2013, 88, 084001
• Taylor et al. (2016) Taylor, S. R., Vallisneri, M., Ellis, J. A., et al. 2016, ApJ, 819, L6
• Taylor et al. (2015) Taylor, S. R., Mingarelli, C. M. F., Gair, J. R., et al. 2015, Physical Review Letters, 115, 041101
• The Astropy Collaboration et al. (2013) The Astropy Collaboration, Robitaille, Thomas P., Tollerud, Erik J., et al. 2013, A&A, 558, A33
• The NANOGrav Collaboration et al. (2015) The NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ, 813, 65
• van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
• Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267
• Wrobel & Nyland (2012) Wrobel, J. M., & Nyland, K. 2012, AJ, 144, 160
• Zahid et al. (2016) Zahid, H. J., Geller, M. J., Fabricant, D. G., & Hwang, H. S. 2016, ApJ, 832, 203
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters