The NANOGrav 11-year Data Set: Pulsar-timing Constraints On The Stochastic Gravitational-wave Background

I. Introduction

The three major collaborations involved in this effort are the North American Nanohertz Observatory for Gravitational-waves (NANOGrav, McLaughlin 2013), the European Pulsar Timing Array (EPTA, Desvignes et al. 2016), and the Parkes Pulsar Timing Array (PPTA, Hobbs 2013). Additionally, the International Pulsar Timing Array (IPTA, Verbiest et al. 2016) exists as an umbrella consortium for data-sharing, coordinated timing campaigns, and joint GW analysis. The increasing sensitivity of PTAs is apparent in the ever-tightening upper limits (van Haasteren et al. 2011; Demorest et al. 2013; Shannon et al. 2013; Lentati et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2016) on the stochastic GWB from the unresolved superposition of SMBHB signals out to redshift $\lesssim 1$.

The road toward detection lies not only through the accumulation of ever longer and more accurate time-of-arrival (TOA) data for larger arrays of monitored pulsars, but also through the development of powerful, robust, and reliable data-analysis methods to demonstrate the presence of GWs in PTA data. In this article, we report substantial advances along both avenues. First, we present our stochastic-GW analysis of NANOGrav’s largest and most sensitive dataset so far, spanning $45$ pulsars and $11.4$ years. See Sec. II and Arzoumanian et al. 2018 for more about this “11-year” dataset. Second, we describe our statistical-inference framework, which was significantly augmented compared to our GW study of the 9-year dataset (Arzoumanian et al. 2016, hereafter 7). Improvements include a practical strategy to isolate the expected signature of stochastic GWs in our data—namely the emergence of a long-timescale noise process that is common to all pulsars, and the positive detection of inter-pulsar spatial correlations with a quadrupolar signature (Hellings & Downs 1983). This strategy is based on Bayesian model selection, and is extensible to large arrays and datasets. Indeed, for the first time with a large pulsar array, we are able to report GW upper limits and GW-vs-noise (“detection”) Bayes factors computed with likelihoods that include spatial correlations – such as the ones predicted by Hellings & Downs 1983 – a goal that had previously proved computationally unfeasible beyond small arrays (Lentati et al. 2015).

This article also features a more robust, Bayesian–frequentist hybrid “optimal-statistic” analysis (Anholm et al. 2009; Demorest et al. 2013; Chamberlin et al. 2015), which complements our primary Bayesian approach. Additionally, we employ a more flexible end-to-end approach for PTA GW searches to constrain astrophysical parameters (characterizing SMBHB populations and environments, as well as cosmic-string properties). This approach uses a set of GW-spectrum simulations that span the parameter-space region of interest, and interpolates them by means of Gaussian processes (GPs) (Williams & Rasmussen 2006; Taylor et al. 2017b), resulting in a flexible new model that is calibrated directly by detailed simulations.

Lastly, but perhaps most importantly, we report on how Solar System ephemeris (SSE) errors can manifest as a false GWB signal in PTA data, for sufficiently long and high-quality datasets. The SSE is used to refer TOA measurements to an inertial frame located at the Solar System barycenter (SSB). Previous GW searches treated the ephemeris as a fixed-parameter model without uncertainties. However, in the course of analyzing the 11-year dataset we discovered that adopting different ephemerides (among the last few published by the Jet Propulsion Laboratory (JPL); see Folkner et al. 2009; Folkner et al. 2014; Folkner et al. 2016; Folkner & Park 2016) leads to significantly different upper-limit and model-comparison statistics. As PTA datasets become larger, longer, and more precise, our GW searches will continue to uncover systematic effects that will limit our sensitivity unless handled appropriately. To this end, we have developed a physical model of ephemeris uncertainties, and we demonstrate that it makes our analysis insensitive to the choice among recent ephemerides.

This paper is laid out as follows: methodological advances are discussed in Sec. III. In Secs. IV and V we report GW upper limits and detection Bayes factors based on the 11-year dataset, as well as new constraints on astrophysical and cosmological sources of low-frequency GWs. In Sec. VI we present our conclusions and discuss prospects for future observations.

For the busy reader, the following summarizes the most consequential results:

Once we take ephemeris uncertainty into account, we find Bayesian model comparison to be inconclusive on the presence of a GWB-like signal in the data (with signal-vs.-noise and spatial-correlation Bayes factors both $\sim 1$). Adopting one of the fixed JPL ephemerides leads to signal-vs.-noise Bayes factors as high as $26\pm 2$ in favor of a GWB-like signal (for JPL ephemeris DE430), suggesting that systematic ephemeris errors can masquerade as GWs—and conversely that modeling these errors can subtract power from a putative GWB signal. This degeneracy will be resolved over the next few years as we collect longer and larger datasets, and as ephemeris accuracy improves with data from current NASA missions.

Using a model of cosmic-string–generated GW spectra that interpolates among extensive string-network simulations (Blanco-Pillado & Olum 2017), we place a 95% upper limit of $5.3(2)\times 10^{-11}$ on the string tension $G\mu/c^{2}$ for a reconnection probability $p=1$. This result is marginalized over ephemeris uncertainties, but neglects inter-pulsar spatial correlations. (Including these is still too taxing computationally; however we argue that our upper limits assuming a variable power-law exponent, described in Sec. IV.1, are affected modestly by correlations, and so should be the cosmic-string result.) Previous studies reported limits of $1.3\times 10^{-10}$ (7) and $8.6\times 10^{-10}$ (Lentati et al. 2015), although different prior assumptions and the lack of ephemeris modeling preclude a direct comparison.

II. The 1111-year Data Set

Our analyses throughout this paper make use of the NANOGrav $11$-year dataset, which consists of the TOA data and pulsar timing models recently presented in 8, and is publicly available online data.nanograv.org. This dataset is derived from timing observations of $45$ millisecond pulsars between July 30th, 2004 to December 31st, 2015. The first five years of data on seventeen pulsars constituted the NANOGrav $5$-year dataset, which we previously published in 22. The $5$-year dataset was augmented with four years of data, reported as the $9$-year dataset in Arzoumanian et al. 2015, which came with the substantial improvements of new broadband instrumentation, a nearly twofold increase in the timing baseline for the original $17$ pulsars, and a more than twofold increase in the total number of observed sources to $37$ pulsars. The present extension of the $9$-year dataset is composed of two years of data that were observed and processed in a nearly identical fashion to the previous augmentation, with the addition of nine pulsars and the removal of one (see 8 for full details). Here we briefly review the instrumentation, observations, and basic data reduction of the entire dataset, referring the reader to 8, 6, and references therein for a thorough description. A sky map of of the pulsars in this data set is shown in Figure 1, with indicators of the time span and data volume for each pulsar.

