Buckets:
THE NANOGRAV NINE-YEAR DATA SET:
LIMITS ON THE ISOTROPIC STOCHASTIC GRAVITATIONAL WAVE BACKGROUND
Z. ARZOUMANIAN1, A. BRAZIER2, S. BURKE-SPOLAOR3,27, S. J. CHAMBERLIN4, S. CHATTERJEE2, B. CHRISTY5, J. M. CORDES2,
N. J. CORNISH6, P. B. DEMOREST3, X. DENG7, T. DOLCH2,25, J. A. ELLIS8,28, R. D. FERDMAN9, E. FONSECA10,
N. GARVER-DANIELS11, F. JENET12, G. JONES13, V. M. KASPI9, M. KOOP7, M. T. LAM2, T. J. W. LAZIO8, L. LEVIN11,
A. N. LOMMEN5, D. R. LORIMER11, J. LUO12, R. S. LYNCH14, D. R. MADISON2, M. A. MC LAUGHLIN11, S. T. McWILLIAMS11,
C. M. F. MINGARELLI15,16,29, D. J. NICE17, N. PALLIYAGURU11, T. T. PENNUCCI18, S. M. RANSOM19, L. SAMPSON6,
S. A. SANIDAS20,21, A. SESANA22, X. SIEMENS4, J. SIMON4, I. H. STAIRS10, D. R. STINEBRING23, K. STOVALL24, J. SWIGGUM11,
S. R. TAYLOR8, M. VALLISNERI8, R. VAN HAASTEREN8,28, Y. WANG26, & W. W. ZHU10,16
*Author order alphabetical by surname
1Center for Research and Exploration in Space Science and Technology and X-Ray Astrophysics Laboratory,
NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA
2Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
3National Radio Astronomy Observatory, 1003 Lopezville Rd., Socorro, NM 87801, USA
4Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA
5Department of Physics and Astronomy, Franklin & Marshall College, P.O. Box 3003, Lancaster, PA 17604, USA
6Department of Physics, Montana State University, Bozeman, MT 59717, USA
7Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802, USA
8Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
9Department of Physics, McGill University, 3600 University St., Montreal, QC H3A 2T8, Canada
10Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
11Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
12Center for Gravitational Wave Astronomy, University of Texas at Brownsville, Brownsville, TX 78520, USA
13Department of Physics, Columbia University, New York, NY 10027, USA
14National Radio Astronomy Observatory, P.O. Box 2, Green Bank, WV 24944, USA
15TAPIR, MC 350-17, California Institute of Technology, Pasadena, CA 91125, USA
16Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, D-53121 Bonn, Germany
17Department of Physics, Lafayette College, Easton, PA 18042, USA
18University of Virginia, Department of Astronomy, P.O. Box 400325, Charlottesville, VA 22904, USA
19National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
20Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands
21Jodrell Bank Centre for Astrophysics, University of Manchester, Manchester, M13 9PL, United Kingdom
22School of Physics & Astronomy, University of Birmingham, Birmingham, B15 2TT, UK
23Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA
24Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87131, USA
25Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, Michigan 49242, USA
26School of Physics, Huazhong University of Science and Technology, 1037 Luoyu Road, Wuhan, Hubei Province 430074, China
27Jansky Fellow
28Einstein Fellow and
29Marie Curie Fellow
Draft version August 13, 2015
ABSTRACT
We compute upper limits on the nanohertz-frequency isotropic stochastic gravitational wave background (GWB) using the 9-year data release from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration. Well-tested Bayesian techniques are used to set upper limits for a GWB from supermassive black hole binaries under power law, broken power law, and free spectral coefficient GW spectrum models. We place a 95% upper limit on the strain amplitude (at a frequency of $\text{yr}^{-1}$ ) in the power law model of $A_{\text{gw}} < 1.5 \times 10^{-15}$ . We find that the upper limit on the individual spectral components is consistent with a white noise spectrum at frequencies greater than the inverse observation time. For a broken power law model, we place priors on the strain amplitude derived from simulations of Sesana (2013) and McWilliams et al. (2014). Using Bayesian model selection we find that the data favor a broken power law to a pure power law with odds ratios of 22 and 2.2 to one for the McWilliams and Sesana prior models, respectively. The McWilliams model is essentially ruled out by the data, and the Sesana model is in tension with the data under the assumption of a pure power law. Using the broken power-law analysis we construct posterior distributions on environmental factors that drive the binary to the GW-driven regime including the stellar mass density for stellar-scattering, mass accretion rate for circumbinary disk interaction, and orbital eccentricity for eccentric binaries, marking the first time that the shape of the GWB spectrum has been used to make astrophysical inferences. We then place the most stringent limits so far on the energy density of relic GWs, $\Omega_{\text{gw}}(f)h^2 < 4.2 \times 10^{-10}$ , yielding a limit on the Hubble parameter during inflation of $H_* = 1.6 \times 10^{-2} m_{\text{Pl}}$ , where $m_{\text{Pl}}$ is the Planck mass. Our limit on the cosmic string GWB, $\Omega_{\text{gw}}(f)h^2 < 2.2 \times 10^{-10}$ , translates to a conservative limit of $G\mu < 3.3 \times 10^{-8}$ – a factor of 4 better than the joint Planck and high- $l$ Cosmic Microwave Background data from other experiments.1. INTRODUCTION
The 1967 discovery of pulsars (Hewish et al. 1968) inaugurated the age of high-energy astrophysics, and in time it led to the indirect but compelling confirmation that compact binaries lose energy to gravitational waves (GWs) in accordance with general-relativistic predictions (Taylor & Weisberg 1989). It is therefore quite fitting that pulsars should now feature prominently in the quest to directly detect low-frequency GWs, and to use them as probes of the unseen, dynamical Universe. In this article, we describe how one can make astrophysical inferences about the underlying GW source population using the shape of the nanoHertz GW spectrum. Using the 9-year NANOGrav data release, Arzoumanian et al. (2015), we provide the first constraints on environmental factors which may be contributing the measured shape of the spectrum. We expect these novel techniques to become standard analysis tools, and as such, provide detailed descriptions of our methods. We also provide a summary of our results here, so that one may get an overview of our new methods and results before delving into the details.
1.1. Gravitational Wave Detection with Pulsar Timing Arrays
Sazhin (1978) and Detweiler (1979) realized that GWs could manifest as otherwise unexplained residuals in the times of arrival (TOAs) of pulsar signals after subtracting a deterministic timing model. This timing model accounts for the intrinsic evolution of pulsar spin, radio frequency-dependent delays to interstellar propagation effects, the astrometric time delays and advances due to the relative motion of the pulsar system with respect to the Earth (indeed, to the observatory), as well as orbital-kinematic and light-propagation effects for pulsars that orbit a binary companion (see, e.g., Lorimer & Kramer 2005). Foster & Backer (1990) pointed out that the timing residuals from an array of pulsars (a pulsar timing array, or PTA) could be analyzed coherently to separate GW-induced residuals, which have distinctive correlations among different pulsars (Hellings & Downs 1983), from other systematic effects, such as clock errors or delays due to light propagation through the interstellar medium.
Today, three international consortia [NANOGrav (McLaughlin 2013), the EPTA (Kramer & Champion 2013), and the PPTA (Hobbs 2013)] are more than ten years into extensive campaigns to search for GWs by timing dozens of individual millisecond pulsars (MSPs), in which the best-timed have rms residuals less than 100 ns (corresponding to GW strain sensitivities $\sim 10^{-15}$ ). The three PTAs collaborate and share data under the aegis of the International Pulsar Timing Array (Hobbs et al. 2010).
In order to robustly detect GWs, one must have a thorough understanding of the underlying noise in the pulsar timing data (see e.g. Cordes 2013; Stinebring 2013, for a detailed review). Template matching errors due to radiometer noise are uncorrelated in both time and frequency, but pulse-jitter noise (Cordes & Shannon 2010) appears to affect all TOAs obtained simultaneously in different frequency channels. Correlated timing noise with a red power spectrum occurs to varying degree in different pulsars. Spin noise (Shannon & Cordes 2010) is achromatic and is much smaller in MSPs compared to objects with stronger magnetic fields and longer spin periods.
Chromatic red noise due to propagation through intervening plasmas (ISM, interplanetary medium, and ionosphere) may also be present if dispersive delays are not removed perfectly or if scattering and refraction effects contribute significantly.
1.2. The stochastic GW background from SMBHBs
PTAs are most sensitive to GWs with frequencies on the order of the inverse timespan of timing observations, where TOA measurement noise averages out most efficiently. The strongest expected sources in this band are supermassive black hole binaries (SMBHBs) with masses of $10^8$ – $10^{10} M_\odot$ , out to $z \simeq 1$ (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003). The binaries form after the hierarchical mergers (Sesana et al. 2004, 2008) of galaxies hosting individual SMBHs (as most galaxies are thought to do, cf. Kormendy & Ho 2013). The most massive and nearest binaries may be detected individually by PTAs through their continuous GW emission. Moreover, the cosmic population of SMBHBs may be observed collectively as a stochastic GW background composed of the incoherent superposition of signals from the binaries. Rosado et al. (2015) discuss which detection is likely to come first.
The main focus of this article is the measurement of an isotropic stochastic GW background from SMBHBs. Anisotropic-background searches based on the formalism and techniques developed by Gair et al. (2014); Taylor & Gair (2013); Mingarelli et al. (2013) are currently underway and will be the subject of a follow-up paper. We briefly consider stochastic backgrounds from relic (or primordial) GWs as well as backgrounds from cosmic (super)strings as complementary studies. Our main focus is to demonstrate how, even in the absence of a positive detection, PTA data can be used to constrain and characterize the astrophysical processes and SMBHB source populations that give rise to the GW background.
The simplest characterization of the stochastic GW background – a power-law Gaussian process with isotropic inter-pulsar correlations – applies if:
- all binaries are assumed to have circular orbits (so each component signal is instantaneously monochromatic);
- all binaries evolve through the PTA band due purely to GW emission, as opposed to environmental effects such as interactions with nearby gas or with stars in the galactic nucleus;
- all binaries are distributed isotropically across the sky in sufficient numbers to fulfill the central limit theorem at all frequencies.
Under these conditions, the observed timing residuals due to the GW background are described fully by the (cross-) power spectral density
where $a$ and $b$ range over the pulsars in the array, $\gamma = 13/3$ for a background composed of SMBHBs, and $\Gamma_{ab}$ is the Hellings–Downs (1983) isotropic correlation coefficient for pulsars $a$and $b$ (a function of the separation angle between their lines of sight, which is normalized so that $\Gamma_{aa} = 1$ ; see also Min-garelli & Sidery 2014). Power-law GW backgrounds are also described (independently of observations) in terms of their characteristic strain
which is related to Eq. (1) by $S_{ij}(f) = \Gamma_{ij} h_c(f)^2 / (12\pi^2 f^3)$ and $\gamma = 3 - 2\alpha$ ( $\alpha = -2/3$ for SMBHBs).
Recent predictions for the value of $A_{\text{gw}}$ , based on models of SMBH–galaxy coevolution and on observational constraints of galaxy assembly and SMBH mass functions, range between $\sim 10^{-15}$ and $10^{-14}$ (McWilliams et al. 2014; Sesana 2013; Ravi et al. 2014—hereafter MOP14, S13, and RWS14). The MOP14 model assumes that all SMBHBs are in circular orbits evolving under GW emission alone, as well as a single black hole–host correlation from McConnell & Ma (2013), yielding an estimate of the stochastic GW background that is roughly four times as large as S13 and RWS14. The S13 model uses a wide variety of galaxy merger rates and empirical black hole–host relations to yield a collection of phenomenological SMBHB merger rates, which are used to compute a distribution of possible GW signals. Note that for this paper we consider distributions from S13 that only use the black hole–host galaxy relations of McConnell & Ma (2013) and Graham & Scott (2013). The RWS14 model also assumes the black-hole–host correlation of McConnell & Ma (2013) but includes the possibility of the SMBHBs evolving in stellar environments, and accounts for non-zero binary eccentricities. Thus, some of these models predict spectral densities that deviate from straight power-law behavior at low frequencies; in that case, we refer the fiducial $A_{\text{gw}}$ to their value at a frequency of $\text{yr}^{-1}$ . Finally, recent results, (Kormendy & Ho 2013), indicate that the black hole–host correlation’s normalization is being revised to larger values with more observations indicating that an even stronger GWB may be expected; however, for this work we use the published results based on McConnell & Ma (2013) to make the most fair comparison among models.
1.3. New results in this paper
Over the last few years, the three PTAs have reported ever improving upper limits on the GW backgrounds of the form (1): $A_{\text{gw}} < 7 \times 10^{-15}$ (NANOGrav, Demorest et al. 2013, hereafter DFG13), $6 \times 10^{-15}$ (EPTA, van Haasteren et al. 2011), $3 \times 10^{-15}$ (EPTA, Lentati et al. 2015, hereafter LTM15), $2.4 \times 10^{-15}$ (PPTA, Shannon et al. 2013), all quoted at 95% confidence and a reference frequency of $\text{yr}^{-1}$ , although differences in the statistical analyses and in the availability and selection of pulsar datasets mean that these numbers are not entirely homogeneous. It is clear that there is significant tension between these observational limits and astrophysical expectations for $A_{\text{gw}}$ . It is important to note that a limit on $A_{\text{gw}}$ does not translate directly to a limit on the SMBHB population because of the finite number of pulsars that contribute to the limit and the stochasticity of the GW signal itself.
This statement can be made more precise. In Sec. 4.2 of this paper we report a new 95% upper limit $A_{\text{gw}} < 1.5 \times 10^{-15}$ , obtained from the Bayesian analysis of NANOGrav’s 9-year, 37-pulsar dataset released in 2015 (Arzoumanian et al. 2015). This limit is five times more constraining than the same analysis applied to NANOGrav’s 5-year dataset (DFG13) (see the
top of Fig. 4 for a comparison of the two posterior probability distributions). Now, following Shannon et al. (2013), we can assess the consistency of our result with astrophysical GW-background models. We find a 0.8% probability that the observed $A_{\text{gw}}$ , as characterized probabilistically by its Bayesian posterior, is drawn from the amplitude distribution developed in MOP14, and a 20% probability that it is drawn from the (very similar) RWS14 and S13 distributions. Correspondingly, the two bottom panels of Fig. 4 show that 9-year observations update significantly the MOP14 and RWS14/S13 amplitude priors, much more so than our 5-year dataset.
In Sec. 4.1 we report also a frequentist, optimal-statistic 95% upper limit $A_{\text{gw}} < 1.3 \times 10^{-15}$ , a fivefold improvement on the analogous result of DFG13; however, the optimal statistic is problematic in the presence of marginally-detectable GW signals, so we offer it only as a proxy for the improving sensitivity of NANOGrav’s observations.
Stochastic GW backgrounds in the PTA band may also originate from quantum fluctuations amplified during inflation (Grishchuk 2005) and from topological broken-symmetry remnants such as (Damour & Vilenkin 2001; Olmez et al. 2010), for which Eq. (1) applies with $\gamma = 5$ (depending on the equation of state $w$ ) and $16/3$ , respectively.
In Sec. 5.2.1 we obtain 95% upper limits $A_{\text{gw}} < 8.1 \times 10^{-16}$ for relic GWs, corresponding to energy-density limits $\Omega_{\text{gw}}(f = 1/\text{yr})h^2 < 4.2 \times 10^{-10}$ , where $h$ parametrizes the Hubble constant $H_0 \equiv h \times 100 \text{ km/s/Mpc}$ . This limit is a factor of 2.9 better than limit reported by the EPTA in Lentati et al. (2015). We then obtain limits on the Hubble parameter during inflation, $H_* = 1.6 \times 10^{-2} m_{\text{Pl}}$ , where $m_{\text{Pl}} \equiv 1/\sqrt{G}$ is the Planck mass, using the method developed by Zhao (2011).
In Sec. 5.2.2 cosmic strings, we find $A_{\text{gw}} < 6 \times 10^{-16}$ and $\Omega_{\text{gw}}(f = 1/\text{yr})h^2 < 2.23 \times 10^{-10}$ , corresponding to a conservative limit on the string tension of $G\mu < 3.3 \times 10^{-8}$ – a factor of 4 better than the joint Planck and high- $l$ Cosmic Microwave Background data from other experiments. If we then restrict ourselves to a GWB produced by the production of large cosmic string loops, as described by Blanco-Pillado et al. (2014), then our string tension limit is much more restrictive: $G\mu < 1.3 \times 10^{-10}$ , a factor of 6.6 times more constraining than an identical analysis performed using the EPTA limit.
1.4. Astrophysical Inference
The mismatch between observations and expectations can be explained in different ways. First, it is possible that the astrophysical models and the three assumptions listed above are correct, but the background amplitude realized in nature lies in the tails of the predicted distributions. This hypothesis obviously wanes as upper limits get more stringent.
Second, it is possible that some of the inputs of the astrophysical models are not estimated correctly; we can then use GW observations to constrain these inputs. For example, in Sec. 5.1.1 we assume measurements of the galaxy mass function and merger rate, and we constrain the scaling between the galaxy bulge mass and the central SMBH mass, which affects the observed $A_{\text{gw}}$ most significantly, through the distribution of binary chirp masses. We find a slight inconsistency between the scaling relation reported by Kormendy & Ho (2013) and our limit, whereas the McConnell & Ma (2013) estimate is consistent within its margin of error. The McConnell & Ma (2013) black hole–stellar velocity dispersion relation underlies the MOP14 predictions, while S13 and RWS14 take intoaccount a variety of alternative black hole–host estimates.
Third, the simple GW-background characterization that yields Eq. (1) may not be realistic, because SMBHs may form with significant eccentricity and retain it into the PTA band, distributing GW emission over a range of frequency harmonics. Furthermore, environmental effects (interactions with stars on centrophilic orbits in galactic nuclei, or with circumbinary gas disks) can accelerate the transit of individual binary systems through the PTA band (see Secs. 5.1.3 and 5.1.2 for details and references). These environmental effects deplete the GW background at low frequencies where PTA measurements are most sensitive (i.e., frequencies $\sim$ the inverse observation timespan), so PTA upper limits may yet be compatible with the MOP14/S13/RWS14 predictions at the higher frequencies where the $\gamma = 13/3$ power law is realized.
To investigate this point, in Sec. 4.2.2 we reanalyze the NANOGrav data using a phenomenological $S_{ij}(f)$ in the form of an inflected power law, parametrized by a turnover frequency $f_{\text{bend}}$ and by a shape parameter $\kappa$ , as proposed by Sampson et al. (2015) [see Eq. (24)]. By combining this enhanced GW-background model with MOP14 and S13/RWS14 amplitude priors, we conclude that the data prefer an inflected spectrum to a moderate degree for MOP14, however this preference is very weak for S13/RWS14 models. Quantitatively, the Bayes factors between enhanced and pure-power-law spectral models for each of the two priors are 22 and 2.2, respectively; graphically, the shading in Fig. 5 represents the frequency-by-frequency posterior probability density for the GW spectrum, which appears significantly inflected for MOP14, and only slightly so for S13/RWS14. The data are not sufficiently informative to constrain the amplitude and shape of the spectrum jointly in the absence of a compact prior, so we cannot produce a unique metric of consistency, as we did above in the case of the simple power-law spectrum.
Beyond this phenomenological characterization, the joint $(A_{\text{gw}}, f_{\text{bend}}, \kappa)$ posteriors can be mapped into constraints for the SMBHB eccentricities (which we do in Sec. 5.1.3) and for the astrophysical variables that govern environmental interactions (Sec. 5.1.2). For the former, we assume for simplicity that all binaries had the same eccentricity $e_0$ when the semi-major axis of their orbits was 0.01 pc and that they evolved purely by GW emission since; we follow Huerta et al. (2015) to construct eccentric binary populations and GW strain spectra. The resulting posteriors on $e_0$ indicate that $e_0 \gtrsim 0.7$ is preferred for the MOP14 prior, and $e_0 \gtrsim 0.5$ is preferred for S13/RWS14 (though still consistent with smaller values). These limits suggest that either SMBHs form with rather high $e_0$ , or that binary eccentricity is not a good explanation for the mismatch between $A_{\text{gw}}$ observations and predictions.
To characterize environmental interactions (see Sec. 5.1.2), we compute the evolution of orbital frequency due to stellar scattering events and to circumbinary gas disk interactions as $df/dt \propto f^{1/3}$ and $f^{4/3}$ , respectively, corresponding to $\kappa = 10/3$ and $7/3$ (since for GW-driven evolution $df/dt \propto f^{11/3}$ ); the frequency $f_{\text{bend}}$ then marks the transition between environmentally- and GW-driven evolution. For the case of stellar scattering, $f_{\text{bend}}$ depends most significantly on the mass density $\rho$ of galactic-core stars. Astrophysical estimates for $\rho$ are quite uncertain with typical values between $10 - 10^4 M_{\odot} \text{ pc}^{-3}$ Dotti et al. (2007) assuming that a majority of our GW sources come from merging elliptical galaxies. Under several simplifying assumptions (e.g., that all binaries have circular orbits, that all galaxies have comparable
densities, and that only a single environmental effect is active), we find that $\rho \gtrsim 10^4 M_{\odot} \text{ pc}^{-3}$ is strongly preferred for the MOP14 amplitude prior, and the data is unconstraining for the S13/RWS14 prior. For the circumbinary disk case, $f_{\text{bend}}$ depends on the accretion rate on to the primary (most massive) BH, $\dot{M}1$ . The accretion rate of the primary BH is a function of $M_1$ , the mass of the primary BH, $\epsilon$ , the radiative efficiency parameter with a canonical value of $\epsilon = 0.1$ and $\kappa{\text{opp}}$ is the disk opacity, $\dot{M}1 \propto M_1 \epsilon^{-1} \kappa{\text{opp}}^{-1}$ . Hence $\dot{M}1$ takes on a range of values, typically $10^{-3} M{\odot} \text{ yr}^{-1} - 1 M_{\odot} \text{ yr}^{-1}$ , see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; Goicovic et al. 2015. In our analysis, we find that the accretion rate $\dot{M}1 \gtrsim 10^{-1} M{\odot} \text{ yr}^{-1}$ is strongly preferred for the MOP14 amplitude prior, and again, the data is unconstraining for the S13/RWS14 prior.
1.5. Plan of the Paper
This paper presents the first analysis that characterizes the spectral amplitude and shape of the GW background; in some cases we find tension between observations and predictions for the GW background from SMBHs. When informed with astrophysical amplitude priors, the data favor phenomenological models that include an inflection point over pure GW-driven power laws. We attempt to interpret the phenomenological posteriors in terms of binary-eccentricity or environmental effects by carrying out our analyses in sequence, investigating different effects separately. Admittedly, we rely on simple analytical models of eccentricity and environments; to some extent, this simplicity becomes necessary if we are to draw any astrophysical conclusion from data that are not yet very informative. Therefore, it will be extremely interesting to improve the sophistication of our analysis, and to apply it to longer, richer datasets, such as the upcoming NANOGrav 11-year dataset, as well as the multiple-PTA datasets assembled by the International Pulsar Timing Array.
The rest of this paper is organized as follows: in Sec. 2 we describe the NANOGrav 9-year pulsar-timing and dataset; in Sec. 3 we discuss our signal and noise models, as well as the statistical framework of our analysis; in Sec. 4 we document important implementation details, and report our results in detail; in Sec. 5 we discuss the astrophysical interpretation of our analysis, and derive limits on astrophysical and cosmological quantities; in Sec. 6 we offer our conclusions.
2. OBSERVATIONS
This paper uses the observations from the NANOGrav 9-year data release, recently presented in Arzoumanian et al. (2015), which contains observations made over a time span from 2004 to 2013. Initially the array consisted of 15 pulsars, and it grew to 37 pulsars over the course of the project. The first five years of data on 17 pulsars constituted the NANOGrav 5-year data release, previously reported by DFG13. In this release, all data have been reprocessed. We give a brief overview of the dataset in this section; for details we refer the reader to Arzoumanian et al. (2015).
2.1. Observatories
The 305 m William E. Gordon Telescope of the Arecibo Observatory (AO) was used to observe pulsars with declinations in the range $0^\circ < \delta < 39^\circ$ ; pulsars outside of this range were observed with the 100 m Robert C. Byrd Green Bank Telescope (GBT) of the National Radio Astronomy Observatory (NRAO). The pulsars PSRs J1713+0747 and B1937+21were observed at both AO and the GBT. The typical observation cadence was about once every month.
At AO, four receivers were used: 327 MHz, 430 MHz, 1400 MHz, and 2100 MHz. Of those, typically two (or more) were used in immediate succession within $\sim 1$ hour per observation session, which is possible at AO since the receiver does not need to be physically replaced for a receiver change. At the GBT, observations were made with two receivers at 800 MHz and 1400 MHz. Observations were only included in the dataset if observations could be made within a time span of 14 days; otherwise the observations were discarded for the lack of information about variations in the interstellar medium dispersion.
Two sets of nearly-identical data acquisition systems were developed specifically for NANOGrav. Early observations through 2012.3 (AO) and 2011.0 (GBT) were recorded by the Astronomical Signal Processor (ASP) and the Green Bank Astronomical Signal Processor (GASP) respectively (Demorest 2007). The later observations beginning at 2012.2 (AO) and 2010.2 (GBT) were recorded using the Puerto Rican Ultimate Pulsar Processing Instrument (PUPPI) and the Green Bank Ultimate Pulsar Processing Instrument (GUPPI) respectively (DuPlain et al. 2008; Ford et al. 2010). These instruments performed real-time coherent dedispersion on the digitized incoming baseband data, and folded the data using a pre-computed ephemeris. After RFI excision, polarization calibration, and flux calibration, the end product of each instrument consisted of total-intensity pulse profiles for a series of frequency channels. These profiles were integrated over the course of an observation, resulting in one or more subintervals of typically 20-30 minutes each.
2.2. Time of Arrival data
TOAs were created using various tools in the PSRCHIVE package: the Fourier-domain algorithm of Taylor (1992) to calculate TOAs, and denoising the pulse profiles via wavelet decomposition. Offsets resulting from latencies between different observing systems were determined as overall timing-model-fit-parameters, which were taken into account when doing the timing noise analysis.
For each pulsar, TOAs were calculated for all frequency channels recorded from a given receiver. The effect of time-varying dispersion was taken into account by including “DMX” parameters in the timing model (DFG13), which essentially allows for an extra delay proportional to $1/\nu^2$ to be fit for, where $\nu$ is radio frequency. Pulse shape evolution with frequency is taken into account differently than in DFG13. Instead of including a phase offset per frequency channel, a heuristic mitigation procedure is used that parameterizes the profile evolution with frequency, the details of which are in Arzoumanian et al. (2015).
In summary, similar to the NANOGrav 5-year release, the NANOGrav 9-year release consists of high-quality, publicly available1 TOAs for 37 pulsars, produced per frequency sub-channel.
3. DATA ANALYSIS METHODS
All the data analysis methods we use in this manuscript (both the Bayesian and frequentist methods) are effectively carried out in the time-domain. We start in Sec. 3.1 by defining the likelihood function and introducing our notation,
which follows that of Arzoumanian et al. (2015). We continue our discussion with the Bayesian approach in Sec. 3.2, and the frequentist approach in the form of the optimal cross-correlation statistic in Sec. 3.3.
3.1. Likelihood
We start our discussion by decomposing our $N_{\text{TOA}}$ timing residual data for a single pulsar $\delta\mathbf{t}$ in its individual constituents as follows:
The term $M\boldsymbol{\epsilon}$ describes inaccuracies in the subtraction of the timing model, where $M$ is called the timing model design matrix, and $\boldsymbol{\epsilon}$ is a vector of small timing model parameter offsets. The term $F\mathbf{a}$ describes all low-frequency signals, including low-frequency (“red”) noise, with a limited number of Fourier coefficients $\mathbf{a}$ . Our harmonics are chosen as integer multiples of the harmonic base frequency $1/T$ , with $T$ the length of our dataset (either of a single pulsar, or the entire array of pulsars). The matrix $F$ then has alternating sine and cosine functions. We note that this is just a particular choice of rank-reduced basis, and we could have chosen many others without influencing our results. The term $U\mathbf{j}$ describes noise that is fully correlated (correlation coefficient of 1) for simultaneous observations at different observing frequencies, but fully uncorrelated in time. The matrix $U$ is an $N_{\text{TOA}} \times N_{\text{epoch}}$ matrix that maps the $N_{\text{epoch}}$ observation sessions to the $N_{\text{TOA}}$ TOAs. The vector $\mathbf{j}$ describes the white noise per observation session that is fully correlated across all observing frequencies. The last term, $\mathbf{n}$ , describes Gaussian white noise that is assumed to be left in the data, after correcting for all known systematics.
The white noise is assumed to be an uncorrelated, Gaussian noise process, with variance-covariance:
where $E_k$ and $Q_k$ are the TEMPO and TEMPO2 “EFAC” and “EQUAD” parameters for each observing backend $k$ , $\delta_{ij}$ is the Kronecker delta function, and $W = \text{diag}{\sigma_i^2}$ are the TOA uncertainties. The matrix elements $(i, j)$ apply to the TOAs corresponding to the observing system labeled by $k$ . In practice we cannot fully separate various contributions to our TOAs, so we have to take into account that our corrections for the various terms of Eq. (3) are inaccurate. We do this by constructing our likelihood from the noise-mitigated timing residuals:
where $\mathbf{r}$ is our best approximation of $\mathbf{n}$ , given our knowledge of all the noise and signal parameters. The likelihood can now be written down as:
where $N = \sum_k N_k$ represents the total effect of all white noise. We have collectively denoted all parameters not directly represented by $\boldsymbol{\epsilon}$ , $\mathbf{a}$ , and $\mathbf{j}$ as $\phi$ . Henceforth, we shall refer to $\boldsymbol{\epsilon}$ as the hyperparameters. We group the reduced-rank signals as follows:
which allows us to elegantly place a Gaussian prior on the coefficients of these random processes. The prior covariance is:
resulting in a prior:
where $\infty$ is a diagonal matrix of infinities, which effectively means we have a uniform unconstrained prior on the timing model parameters $\epsilon$ . As described in Arzoumanian et al. (2015), this representation allows us to analytically marginalize Eq. (6) times Eq. (9) over the waveform coefficients $\mathbf{b}$ of the noise, resulting in a drastic reduction of the dimensionality of our posterior (Lentati et al. 2013; van Haasteren & Vallisneri 2014, 2015):
with $C = N + TBT^T$ . The Woodbury matrix identity (Woodbury 1950) can be used to evaluate Eq. (10) efficiently.
The parameters that describe $B$ are the hyperparameters $\phi$ . The hyperparameters of the diagonal matrix $\mathcal{J}$ are the per-backend TEMPO2 ECORR parameters. The matrix $\varphi$ represents the spectrum of the low-frequency noise and the stochastic gravitational waves, and it therefore contains terms correlated between pulsars. Denoting pulsar number with $(a, b)$ , and frequency bin with $(i, j)$ , we can write:
where $\eta_a$ is the spectrum of the low-frequency noise of pulsar $a$ , $\rho$ is the spectrum of the GWB, and $\Gamma_{ab}$ is the signal correlation matrix. The elements of the signal correlation matrix consist of the overlap reduction function for a GWB signal, which is a dimensionless function that quantifies the correlated response of the pulsars to a stochastic GWB (Mingarelli & Sidery 2014). The quantities $\rho_i$ and $\eta_{ai}$ can be either modeled as independent model parameters (i.e., per frequency), or as modeled spectral density with a specific shape (e.g., a power-law model). We note that $\varphi$ can be represented by a block-diagonal matrix, where each block corresponds to a specific frequency bin; all frequencies are theoretically independent degrees of freedom.
3.2. Bayesian analysis
As an alternative to the orthodox frequentist approach to data analysis, Bayesian inference is a method of statistical inference in which Bayes rule of conditional probabilities is used to update one's knowledge as observations are acquired. Given a model $\mathcal{H}$ , model parameters $\Theta$ , and observations $\mathcal{D}$ , we write Bayes rule as:
where $Pr(\Theta|\mathcal{D}, \mathcal{H}) \equiv Pr(\Theta)$ is the posterior (probability distribution) of the parameters, $Pr(\mathcal{D}|\Theta, \mathcal{H}) \equiv L(\Theta)$ is the likelihood, $Pr(\Theta|\mathcal{H}) \equiv \pi(\Theta)$ is the prior (probability distribution), and $Pr(\mathcal{D}|\mathcal{H}) \equiv \mathcal{Z}$ is the marginal likelihood or evidence.
The left-hand side of Eq. (12) can be regarded as the “output” of the Bayesian analysis, and the right-hand side is the
“input”. Indeed, provided we have a generative model of our observations (meaning we can simulate data, given the model parameters), we know the likelihood and prior. However, for parameter estimation we would like to know the posterior, and for model selection we need the evidence.
For parameter estimation, the evidence $\mathcal{Z}$ is usually ignored, and one can use $L(\Theta)\pi(\Theta)$ directly to estimate confidence intervals. Typically one provides confidence intervals for single components and pairs of elements of $\Theta$ . This involves an integral over $Pr(\Theta)$ over all but one or two parameters, a process called marginalization. When $\Theta$ is higher-dimensional, Monte-Carlo sampling methods are typically used to perform this multi-dimensional integral. We use Markov Chain Monte Carlo methods in this work to sample the posterior distribution.
Model selection between two models $\mathcal{H}_0$ and $\mathcal{H}_1$ can be carried by calculating the “Bayes factor”: the ratio between the evidence for the two models. Assuming we have a prior degree of belief of how likely the two model are ( $Pr(\mathcal{H}_0)$ and $Pr(\mathcal{H}_1)$ ), we can write:
where $\mathcal{B}_{10}(\mathcal{D}) \equiv \mathcal{Z}_1/\mathcal{Z}_0$ is the Bayes factor, and $O$ is the odds ratio. The odds ratio can be obtained by calculating the evidence $\mathcal{Z}$ for each model separately (e.g. with Nested Sampling or thermodynamic integration), or by calculating the Bayes factor $\mathcal{B}$ between two models directly (e.g. with trans-dimensional markov chain Monte Carlo methods).
3.3. Optimal cross-correlation statistic
The optimal statistic (Anholm et al. 2009; Demorest et al. 2013; Chamberlin et al. 2015) is a frequentist estimator of the isotropic GWB strain-spectrum amplitude in the weak-signal regime. The estimator maximizes the likelihood of Eq. (10), and it can be written as:
where $\mathbf{P}a = \langle \delta\mathbf{t}a \delta\mathbf{t}a^T \rangle$ is auto-covariance matrix of the residuals of a single pulsar. This is $C$ of Eq. (10) for a noise-only model. The term $A{\text{gw}}^2 \tilde{\mathbf{S}}{ab} = \mathbf{S}{ab} = \langle \delta\mathbf{t}a \delta\mathbf{t}b^T \rangle$ represents the signal covariance between pulsar $a$ and pulsar $b$ . As described in the previous section, our model contains no other signals with a non-zero correlation coefficient between different pulsars. The normalization of the optimal statistic is such that $\langle \hat{A}{\text{gw}}^2 \rangle = A{\text{gw}}^2$ .
Following Chamberlin et al. (2015), we also quote expressions for the variance, and the signal-to-noise ratio (S/N, in power) of the optimal statistic. The variance of the optimal statistic in the absence of a GWB signal is given by:
Although this expression is not valid in general, in the weak-signal regime, in which the cross correlated power is ignored from the variance, this can be used as an approximation of the variance in $\hat{A}_{\text{gw}}^2$ . The S/N for a given signal and noiserealization is:
which has an expectation value of
These expressions generalize the detection significance estimator provided in Jenet et al. (2005), properly taking into account the spectrum of the signal and the noise, as well as details such as the irregular sampling. The S/N here has the same meaning, which if interpreted as a zero-mean unit-variance normal distribution, can be used to place upper-limits on the GWB amplitude. Clearly, this interpretation is only appropriate in the weak-signal regime, but it serves as an independent sanity check for our other methods. Setting our statistical significance threshold $\alpha = 0.95$ , we can place a one-sided upper-limit as:
We note that all the usual caveats of frequentist-type upper-limits apply to this methodology as well, as no prior information is used. For instance, it is possible to set an upper-limit of $A_{\text{ul}} = 0$ in a dataset without a (detectable) signal, which theoretically happens with probability $\alpha$ .
It was shown in Chamberlin et al. (2015) that the optimal statistic is identical to the cross-correlation statistic used by DFG13. This alternative interpretation of the optimal statistic allows us to obtain a measure of the cross power between pulsars. The cross power is the amount of correlated power between the timing residuals of different pulsars, and one would expect this cross power to follow the Hellings and Downs cross-correlation signature for a detectable GWB signal. The cross power and the uncertainty estimates are given by
which is independent of any specific overlap reduction function and $\hat{A}{\text{gw}}^2 \Gamma{ab} \hat{\mathbf{S}}{ab} = \mathbf{S}{ab}$ . One can then fit these correlation coefficients assuming a particular overlap reduction function by minimizing
where $\Gamma_{ab}$ are the cross-correlation coefficients given by the overlap reduction function for an isotropic GWB. It can easily be shown that the best fit value of $A_{\text{gw}}$ is then the optimal statistic value of Eq. (14).
4. RESULTS
4.1. Optimal Statistic Analysis
The optimal cross correlation statistic was applied to the full 9-year NANOGrav dataset. The maximum-likelihood single-pulsar noise values were obtained by independent noise analyses on each pulsar. The maximum-likelihood amplitude
Figure 1. (Top Panel): Cross-correlated power vs. angular separation obtained through an Optimal Statistic Analysis. The dashed red curve shows the maximum-likelihood amplitude mapped onto the Hellings and Downs coefficients. (Bottom Panel): Histogram of the number of pulsars in each bin (red right axis) and a weighted histogram (blue left axis) using the $1\text{-}\sigma$ uncertainties in the cross correlation as weights.
and SNR of this search are $\hat{A}{\text{gw}} = 8.9 \times 10^{-16}$ and $\rho = 0.8$ , respectively, indicating little evidence for the expected GWB cross correlations. The resulting upper limit using this method is $A{\text{gw}}^{95%} = 1.3 \times 10^{-15}$ at a reference frequency of $\text{yr}^{-1}$ , which is 5.4 and 1.5 times more stringent than the limits using this method presented in DFG13 and LTM15, respectively. The corresponding signal-to-noise ratio is 0.8, indicating little evidence for the expected cross correlations from a GWB. It should be noted that the limiting technique presented in Anholm et al. (2009) and DFG13 does not strictly have proper frequentist coverage in the presence of any measurable GWB signal2; therefore, this limit will serve as a proxy to our improved sensitivity and not a strict upper limit.
As a by-product of the optimal statistic analysis, we can obtain the maximum-likelihood cross-correlation values for each pulsar pair in the analysis. In Figure 1 we plot the cross-correlated power vs. angular separation in the top panel. We have binned the values into 10 degree bins using a weighted averaging technique. The red curve shows the maximum-likelihood correlated power fit. It is clear that the cross-correlated power is still consistent with zero signal and that we have much more sensitivity at values of $\zeta < 100^\circ$ , which is a by-product of the fact that our best timed pulsars all lie at a similar position on the sky. This is illustrated in the bottom panel of Figure 1, where we plot the histogram of the number of pulsars in each ten-degree bin (red) as well as the most significant bins (blue). To calculate the most significant bins we use the one-sigma uncertainties on each maximum-likelihood cross-correlation value as weights in the histogram.
2 In the limit of zero GWB signal, this method actually over covers in the frequentist sense. We have found by comparison with other frequentist bounding techniques and Bayesian upper limits that the 5-year limit published in DFG13 is indeed robust and does not suffer to under coverage.#### 4.2. Bayesian Analysis
We have used the Bayesian data analysis techniques described in Section 3 to place upper limits on the GWB parameterized using the standard power law, spectral components and a broken-power-law. Bayesian data analysis of modern PTA data sets is quite difficult to perform in general due to the large parameter spaces ( $\sim 500$ for the current NANOGrav dataset assuming a power law red noise and GWB spectrum) and likelihood evaluation time. Even with the efficient methods described in Section 3 we are still prohibited from carrying out a full Bayesian search where we vary all noise parameters and GW parameters for all pulsars simultaneously.
In this work we ameliorate these two problems by using a two-tiered approach to noise modeling and GWB analysis. To avoid including all noise parameters in our GWB analysis, we first carry out single-pulsar noise analyses on each pulsar in which we include EFAC, EQUAD, and ECORR for each backend/receiver system and red noise parameterized as a power law. We then perform the GWB analysis while holding all EFAC, EQUAD, and ECORR values fixed to their mean value from the single-pulsar analysis but allowing the red noise and GWB parameters to vary. By holding these white noise parameters fixed, we reduce the computational burden since large matrix products can be pre-computed. Furthermore, by holding all white noise values fixed we reduce the number of free parameters drastically from $\sim 500$ to $2N_{\text{psr}} + N_{\text{gw}}$ , where $N_{\text{psr}}$ and $N_{\text{gw}}$ are the number of pulsars in the array and number of parameters describing the GWB, respectively.
We further reduce the computational burden in two ways: by ignoring cross correlations and by only using a subset of pulsars. We choose to ignore cross correlations in this work as they do not contribute to the upper limit in the absence of a signal and only serve to greatly increase our computational burden by requiring that we invert a large dense matrix for every iteration in the analysis. That the cross correlations have no bearing on the upper limit has been argued in Shannon et al. (2013) and further shown in LTM15. A real GWB signal will reveal itself as a strong common red noise term well before the cross correlations become detectable, and since we see no evidence for a common red noise term, we are justified in dropping these terms. Lastly, we choose to only include a subset of pulsars in our analysis as not all pulsars contribute to our upper limit either due to short timing baselines or large measurement errors. To choose this subset of pulsars, we first carry out our single pulsar analyses mentioned above but now include an extra red noise process with power spectral index fixed to $13/3$ (i.e., SMBHBs). The pulsars are then sorted based on their single pulsar upper limits on the GWB. We then carry out the GWB analysis mentioned above to compute an upper limit using an increasing number of pulsars in the sorted list. This process is continued until the upper limit saturates. In other words, including more pulsars beyond this point does not change the upper limit. Using this method we choose to use the 18 pulsars shown in Table 1.
To compute the posterior probability of Eq. (10) and to map out the multi-dimensional parameter space we make use of the pulsar timing data analysis suite PAL23 in conjunction with the parallel-tempering Markov Chain Monte-Carlo (PTMCMC) code4. The details of the PTMCMC sampler can
3 https://github.com/jellis18/PAL2
4 https://github.com/jellis18/PTMCMCSampler
Table 1
Pulsars used in GWB analysis.
| PSR Name | [] | RMSa [] | Timing Baseline [yr] |
|---|---|---|---|
| J1713+0747 | 1.96 | 0.116 | 8.76 |
| J1909–3744 | 4.5 | 0.081 | 9.04 |
| J1640+2224 | 11.8 | 0.158 | 8.9 |
| J1600–3053 | 12.3 | 0.197 | 5.97 |
| J2317+1439 | 13.6 | 0.267 | 8.87 |
| J1918–0642 | 16.0 | 0.34 | 9.01 |
| J1744–1134 | 16.1 | 0.334 | 9.21 |
| J0030+0451 | 19.8 | 0.723 | 8.76 |
| J0613–0200 | 23.4 | 0.592 | 8.58 |
| J1614–2230 | 23.9 | 0.189 | 5.09 |
| B1855+09 | 26.6 | 1.338 | 8.86 |
| J1853+1303 | 31.0 | 0.235 | 5.6 |
| J2145–0750 | 33.0 | 0.37 | 9.07 |
| J1455–3330 | 37.9 | 0.694 | 9.21 |
| J1012+5307 | 43.3 | 1.197 | 9.21 |
| J1741+1351 | 56.8 | 0.103 | 4.24 |
| J2010–1323 | 83.3 | 0.312 | 4.08 |
| J1024–0719 | 92.9 | 0.28 | 4.01 |
aWeighted root-mean-square of epoch-averaged post-fit timing residuals. See Arzoumanian et al. (2015) for more details.
be found in Appendix C of Arzoumanian et al. (2014). The parameters and prior ranges used in the analysis are shown in Table 2. As shown in the table, we have chosen to use uniform priors on both the red noise and GWB amplitude parameters. The GWB amplitude prior is chosen so that it is proper at $A_{\text{gw}} = 0$ , that is to say, it must converge to a finite value in the limit of zero amplitude, otherwise the upper limit would depend on the lower bound of the prior which is undesirable. However, for the red noise amplitude we have no such restriction. A log-uniform prior would result in the most unbiased parameter estimation and a uniform prior would bias the parameter estimation towards higher red noise amplitude and lower red noise spectral indices (since there is a strong covariance between these two parameters); however, the log uniform prior would cause intrinsic red noise to be modeled by the GWB amplitude which is also undesirable. Because of these considerations we choose uniform priors on both as it gives equal weight to both red noise and the GWB. We also find that this prior results in a much more robust upper limit when using different numbers of frequencies in the rank-reduced approximations of Eq. (9).
4.2.1. Power-law and Spectral Limits
When computing upper limits on the dimensionless strain amplitude we fix the spectral index ( $13/3$ in the case of SMBMBs) and adopt a uniform prior on $A_{\text{gw}}$ . When performing a spectral analysis we again use priors that correspond to a uniform prior on $A_{\text{gw}}$ , this results in a prior that is uniform in the square root of the power spectrum coefficients.
Figure 2 show the results of the power law and spectral analysis along with relevant astrophysical model predictions. The solid black and long dashed black lines are the 95% upper limits from the spectral and power-law analyses, respectively. The green, blue, and red shaded regions are the one-sigma prediction on the strain spectra from MOP14, RWS14, and S13, respectively. We find an upper limit on the dimensionless strain amplitude of $A_{\text{gw}}^{95%} = 1.5 \times 10^{-15}$ , a factor of 1.6 better than the most constraining published upper limit to date (Shannon et al. 2013), a factor of 2 more constrainingTable 2
Summary model parameters and prior ranges.
| Parameter | Description | Prior | Comments |
|---|---|---|---|
| White Noise | |||
| EFAC per backend/receiver system | uniform in [0, 10] | Only used in single pulsar analysis | |
| EQUAD per backend/receiver system | uniform in logarithm [-8.5, -5] | Only used in single pulsar analysis | |
| ECORR per backend/receiver system | uniform in logarithm [-8.5, -5] | Only used in single pulsar analysis | |
| Red Noise | |||
| Red noise power law amplitude | uniform in [, ] | 1 parameter per pulsar | |
| Red noise power law spectral index | uniform in [0, 7] | 1 parameter per pulsar | |
| GWB | |||
| GWB power law amplitude | uniform in [, ] | 1 parameter for PTA for power-law models | |
| GWB power law spectral index | delta function | Fixed to different values depending on analysis | |
| GWB power spectrum coefficients at frequency | uniform in [, ]a | 1 parameter per frequency | |
| GWB broken power-law amplitude | log-normalb for models A(B) |
1 parameter for PTA for broken power law models | |
| GWB broken power-law low-frequency spectral index | uniform in [0, 7] | 1 parameter for PTA for broken power law models | |
| GWB broken power-law bend frequency | uniform in logarithm [-9, -7]c | 1 parameter for PTA for broken power law models | |
a The prior uniform in $\rho_i^{1/2}$ is chosen to be consistent with a uniform prior in $A_{\text{gw}}$ for the power law model since $\varphi_i \propto A_{\text{gw}}^2$ .
b These values are quoted in log base 10 and are obtained from MOP14 and S13.
c We choose different prior values on $f_{\text{bend}}$ when mapping to astrophysical model parameters as described in Section 5.
Figure 2. Strain amplitude vs. GW frequency. The solid black and long dashed black lines are the 95% upper limits from our spectral and power-law analyses. The red, blue and green shaded regions are the one-sigma predictions from the models of S13, RWS14, and MOP14. The green shaded region uses the simulation results from MOP14, but replaces the fit to the GWB predictions used in that paper with the functional form given by Eq. (24).
than the recent EPTA upper limit of LTM15, and a factor of 5 more constraining than the 5 year data release upper limit when applying the same Bayesian analysis. Furthermore, we find a slightly less constraining upper limit when using the free spectrum model (power-law equivalent upper limit of $2 \times 10^{-15}$ ). This is to be expected since the free spectrum model has many more degrees of freedom (we use 50 free amplitudes for each of the 50 frequencies in this case) over the power law parameterization (1 degree of freedom). Thus, since the power law model can leverage extra information at all frequencies, as opposed to the spectrum model where each frequency is independent of the others, more constraining upper limits are expected from a power law model. We also find that the upper limit on the strain spectrum from the spectrum analysis is consistent with white noise (i.e., $h_c^{\text{white}}(f) \propto f^{3/2}$ ) at
frequencies $\gtrsim 3/T$ , where $T$ is the length of the longest set of residuals in the data set, which indicates that our GWB upper limits are coming from the three lowest frequency bins. This behavior is to be expected since we have several well timed pulsars that do not span the full 9-year baseline (see Table 1) and thus will have peak sensitivity at frequencies greater than $1/T$ .
From inspection of Figure 2 we see that our 95% upper limit is within at least the 2-sigma confidence region of all three astrophysical models and is sensitive to a potential turnover in the spectrum due to environmental coupling factors. We wish to determine the level of consistency between our data and the power-law models displayed in Figure 2. To accomplish this we follow the method applied in Shannon et al. (2013). Given that we have a model $M$ for the value of the GW amplitude $A_{\text{gw}}$ whose probability distribution function is denoted $p(A_{\text{gw}}|M)$ and that we have a probability distribution function for $A_{\text{gw}}$ given the data, denoted $p(A_{\text{gw}}|d)$ , where $d$ represents the data, the probability that we measure a value of $\hat{A}{\text{gw}}$ greater than that predicted by the model, $A{\text{gw}}^M$ , is given by the law of total probability
Therefore, low values of $P(\hat{A}{\text{gw}} > A{\text{gw}}^M)$ indicate that the range of $A_{\text{gw}}$ that is consistent with our data is inconsistent with the model $M$ , and vice versa. To carry out this procedure the distribution $p(A_{\text{gw}}|d)$ is simply the marginalized posterior distribution when using the uniform prior on $A_{\text{gw}}$ . We use log-normal distributions to model the MOP14, S13, and RWS14, models. Since the models of RWS14, and S13 predict nearly the same GWB amplitude distribution (assuming a power-law only) we make no distinction between these two models. Furthermore, the model distributions on $A_{\text{gw}}$ , given by log-normal distributions have mean and standard deviations of (-14.4, -15) and (0.26, 0.22) for the MOP14 (hereafter model A) and S13/RWS14 (hereafter model B) mod-Figure 3. Upper limit on the GWB as a function of power spectral index.
els, respectively. Using the aforementioned distributions and Eq. (21) we find that our data are 0.8% and 20% consistent with models A and B, respectively, under the assumption of a power-law. This indicates that either the assumptions that go into these models are incorrect, our universe is a realization of the GWB that has an amplitude in the tail of the probability distributions mentioned above, or that environmental effects are depleting SMBHB sources at low frequencies making the power-law assumption faulty. The implications of this last point are discussed in Section 5.
In addition to our power law limits on the stochastic background (i.e., strain spectral index $-2/3$ ), we have also computed the upper limit on the GWB for a range of different spectral indices. In Figure 3 we plot the upper limits obtained at varying spectral indices (red points) vs. power spectral index. We also provide the best fit model for the upper limit as a function of power spectral index where we find $A_{\text{gw}}^{95%} \propto 10^{-0.4\gamma} \propto T^{0.83\alpha}$ . This differs from DFG13 where they find $A_{\text{gw}}^{95%} \propto T^\alpha$ , arguing that this is due to the fact that the sensitivity is dominated by the lowest frequency of $1/T$ . Our fit, giving a slightly weaker dependence on $\alpha$ is consistent with what we have seen above, namely that our limits are not completely dominated by the lowest frequency.
In a Bayesian analysis, the posterior distribution is the prior distribution updated by the data. Here we illustrate this by comparing our power-law upper limits, using identical methods, on the 5-year (DFG13) and 9-year (Arzoumanian et al. 2015) NANOGrav data releases. The results of this comparison are shown in Figure 4 where we plot the marginalized posterior distributions of $\log_{10} A_{\text{gw}}$ for the 9- and 5-year data releases in blue and red, respectively. The gaussian prior distributions described above are shown in green for model A and model B. For the uniform prior case we see quite a dramatic improvement (i.e., the factor of 5 mentioned above) in the upper limits. For model A; the 5-year dataset does somewhat inform the prior, whereas the 9-year data set results in a posterior that is largely inconsistent with the prior distribution. For model B the 5-year data set do not inform the prior at all, whereas the 9-year data set does indeed update the prior.
4.2.2. Broken Power-law Limits
We place constraints on the strength of environmental coupling effects that will likely affect our GWB signal at low frequencies (i.e., large orbital separations) via a simple parameterization of the GWB spectrum that allows for a “bend” frequency at which there is a transition from environmentally-driven evolution to GW-driven evolution. The following dis-
Figure 4. Marginalized posterior probability density of $\log_{10} A_{\text{gw}}$ computed using the nine (blue) and five (red) year NANOGrav data releases for uniform, MOP14 model gaussian, and S13/RWS14 model gaussian prior distributions. The gaussian priors are shown in green.
ussion and analysis techniques are based on Sampson et al. (2015). Here we give a brief overview of this more generalized GWB spectrum.
The characteristic amplitude of a stochastic background from an ensemble of SMBHBs in circular orbits is (Phinney 2001; Sesana et al. 2008; McWilliams et al. 2014)
where $d^3N/(dz d\mathcal{M} dt)$ is the differential number of inspiraling binaries per unit $z$ , $\mathcal{M}$ and $t$ , where $z$ is the redshift, $\mathcal{M} = (m_1 m_2)^{3/5} / (m_1 + m_2)^{1/5}$ is the chirp mass of the binary, and $t$ is the time measured in the binary rest frame. The $dt/d\ln f$ term describes the frequency evolution of the binary system, and $h(f)$ is the strain spectrum emitted by a single circular binary with orbital frequency $f/2$ . Typically, it has been assumed that the binary is purely GW-driven which results in our usual expression for $h_c(f)$ given in Eq. (2); however, physical mechanisms other than GW radiation that are necessary to drive the binary to coalescence (Milosavljević & Merritt 2003) will change the frequency dependence (i.e., the $dt/d\ln f$ term) of this equation (see Colpi 2014, for a review of SMBHB coalescence). Following Sampson et al. (2015) we can generalize the frequency dependence of the strain spectrum to
where $i$ ranges over many physical processes that are driving the binary to coalescence. If we restrict this sum to GW-driven evolution and an unspecified physical process then the strain spectrum is now
where $f_{\text{bend}}$ and $\kappa$ are the parameters that encode information about the physical processes (other than GW radiation) driving the binary evolution. As mentioned above, there could be many physical processes contributing to the frequency evolution of the SMBHB system; however, at current sensitivity it is very unlikely that our data can distinguish them. Thus we follow Sampson et al. (2015) and adopt this slightly simplified spectrum to place constraints on possible environmental coupling mechanisms.
The above discussion has focused on the frequency evolution of SMBHB. The other piece of the equation is the merger rate of SMBHs, which will set the overall amplitude scale for the GWB somewhat independently (assuming the last parsec problem is solved) of the physical mechanisms that drive the system to merger. Here we use the same log-normal distributions introduced in Section 4.2.1 for models A and B as prior probability distributions for the GWB amplitude $A$ in Eq. (24). In the following analysis we also fix $\alpha = -2/3$ and use uniform priors on $\log_{10} f_{\text{bend}} \in [-9, -7]$ and $\kappa \in [0, 7]$ unless stated otherwise.
Here, we have run an identical analysis to that of Section 4.2.1 except that the GWB spectrum model is now that of Equation (24) and we adopt the aforementioned priors on $A$ , $f_{\text{bend}}$ , and $\kappa$ . Figures 5 and 6 show the results of this analysis. Figure 5 shows the posterior probability density of the GWB spectrum defined in Equation (24) with model A on the left and model B on the right. This probability density is constructed by drawing values of $A$ , $f_{\text{bend}}$ , and $\kappa$ from the joint probability distribution shown in Figure 6, constructing the spectrum at each frequency via Equation (24) and then producing a histogram of the spectral power at each frequency. The solid black lines in Figure 5 represent the 95% credible region and the median of the GWB spectrum. The dashed line is the upper limit on $A_{\text{gw}}$ using the purely GW-driven spectrum (i.e., no transition frequency) and the gaussian amplitude priors from models A and B, respectively. Lastly the thin solid black line is the 95% upper limit on the GWB spectrum from the spectral analysis presented in Section 4.2.1. By inspecting the inferred GWB spectrum one can determine that the data prefer a GWB spectrum that has a definitive transition from GW-driven to environmentally-driven within the pulsar timing frequency band for Model A, whereas the data does not significantly constrain the shape of the spectrum for model B.
This can be seen further in the joint posterior distributions in Figure 6 in which the probability distributions (blue) for the bend frequency parameter, $f_{\text{bend}}$ , and spectral index parameter, $\kappa$ , are significantly different from the prior distribution (green) for model A and not significantly different for model B. When this same analysis is carried out on the 5-year data release we find that the data can only slightly inform the prior on $f_{\text{bend}}$ for model A and gives no information on the other parameters in either model A or B. This, again, indicates that the 9-year data release provides us with much more information about the shape of the GWB strain spectrum.
Finally, we can be more quantitative and apply Bayesian model selection to this problem by computing the Bayes factor between the broken-power-law (Equation (24)) and pure-power-law (Equation (2)) parameterizations for both models A and B. When this analysis is carried out we arrive at Bayes factors of $22.2 \pm 1.1$ and $2.23 \pm 0.15$ for models A and B, respectively. These Bayes factors were computed using parallel tempering and a custom thermodynamic integration implementation (See Sec. 6.1.2 of Cornish & Littenberg 2015).
This analysis shows, for the first time, that PTAs are en-
tering a regime where even in the case of a non-detection meaningful constraints can be placed on the dynamical history of the SMBHB population. Furthermore, this analysis shows that when placing upper limits to make statements about the full range of astrophysical merger scenarios it is no longer valid to consider only the classic strain amplitude, but one must instead frame the question in terms of measuring the amplitude and shape of the GWB spectrum. As we have seen in the above analysis and as can be seen clearly in Figure 6, given a model for the SMBHB merger physics (i.e., a prior on $A$ ) and discarding the assumption of a purely GW-driven signal (i.e., a broken-power-law model), it is very difficult to rule out any of the GWB amplitude parameter space with any certainty unless one has strong a-priori knowledge on the shape of the spectrum. However, we can begin to place constraints on the environmental coupling effects that drive the system to the GW-dominated regime.
5. DISCUSSION
5.1. Astrophysical Model Limits
While the parameter estimation of the previous section is interesting in and of itself, some of the most interesting science available from the NANOGrav data is accessible only by relating these parameters to properties of the source populations. Here we attempt to interpret the phenomenological posteriors on the GW spectrum in terms of black hole-host galaxy relations, environmental effects, or binary-eccentricity by carrying out our analyses in sequence, investigating different effects separately. In Section 5.1.1 we use the results of our power-law analyses of Section 4.2.1 to provide constraints on the parameters of scaling relations between host galaxies and black holes. We then go beyond the assumption of a power-law spectrum in Section 5.1.2 to investigate how our constraints on the shape of the characteristic strain spectrum from Section 4.2.2 map to constraints on the environment of SMBH binaries. Finally in Section 5.1.3, we probe the eccentricity of binaries before they entered the PTA band.
5.1.1. Constraints on Host Galaxy – Black Hole Scaling Relations
If the gravitational wave spectrum is assumed to be created by an ensemble of binary SMBHs that are formed following galaxy mergers (spectral index of $-2/3$ ), and whose evolution is assumed to be dominated by GW emission throughout the PTA-sensitive waveband, then we can trace the expected binary SMBH population using observations of galaxy merger rates, the galaxy stellar mass function, and the black hole-host galaxy relation. This is the approach taken in S13, and Ravi et al. (2015). Assuming equal fractional uncertainties in these parameters, the black hole-host galaxy relation will have the largest impact on the predicted level of the GW background. This is due to the much stronger dependence of the GW background on the chirp mass of each source than on the number of sources.
Middleton et al. (2015) shows that it is difficult to extract information from PTA limits without making significant assumptions about the shape of the black hole merger rate density. If instead a galaxy merger rate density calculated from observed galaxy parameters is used as a proxy for the black hole merger rate density, then a limit on the GW background can be directly translated into a limit on the input galaxy parameter spaces (Simon & Burke-Spolaor 2015). For this paper, we focus specifically on the scaling relation between host galaxy properties and black hole mass (e.g. Häring & RixFigure 5. Probability density plots of the recovered GWB spectra for models A and B using the broken-power-law model parameterized by ( $A_{\text{gw}}$ , $f_{\text{bend}}$ , and $\kappa$ ) as discussed in the text. The thick black lines indicate the 95% credible region and median of the GWB spectrum. The dashed line shows the 95% upper limit on the amplitude of purely GW-driven spectrum using the Gaussian priors on the amplitude from models A and B, respectively. The thin black curve shows the 95% upper limit on the GWB spectrum from the spectral analysis.
Figure 6. One- and two-dimensional posterior probability density plots of the spectrum model parameters $A_{\text{gw}}$ , $f_{\text{bend}}$ , and $\kappa$ . In the one-dimensional plots, we show the posterior probability from the 9-year data set (blue), the 5-year dataset (dashed red) and the prior distribution used in both analyses (green). In the two dimensional plots we show a heat map along with the one (solid), two (dashed), and three (dash-dotted) sigma credible regions. model A is on the left and model B is on the right.
(2004), McConnell & Ma (2013)) as it is the observed parameter that is most easily constrained by NANOGrav data. Specifically, we constrain the $M_{\bullet} - M_{\text{bulge}}$ relation:
In addition to $\alpha$ and $\beta$ , observational measurements of this relation also fit for $\epsilon$ , the intrinsic scatter of individual galaxy measurements around the common $\alpha$ , $\beta$ trend line. In practice, $\alpha$ and $\epsilon$ have the greatest impact on predictions of $A_{\text{gw}}$ , and all observational measurements agree with $\beta \approx 1$ .
PTAs are most sensitive to binary SMBHs where both black holes are $\gtrsim 10^8 M_{\odot}$ (e.g. Sesana et al. (2008)). Therefore $M_{\bullet} - M_{\text{bulge}}$ relations that are derived including the most massive systems are the most relevant to understanding the population in the PTA band. Several recent measurements of the $M_{\bullet} - M_{\text{bulge}}$ relation specifically include high-galaxy-mass measurements, e.g. those from Brightest Cluster Galaxies (BCGs). As these fits include the high-mass black holes that we expect to dominate the PTA signals, we take these as the “gold standard” for comparison with PTA limits (KormendyFigure 7. The above plot shows the translation of the marginalized posterior distribution of $A_{\text{gw}}$ , Fig. 4, into the black hole-host galaxy parameter space, which is characterized by an intercept $\alpha$ , a slope $\beta$ , and an intrinsic scatter $\epsilon$ . Flat priors are used for $\alpha$ , $\beta$ , and $\epsilon$ . $\beta$ is not informed by the distribution of $A_{\text{gw}}$ , while both $\alpha$ and $\epsilon$ are, with a limit on $\alpha$ being more strongly set. The curves show the 1, 2, and $3\sigma$ contours. Relevant observational measurements are also shown, with McConnell & Ma (2013) in blue and Kormendy & Ho (2013) in magenta. Since $\beta$ is not strongly informed by the upper limit, we can set an upper limit in $\alpha$ - $\epsilon$ space by marginalizing over $\beta$ . That upper limit is shown in Fig. 8.
& Ho 2013; McConnell & Ma 2013, e.g.).
The translation of an upper limit on $A_{\text{gw}}$ to the black hole-host galaxy parameter space is calculated as follows:
where the posterior of $A_{\text{gw}}$ , $p(A_{\text{gw}}(\alpha, \beta, \epsilon, \theta) | \text{PTA})$ , is the marginalized posterior distribution of $A_{\text{gw}}$ , which is shown in Fig. 4; $A_{\text{gw}}(\alpha, \beta, \epsilon, \theta)$ is the prediction of $A_{\text{gw}}$ calculated from models similar to S13; $\theta$ represents the galaxy stellar mass function and the galaxy merger rate; and $p(\alpha, \beta, \epsilon | \text{PTA})$ is the marginalized posterior distribution of the black hole-host galaxy relation, which is shown in Fig. 7. For this analysis, we use two leading measurements of the galaxy stellar mass function, Ilbert et al. (2013) and Tomczak et al. (2014), and two measurements of the galaxy merger rate, Robotham et al. (2014) and Keenan et al. (2014), as the basis for simulating a local population of binary SMBHs. A flat prior is used for $\alpha$ , $\beta$ , and $\epsilon$ , and the posterior on $A_{\text{gw}}$ using a uniform prior, as seen in Fig. 4, is directly translated into this parameter space. The result of which is shown in Fig. 7. $\beta$ is clearly not informed by a PTA posterior, but the combination of $\alpha$ and $\epsilon$ are, with the strongest limit being set on $\alpha$ .
Fig. 8 shows the translation of our posterior on $A_{\text{gw}}$ into $\alpha$ - $\epsilon$ parameter space with observational measurements of the parameters from Kormendy & Ho (2013) and McConnell & Ma (2013). Assuming a power-law analysis of the S13 model, as described in Simon & Burke-Spolaor (2015), there is a slight inconsistency between our upper limit and the Kormendy & Ho (2013) measurement. Potential solutions to an inconsistency include: a significant number of black hole binaries do
Figure 8. The above plot shows the translation of the 95% upper limit on $A_{\text{gw}}$ , Fig. 4, into the parameter space $\alpha$ - $\epsilon$ , which characterizes the black hole-host galaxy relation as described in Equation. (25). The parameter space above the line is inconsistent with the power-law analysis of the S13 model, as described in Simon & Burke-Spolaor (2015). Observational measurements of this parameter space are shown with errorbars.
Figure 9. By introducing the parameter $T_{\text{stall}}$ , as described in Simon & Burke-Spolaor (2015), we can start to explore the inconsistency of our upper limit with power-law models for the GW background. In the above plot, we allow $T_{\text{stall}}$ to vary while using the $M_{\bullet} - M_{\text{bulge}}$ relation Kormendy & Ho (2013). The probability of $T_{\text{stall}}$ is a direct translation of the posterior on $A_{\text{gw}}$ from Fig. 4. The blue line is the 95% lower limit on $T_{\text{stall}}$ , which we set at 0.73 Gyr. While this is not sufficiently constraining to make meaningful astrophysical statements, this parameter may be useful for future PTA upper limits.
not reach the GW-dominant regime in our assumed timescale (they ‘stall’); the ‘classical’ assumption of a power-law strain spectrum in the PTA band is incorrect and in fact there is a turn-over in the strain spectrum at lower frequencies (see Sec. 4.2.2); or that the measured astronomical parameters are not correct for the population of binary SMBHs in the PTA band.
As the possibility for a different strain spectrum curve is discussed in Sec. 4.2.2, let us explore the potential for ‘stalling’ within the model described so far in this section. Using the galaxy merger rate density as a proxy for the black hole merger rate density implies as assumption that the events occur at a similar cosmological time. If there was significant stalling in the binary black hole population, then these events would be offset in cosmological time by some ‘stalling timescale’. There is then nothing in the model to allow for anything other than an efficient solution to the ‘last-parsec’ problem Begelman et al. (1980). As done in Simon & Burke-Spolaor (2015), we introduce a variable, $T_{\text{stall}}$ , which is a mea-sure of this offset in time between the assumed galaxy merger rate density and the black hole merger rate density used to model a GW background. Fig. 9 shows the translation of the posterior distribution of $A_{\text{gw}}$ , Fig. 4, into a probability distribution of $T_{\text{stall}}$ using the Kormendy & Ho (2013) measurements of the $M_{\bullet} - M_{\text{bulge}}$ relation. Using this we set a 95% lower limit on $T_{\text{stall}}$ of 0.73 Gyr, which is not sufficiently constraining to indicate which of the many suggested ‘solutions’ to the ‘last-parsec’ problem are not occurring for the systems in the PTA band. However, this parameter may be useful for future PTA upper limits as a probe of the level of ‘stalling’ in the binary black hole population.
5.1.2. Constraints on binary environmental influences
The cores of galactic merger remnants can harbor stars with little angular momentum and almost radial trajectories which intersect the central galactic region (centrophilic orbits). These stars can undergo three-body interactions with the resident supermassive black-hole binary, causing the stars to be ejected, which results in energy and angular momentum being extracted from the black hole system, and leading to binary hardening (Quinlan 1996).5 Additionally, the formation of circumbinary gaseous disks can lead to interactions which extract energy and angular momentum from the binary orbit, driving it towards smaller orbital separations (Ivanov et al. 1999). We expect that, in the type of gas-poor galaxies which dominate the nanohertz GW signal, hardening from stellar scattering will dominate over circumbinary interactions, but we consider both in the following. We begin with a discussion of how these environmental mechanisms drive the evolution of the binary, then provide constraints on the frequency at which the characteristic strain spectrum exhibits a turnover from the analysis in Sec. 4.2.2. We finish by mapping these frequencies to constraints on the astrophysical environment of SMBH binaries emitting GWs in the nanohertz band.
Environmentally-driven orbital evolution — We use the formalism of Quinlan (1996) to define the evolution of a (circular) binary due to three-body stellar scattering events, where
where $a$ is the binary orbital semi-major axis, $\rho$ is the mass density of galactic core stars, $H$ is a dimensionless hardening rate which takes a value of $\sim 15$ , and $\sigma$ is the velocity dispersion of core stars. Using Kepler’s third law, we can rearrange this equation to solve for the rate of orbital frequency evolution, which gives $df/dt \propto f^{1/3}$ . Since the binary orbital evolution will be due to a combination of environmental influences and GW emission, $df/dt$ is actually a sum over all mechanisms, as in Eq. (23). We know that the rate of frequency evolution due to GW emission is $df/dt \propto f^{11/3}$ (Peters & Mathews 1963), hence, in the language of the parametrized spectrum model in this paper given in Eq. (24), $\kappa = 10/3$ for binary hardening by three-body stellar scattering.
Likewise, the evolution of a circular binary due to circumbinary disk interaction is modeled within the $\alpha$ -disk formalism (Ivanov et al. 1999; Haiman et al. 2009; Sesana 2013)
5 We assume that all galactic merger remnants maintain the same mass density of core stars throughout the binary merger. The subtleties of loss-cone replenishing impact the evolution of the binary and of the central density profile within a factor of $\sim 2$ (a few at most), as shown by Sesana & Khan (2015) and Vasiliev et al. (2015).
as (Sesana 2013)
where $\dot{M}_1$ is the accretion rate onto the primary black hole, $\mu$ is the binary reduced mass, and $a_0$ is a characteristic orbital separation at which the enclosed disk mass equals the mass of the secondary black hole. The latter can be expressed as (Ivanov et al. 1999)
where $\alpha$ is a disk viscosity parameter, $M_{1,2}$ are the binary black hole masses; $\dot{M}_E = 4\pi GM_1/c\kappa_T$ is the Eddington accretion rate of the primary ( $\kappa_T$ is the Thompson opacity coefficient); and $r_g = 2GM_1/c^2$ is the Schwarzschild radius of the primary.
As in the stellar scattering case, we can rearrange this equation to determine the orbital frequency evolution. In this model, $df/dt \propto f^{4/3}$ , and so $\kappa = 7/3$ for $\alpha$ -disk binary interactions.
Constraints on spectral turnover frequency — We define the spectral turnover frequency in the obvious way to mean the frequency at which the characteristic strain spectrum exhibits a change in slope. If the low-frequency slope is positive, this will correspond to the point at which the spectrum is maximized. One must note that setting $f = f_{\text{bend}}$ in Eq. (24) does not maximize the characteristic strain spectrum. Rather the turnover frequency will be a function of both $f_{\text{bend}}$ and $\kappa$ :
We can combine our measurements of $f_{\text{bend}}$ and $\kappa$ from the analysis in Sec. 4.2.2 to compute the probability distribution of spectral turnover frequencies. We find that placing numerical constraints on $f_{\text{bend}}$ is difficult as the posterior is heavily influenced by the prior, namely the upper and lower bounds for the uniform priors used in this analysis. In the following analysis we set the lower bound on $f_{\text{bend}}$ by requiring that the power at $f = 1/T_{\text{span}}$ differs by no more than 10% from a pure-power-law for any value of $\kappa$ . By using this prior, we ensure that we can recover a pure-power law spectrum (in our frequency range) for any value of $\kappa$ . The upper bound of $f_{\text{bend}}$ is chosen based on the specific environmental coupling mechanism we are considering.
We emphasize that the probability distributions of $f_{\text{bend}}$ , $\kappa$ , and $f_{\text{turn}}$ are not distributions of these parameters over the SMBH binary population. Rather our posterior distribution is illustrating the spread of our beliefs in the measurement of the single $f_{\text{bend}}$ and $\kappa$ model parameters.
Constraints on environmental parameters — We can now extract astrophysical constraints from our constraints on the transition and spectral turnover frequencies. By equating the rate of orbital evolution, $da/dt$ , due to environmental mechanisms and GW emission, we can deduce the characteristic transition frequency, $f_{\text{bend}}$ , between these influences. We firstly considerstellar-scattering, for which the transition frequency is given by
where $M$ is the total binary mass, $M_8 \equiv M/(10^8 M_\odot)$ , $q = M_2/M_1$ , $q_r = q/(1+q)^2$ , $\sigma_{200} \equiv \sigma/(200 \text{ km/s})$ , $\rho_3 \equiv \rho/(10^3 M_\odot \text{ pc}^{-3})$ , and $H_{15} \equiv H/15$ . In the second line we have used the $M-\sigma$ relation (McConnell et al. 2011) to replace velocity dispersion, where $M_8 \approx 1.9 \sigma_{200}^5$ . It follows from the McConnell et al. (2011) $M-\sigma$ relation that the velocity dispersion term on the first line of Eq. (31) has the mass scaling $\sigma_{200}^{-3/10} \propto M_8^{-3/50}$ , which is a small modification to the exponent of the other mass term, $M_8^{-2/5}$ . Hence, $f_{\text{bend}}$ has a relatively weak dependence on the mass scaling of $\sigma$ . Nevertheless we can express the transition frequency in terms of variables of a parametrized $M-\sigma$ relation, where $\log_{10}(M/M_\odot) = a + b \log_{10} \sigma_{200}$ , such that $M_8 = 10^{a-8} \sigma_{200}^b$ , and
Finally, the weak scaling of $f_{\text{bend}}$ with $H$ , and the $\lesssim 20%$ deviations of this parameter away from 15 seen in numerical scattering experiments (Sesana & Khan 2015) justifies our keeping this parameter fixed at its fiducial value of $H_{15} = 1$ . Astrophysical estimates on $\rho$ are quite uncertain with estimated values around $10-10^4 M_\odot \text{ pc}^{-3}$ for typical environments (Dotti et al. 2007). The variation of estimates over several orders of magnitude is why we choose to investigate only $\rho$ in our stellar-scattering constraints.
The equivalent transition frequency for $\alpha$ -disk interaction is
where $a_0$ is as previously defined. We adopt the fiducial value of the disk viscosity parameter $\alpha \sim 10^{-2}$ used in Ivanov et al. (1999). The very weak dependence of the bend frequency on this parameter, $f_{\text{bend}} \propto \alpha^{6/49}$ will significantly dampen the influence of any deviations from this fixed value. Hence, in our constraints on the influence of circumbinary disk interactions, we only vary the accretion rate of gas onto the primary black-hole, $\dot{M}_1$ , of which estimates in the literature vary over several orders of magnitude – typically $10^{-3} M_\odot \text{ yr}^{-1} - 1 M_\odot \text{ yr}^{-1}$ , see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; Goicovic et al. 2015.
Equations (31)-(34) indicate the GW frequency at which a single circular binary will transition between being environmentally driven and being GW driven. Of course, $f_{\text{bend}}$ is not the quantity that can be extracted from a spectral analysis – $f_{\text{turn}}$ is what we can measure. Our analysis of the stochastic GW background has provided us with constraints on the characteristic transition frequency of the entire population, and thus the turnover of the spectrum. Hence, our path to providing constraints on environmental parameters requires us to numerically construct characteristic strain spectra for populations of SMBH binaries in contact with their environment. We can then construct numerical mappings between the environmental parameters of interest (core stellar mass density, $\rho$ , and primary black hole accretion rate, $\dot{M}_1$ ) and the turnover of the spectrum. We use the formalism of Phinney (2001) via
Eq. (22), where the differential comoving number density of merging binaries per redshift and component masses is constructed as in McWilliams et al. (2014):
where $P(z)$ encapsulates redshift dependent factors, and $\phi(M_{{1,2}})$ is the number density of black holes of a certain mass, which is given by a (redshift dependent) Schechter function modified at the high-mass end by a lognormal component to accommodate recent high mass BCG (brightest cluster galaxy) discoveries (Lin et al. 2010).
The combined influence on the binary orbital evolution of GW emission and environmental couplings is modeled with the sum in Eq. (23), where either stellar scattering or disk interactions are included i.e. we consider one mechanism at a time. By considering all binary environments to have the same $\rho$ or $\dot{M}_1$ , we can determine the spectral turnovers from the numerically computed strain spectra, iterating over many values of these environmental parameters to deduce a mapping. We draw binary systems from the ranges $z \in [0, 1]$ , $M_1 \in [10^7, 10^{10}] M_\odot$ , and $q \in [0.1, 1]$ . The results for our fiducial assumptions are shown in the top panel of Fig. 10, with the stellar density required to give a certain turnover frequency shown in the left panel, and the primary accretion rate required to give a certain turnover frequency shown in the right panel.
In the lower panels of Fig. 10 we plot the posterior distributions of the stellar density, $\rho$ , for stellar hardening, and mass accretion rate, $\dot{M}1$ for circumbinary disk interaction. In this analysis we perform the Bayesian parameter estimation for fixed values of $\kappa$ corresponding to the appropriate values for stellar hardening ( $\kappa = 10/3$ ) and circumbinary disk interaction ( $\kappa = 7/3$ ). These posteriors are constructed by first converting the marginalized distributions on $f{\text{bend}}$ to a distribution for $f_{\text{turn}}$ via Eq. (30) and then using the empirical mapping to convert $f_{\text{turn}}$ to the appropriate astrophysical parameter. Again, we do not place numerical confidence limits on $\rho$ or $\dot{M}_1$ since the data does not constrain the prior distribution at large values. Nonetheless, from inspection of Fig. 10 we see that the MOP14 model heavily prefers $\rho \gtrsim 10^4 M_\odot \text{ pc}^{-3}$ and $\dot{M}_1 \gtrsim 10^{-1} M_\odot \text{ yr}^{-1}$ , while the S13 model is largely unconstraining for both mechanisms. Typical densities of massive elliptical galaxies at the MBH influence radius is $\sim 10^3 M_\odot \text{ pc}^{-3}$ , making the MOP14 model hard to reconcile with observations, even if we consider that massive ellipticals were likely factor 2–3 more compact at $z = 1$ (Sesana, unpublished). Our results approach the upper end (for the MOP14 prior) of the expected range of $\dot{M}$ , $10^{-3} M_\odot \text{ yr}^{-1} - 1 M_\odot \text{ yr}^{-1}$ , observed in the local universe and predicted via simulations, see e.g. Di Matteo et al. 2001; Armitage & Natarajan 2002; Goicovic et al. 2015. Furthermore, Dotti et al. (2015) predict that $\dot{M}_1 \ll 10^{-1} M_\odot \text{ yr}^{-1}$ for BH masses of $10^9 M_\odot$ and redshifts $z < 1$ ; however, these are average accretion rates, and short, episodic accretion triggered by galaxy mergers could occur at a higher rate.
We go beyond the fiducial assumptions for the case of stellar hardening since it is the most likely environmental coupling mechanism for SMBH binaries. When we increase the low-mass cutoff of systems which contribute to the characteristic strain budget this further constrains the stellar mass density. This is seen most easily in Eqn. (31), where one mustFigure 10. (top): Empirical mapping from $f_{\text{turn}}$ to $\rho$ (left) and $\dot{M}1$ (right). (bottom): Posterior distributions for the mass density of stars in the galactic core (left) and the accretion rate of the primary black hole from a circumbinary disk (right). These distributions are constructed by first converting the marginalized distribution of $f{\text{bend}}$ to a distribution of $f_{\text{turn}}$ via Eq. (30), and then using the empirical mapping described in the text to convert from $f_{\text{turn}}$ to the astrophysical quantities $\rho$ and $\dot{M}_1$ , respectively.
raise the stellar mass density to match a corresponding increase in binary mass so that the transition frequency is maintained. Furthermore, modeling the distribution of black holes masses in Eq. (35) without the lognormal component or redshift evolution will increase the contribution of lower mass binaries to the GW strain budget, leading to smaller stellar mass density constraints than reported in Fig. 10. Varying the normalization, $a$ , and exponent, $b$ , of the $M-\sigma$ relation such that $a \in [7, 9]$ and $b \in [4, 6]$ has very little impact on the environmental constraints.
5.1.3. Constraints on SMBH binary population eccentricity
It is not only the astrophysical environment of SMBH binaries that can induce a bend in the characteristic strain spectrum. Binaries with non-zero eccentricity emit GWs at a spectrum of harmonics of the orbital frequency. The cumulative effect over the entire population can lead to a depletion of the low frequency strain spectrum (Enoki et al. 2007; Sesana 2013; Ravi et al. 2014; Huerta et al. 2015), and a turnover whose shape can be captured with the parametrized spectrum model employed in this paper. Hence, we can use our $f_{\text{turn}}$ posterior from the marginalization of $f_{\text{bend}}$ over all $\kappa$ to deduce constraints on the eccentricity of binaries at some reference orbital separation. Our approach follows from the previous astrophysics constraints, where we build populations and strain spectra which have varying eccentricities at a fixed semi-major axis of 0.01 pc, then construct a relationship between this eccentricity and the spectral turnover. An important modeling assumption we make here is that binaries are (and always have been) driven entirely by GW emission. This factors into how we model $df_r/dt_r$ and how we evolve the binary eccentricity, where we assume binaries could have eccentricities arbitrarily close to 1 in the past.
We construct eccentric populations and corresponding strain spectra using the formalism of Huerta et al. (2015). The resulting relationship between the spectral turnover frequency and the eccentricity of all binaries at a semi-major axis of 0.01 pc is shown in the top panel of Fig. 11, along
Figure 11. Same as Figure 10 except now we display the empirical mapping (top) and posterior distribution (bottom) for the eccentricity of SMBH binaries when they had a semi-major axis of 0.01 pc.
with the corresponding eccentricity posteriors from each amplitude prior in the bottom panel. The high turnover frequency obtained with the MOP14 prior leads to an eccentricity posterior distribution that largely favors $e_0 \gtrsim 0.7$ while the S13 prior leads to an eccentricity posterior that is consistent with smaller eccentricities, more weakly favoring $e_0 \gtrsim 0.5$ . We emphasize that, whilst these eccentricities seem rather large, it is well established that binaries evolving in stellar environments tend to increase their eccentricity (Quinlan 1996). It is therefore likely that most binaries can get to $e \sim 0.5-0.7$ along their evolution (see tracks in Sesana 2010). The eccentricity growth rate is generally larger for smaller binary mass ratios, and for larger initial eccentricities. The latter is indeeda crucial parameter; if, following galaxy mergers, the MBHB already has a significant ( $e \gtrsim 0.5$ ) eccentricity at the moment of formation, the subsequent evolution will almost certainly drive it to $e > 0.9$ . Given that the MBHB eccentricity at formation is hard to determine (Aarseth 2003; Hemsendorf et al. 2002; Berentzen et al. 2009; Amaro-Seoane et al. 2009), it is impossible to draw astrophysical conclusions from the constraints above. Furthermore, in reality the history of a binary's eccentricity will see phases of growth and circularization depending upon the interplay of environmental factors and GW emission (e.g., Sesana 2010; Kocsis & Sesana 2011). This should be considered in future analyses.
5.2. Cosmological Model Limits
SMBH binary mergers are not the only possible source of a GWB in the pulsar timing band. In this subsection, we use our power-law spectrum results, shown in Fig. 3, to constrain the fractional energy density of the Universe in relic GWs, along with stringent limits on the tension of cosmic strings.
5.2.1. Relic GWs
Quantum fluctuations of the gravitational field in the early universe, amplified by an inflationary phase, are expected to produce a stochastic relic GW background, see e.g. Grishchuk (1976, 1977); Starobinsky (1980); Linde (1982). Observations of this background would provide a unique insight into the highly curved spacetime of the early universe at less than $10^{-32}$ s after the Big Bang and at energy scales of $10^{16}$ GeV, when quantum mechanics and general relativity should reconcile, BICEP2/Keck et al. (2015); Ade et al. (2014). This background is expected to produce a characteristic signature in the polarization of the Cosmic Microwave Background (CMB) radiation, as well as CMB temperature anisotropies Grishchuk (2005). In the context of PTAs, the amplitude of the relic GW background is of astrophysical and cosmological interest due to the fact that it intrinsically depends on the equation of state of the early universe, $w$ , and thus the Hubble constant in the inflationary stage $H_*$ , as well as the tensor-to-scalar ratio $r$ , see e.g. Zhao (2011); Zhao et al. (2013).
Specifically, the spectral index of the stochastic GWB, $\alpha$ , is related to the equation of state of the early universe, $w$ , via
where $n_t$ is the spectral index of the primordial power spectrum, usually set to 0. In current hot big bang models, $w = 1/3$ and $n_t = 0$ , thus $\alpha = -1$ . We therefore fix the spectral index to -1 and apply the Bayesian analysis methods described in Sec 3.2 to the 9-year NANOGrav dataset. Using analysis methods described in 4.2, we obtain a 95% upper limit of on the amplitude of the relic GW background of
assuming a power spectrum for the characteristic strain with $\alpha = -1$ at a reference frequency of $f = \text{yr}^{-1}$ . This then constrains the GW energy density content per unit logarithmic frequency, divided by the critical energy density, $\rho_c$ , to close the Universe:
Table 3
Values for $\Omega_{\text{gw}}(f)h^2$ reported at 1/yr. Cosmological parameters used for $H_*$ are $h = 0.702$ , $T_{\text{CMB}} = 0.276$ K, $\Omega_\Lambda = 0.725$ , $\Omega_m = 0.275$ , and $z_{\text{eq}} = 3454$ . Values of $H_*$ are given in multiples of the Planck mass $m_{\text{pl}} \equiv 1/\sqrt{G}$ .
| EoS, | Spectral index, | |||
|---|---|---|---|---|
| 0 | -2 | 5.4 | ||
| 1/3 | -1 | |||
| 1 | -1/2 |
where $f$ is the frequency, $\rho_c = 8\pi/(3H_0^2)$ is the critical energy density required to close the Universe, $H_0$ is the Hubble expansion rate, and $\rho_{\text{gw}}$ is the total energy density in GWs (Allen & Romano 1999; Maggiore 2000). The NANOGrav limit is therefore
in a radiation-dominated universe with equation of state of $w = 1/3$ . This new limit is a factor of 2.9 better than the previous best limit of $\Omega_{\text{gw}}(f)h^2 = 1.2 \times 10^{-9}$ from LTM15.
Although a radiation-dominated era is usually assumed to follow inflation in the hot big bang paradigm, there is currently no evidence to show this held before Big Bang Nucleosynthesis (BBN). In fact, the existence of a reheating stage or the existence of a cosmic phase transition both violate this assumption (Grishchuk 2001; Watanabe & Komatsu 2006; Zhao 2011). For completeness, we now allow for other equations of state of the early universe before the BBN stage. The same analysis can be repeated assuming different equations of state for the early universe: a matter-dominated universe would have $w = 0$ and therefore by Eq. (36), $\alpha = -2$ . If the universe were instead dominated by the kinetic energy of the inflaton, then $w = 1$ and $\alpha = -1/2$ .
Finally, we place limits on the Hubble parameter, $H_*$ , in the inflationary stage using methods developed in Zhao (2011). There, they introduced a way of translating the upper limit on a primordial GW background to a constraint on $H_*$ . Using typical cosmological parameters $h = 0.702$ , $T_{\text{CMB}} = 0.276$ K, $\Omega_\Lambda = 0.725$ , $\Omega_m = 0.275$ , and $z_{\text{eq}} = 3454$ , they obtain the following relation:
where $k_* = 2\pi f_*$ and is reported at a reference frequency of $f_* = 1/\text{yr}$ and $m_{\text{pl}} \equiv 1/\sqrt{G}$ is the Planck mass. Using Eq. (40), we can then limit $H_*$ for a fixed equation of state. For example, using the limit on $\Omega_{\text{gw}}h^2 = 4.2 \times 10^{-10}$ for $w = 1/3$ , see Table 3, one can place a limit $H_* = 1.6 \times 10^{-2} m_{\text{pl}}$ . Results for $w = 0, 1$ are in Table 3.
In LTM15, the limit on a primordial GW background with $w = 1/3$ and $n_t = 0$ is $\Omega_{\text{gw}}(f)h^2 = 1.2 \times 10^{-9}$ , resulting in $H_* = 2.74 \times 10^{-2} m_{\text{pl}}$ . One can see that the NANOGrav limit on $H_*$ is an improvement of 1.7 over the EPTA, and the most stringent limit to date.
5.2.2. GW Background from Cosmic (Super)strings
Cosmic strings are linear topological defects that can be produced in the early universe via phase transitions (Kibble 1976; Vilenkin 1981; Vilenkin & Shellard 2000). Cosmic superstrings are fundamental strings stretched to cosmological scales and arise in a wide range of string-theory-inspired cosmological scenarios (Sarangi & Tye 2002; JonesFigure 12. (left): Cosmic string constraints in terms of string tension $G\mu$ and reconnection probability $p$ using the results of recent cosmic string simulations described in Blanco-Pillado et al. (2014). (right): Cosmic string tension $G\mu$ vs loop size parameterized by $\alpha_{cs}$ using the model described in LTM15. The shaded area is ruled out by our GW upper limit in both panels.
et al. 2003; Copeland et al. 2004). Cosmic (super)strings produce a stochastic background of GWs as well as individual bursts (Damour & Vilenkin 2001; Damour & Vilenkin 2005; Siemens et al. 2006, 2007; Ölmez et al. 2010).
Our limits on the amplitude of the stochastic background can also be used to constrain the parameter space of cosmic (super)strings. Recent simulations (Blanco-Pillado et al. 2014) have shown that cosmic (super)string loop densities are dominated by loops that formed at scales comparable to the Hubble size at the time of formation, even though only about 10% of loops are formed with such large sizes. We use the loop distributions derived by Blanco-Pillado et al. (2014), specifically Eqs. (63), (65), and (67) of that reference with loop size $\alpha_{cs} = 0.05$ , together with the techniques described in Ölmez et al. (2010) to compute the stochastic background produced by cosmic string cusps. The cosmological parameters we used are taken from the Planck 2015 data (Ade et al. 2015). In this case the relevant parameters are the string tension $G\mu$ and the reconnection probability $p$ . We explore this parameter space and exclude regions where the cosmic (super)string network would have resulted in a stochastic background amplitude larger than that ruled out by our measurements. The left panel of Figure 12 shows the results of our analysis. On the y-axis we show the reconnection probability and on the x-axis the string tension. The gray shaded area shows the region of cosmic string parameter space that is ruled out. Note that for $p = 1$ our data only allow for cosmic (super)strings with tensions $G\mu < 1.3 \times 10^{-10}$ .
Recently LTM15 presented a comprehensive and fully general overview of cosmic string limits from the EPTA, and found a conservative limit on the string tension to be $G\mu < 1.3 \times 10^{-7}$ . The limit is conservative in the sense that it is found by considering a wide range of loop sizes and taking the upper limit to be the largest possible value of $G\mu$ consistent with the data. The limit was identical to that set by the Planck Collaboration, combining Planck and high- $l$ Cosmic Microwave Background data with Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT), cf. Planck Collaboration et al. (2014). While the calculation in LTM15 was not carried out explicitly for the Blanco-Pillado et al. (2014) simulations we can use their published limit on $\Omega_{gw}(f)h^2 = 1.2 \times 10^{-9}$ for cosmic strings to place a limit of $G\mu < 8.6 \times 10^{-10}$ . Our limit for this model is therefore
roughly a factor of 6.6 times more constraining than the inferred previous limit. Using the same analysis developed by the EPTA, (Sanidas et al. 2012, 2013), we compute the upper limit on the string tension $G\mu$ as a function of loop size $\alpha_{cs}$ as shown in the right panel of Figure 12. Our conservative limit on cosmic string tension using this range of cosmic string models is $G\mu < 3.3 \times 10^{-8}$ , a factor of 4 better than both the combined Planck, ACT, SPT limit and the EPTA limit.
6. SUMMARY AND CONCLUSIONS
This paper reports on the search for an isotropic stochastic GW background in NANOGrav’s 9-year dataset. We do not find positive statistical evidence for the presence of such a signal. Following up on a series of earlier results by the three PTAs, we report new upper limits on the amplitude of backgrounds described by power-law spectra:
• For an astrophysical background of SMBH binaries (corresponding to a timing-residual spectral density with exponent $\gamma = 13/3$ ), we find a 95% confidence limit $A_{gw} < 1.5 \times 10^{-15}$ , five times more constraining than the analogous limit for NANOGrav’s 5-year dataset (DFG13). Under the assumption of purely GW-driven evolution, leading to an unbroken $\gamma = 13/3$ power law, we compute the probability that our constraint is consistent with the MOP14 and S13/RWS14 theoretical predictions for $A_{gw}$ as 0.8% and 20%, respectively, essentially ruling out the MOP14 model and placing the S13/RWS14 model in tension with our data. [Sec. 4.2.1.]
• We verify the consistency of our limit with previously reported scaling relations between SMBH mass and galactic bulge mass, adopting fiducial estimates for galaxy merger rates and the stellar mass function. Under the assumption of circular GW-driven binaries, our limit is slightly inconsistent with the Kormendy & Ho (2013) relation, and consistent within the error margin for the McConnell & Ma (2013) relation. [Sec. 5.1.1.]
• We also perform an optimal-statistic (cross-correlation) analysis, and find limits that are 5.4 and 1.5 times more constraining than the analogous DFG13 and LTM15 results. The cross-correlation SNR is 0.8, indicating thatthere is little evidence for inter-pulsar residual correlations induced by GWs. [Sec. 4.1.]
• We extend the power-law background search to generic spectral indices, and place the most stringent limits so far on the energy density of relic GWs, $\Omega_{\text{gw}}(f = 1/\text{yr})h^2 < 4.2 \times 10^{-10}$ , for a $w = 1/3$ early-Universe equation of state. From this we obtain limits on the Hubble parameter during inflation, $H_* = 1.6 \times 10^{-2} m_{\text{pl}}$ . [Secs. 4.2.1, 5.2.1]
• We place the most stringent limits to date on a GW background generated by a network of cosmic strings, $\Omega_{\text{gw}}(f = 1/\text{yr})h^2 < 2.2 \times 10^{-10}$ , which translates into a conservative upper limit on cosmic string tension $G\mu < 3.3 \times 10^{-8}$ , using the model presented in LTM15. This is a factor of 4 better than both the combined Planck, ACT, SPT limit and the EPTA limit. Using the recent models of Blanco-Pillado et al. (2014) we find $G\mu < 1.3 \times 10^{-10}$ , a factor of 6.6 times more constraining than an identical analysis using the EPTA limit. [Secs. 4.2.1, and 5.2.2.]
We further probe the interface between PTA observations and SMBH-binary population estimates by analyzing the 9-year dataset in terms of a GW background described by an inflected power law [Eq. (24)]. We derive joint posteriors for the spectral parameters (the amplitude $A_{\text{gw}}$ , the inflection frequency $f_{\text{bend}}$ , and the shape exponent $\gamma$ ) assuming $A_{\text{gw}}$ priors from MOP14 and S13/RWS14 (see Fig. 5). For both priors (but more so for MOP14), the inflected power-law model is preferred to an unbroken power law. The $f_{\text{bend}}$ posterior can be used to infer astrophysical information about the effects that may reduce GW power at lower frequencies, such as the initial eccentricity of SMBH binaries, and the environmental influences of stars and gas in galactic nuclei. To summarize:
- • We find that the data prefer an inflected spectrum over a pure power law, with Bayes factors $\sim 22$ and 2.2 for the MOP14 and S13 amplitude priors, respectively. The same analysis, run on the 5-year DFG13 dataset, provides little to no information about the shape of the GWB spectrum. [Sec. 4.2.1.]
- • Under several simplifying assumptions, we map the posterior distribution of $f_{\text{bend}}$ into posterior distributions for the nuclear stellar mass density $\rho$ (which modulates the strength of binary frequency evolution by stellar scattering, cf. Sec. 5.1.2), the SMBH mass accretion rate from circumbinary disks $\dot{M}_1$ (which can be linked to binary frequency evolution by interactions with gas, cf. Sec. 5.1.2), and the initial value of SMBH-binary orbital eccentricity $e_0$ (which distributes GW power to higher frequency harmonics, cf. Sec. 5.1.3).
For the last decade (and longer), the three PTA consortia have been engaged in an accelerating race toward higher GW sensitivities—the analysis presented in this paper represents the latest stage of the race, but not the last. While our investigation cannot claim the ultimate prize, a positive detection, it is the first to use information about the amplitude and shape of GW background to make concrete (if limited) astrophysical statements about the dynamics and environments of SMBH binaries. The era of GW astrophysics is truly upon us.
Author contributions. An alphabetical-order author list was used for this paper in recognition of the fact that a large, decade-timescale project such as NANOGrav is necessarily the result of the work of many people. All authors contributed to the activities of the NANOGrav collaboration leading to the work presented here, and reviewed the manuscript, text, and figures prior to the paper’s submission. Additional specific contributions to this paper are as follows. ZA, KC, PBD, TD, RDF, EF, MEG, GJ, MJ, MTL, LL, MAM, DJN, TTP, SMR, IHS, KS, JKS, and WWZ made observations for this project and developed timing models. JAE coordinated the writing of the paper. JAE, RvH, and CMFM led this search by directly running the analysis pipelines. JAE, RvH, LS, NJC, SRT, and MV designed and tuned the statistical analysis. SRT, CMFM, STM, AS, JS, SB-S, XS, SAS, LS, MV, NJC, and JAE developed the interpretation of astrophysical results. CMFM developed and interpreted the relic GW results. XS, SAS, CMFM developed and interpreted the cosmic strings results. JAE, SRT, CMFM, JS, SB-S, MV, RvH, XS, STM, and AS wrote the paper, collected the bibliography, prepared figures and tables.
Acknowledgments. The work of ZA, AB, SB-S, SJC, SC, BC, NJC, JMC, PBD, TD, JAE, NG-D, FJ, GJ, MTL, TJWL, LL, ANL, DRL, JL, RSL, DRM, MAM, STM, DJN, NP, TTP, SMR, XS, DRS, KS, JS, MV, RvH and YW was partially supported through the National Science Foundation (NSF) PIRE program award number 0968296. NANOGrav research at UBC is supported by an NSERC Discovery Grant and Discovery Accelerator Supplement and by the Canadian Institute for Advanced Research. D.R.M. acknowledges partial support through the New York Space Grant Consortium. JAE and RvH acknowledge support by NASA through Einstein Fellowship grants PF4-150120 and PF3-140116, respectively. MV acknowledges support from the JPL RTD program. Portions of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. CMFM was supported by a Marie Curie International Outgoing Fellowship within the 7th European Community Framework Programme. SRT was supported by appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities through a contract with NASA. SAS acknowledges funding from an NWO Vidi fellowship (PI: JTW Hessels). AS is supported by a University Research Fellowship of the Royal Society. This work was supported in part by National Science Foundation Grant No. PHYS-1066293 and by the hospitality of the Aspen Center for Physics. This research was performed in part using the Zwicky computer cluster at Caltech supported by NSF under MRI-R2 award no. PHY-0960291 and by the Sherman Fairchild Foundation. A majority of the computational work was performed on the Nemo cluster at UWM supported by NSF grant No. 0923409. Parts of the analysis in this work were carried out on the Nimrod cluster made available by S.M.R. Data for this project were collected using the facilities of the National Radio Astronomy Observatory and the Arecibo Observatory. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is operated by SRI International under a cooperative agreement with the NSF (AST-1100968), and in alliance with Ana G. Méndez-Universidad Metropolitana and the Universities Space Research Association.REFERENCES
Aarseth, S. J. 2003, Ap&SS, 285, 367
Ade, P. A. R., Aikin, R. W., Barkats, D., Benton, S. J., Bischoff, C. A., Bock, J. J., Brevik, J. A., Buder, I., Bullock, E., Dowell, C. D., Duband, L., Filippini, J. P., Fliescher, S., Golwala, S. R., Halpern, M., Hasselfield, M., Hildebrandt, S. R., Hilton, G. C., Hristov, V. V., Irwin, K. D., Karkare, K. S., Kaufman, J. P., Keating, B. G., Kernasovskiy, S. A., Kovac, J. M., Kuo, C. L., Leitch, E. M., Lueker, M., Mason, P., Netterfield, C. B., Nguyen, H. T., O'Brien, R., Ogburn, R. W., Orlando, A., Pryke, C., Reintsema, C. D., Richter, S., Schwarz, R., Sheehy, C. D., Staniszewski, Z. K., Sudiwala, R. V., Teply, G. P., Tolan, J. E., Turner, A. D., Vieregg, A. G., Wong, C. L., & Yoon, K. W. 2014, Phys. Rev. Lett., 112, 241101
Ade, P. A. R. et al. 2015
Allen, B. & Romano, J. D. 1999, 59, 102001
Amaro-Seoane, P., Miller, M. C., & Freitag, M. 2009, ApJ, 692, L50
Anholm, M., Ballmer, S., Creighton, J. D. E., Price, L. R., & Siemens, X. 2009, Phys. Rev. D, 79, 084030
Armitage, P. J. & Natarajan, P. 2002, ApJ, 567, L9
Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., Chamberlin, S. J., Chatterjee, S., Cordes, J. M., Demorest, P. B., Deng, X., Dolch, T., Ellis, J. A., Ferdman, R. D., Garver-Daniels, N., Jenet, F., Jones, G., Kaspi, V. M., Koop, M., Lam, M. T., Lazio, T. J. W., Lommen, A. N., Lorimer, D. R., Luo, J., Lynch, R. S., Madison, D. R., McLaughlin, M. A., McWilliams, S. T., Nice, D. J., Palliyaguru, N., Pennucci, T. T., Ransom, S. M., Sesana, A., Siemens, X., Stairs, I. H., Stinebring, D. R., Stovall, K., Swiggum, J., Vallisneri, M., van Haasteren, R., Wang, Y., Zhu, W. W., & NANOGrav Collaboration. 2015, submitted to ApJ
—. 2014, ApJ, 794, 141
Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
Berentzen, I., Preto, M., Berczik, P., Merritt, D., & Spurzem, R. 2009, ApJ, 695, 455
BICEP2/Keck, Planck Collaborations, :, Ade, P. A. R., Aghanim, N., & et al. 2015, ArXiv e-prints
Blanco-Pillado, J. J., Olum, K. D., & Shlaer, B. 2014, Phys. Rev. D, D89, 023512
Chamberlin, S. J., Creighton, J. D. E., Siemens, X., Demorest, P., Ellis, J., Price, L. R., & Romano, J. D. 2015, Phys. Rev. D, 91, 044048
Colpi, M. 2014, Space Sci. Rev., 183, 189
Copeland, E. J., Myers, R. C., & Polchinski, J. 2004, JHEP, 06, 013
Cordes, J. M. 2013, Classical and Quantum Gravity, 30, 224002
Cordes, J. M. & Shannon, R. M. 2010
Cornish, N. J. & Littenberg, T. B. 2015, Classical and Quantum Gravity, 32, 135012
Damour, T. & Vilenkin, A. 2001, 64, 064008
Damour, T. & Vilenkin, A. 2005, Phys. Rev. D, D71, 063510
Demorest, P. B. 2007, PhD thesis, University of California, Berkeley
Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., Nice, D., Ransom, S., Stairs, I. H., Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., Chamberlin, S. J., Cordes, J. M., Ellis, J., Finn, L. S., Freire, P., Giampanis, S., Jenet, F., Kaspi, V. M., Lazio, T. J., Lommen, A. N., McLaughlin, M., Palliyaguru, N., Perrodin, D., Shannon, R. M., Siemens, X., Stinebring, D., Swiggum, J., & Zhu, W. W. 2013, 762, 94
Detweiler, S. 1979, 234, 1100
Di Matteo, T., Carilli, C. L., & Fabian, A. C. 2001, ApJ, 547, 731
Dotti, M., Colpi, M., Haardt, F., & Mayer, L. 2007, MNRAS, 379, 956
Dotti, M., Merloni, A., & Montuori, C. 2015, MNRAS, 448, 3603
DuPlain, R., Ransom, S., Demorest, P., Brandt, P., Ford, J., & Shelton, A. L. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7019, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
Enoki, M., Enoki, M., Nagashima, M., & Nagashima, M. 2007, Progress of Theoretical Physics, 117, 241
Ford, J. M., Demorest, P., & Ransom, S. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7740, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 0
Foster, R. S. & Backer, D. C. 1990, 361, 300
Gair, J., Romano, J. D., Taylor, S., & Mingarelli, C. M. F. 2014, Phys. Rev. D, 90, 082001
Goicovic, F. G., Cuadra, J., Sesana, A., Stasyszyn, F., Amaro-Seoane, P., & Tanaka, T. L. 2015, ArXiv e-prints
Graham, A. W. & Scott, N. 2013, ApJ, 764, 151
Grishchuk, L. 2001, in Lecture Notes in Physics, Vol. 562, Gyros, Clocks, Interferometers...: Testing Relativistic Gravity in Space, ed. C. Löffler, J. Zimmerzahl, C. Everitt, & F. Hehl (Springer Berlin Heidelberg), 167–192
Grishchuk, L. P. 1976, Soviet Journal of Experimental and Theoretical Physics Letters, 23, 293
Grishchuk, L. P. 1977, in Annals of the New York Academy of Sciences, Vol. 302, Eighth Texas Symposium on Relativistic Astrophysics, ed. M. D. Papagiannis, 439
—. 2005, Physics Uspekhi, 48, 1235
Haiman, Z., Kocsis, B., & Menou, K. 2009, The Astrophysical Journal, 700, 1952
Häring, N. & Rix, H.-W. 2004, ApJ, 604, L89
Hellings, R. W. & Downs, G. S. 1983, 265, L39
Hemsendorf, M., Sigurdsson, S., & Spurzem, R. 2002, ApJ, 581, 1256
Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins, R. A. 1968, 217, 709
Hobbs, G. 2013, Classical and Quantum Gravity, 30, 224007
Hobbs, G., Archibald, A., Arzoumanian, Z., Backer, D., Bailes, M., Bhat, N. D. R., Burgay, M., Burke-Spolaor, S., Champion, D., Cognard, I., Coles, W., Cordes, J., Demorest, P., Desvignes, G., Ferdman, R. D., Finn, L., Freire, P., Gonzalez, M., Hessels, J., Hotan, A., Janssen, G., Jenet, F., Jessner, A., Jordan, C., Kaspi, V., Kramer, M., Kondratiev, V., Lazio, J., Lazaridis, K., Lee, K. J., Levin, Y., Lommen, A., Lorimer, D., Lynch, R., Lyne, A., Manchester, R., McLaughlin, M., Nice, D., Oslowski, S., Pilia, M., Possenti, A., Purver, M., Ransom, S., Reynolds, J., Sanidas, S., Sarkissian, J., Sesana, A., Shannon, R., Siemens, X., Stairs, I., Stappers, B., Stinebring, D., Theureau, G., van Haasteren, R., van Straten, W., Verbiest, J. P. W., Yardley, D. R. B., & You, X. P. 2010, Classical and Quantum Gravity, 27, 084013
Huerta, E. A., McWilliams, S. T., Gair, J. R., & Taylor, S. R. 2015, ArXiv e-prints
Ilbert, O., McCracken, H. J., Le Fèvre, O., Capak, P., Dunlop, J., Karim, A., Renzini, M. A., Caputi, K., Boissier, S., Arnouts, S., Aussel, H., Comparat, J., Guo, Q., Hudelot, P., Kartaltepe, J., Kneib, J. P., Krogager, J. K., Le Floch, E., Lilly, S., Mellier, Y., Milvang-Jensen, B., Moutard, T., Onodera, M., Richard, J., Salvato, M., Sanders, D. B., Scoville, N., Silverman, J. D., Taniguchi, Y., Tasca, L., Thomas, R., Toft, S., Tresse, L., Vergani, D., Wolk, M., & Zirm, A. 2013, A&A, 556, A55
Ivanov, P. B., Papaloizou, J. C. B., & Polnarev, A. G. 1999, Monthly Notices of the Royal Astronomical Society, 307, 79
Jaffe, A. H. & Backer, D. C. 2003, 583, 616
Jenet, F. A., Hobbs, G. B., Lee, K. J., & Manchester, R. N. 2005, 625, L123
Jones, N. T., Stoica, H., & Tye, S.-H. H. 2003, Physics Letters B, 563, 6
Keenan, R. C., Foucaud, S., De Propis, R., Hsieh, B. C., Lin, L., Chou, R. C. Y., Huang, S., Lin, J. H., & Chang, K. H. 2014, ApJ, 795, 157
Kibble, T. W. B. 1976, Journal of Physics A Mathematical General, 9, 1387
Kocsis, B. & Sesana, A. 2011, Monthly Notices of the Royal Astronomical Society, 411, 1467
Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 111
Kramer, M. & Champion, D. J. 2013, Classical and Quantum Gravity, 30, 224009
Lentati, L., Alexander, P., Hobson, M. P., Taylor, S., Gair, J., Balan, S. T., & van Haasteren, R. 2013, 87, 104021
Lentati, L., Taylor, S. R., Mingarelli, C. M. F., Sesana, A., Sanidas, S. A., Vecchio, A., Caballero, R. N., Lee, K. J., van Haasteren, R., Babak, S., Bassa, C. G., Brem, P., Burgay, M., Champion, D. J., Cognard, I., Desvignes, G., Gair, J. R., Guillemot, L., Hessels, J. W. T., Janssen, G. H., Karuppusamy, R., Kramer, M., Lassus, A., Lazarus, P., Liu, K., Oslowski, S., Perrodin, D., Petiteau, A., Possenti, A., Purver, M. B., Rosado, P. A., Smits, R., Stappers, B., Theureau, G., Tiburzi, C., & Verbiest, J. P. W. 2015, ArXiv e-prints
Lin, Y.-T., Ostriker, J. P., & Miller, C. J. 2010, The Astrophysical Journal, 715, 1486
Linde, A. D. 1982, Physics Letters B, 108, 389
Lorimer, D. & Kramer, M. 2005, Cambridge Observing Handbooks for Research Astronomers, Vol. 4, Handbook of Pulsar Astronomy, 1st edn. (Cambridge, U.K.; New York, U.S.A.: Cambridge University Press)
Maggiore, M. 2000, Phys. Rep., 331, 283
McConnell, N. J. & Ma, C.-P. 2013, 764, 184
McConnell, N. J., Ma, C.-P., Gebhardt, K., Wright, S. A., Murphy, J. D., Lauer, T. R., Graham, J. R., & Richstone, D. O. 2011, Nature, 480, 215
McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008
McWilliams, S. T., Ostriker, J. P., & Pretorius, F. 2014, ApJ, 789, 156
Middleton, H., Del Pozzo, W., Farr, W. M., Sesana, A., & Vecchio, A. 2015, ArXiv e-prints
Milosavljević, M. & Merritt, D. 2003, ApJ, 596, 860
Mingarelli, C. M. F. & Sidery, T. 2014, Phys. Rev. D, 90, 062011
Mingarelli, C. M. F., Sidery, T., Mandel, I., & Vecchio, A. 2013, 88, 062005
Olmez, S., Mandic, V., & Siemens, X. 2010, Phys. Rev. D, D81, 104028
Ölmez, S., Mandic, V., & Siemens, X. 2010, Phys. Rev. D, 81, 104028Peters, P. C. & Mathews, J. 1963, Physical Review, 131, 435
Phinney, E. S. 2001, ArXiv Astrophysics e-prints
Planck Collaboration, Ade, P. A. R., Aghanim, N., Armitage-Caplan, C., Arnaud, M., & et al. 2014, A&A, 571, A25
Quinlan, G. D. 1996, New Astronomy, 1, 35
Rajagopal, M. & Romani, R. W. 1995, ApJ, 446, 543
Ravi, V., Wyithe, J. S. B., Shannon, R. M., & Hobbs, G. 2015, MNRAS, 447, 2772
Ravi, V., Wyithe, J. S. B., Shannon, R. M., Hobbs, G., & Manchester, R. N. 2014, MNRAS, 442, 56
Robotham, A. S. G., Driver, S. P., Davies, L. J. M., Hopkins, A. M., Baldry, I. K., Agius, N. K., Bauer, A. E., Bland-Hawthorn, J., Brough, S., Brown, M. J. I., Cluver, M., De Propriis, R., Drinkwater, M. J., Holwerda, B. W., Kelvin, L. S., Lara-Lopez, M. A., Liske, J., López-Sánchez, Á. R., Loveday, J., Mahajan, S., McNaught-Roberts, T., Moffett, A., Norberg, P., Obreschkow, D., Owers, M. S., Penny, S. J., Pimblet, K., Prescott, M., Taylor, E. N., van Kampen, E., & Wilkins, S. M. 2014, MNRAS, 444, 3986
Rosado, P. A., Sesana, A., & Gair, J. 2015, MNRAS, 451, 2417
Sampson, L., Cornish, N. J., & McWilliams, S. T. 2015, Phys. Rev. D, 91, 084055
Sanidas, S. A., Battye, R. A., & Stappers, B. W. 2012, Phys. Rev. D, 85, 122003
—. 2013, ApJ, 764, 108
Sarangi, S. & Tye, S.-H. H. 2002, Physics Letters B, 536, 185
Sazhin, M. V. 1978, 22, 36
Sesana, A. 2010, 719, 851
Sesana, A. 2013, Classical and Quantum Gravity, 30, 224014
Sesana, A. 2013, MNRAS, 433, L1
Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623
Sesana, A. & Khan, F. M. 2015, ArXiv e-prints
Sesana, A., Vecchio, A., & Colacino, C. N. 2008, 390, 192
Shannon, R. M. & Cordes, J. M. 2010, 725, 1607
Shannon, R. M., Ravi, V., Coles, W. A., Hobbs, G., Keith, M. J., Manchester, R. N., Wyithe, J. S. B., Bailes, M., Bhat, N. D. R., Burke-Spolaor, S., Khoo, J., Levin, Y., Oslowski, S., Sarkissian, J. M., van Straten, W., Verbiest, J. P. W., & Want, J.-B. 2013, Science, 342, 334
Siemens, X., Creighton, J., Maor, I., Ray Majumder, S., Cannon, K., & Read, J. 2006, Phys. Rev., D73, 105001
Siemens, X., Mandic, V., & Creighton, J. 2007, Phys. Rev. Lett., 98, 111101
Simon, J. & Burke-Spolaor, S. 2015, in preparation
Starobinsky, A. A. 1980, Physics Letters B, 91, 99
Stinebring, D. 2013, Classical and Quantum Gravity, 30, 224006
Taylor, J. H. 1992, 341, 117
Taylor, J. H. & Weisberg, J. M. 1989, ApJ, 345, 434
Taylor, S. R. & Gair, J. R. 2013, 88, 084001
Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., Labbé, I., Straatman, C. M. S., Papovich, C., Glazebrook, K., Allen, R., Brammer, G. B., Kacprzak, G. G., Kawinwanichakij, L., Kelson, D. D., McCarthy, P. J., Mehrtens, N., Monson, A. J., Persson, S. E., Spitler, L. R., Tilvi, V., & van Dokkum, P. 2014, ApJ, 783, 85
van Haasteren, R., Levin, Y., Janssen, G. H., Lazaridis, K., Kramer, M., Stappers, B. W., Desvignes, G., Purver, M. B., Lyne, A. G., Ferdman, R. D., Jessner, A., Cognard, I., Theureau, G., D’Amico, N., Possenti, A., Burgay, M., Corongiu, A., Hessels, J. W. T., Smits, R., & Verbiest, J. P. W. 2011, 414, 3117
van Haasteren, R. & Vallisneri, M. 2014, Phys. Rev. D, 90, 104012
—. 2015, MNRAS, 446, 1170
Vasiliev, E., Antonini, F., & Merritt, D. 2015, ArXiv e-prints
Vilenkin, A. 1981, Physics Letters B, 107, 47
Vilenkin, A. & Shellard, E. P. S. 2000, Cosmic Strings and Other Topological Defects (Cambridge University Press)
Watanabe, Y. & Komatsu, E. 2006, Phys. Rev. D, 73, 123515
Woodbury, M. A. 1950, Memorandum report, 42, 106
Wyithe, J. S. B. & Loeb, A. 2003, 590, 691
Zhao, W. 2011, Phys. Rev. D, 83, 104021
Zhao, W., Zhang, Y., You, X.-P., & Zhu, Z.-H. 2013, Phys. Rev. D, 87, 124012
Xet Storage Details
- Size:
- 135 kB
- Xet hash:
- d22b6a458258b93ed6e89ab1a36d823a040c14f4936bac208e9f4fd81c809e8b
Xet efficiently stores files, intelligently splitting them into unique chunks and accelerating uploads and downloads. More info.