We obtained all data using the 100-m Robert C. Byrd Green Bank Telescope (GBT) of the Green Bank Observatory greenbankobservatory.org/telescopes/gbt/ and the 305-m William E. Gordon Telescope (Arecibo) of Arecibo Observatory outreach.naic.edu/ao/. Sources within Arecibo’s declination range ($0^{\circ}<\delta<39^{\circ}$) were observed there due to its superior sensitivity, and only pulsars J1713+0747 and B1937+21 were observed at both telescopes. Excluding early portions from the $5$-year dataset, we observed each source roughly once a month for the entire dataset. In addition, some pulsars have been observed weekly in a campaign to increase our sensitivity to individual sources of GWs (Arzoumanian et al. 2014, Section 6.1). Specifically, two pulsars have been observed weekly at the GBT since $2013$ (PSRs J1737$+$0747 and J1909$-$3744), and five pulsars have been observed weekly at Arecibo since $2015$ (PSRs J0030$+$0451, J1640$+$2224, J1713$+$0747, J2043$+$1711, and J2317$+$1439).

During most epochs, This excludes the weekly observations, which were performed at 1.4 GHz only, as well as epochs for which receivers were unavailable for technical reasons. we observed sources in two widely separated frequency bands, in order to accurately remove the frequency-dependent dispersion delay introduced by the ionized interstellar medium (ISM). At the GBT, we used the 820-MHz and 1.4-GHz receivers for all observations. Since mechanical and time constraints prohibit alternating continually between the two receivers, observations in the two bands were always separated, typically by several days. At Arecibo, we observed all pulsars at 1.4 GHz, plus a second frequency band (centered at either 430 MHz or 2.3 GHz) chosen depending on the spectrum and ISM characteristics of each pulsar Pulsar J2317+1439 was originally observed with the 327 and 430-MHz receivers, but in 2014 we replaced the former with the 1.4-GHz receiver.. Pulsars observed at Arecibo are always observed in the two frequency bands one after another, separated by a few minutes.

For approximately the first six years, data were acquired with an identical pair of backend instruments that have since been decommissioned (GASP at the GBT, ASP at Arecibo). Since 2010 and 2012, respectively, the broadband-capable backend clones GUPPI (at the GBT) and PUPPI (at Arecibo) have been used for taking data.

II.2. Processing & Time-Of-Arrival Data

The $5$-year TOA dataset was left mostly untouched as a subset of the $11$-year dataset, except for reprocessing under DE$436$. All of the GUPPI and PUPPI profile data, however, were reprocessed from scratch to make a consistent set of TOAs. The TOAs were generated using standard template-matching cross-correlation methods, using only the total intensity profiles, producing one TOA per frequency channel per temporal subintegration. Existing template profiles were reused for pulsars that were part of the $9$-year dataset, and created for new pulsars.

An additional set of procedures culled “outlier”, low signal-to-noise (non-Gaussian distributed), or otherwise corrupt TOAs from the dataset using methods described in Vallisneri & van Haasteren 2017. The $11$-year dataset comprises a total of $309,201$ TOAs. All data reduction was completed using PSRCHIVE psrchive.sourceforge.net (Hotan et al. 2004) and custom NANOGrav processing scripts github.com/demorest/nanopipe.

II.3. Timing Models & Noise Analysis

Timing models from the $9$-year dataset were refit to the extended dataset and updated to include new parameters when deemed necessary on the basis of statistical significance tests. We fit timing models for newly-added pulsars using a procedure similar to that described in 6. All timing models were created or updated using the standard timing software TEMPO tempo.sourceforge.net and TEMPO2 bitbucket.org/psrsoft/tempo2.git (Hobbs et al. 2006; Edwards et al. 2006), and crosschecked for consistency.

A standard noise model was also fit simultaneously with the timing model as described in 6 and 8. Each pulsar’s white noise model includes a scale parameter on the TOA uncertainties (EFAC), an added variance (EQUAD), and a per-epoch variance (ECORR) for each observing system (i.e. a unique combination of backend and receiver). In addition, a red noise process for each pulsar was modeled by a power-law spectral density described by an amplitude and spectral index. The inclusion of a red process in the noise model was not favored by all pulsars, but we include it in all subsequent analyses since this does not affect parameter constraints. In the analyses described in the subsequent sections, we vary the pulsars’ red noise parameters and the parameters of the gravitational wave background, but fix the white noise parameters. Allowing the white noise parameters to vary does not alter the results, but significantly increases the computation time.

The SSE model used for the original analysis of the $5$-year dataset (Demorest et al. 2013) was DE$405$ (Standish 2004), while for the $9$-year dataset (Arzoumanian et al. 2015) all data (whether new, or from the $5$-year dataset) was modeled with DE$421$ (Folkner et al. 2009). For the $11$-year dataset we use DE$436$ (Folkner & Park 2016) as the fiducial SSE under which the data is processed and released. We do not need separate dataset releases for the different SSEs that we investigate in the following, since our GWB analysis incorporates marginalization over all affected processes, such as the individual timing and red-noise models.

III. Data Analysis Methods

Characterizing all deterministic and noise processes in each pulsar, as well as teasing out a putative GWB signature from the cross-correlation of large datasets, requires a robust and sophisticated statistical framework. In the following we describe the major new features of the NANOGrav PTA analysis framework, as updated from 7. Sec. III.1 describes our use of Bayesian inference as it pertains to computing GWB upper limits and detection statistics. Sec. III.2 outlines how a GWB manifests in our data as a long-timescale stochastic process with a distinctive correlation signature between pulsars. In Sec. III.3 we describe how the Solar System ephemeris model appears in our pulsar-timing analysis, and our new Bayesian scheme to mitigate its uncertainties. The structure of our generative signal and noise model is outlined in Sec. III.4, followed in Sec. III.5 by the definition of our frequentist estimator for the GWB amplitude and significance. Finally, in Sec. III.6 we list and provide links for all open-source software used in our GWB analysis.

We primarily employ Bayesian inference (see, e.g., Gregory 2005) to extract physical information from our data, deriving marginalized posterior distributions and credible regions, basing upper limits on credible intervals, and relying on ratios of evidences (a.k.a. Bayes factors) to compare models with different assumptions and parametrizations. We explore our high-dimensional parameter space stochastically, using the parallel-tempering Markov Chain Monte Carlo (MCMC) sampler (Ellis & van Haasteren 2017b) described in the appendices of Arzoumanian et al. 2014.

with $x=0.95$ and $N$ the number of (quasi-)independent samples Quasi-independence here refers to samples separated by one auto-correlation chain length. in the chain.

As our PTA dataset becomes longer and more sensitive, we expect that evidence for the presence of GWs will emerge in two phases: first, as red-spectrum processes with the same amplitude in each pulsar, and with spectral slope consistent with an SMBHB population; later (perhaps several years), and conclusively, as Hellings–Downs spatial correlations predicted for an isotropic GWB. We note that anisotropic GWBs will have different (but predictable) spatial correlations (Mingarelli et al. 2013; Taylor & Gair 2013; Mingarelli & Sidery 2014; Gair et al. 2014).

Correspondingly, we characterize evidence for a GWB in the 11-year dataset in two steps. We first obtain the Bayes factor for a dataset model that includes a red-spectrum process with common statistical properties in all pulsars (but is uncorrelated between them), against a model with only per-pulsar noise processes. This is signal-vs.-noise model selection. We then obtain the Bayes factor for Hellings–Downs inter-pulsar spatial correlations vs. no correlations at all. This is spatial-correlation model selection, which we consider the definitive scheme for GWB detection. We also perform variants of these comparisons—for instance, we compare the Hellings–Downs and uncorrelated process against processes with monopolar (akin to long-timescale clock errors) and dipolar (akin to SSE errors) spatial correlations.

For nested models (in our case, a signal-plus-noise model $\mathcal{H}_{1}$ and a noise-only model $\mathcal{H}_{0}$ obtained by fixing the GW amplitude to 0) we employ the Savage–Dickey formula (Dickey 1971)

For disjoint models (in our case, a model consisting of a Hellings–Downs-correlated red process plus pulsar noise, vs. a model consisting of a common-amplitude, spatially-uncorrelated red process plus pulsar noise) we use a product-space method (Carlin & Chib 1995; Godsill 2001; Hee et al. 2016). In this method we define a super-model that contains all parameters from all models under consideration, as well as an additional model-indexing variable that determines which model is ‘‘active’’ and used to evaluate the likelihood. This variable is technically discrete, but it can be sampled continuously and cast to an integer to choose the active model. (In our example, where the parameters are actually the same in both models, the index variable would simply toggle Hellings–Downs correlations in the evaluation of the likelihood.) The ratios of posterior probabilities for two model indices approximate the corresponding Bayes factor. We follow Cornish & Littenberg 2015 to estimate Bayes-factor uncertainties.

Evaluating the multi-pulsar likelihood is very computationally expensive when we account for inter-pulsar spatial correlations. In that case, we accelerate inference by running at least ten parallel copies of each spatially correlated analysis. These subchains can then be concatenated to form a much larger chain. Each subchain is analyzed to determine that it has ‘‘burned in’’ In MCMC analysis, some early sampled points must be disregarded before the chain can be considered to be sampling from the true posterior probability distribution. The disregarded early portion of the chain is called the “burn in” stage. before combining it with others. To derive upper limits and Savage–Dickey Bayes factors, we simply append the subchains together and proceed as described above. For product-space Bayes factors, we obtain the factor itself from the combined subchains, but we estimate uncertainties in each subchain separately, then add them in quadrature (Cornish & Littenberg 2015).

Arbitrary rules of thumb have been given to interpret the statistical significance of Bayes factors of different magnitudes (see, e.g., Jeffreys 1961; Kass & Raftery 1995), but it is hard to find agreement beyond the trivial statement that factors $\sim 1$ are inconclusive, while very large or small factors point to a strong preference for either model. In the context of a detection scheme, it seems appropriate to examine the frequentist distribution of the Bayes factor, and to set detection thresholds as a function of false-alarm probability (Vallisneri 2012). The sky-scramble and phase shifts methods (Taylor et al. 2017a; Cornish & Sampson 2016) have been proposed to produce a background distribution of the Bayes factor for which spatial correlations are effectively removed from the data. By contrast, we currently lack a practical approach to establish the significance of a common uncorrelated process; such an approach would likely involve a combination of inference runs on simulated data and cross-validation experiments, such as comparing results for subsets of the dataset. As we shall see, all the ephemeris-marginalized Bayes factors obtained in this paper are close to unity, and can be deemed inconclusive without a frequentist analysis.

III.2. Gravitational-wave strain spectrum

The observed timing residuals due to a GWB with characteristic strain $h_{c}(f)$ are described by the cross-power spectral density

where $\Gamma_{ab}$ is the overlap reduction function (ORF), which describes correlations between pulsars $a$ and $b$ in the array. In the case of an isotropic background from SMBHBs the ORF is given by Hellings & Downs 1983 (hereafter referred to as H.–D. correlations). Other correlated effects such as systematic errors in the Solar System ephemeris or clocks can also be described by a timing-residual spectrum that includes a different ORF.

In this paper we consider four models of the GWB spectrum:

A population of inspiraling SMBHBs in circular orbits, evolving by GW emission alone produces a characteristic GW-strain spectrum, expressed as

with $\alpha=-2/3$ (Phinney 2001). Different spectral slopes can be used to model relic radiation from the early Universe, under different assumptions for the equation of state of the Universe post-inflation/pre–Big-Bang-Nucleosynthesis (see Sec. V.3). We find it expedient to perform our analysis in terms of the timing-residual spectral index $\gamma=3-2\alpha$, such that

The fiducial SMBHB $\alpha=-2/3$ then corresponds to $\gamma=13/3$.

If SMBHBs remain coupled to the dynamics of their galactic environments as they evolve into the nanohertz band, the nanohertz GW strain spectrum will be more complex than described by Eq. (4). This may be the case if three-body scattering of stars from the galactic-center loss cone (Quinlan 1996; Sesana et al. 2006, e.g.) or interaction with a viscous circumbinary disk (Kocsis & Sesana 2011; Haiman et al. 2009, e.g.) are a stronger dynamical influence than GW emission at wide orbital separations. When the binary reaches milliparsec separations, GW emission will always be dominant. Sampson et al. 2015 introduced a broken power-law model,

to model such spectra, where the slope transitions from positive at low frequencies to the canonical $-2/3$ at higher frequencies. The frequency at which the transition occurs encodes information about the typical binary’s orbital evolution and astrophysical environment.

To characterize the GW-strain sensitivity of our dataset as a function of frequency, we adopt independent uniform priors for the dimensionless-strain amplitudes of each sine–cosine pair of red-process Fourier components (see Sec. III.4), corresponding to frequencies $k/T$, with $k=1,\ldots,N$, where $T$ is the longest timespan in the combined dataset, and $N$ (set to 50 in this paper) is the number of Fourier component pairs. We then derive a joint posterior for all amplitudes.

This model was introduced by Taylor et al. 2017b as a way to perform searches that are directly informed by detailed source-astrophysics simulations, and to sample the posteriors of the binary environment and dynamics parameters that affect the GW spectrum without generating a new simulation for each likelihood evaluation. In practice, we perform simulations over a grid in the parameter space of interest, and for each simulation we compute the GW characteristic strain spectrum. We then train a Gaussian process (Williams & Rasmussen 2006) to interpolate over all spectra in parameter space, allowing spectral amplitudes to be predicted at any other point with an associated normal uncertainty. We then use these predictions and uncertainties as priors on the strain amplitude at each frequency within the free-spectrum model.

III.3. Solar System ephemeris errors and uncertainties

A Solar System ephemeris is used in pulsar timing to convert observatory TOAs to an inertial frame centered at the Solar System barycenter, factoring out all effects due to Earth’s motion. The dominant correction to the TOAs is the Roemer delay—the classical light-travel time between the geocenter and the Solar System barycenter. Pulsar-timing studies have typically relied on the latest SSE released by JPL, adopting it as a model with fixed parameters—that is, without including any SSE parameter uncertainties or corrections in timing-model fits. In the early stages of our analysis of the NANOGrav 11-year dataset, we became aware that the choice of SSE among the latest few released by JPL has a measurable impact on our GWB upper limits and model-comparison Bayes factors. Indeed, the abundance and precision of NANOGrav’s measurements are now such that the accuracy to which we can estimate the Earth’s orbit around the SSB limits our sensitivity to GWs. SSE errors have been speculated on as a source of potential bias in PTA GW detection efforts (Tiburzi et al. 2016), but this paper marks the first time that this effect has been rigorously studied with real datasets.

The JPL SSEs https://ssd.jpl.nasa.gov/?ephemerides, as well as the French INPOP https://www.imcce.fr/inpop, fit the orbits and masses of a large set of Solar System bodies to a heterogenous dataset collected over the last few decades, using spacecraft ranging, direct planetary radar ranging, spacecraft VLBI, and (for the Moon) laser-ranging of retroreflectors left by the Apollo missions. The orbits are integrated numerically from initial conditions (“epoch” positions and velocities), which are the parameters that are fit for, together with other quantities such as the masses of minor Solar System bodies [but not planet masses, which are estimated separately from observed motions in planetary systems (Folkner et al. 2009)]. The resulting SSEs are distributed as Chebyshev polynomials over a range of dates; notably, they do not include estimates of orbit uncertainties and of possible systematics.

To investigate the effects of SSE errors, we repeated all upper-limit and model-comparison analyses in this paper using the four most recent JPL SSEs [DE421, released in 2008 (Folkner et al. 2009); DE430 (Folkner et al. 2014); DE435 (Folkner et al. 2016); DE436 (Folkner & Park 2016)]; for the simplest analysis, we used also the French INPOP$13$c (Fienga et al. 2014). The orbit of Earth relative to the Sun is consistent at the $10$-m level across these ephemerides, after accounting for an overall rotation w.r.t. the International Celestial Reference Frame, which originates from updated very-long-baseline-interferometry observations of spacecraft at Mars. However, the orbit of the Sun w.r.t. the SSB and (therefore) the orbit of Earth w.r.t. the SSB match only at the $100$-m level. This discrepancy is attributed to differences in the estimated masses and positions of Jupiter, Uranus, and Neptune. Hence, our GW analysis shows significant systematic differences among the upper limits and Bayes factors computed using different ephemerides. Near-future efforts may lead to improvements in the ephemeris accuracy that are appropriate for pulsar timing, namely: $(i)$ estimates of Jupiter’s orbit will be improved by including Juno spacecraft data in the SSE fit; $(ii)$ ranging data from Cassini may better estimate the mass of Uranus; $(iii)$ Gaia data may improve orbit estimates for Uranus and Neptune; and finally $(iv)$ pulsar-timing data may be used to improve the estimate of Neptune’s mass.

Thus, we present GW upper limits and model-comparison Bayes factors that are marginalized over these SSE uncertainty parameters. We regard these BayesEphem limits and Bayes factors as our fiducial results in this paper. To derive them, we constrain the outer-planet masses using the current IAU best estimates (IAU 2017; Jacobson et al. 2000; Jacobson et al. 2006; Jacobson 2014; Jacobson 2009), and use IAU uncertainties to set Gaussian priors. The rate of rotation about the ecliptic pole is left unconstrained. We experimented with setting priors for Jupiter’s orbital elements using estimated uncertainties, Folkner & Park 2017 estimate uncertainties in Jupiter and Saturn orbits by comparing fits that use independent subsets of the data for each planet. but we find better results using uninformative priors. This is not surprising, because Jupiter’s orbital elements are highly correlated with those of the other planets, and our linearized correction of Jupiter’s orbit cannot account for those correlations. Nevertheless, the resulting variations of Earth’s orbit are comparable with the systematic differences that we observe across JPL SSEs, which we take as evidence that the BayesEphem uncertainty parameters are representative of true SSE uncertainties.

III.4. Data model and likelihood

Except for Gaussian-process spectrum emulation and for the treatment of SSE errors, the data model used in this paper matches that of 7 closely, so we refer the reader to that publication for an overview of noise modeling, marginalization over timing-model parameters, our rank-reduced formalism for time-correlated processes (e.g., timing noise or GWB), and the PTA likelihood.

The rank-reduced formalism refers to the expansion of processes on a sine–cosine Fourier basis with frequencies $k/T$, where $T$ is the span between the minimum and maximum TOA in the array. The number of basis vectors is chosen to be high enough that inference results are insensitive to adding more: we use 30 for all applications except for the free-spectrum GWB model, for which we use 50.

As for the PTA likelihood, we introduced a significant change compared to 7. “ECORR” (jitter-like) noise is fully correlated for simultaneous observations at different observing frequencies, but fully uncorrelated in time. In 7, we treated ECORR degrees of freedom by assigning them “exploder” basis vectors, and then analytically marginalizing their coefficients simultaneously with timing-model, red-noise, and GWB coefficients. Doing so becomes computationally prohibitive when H.–D. correlations are included. In this paper, we include ECORR noise as block-diagonal entries (one block per epoch per backend–receiver system) in the otherwise diagonal white-noise covariance matrix, and invert the matrix using the fast Sherman & Morrison 1950 formula. Doing so eliminates a significant computational bottleneck.

As in 7, computational efficiency is also helped by fixing all white-noise parameters to their 1D maximum a posteriori values from single-pulsar noise studies. This choice is justified empirically by the very small variance of white-noise parameters.

Our upper-limit and model-comparison studies are performed under a variety of assumptions about the presence of red-spectrum processes: in addition to individual red-spectrum timing noise for each pulsar, we model the GWB as a spatially uncorrelated common process (a computational simplification appropriate in the weak-GWB limit, used in 7) and as a Hellings–Downs-correlated common process (our fiducial GWB model); we also consider common processes with different correlations (dipolar, as appropriate for SSE errors, and monopolar, as appropriate for long-timescale clock errors). Table 1 describes the nine models used in this paper, which are labeled 1, 2A–D, and 3A–D. In model-class $1$ only intrinsic pulsar noise processes are included; in model-class $2$ there are intrinsic pulsar noise processes, as well as non-GW noise processes that induce inter-pulsar spatial correlations (such as clock and SSE errors); in model-class $3$ we include a GWB signal. The roman characters given after the model-class number indicate the specific combination of noise and signal processes forming the model.

We perform each analysis by adopting each of the DE421, DE430, DE435, and DE436 (and occasionally INPOP$13$c) ephemerides as fixed-parameter models, and by marginalizing over SSE uncertainties using BayesEphem. Our Bayesian priors for all parameters are described in Table 2.

III.5. Optimal Statistic

which is a measure of the significance of inter-pulsar spatial correlations. When drawing comparisons between results produced using this frequentist technique and our Bayesian techniques, the relevant model selection is between models $3$A and $2$A.

III.6. Software

We generated most of the results in this paper using the open-source software package NX01 https://github.com/stevertaylor/NX01 (Taylor 2017), which implements the PTA likelihood and priors. NX01 was validated on a wide range of problems, including several 11-year analyses, by cross-comparison with the well-established PAL2 https://github.com/jellis18/PAL2 (Ellis & van Haasteren 2017a) and with NANOGrav’s new flagship package, enterprise https://github.com/nanograv/enterprise (Ellis et al. 2017). We perform MCMC using PTMCMCSampler https://github.com/jellis18/PTMCMCSampler (Ellis & van Haasteren 2017b), which implements a variety of proposal schemes (adaptive Metropolis, differential evolution, parallel tempering, etc.), which can be used together in the same run.

As a companion to this paper, we are releasing a Docker https://github.com/nanograv/11yr_stochastic_analysis image that contains a full stack of our software (including all required libraries), and that can be used to reproduce the upper limits, Bayes factors, as well as many of the figures of this paper, using enterprise.

IV. Results

As discussed in Sec. III.4, we perform analyses for variants of our data model that reflect different assumptions about common red-spectrum processes, as listed in Table 1, and under four JPL ephemerides as well as BayesEphem (in select cases we include also the French INPOP13, which yields results broadly similar to DE430).

Following 7, we present upper limits on the strain amplitude of a GWB modeled as a power law and as a free spectrum (see Sec. III.2).

Comparing the columns of Table 4 shows how the upper limits vary under different assumptions on the presence of spatially correlated common processes in the data. The limits are slightly more stringent if we model the GWB as a spatially uncorrelated common process (model 2A in the second column), indicating that Hellings–Downs correlations help the likelihood isolate a GW-like signal (whether real, or due to random noise fluctuations). Introducing additional spatially correlated processes (with ephemeris-error–like dipolar correlations, clock-error–like monopolar correlations, or both, corresponding to models 3B, 3D, and 3C) reduces upper limits for the individual ephemerides but not for BayesEphem, suggesting that the same realization of inter-pulsar signal correlations can be picked up by different ORFs, and that dipole and monopole processes can absorb some, but not all, of the systematic bias caused by ephemeris error.

In Figure 2 we show the 95% upper limit for the amplitude of an uncorrelated common process (model 2A) as a function of $\gamma$. In the absence of red noise, and if the lowest sampling frequency ($1/T$) dominated our sensitivity, we would expect these constraints to scale as $\propto T^{-\gamma/2}$, where $T$ is the longest timing baseline across the entire PTA. We find the actual scaling to be closer to $\propto T^{-0.4\gamma}$, indicating that red noise is present and that more than one frequency component contributes to the likelihood.

IV.2. Bayesian model-comparison evidence for GWs

In Tables 5 and 6 and in Fig. 4, we show Bayes factors for two sets of model comparisons performed on the 11-year dataset to quantify the statistical evidence for a stochastic GWB and for coherent sources of systematic errors that lead to spatially correlated residuals. The first four columns of Table 5 and the graph on the left of Fig. 4 are diagnostic of the multilevel decision scheme outlined above in Sec. III.1. Adopting the JPL ephemerides as fixed-parameter models, the data favor the presence of a common uncorrelated process in all pulsars, to various degrees and especially so for DE430, and they favor slightly the presence of Hellings–Downs inter-pulsar correlations. However, this preference disappears if we marginalize over the ephemeris uncertainties.

The convergence of the solid lines to a flatter common shape demonstrates that our modeling of ephemeris uncertainties bridges the four ephemerides successfully, removing spurious evidence for GWs, or potentially absorbing a true GW signal. However, if a true GW signal is present, it happens to be significantly covariant with the systematic differences in the Roemer delays induced by the last few ephemerides; furthermore, the signal appears to weaken as we shift from older (DE421, DE430) to newer, plausibly more accurate ephemerides (DE435, DE436), although this trend is not entirely consistent. In this paper, we do not attempt to quantify whether these circumstances are realized often in the ensemble of possible datasets similar to ours; nevertheless, these circumstances motivate our choice of marginalizing over ephemeris uncertainties as the principled Bayesian strategy for our analysis.

The six rightmost columns of Table 5, as well as Table 6 and the graph on the right of Figure 4, document the degree to which the data favor the presence of timing-residual components with different spatial correlations. Components with both dipolar (ephemeris-error–like) and monopolar (clock-error–like) correlations are disfavored, although this conclusion is significantly weakened if we marginalize over ephemeris uncertainties. At the same time, the evidence for quadrupolar (GWB-like) correlations is weakened when the model allows for other spatially correlated processes. This is not unexpected, since spatial correlations with different multipolar structures only become truly orthogonal in the limit of many equally low-noise pulsars.

Indeed, discrimination of monopolar, dipolar, and quadrupolar correlation signatures will improve as our datasets gain more and more pairs of high–timing-precision pulsars with a broad distribution of angular separations. We plan to characterize discrimination requirements (on pulsar number, timing quality, and sky position) in our upcoming paper on SSE error modeling.

We performed a small number of simulations to test the impact of BayesEphem on our GWB detection prospects over the next few years. To this end, we produced realistic $15$-yr datasets To produce the datasets, we used actual observation epochs for the $34$ NANOGrav pulsars, and set residuals equal to white measurement noise plus red-spectrum intrinsic noise, at levels consistent with those estimated for the actual data (8). We rescaled TOA uncertainties by a factor $1.5$, which calibrates the noise-only simulated dataset so that its $11.4$ year “slice” has the same (DE$436$, model $2$A) GWB upper limit as the real data. We extended the dataset baseline to $15$ years by drawing observation epochs and TOA measurement errors from distributions of these quantities over the last $3$ years of real data. using DE$436$ and injecting GWBs of various amplitudes, and we analyzed the full datasets, as well as their $11.4$ yr “slices,” using DE$430$ and BayesEphem. We chose DE$430$ because it led to the highest signal-vs.-noise Bayes factor (model $2$A-vs.-$1$) and upper limits for the actual data.

For a noise-only simulation, we find that unmodeled systematic offsets between DE$436$ and DE$430$ are interpreted as a common red-spectrum process with a signal-vs.-noise Bayes factor (model $2$A-vs.-$1$) of $\sim 2$ in $11.4$ years of data, and $\sim 20$ in $15$ years of data. By contrast, BayesEphem is able to account for the offsets, reducing Bayes factors to levels consistent with noise fluctuations. As we increase the injected GWB amplitude, model $2$A-vs.-$1$ Bayes factors remain low for $11.4$ yrs of data, even for amplitudes comparable to our fiducial upper limits. The same is true for model $3$A-vs.-$2$A Bayes factors (the definitive spatial-correlation test for GWBs), which are plotted in Figure 6.

For $15$ years of data, the scaling of Bayes factors with injected GWB amplitude is comparable for both DE$430$ and BayesEphem. Remarkably, the potential covariance of BayesEphem parameters with GWB amplitude does not inhibit signal detection in the near future, even at astrophysically-pessimistic levels ($\sim 5\times 10^{-16}$, consistent with Sesana et al. 2016). Thus, while SSE errors may spuriously produce early signs of a GWB (i.e. a common red-spectrum process), their mitigation with BayesEphem will not impair prospects for near-future GWB detection. We regard our simulations as conservative, since additional pulsars, as well as improved timing precision and SSE accuracy, will accelerate progress toward detection.

IV.3. Optimal statistic

Table 8compares the noise-marginalized optimal statistic computed for Hellings–Downs spatial correlations with variants of the statistic that model dipolar and monopolar correlations. In addition to computing the optimal statistic using individual ephemerides, we also use BayesEphem to marginalize over the ephemeris uncertainty. For all of these analyses, we find no evidence for a common process with either Hellings–Downs, monopolar, or dipolar spatial correlations.

The upper half of Figure 7 shows the mean noise-marginalized cross-correlated power between pulsar pairs as a function of angular distribution, averaged into 10 degree bins. There is no evidence of the Hellings–Downs correlations characteristic of isotropic GWBs. The lower half of the plot shows a histogram of angular separations for the pulsar pairs in our dataset: NANOGrav is currently most sensitive to angular separations between $30^{\circ}$ and $60^{\circ}$, which correspond to the smallest errors in the cross-correlation plot.

IV.4. Comparison of 9-year and 11-year results

The 9-year analysis of 7 adopted DE421 as a fixed-parameter model without uncertainties, and did not include Hellings–Downs correlations. Thus, a straight comparison can be made with the 11-year DE421 model-2A results: the $\gamma=13/3$ upper limit remains at $1.5\times 10^{-15}$, while the $\gamma=13/3$ Bayes factor vs. pulsar noise changes from $0.81$ to $8.3$; however, this comparison is not very significant given what we have learned about ephemeris errors.

The model-2A Bayes factors vs. pulsar noise are $0.910(7)$ for $\gamma=13/3$ and $1.210(4)$ for $\gamma\in$, while they are $1.27(1)$ and $2.29(3)$ for model 3A. All Bayes factors under BayesEphem are comparably uninformative for the 9-year and 11-year datasets.

V. Limits on Astrophysical Models

Some of the most exciting science made possible by the NANOGrav data is realized when we use the GWB constraints to confront the astrophysics of various source populations. Now, the most likely source population for PTAs is SMBHBs. In 7 we introduced simple PTA constraints on SMBHB population parameters, but due to methodological limitations we were unable to deliver a realistic analysis, i.e. we derived constraints for the parameters describing a broken–power-law spectrum, and then reinterpreted those constraints in terms of SMBHB effects that could alter the spectrum, taken one at a time. In this paper we adopt the modeling framework developed by Taylor et al. 2017b to go much further: we use a set of population-synthesis simulations to explore the effects of population parameters on the GWB spectrum, then constrain those population parameters directly from the data. We apply the same method also to the most recent cosmic-string models.

The eccentricity evolution in this model follows the prescription first derived in Quinlan 1996, and later expanded upon in Sesana 2010. However, recent work in Rasskazov & Merritt 2017a (Sesana et al. 2011; Gualandris et al. 2012; Mirza et al. 2017, see also) has shown that eccentricity evolution can be damped by the rotation of the central stellar bulge, which would lessen the effect of extreme initial eccentricities.

V.2. Cosmic strings

Cosmic strings are linear topological defects that can form in the early Universe as a result of symmetry-breaking phase transitions (Kibble 1976; Vilenkin 1981; Vilenkin 1985; Vilenkin & Shellard 2000). Strings that form with lengths greater than the horizon are known as “long” or “infinite” strings, while smaller strings form loops. If two strings meet one another they can exchange partners, and small portions of string can be chopped off with a reconnection probability $p$. For classical strings $p=1$, but String-Theory–inspired models may have $p<1$. This is due to the fact that fundamental strings interact probabilistically, and also that in these models an intersection occurring in the usual three spatial dimensions need not occur in higher compactified dimensions. Cosmic string networks evolve toward an attractor solution known as the “scaling regime” in which the statistical properties of the system (such as the average size of loops or the distance between long strings) scale with the cosmic time, and the energy density of the string network is a small constant fraction of the radiation or matter density. Cosmic strings have tensions equal to their mass per unit length, $\mu$. This tension is so high that strings oscillate relativistically under their own tension, decaying solely through the emission of GWs, and shrinking in size. The formation of loops and their subsequent decay by GW emission is the mechanism by which the string network loses energy and reaches the scaling regime. The GW spectrum from cosmic string networks is exceptionally broadband, covering all regions of LIGO, LISA, and PTA sensitivity. For our purposes, we describe the parameter space of cosmic strings in terms of their dimensionless tension, $G\mu/c^{2}$, and their reconnection probability, $p$.

We take a more self-consistent approach than previous PTA analyses. Rather than re-fit posterior samples (from power-law or free spectrum searches) to cosmic string models (Arzoumanian et al. 2016; Lentati et al. 2015), we train a GP interpolant on output from the most up-to-date string population simulations. Blanco-Pillado & Olum 2017 and Blanco-Pillado et al. 2018 performed a complete end-to-end calculation of the stochastic GW background expected from a network of cosmic strings, namely: $(i)$ simulation of the long-string network to find a representative sample of loop sizes and shapes; $(ii)$ modeling of loop shape deformations due to gravitational back-reaction; $(iii)$ GW spectrum computed for each loop; $(iv)$ evaporation and production modeled to find the distribution of loops over $z$; $(v)$ integration of the GW spectrum of each loop over the redshift-dependent loop distribution; and finally $(vi)$ integration over cosmological time to find the present-day GW background.

The output from these simulations corresponds to GW energy density spectra at a range of string tension values, $G\mu/c^{2}$, over $25$ orders of magnitude in frequency and has been made publicly available. http://cosmos.phy.tufts.edu/cosmic-string-spectra/ We convert these to characteristic strain, then at each frequency-bin in our PTA analysis we train a GP to emulate the strain as a function of string tension. We expand our model to include reconnection probability, $p$, by analytically scaling the fiducial $p=1$ strain spectrum by $(1/p)^{1/2}$ (Sakellariadou 2005). We then use this model (with all features of the cosmic-string spectrum included) to analyze the NANOGrav 11-year dataset. We do not model signal finiteness or anisotropy due to bright resolvable cosmic-string bursts, since this is only expected when initial loop sizes are very small ($\lesssim 10^{-8}$) (Kuroyanagi et al. 2017).

Figure 12shows the $95\%$ upper limit on string tension as a function of reconnection probability. The shaded region enclosed by the solid black line indicates parameter space that is excluded by the NANOGrav 11-year dataset under the assumptions of the Blanco-Pillado & Olum 2017 cosmic string simulations. For $p=1$ the string tension is constrained to be $G\mu/c^{2}<5.3(2)\times 10^{-11}$. At this level we would not expect any measurable effects in the CMB power spectrum, nor through gravitational lensing (Blanco-Pillado et al. 2018). PTAs are currently the best experiment with which to detect cosmic strings, and to place stringent limits on the string parameter space.

By contrast, the NANOGrav 9-year dataset (7) constraints on string tension (shown as an excluded region with a dashed line boundary) were computed under the assumptions of older string simulations (Blanco-Pillado et al. 2014), and were obtained by re-sampling the posterior distribution of a power-law GWB spectrum. For $p=1$ the string tension was constrained to be $G\mu/c^{2}<1.3\times 10^{-10}$. Finally, even though the most recent EPTA constraints on cosmic strings (Lentati et al. 2015) were not computed under the assumptions of the Blanco-Pillado et al. 2014 simulations, in 7 the constraints were converted to get a corresponding limit on the string tension of $G\mu/c^{2}<8.6\times 10^{-10}$. Thus, the constraints on cosmic string tension from the NANOGrav 11-year dataset are $2.5$ times better than the NANOGrav 9-year dataset, and $16.2$ times better than the most recent EPTA analysis. The $9$- to $11$-year improvement is to be expected, since BayesEphem analyses of the $11$-year dataset give consistently more constraining GWB limits than DE$421$ analyses of the $9$-year dataset. There are a few other notable caveats to these comparisons; $(i)$ the NANOGrav $9$-year and EPTA analyses were performed under a fixed JPL SSE model, while the NANOGrav $11$-year analysis uses BayesEphem; $(ii)$ the simulation advances of Blanco-Pillado & Olum 2017 with respect to Blanco-Pillado et al. 2014 impede a direct comparison. However, the additional $\sim 2$ years of data in the new NANOGrav dataset, the new SSE uncertainty modeling, and the improved end-to-end analysis with simulated cosmic-string spectra all combine to increase NANOGrav’s sensitivity to the cosmic-string parameter space.

V.3. Primordial gravitational-waves

According to the theory of inflation, quantum fluctuations in the spacetime geometry of the early Universe are amplified to cosmological scales. Inflation leaves a background of relic primordial GWs that may be observable today (Grishchuk 1976; Grishchuk 1977; Starobinsky 1980; Linde 1982; Fabbri & Pollock 1983). Studies of the cosmic microwave background (CMB) that attempt to observe these GWs indirectly through their imprint of tensor-mode CMB polarizations are limited to probing the surface of last scattering, roughly 300,000 years after the Big Bang (Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997; BICEP2/Keck et al. 2015). By contrast, GW observations can in principle observe a much earlier epoch in the history of the Universe, extending back to as little as $10^{-32}$ s post Big Bang. Indeed, the spectral index of the primordial GWB is determined by the equation-of-state parameter $w$ in the immediate post-inflation, pre–Big-Bang-Nucleosynthesis Universe, and by the tensor index $n_{t}$, which depends on the detailed dynamics of inflation (see Grishchuk 2005 and references therein). The primordial spectral dependencies are typically stated in terms of GWB-$\alpha$ as in Lasky et al. 2016 and 7. We can express GWB-$\gamma$ (see Eq. 4) for a primordial spectrum as

These limits constrain the energy density spectrum of the primordial GWB by way of

where $h$ is the dimensionless Hubble parameter, $H_{0}=100$ km s-1 Mpc-1, and $h_{c}$ is the characteristic GW strain. For a radiation-dominated post-inflationary Universe, we obtain

after marginalizing over SSE uncertainties. This is a 20% improvement over the result quoted in 7; that number, however, should be revised upward significantly due to SSE bias. Referring back to the bottom panel of Figure 3, we see that the energy-density sensitivity of our PTA dataset is dominated by the lowest few frequencies, which individually have $95\%$ upper limit values of $\sim 10^{-9}$, but which in combination beat the limit down to the value quoted in Equation 11.

VI. Summary and Conclusions

This paper reports on the search for an isotropic stochastic GW background (GWB) in NANOGrav’s $11$-year dataset. We targeted a GW signal with predominantly low-frequency power, and so analyzed only those pulsars that have greater than $3$ years of observations, corresponding to $34$ out of the $45$ in the data release. Our investigations encompassed different models of the GWB strain spectrum, spatial correlations between pulsars, and Solar System ephemeris (SSE). The latter influence was rigorously studied, and led to the major discovery of this paper:

We found significant variations in GW upper-limits and detection statistics when the dataset was analyzed under different published models of the SSE. These models are primarily from the Jet Propulsion Laboratory (JPL), ranging from DE$421$ to DE$436$. We also performed a limited analysis with INPOP$13$c.

The ratio of Bayesian evidences between models that include a GWB versus only intrinsic pulsar noise processes varies between $\sim 2$ and $\sim 26$ in favor of a GWB, while the odds favoring GW-induced spatial correlation between pulsars vary between $1.18:1$ and $1.63:1$. The frequentist analog to the Bayesian odds-ratio (known as the “optimal-statistic”) gives a signal-to-noise ratio for GW-induced spatial correlations that varies between $0.57$ and $0.87$.

This discovery has major ramifications on how we interpret previous PTA results, and also how our analysis methodology must be revised for future searches.

We formulated a perturbative model (“BayesEphem”) that acts to bridge the systematic offsets in the various published models of the SSE, resulting in the first pulsar-timing constraints on GWs that are robust against Solar System uncertainties. This model corrects for coordinate-frame drift, uncertainties in gas-giant masses, and uncertainties in Jupiter’s orbital elements.

Under this new model, the upper limit on the strain amplitude becomes $1.34\times 10^{-15}$ for a common red-spectrum process, and $1.45\times 10^{-15}$ for a GWB. Adding further spatially-correlated processes in the model served to worsen these limits only slightly.

The evidence ratio for models that include a GWB versus only intrinsic pulsar noise processes is $1$ for a GWB with fixed spectral slope, and $0.70$ if the spectral slope is varied. The odds ratio favoring GW-induced spatial correlations between pulsars is $1.08:1$ if the spectral slope is fixed, or $1.15:1$ if the slope is varied. The frequentist optimal-statistic gives a signal-to-noise ratio for GW-induced spatial correlations of $0.09$, where the spectral slope is necessarily fixed at the fiducial value of $-2/3$. Both the Bayesian and frequentist analysis show inconclusive evidence for a GW-like red-spectrum process and quadrupolar inter-pulsar spatial correlations.

We also performed a systematic study of spatially-correlated processes in the PTA dataset under different ephemerides, tabulating upper limits and evidence ratios for various combinations of a common red-spectrum process, GWB, stochastic clock error, and stochastic SSE uncertainty. With BayesEphem the presence of these additional spatially-correlated processes slightly worsens the GW upper limits, but all remain broadly consistent within uncertainties. Dipole spatial correlations between pulsars seem most disfavored under BayesEphem, likely because we have dealt with the most plausible source of such correlations with our deterministic SSE-uncertainty modeling. Uncertainties in the evidence and odds ratios (in addition to their absolute values being around unity) prevent us being able to make strong statements. The NANOGrav $11$-year dataset is only weakly informative of spatial correlations between pulsars.

Over the last few years, the PTA community has made great strides in gathering ever larger, higher-quality datasets, and in developing sophisticated analysis methods that can deal with the complex noise budgets and subtle systematics typical of pulsar timing, while interfacing ever more closely and robustly with the astrophysics of GW sources. The sequence of recent stochastic-GW papers [for NANOGrav, Demorest et al. 2013, Arzoumanian et al. 2016, this paper] is a fitting witness to this growth. We expect this effort to be rewarded by nanohertz GW detection within the next several years (Taylor et al. 2016), if the steadfast pursuit of methodological rigor and physical insight remains our cynosure.

References