The NANOGrav 15-year Data Set: Evidence for a Gravitational-Wave Background

The NANOGrav Collaboration

Corresponding author: The NANOGrav Collaboration [email protected] $^{21}$ National Radio Astronomy Observatory, 1003 Lopezville Rd., Socorro, NM 87801, USA

$^{22}$ Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA

$^{23}$ Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA

$^{24}$ School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA

$^{25}$ Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA

$^{26}$ Bionic Health, 800 Park Offices Drive, Research Triangle Park, NC 27709

$^{27}$ Department of Astronomy, University of Maryland, College Park, MD 20742

$^{28}$ Center for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771

$^{29}$ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

$^{30}$ Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA

$^{31}$ Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

$^{32}$ Department of Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI 48109, USA

$^{33}$ Department of Astronomy, University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA

$^{34}$ Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

$^{35}$ University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA

$^{36}$ Department of Physics, State University of New York at Oswego, Oswego, NY, 13126, USA

$^{37}$ NASA Marshall Space Flight Center, Huntsville, AL 35812, USA

$^{38}$ Department of Physics and Astronomy, Haverford College, Haverford, PA 19041, USA

$^{39}$ Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada

$^{40}$ Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA

$^{41}$ Department of Physics, University of California, Berkeley, CA 94720, USA

$^{42}$ Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA

$^{43}$ E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

$^{44}$ Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK

$^{45}$ International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia

$^{46}$ Department of Physics, Yale University, New Haven, CT 06520, USA

$^{47}$ Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

$^{48}$ Department of Astronomy, Yale University, 52 Hillhouse Ave, New Haven, CT 06511

$^{49}$ Black Hole Initiative, Harvard University, 20 Garden Street, Cambridge, MA 02138

$^{50}$ Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George St., Toronto, ON M5S 3H4, Canada

$^{51}$ Department of Physics, Lafayette College, Easton, PA 18042, USA

$^{52}$ Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA

$^{53}$ Institute of Physics and Astronomy, Eötvös Loránd University, Pázmány P. s. 1/A, 1117 Budapest, Hungary

$^{54}$ Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA

$^{55}$ Department of Physics, University of Puerto Rico, Mayagüez, PR 00681, USA

$^{56}$ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

$^{57}$ Department of Physics, Texas Tech University, Box 41051, Lubbock, TX 79409, USA

$^{58}$ Department of Physics, Penn State Abington, Abington, PA 19001, USA

$^{59}$ Institute for Theoretical Physics, University of Münster, 48149 Münster, Germany

$^{60}$ Giant Army, 915A 17th Ave, Seattle, WA 98122

$^{61}$ Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA

$^{62}$ Center for Astrophysics, Harvard University, 60 Garden St, Cambridge, MA 02138

$^{63}$ Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA

$^{64}$ Department of Physics, Ben-Gurion University of the Negev, Be’er Sheva 84105, Israel

$^{65}$ Feza Gürsey Institute, Bogazici University, Kandilli, 34684, Istanbul, Turkey

$^{66}$ Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Callinstrasse 38, D-30167, Hannover, Germany

$^{67}$ Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208

$^{68}$ Adler Planetarium, 1300 S. DuSable Lake Shore Dr., Chicago, IL 60605, USA

Introduction

Almost a century had to elapse between Einstein’s prediction of gravitational waves (GWs, Einstein [45]) and their measurement from a coalescing binary of stellar-mass black holes (Abbott et al. [1]). However, their existence had been confirmed in the late 1970s through measurements of the orbital decay of the Hulse–Taylor binary pulsar (Hulse & Taylor [76]; Taylor et al. [133]). Today, pulsars are again at the forefront of the quest to detect GWs, this time from binary systems of central galactic black holes.

Black holes with masses of $10^{5}$–$10^{10}\ M_{\odot}$ exist at the center of most galaxies and are closely correlated with the global properties of the host, suggesting a symbiotic evolution (Magorrian et al. [94]; McConnell & Ma [96]). Galaxy mergers are the main drivers of hierarchical structure formation over cosmic time (Blumenthal et al. [27]) and lead to the formation of close massive–black-hole binaries long after the mergers (Begelman et al. [26]; Milosavljević & Merritt [101]). The most massive of these (supermassive black-hole binaries, SMBHBs, with masses $10^{8}$–$10^{10}\ M_{\odot}$) emit GWs with slowly evolving frequencies, contributing to a noise-like broadband signal in the nHz range (the GW background, GWB; Rajagopal & Romani [116]; Jaffe & Backer [78]; Sesana et al. [128]; McWilliams et al. [98]; Burke-Spolaor et al. [28]). If all contributing SMBHBs evolve purely by loss of circular orbital energy to gravitational radiation, the resultant GWB spectrum is well described by a simple $f^{-2/3}$ characteristic-strain power law (Phinney [111]).

However, GWB signals that are not produced by populations of inspiraling black holes may also lie within the nHz band; these include primordial GWs from inflation, scalar-induced GWs, and GW signals from multiple processes arising due to cosmological phase transitions, such as collisions of bubbles of the post-transition vacuum state, sound waves, turbulence, and the decay of any defects such as cosmic strings or domain walls that may have formed (see, e.g., Guzzetti et al. [61]; Caprini & Figueroa [29]; Domènech [43], and references therein).

The detection of nHz GWs follows the template outlined by Pirani [112, 113], whereby we time the propagation of light to measure modulations in the distance between freely falling reference masses. Estabrook & Wahlquist [49] derived the GW response of electromagnetic signals traveling between Earth and distant spacecraft, sparking interest in low-frequency GW detection. Sazhin [126] and Detweiler [41] described nHz GW detection using Galactic pulsars and (effectively) the solar system barycenter as references, relying on the regularity of pulsar emission and planetary motions to highlight GW effects. The fact that pulsars are such accurate clocks enables precise measurements of their rotational, astrometric, and binary parameters (and more) from the times-of-arrival of their pulses, which are used to develop ever-refining end-to-end timing models. Hellings & Downs [70] made the crucial suggestion that the correlations between the time-of-arrival perturbations of multiple pulsars could reveal a GW signal buried in pulsar noise; Romani [121] and Foster & Backer [51] proposed that a pulsar timing array (PTA) of highly stable millisecond pulsars (Backer et al. [23]) could be used to search for a GWB. Nevertheless, the first multi-pulsar, long-term GWB limits were obtained by analyzing millisecond-pulsar residuals independently, rather than as an array (Stinebring et al. [132]; Kaspi et al. [84]).

From a statistical-inference standpoint, the problem of detecting nHz GWs in PTA data is analogous to GW searches with terrestrial and future space-borne experiments, in which the propagation of light between reference masses is modeled with physical and phenomenological descriptions of signal and noise processes. It is distinguished by the irregular observation times, which encourage a time- rather than Fourier-domain formulation, and by noise sources (intrinsic pulsar noise, interstellar-medium–induced radio-frequency–dependent fluctuations, and timing-model errors) that are correlated on timescales common to the GWs of interest. This requires the joint estimation of GW signals and noise, which is similar to the kinds of global fitting procedures already used in terrestrial GW experiments, and proposed for space-borne experiments. GW analysts have therefore converged on a Bayesian framework that represents all noise sources as Gaussian processes [144, 145], and relies on model comparison (i.e., Bayes factors, which are ratios of fully marginalized likelihoods) to define detection (see, e.g., Taylor [134]). This Bayesian approach is nevertheless complemented by null hypothesis testing, using a frequentist detection statistic⁵ (the “optimal statistic” of Anholm et al. [12]; Demorest et al. [39]; Chamberlin et al. [32]) averaged over Bayesian posteriors of the noise parameters (Vigeland et al. [147]).

The GWB—rather than GW signals from individually resolved binary systems—is likely to become the first nHz source accessible to PTA observations (Rosado et al. [123]). Because of its stochastic nature, the GWB cannot be identified as a distinctive phase-coherent signal in the way of individual compact-binary-coalescence GWs. Rather, as PTA data sets grow in extent and sensitivity one expects to first observe the GWB as excess low-frequency residual power of consistent amplitude and spectral shape across multiple pulsars ([122, 115]). An observation following this behavior was reported in 2020 (Arzoumanian et al. [19], henceforth NG12gwb) for the 12.5-year data set collected by the North American Nanohertz Observatory for Gravitational waves (NANOGrav, McLaughlin [97]; Ransom et al. [117]), and then confirmed ([58, 33]) by the Parkes Pulsar Timing Array (PPTA, Manchester et al. [95]) and the European Pulsar Timing Array (EPTA, Desvignes et al. [40]), following many years of null results and steadily decreasing upper limits on the GWB amplitude.

A combined International Pulsar Timing Array (IPTA, Perera et al. [109]) data release consisting of older data sets from the constituent PTAs also confirmed this observation (Antoniadis et al. [13]). Nevertheless, the finding of excess power cannot be attributed to a GWB origin merely by the consistency of amplitude and spectral shape, which could arise from intrinsic pulsar processes of similar magnitude ([60, 152]), or from a common systematic noise such as clock errors (Tiburzi et al. [140]). Instead, definitive proof of GW origin is sought by establishing the presence of phase-coherent inter-pulsar correlations with the characteristic spatial pattern derived by Hellings and Downs (1983, henceforth HD): for an isotropic GWB, the correlation between the GW-induced timing delays observed at Earth for any pair of pulsars is a universal, quasi-quadrupolar function of their angular separation in the sky. Even though this correlation pattern is modified if there is anisotropy in the GWB—which may be the case for a GWB generated by a SMBHB population ([104, 136, 37, 103, 105, 120]—the HD template is effective for detecting even anisotropic GWBs in all but the most extreme scenarios ([37, 36, 139, 25, 10]).

In this letter we present multiple lines of evidence for an excess low-frequency signal with HD correlations in the NANOGrav 15-year data set (Figure 1). Our key results are as follows. The Bayes factor between an HD-correlated, power-law GWB model and a spatially uncorrelated common-spectrum power-law model ranges from 200 to 1,000, depending on modeling choices (Figure 2). The noise-marginalized optimal statistic, which is constructed to be selectively sensitive to HD-correlated power, achieves a signal-to-noise ratio of $\sim5$ (Figure 3 and Figure 4). We calibrated these detection statistics by removing correlations from the 15-year data set using the phase-shift technique, which removes inter-pulsar correlations by adding random phase shifts to the Fourier components of the common process (Taylor et al. [137]). We find false-alarm probabilities of $p = 10^{-3}$ and $p = 5 \times10^{-5}$ for the observed Bayes factor and optimal statistic, respectively (see Figure 3).

For our fiducial power-law model ($f^{-2/3}$ for characteristic strain and $f^{-13/3}$ for timing residuals) and a log-uniform amplitude prior, the Bayesian posterior of GWB amplitude at the customary reference frequency $1\ \mathrm{yr}^{-1}$ is $A_{\mathrm{GWB}} = 2.4^{+0.7}_{-0.6} \times10^{-15}$ (median with 90% credible interval), which is compatible with current astrophysical estimates for the GWB from SMBHBs (e.g., Burke-Spolaor et al. [28]; Agazie et al. [3]. This corresponds to a total integrated energy density of $\Omega_{\mathrm{gw}} = 9.3^{+5.8}_{-4.0} \times10^{-9}$ or $\rho_{\mathrm{gw}} = 7.7^{+4.8}_{-3.3} \times10^{-17}$ ergs cm$^{-3}$ (assuming $H_0 = 70$ km/s/Mpc) in our sensitive frequency band. For a more general model of the timing-residual power spectral density with variable power-law exponent $-\gamma$, we find $A_{\mathrm{GWB}} = 6.4^{+4.2}_{-2.7} \times10^{-15}$, and $\gamma= 3.2^{+0.6}_{-0.6}$. See panel (b) of Figure 1 for $A_{\mathrm{GWB}}$ and $\gamma$ posteriors. The posterior for $\gamma$ is consistent with the value of 13/3 predicted for a population of SMBHBs evolving by GW emission, although smaller values of $\gamma$ are preferred; however, the recovered posteriors are consistent with predictions from astrophysical models (see Agazie et al. [3]). We also note that, unlike our detection statistics (which are calibrated under our modeling assumptions), the estimation of $\gamma$ is very sensitive to minor details in the data model of a few pulsars.

The rest of this paper is organized as follows. We briefly describe our data set and data model in § 2. Our main results are discussed in detail in § 3 and § 4; they are supported by a variety of robustness and validation studies, including a spectral analysis of the excess signal (§ 5.2), a correlation analysis that finds no significant evidence for additional spatially correlated processes (§ 5.3), and cross-validation studies with single-telescope data sets and leave-one-pulsar-out techniques (§ 5.4). In the past two years we have performed an end-to-end review of the NANOGrav experiment, to identify and mitigate possible sources of systematic error or data set contamination: our improvements and considerations are partly described in a set of companion papers: on the NANOGrav statistical analysis as implemented in software (Johnson et al. [80]), on the 15-year data set (Agazie et al. [3], hereafter NG15), and on pulsar models (Agazie et al. [3], hereafter NG15detchar). More companion papers address the possible SMBHB (Agazie et al. [3]) and cosmological (Afzal et al. [2]) interpretations of our results, with several more GW searches and signal studies in preparation. We look forward to the cross-validation analysis that will become possible with the independent data sets collected by other IPTA members.

The 15-year Data Set and Data Model

The NANOGrav 15-year data set⁷ (NG15) contains observations of 68 pulsars obtained between July 2004 and August 2020 with the Arecibo Observatory (Arecibo), the Green Bank Telescope (GBT), and the Very Large Array (VLA), augmenting the 12.5-year data set (Alam et al. [8]) with 2.9 years of timing data for the 47 pulsars in the previous data set, and with 21 new pulsars.⁸ For this paper we analyze narrow-band times of arrivals (TOAs), which are computed separately for sub-bands of each receiver, and focus on the 67 pulsars with a timing baseline $\ge3$ years. We adopt the TT(BIPM2019) timescale and the JPL DE440 ephemeris (Park et al. [108]), which improves Jupiter’s orbit with ranging and VLBI observations of the Juno spacecraft. Uncertainties in the Jovian orbit impacted NANOGrav’s 11-year GWB search (Arzoumanian et al. [18]; Vallisneri et al. [143]), but they are now negligible.

For each pulsar, we fit the TOAs to a timing model that includes pulsar spin period, spin period derivative, sky location, proper motion, and parallax. While not all pulsars have measurable parallax and proper motion, we always include these parameters because they induce delays with the same frequencies for all pulsars ($f = 0.5\,\mathrm{yr}^{-1}$ for parallax and $f = \mathrm{yr}^{-1}$ plus a linear envelope for proper motion), so there is a risk that a parallax or proper motion signal could be misidentified as a GW signal. Fitting for these parameters in all pulsars reduces our sensitivity to GWs at those frequencies; however, this effect is minimal for GWB searches since these frequencies are much higher than the frequencies at which we expect the GWB to be significant. For binary pulsars, the timing model includes also five orbital elements for binary pulsars and additional non-Keplerian parameters when these improve the fit as determined by an $F$ test. We fit variations in dispersion measure with a piecewise constant “DMX” function (Arzoumanian et al. [16]; Jones et al. [81]). The individual analysis of each pulsar provides best-fit estimates of the timing residuals $\delta t$, of white measurement noise, and of intrinsic red noise, modeled as a power law (Cordes [34]; Lam et al. [88]; Jones et al. [81]).⁹ White measurement noise is described by three parameters: a linear scaling of TOA uncertainties (“EFAC”), white noise added to the TOA uncertainties in quadrature (“EQUAD”), and noise common to all sub-bands at the same epoch (“ECORR”), with independent parameters for every receiver/backend combination (see NG15detchar). We summarize white noise by its maximum a posteriori (MAP) covariance $C$. See App. A for more details of our instruments, observations, and datareduction pipeline: a complete discussion of the data set can be found in NG15.

In our Bayesian GWB analysis, we model $\boldsymbol{\delta t}$ as a finite Gaussian process consisting of time-correlated fluctuations that include intrinsic red pulsar noise and (potentially) a GW signal, along with timing-model uncertainties [144, 145, 134]. The red noise is modeled with Fourier basis $\mathbf{F}$ and amplitudes $\mathbf{c}$ [92]. All Fourier bases (the columns of $\mathbf{F}$) are sines and cosines computed on the TOAs with frequencies $f_i=i/T$, where $T=16.03$ yr is the TOA extent. The timing-model uncertainties are modeled with design-matrix basis $\mathbf{M}$ and coefficients $\boldsymbol{\epsilon}$. The single-pulsar log likelihood is then

with

The prior for the $\boldsymbol{\epsilon}$ is taken to be uniform with infinite extent, so the posterior is driven entirely by the likelihood. The set of the $\{\mathbf{c}\}$ for all pulsars take a joint normal prior with zero mean and covariance

here $a,b$ range over pulsars and $i,j$ over Fourier components; $\delta_{ij}$ is Kronecker’s delta. The term $\varphi_{ai}$ describes the spectrum of intrinsic red noise in pulsar $a$, while $\Phi_{ab,i}$ describes processes with common spectrum across all pulsars and (potentially) phase-coherent inter-pulsar correlations. The $\{\mathbf{c}\}$ prior ties together the single-pulsar likelihoods (Equation (1)) into a joint posterior, $p(\mathbf{c},\boldsymbol{\epsilon},\boldsymbol{\eta}|\boldsymbol{\delta t})\propto p(\boldsymbol{\delta t}|\mathbf{c},\boldsymbol{\epsilon})p(\mathbf{c},\boldsymbol{\epsilon}|\boldsymbol{\eta})p(\boldsymbol{\eta})$, where we have dropped subscripts to denote the concatenation of vectors for all pulsars, and where $\boldsymbol{\eta}$ denotes all the hyperparameters (such as red-noise and GWB power spectrum amplitudes) that determine the covariances. We marginalize over $\mathbf{c}$ and $\boldsymbol{\epsilon}$ analytically, and use Markov chain Monte Carlo techniques (see App. B) to estimate $p(\boldsymbol{\eta}|\boldsymbol{\delta t})$ for different models of the intrinsic red noise and common spectrum.

The data-model variants adopted in this paper all share this probabilistic setup, but differ in the structure and parametrization of $\Phi_{ab,i}$. For a model with intrinsic red noise only (henceforth IRN), $\Phi_{ab,i}=0$; for common-spectrum spatially-uncorrelated red noise (CURN), $\Phi_{ab,i}=\delta_{ab}\Phi_{\mathrm{CURN},i}$; for an isotropic GWB with Hellings–Downs correlations (HD), $\Phi_{ab,i}=\Gamma(\xi_{ab})\Phi_{\mathrm{HD},i}$, with $\Gamma$ the Hellings–Downs function of pulsar angular separations $\xi_{ab}$

In NG12gwb we established strong Bayesian evidence for CURN over IRN; finding that HD is preferred over CURN would point to the GWB origin of the common-spectrum signal. We also investigate other spatial correlation patterns, e.g., monopole or dipole, introduced in §5.3.

Throughout this paper, we set the spectral components $\varphi_{ai}$ of intrinsic pulsar noise (which have units of $\mathrm{s}^{2}$, as appropriate for the variance of timing residuals) to a power law,

introducing two dimensionless hyperparameters for each pulsar: the intrinsic-noise amplitude $A_a$ and spectral index $\gamma_a$. We use log-uniform and uniform priors, respectively, on these hyperparameters; their bounds are described in App. B. More sophisticated intrinsic-noise models are discussed in §5.1 and NG15detchar. In models CURN$^\gamma$ and HD$^\gamma$, the common spectra $\Phi_{\mathrm{CURN},i}$ and $\Phi_{\mathrm{HD},i}$ follow Equation (6),

introducing hyperparameters $A_{\mathrm{CURN}},\gamma_{\mathrm{CURN}}$ and $A_{\mathrm{HD}},\gamma_{\mathrm{HD}}$ respectively. However, we set $\gamma_{\mathrm{HD}}=13/3$ for the GWB from a stationary ensemble of inspiraling binaries, and refer to that fiducial model as HD$^{13/3}$. For specific “free spectrum” studies we will instead model the individual $\Phi_{\mathrm{CURN},i}$ or $\Phi_{\mathrm{HD},i}$ elements, and refer to models CURN$^{\mathrm{free}}$ and HD$^{\mathrm{free}}$. Throughout this article we use frequencies $f_i=i/T$ with $i=1$–$30$ for intrinsic noise ($f=2$–$59$ nHz), covering a frequency range over which pulsar noise transitions from red-noise-dominated to white-noise-dominated. For common-spectrum noise, we limit the frequency range in order to reduce correlations with excess white noise at higher frequencies. Following NG12gwb, we fit a CURN$^\gamma$ model enhanced with a power-law break to our data, and limit frequencies to the MAP break frequencies ($i=1$–$14$ or $f=2$–$28$ nHz; see App. C).

Bayesian analysis

When fit to the 15-year data set, the CURN$^\gamma$ and HD$^\gamma$ models agree on the presence of a loud time-correlated stochastic signal with common amplitude and spectrum across pulsars¹⁰. The joint $A_{\mathrm{HD}}$–$\gamma_{\mathrm{HD}}$ Bayesian posterior is shown in panel (b) of Figure 1, with 1-D marginal posteriors in the horizontal and vertical sub-plots. The posterior medians and 5–95% quantiles are $A_{\mathrm{HD}} = 6.4^{+4.2}_{-2.7} \times10^{-15}$ and $\gamma_{\mathrm{HD}} = 3.2^{+0.6}_{-0.6}$. The thicker curve in the vertical subplot is the $A_{\mathrm{HD}}$ posterior for the $\mathrm{HD}^{13/3}$ model, for which $A_{\mathrm{HD},13/3} = 2.4^{+0.7}_{-0.6} \times10^{-15}$. These amplitudes are compatible with astrophysical expectations of a GWB from inspiraling SMBHBs (see § 6). The $A_{\mathrm{HD}}$ posterior has essentially no support below $10^{-15}$.

The strong $A_{\mathrm{HD}}$–$\gamma_{\mathrm{HD}}$ correlation is an artifact of using the conventional frequency $f_{\mathrm{ref}} = 1\,\mathrm{yr}^{-1}$ in Equation (6), and it largely disappears when $f_{\mathrm{ref}}$ is moved to the band of greatest PTA sensitivity; see the dashed contours in panel (b) of Figure 1 for $f_{\mathrm{ref}} = (10\,\mathrm{yr})^{-1}$. The $\gamma_{\mathrm{HD}}$ posterior is in moderate tension with the theoretical universal binary-inspiral value $\gamma_{\mathrm{HD}} = 13/3$, which lies at the 99% credible boundary: smaller values of $\gamma_{\mathrm{HD}}$ could be an indication that astrophysical effects, such as stellar scattering and gas dynamics, play a role in the evolution of SMBHBs emitting GWs in this frequency range (see § 6 and Agazie et al. [3]). This highlights the importance of measuring this parameter. Furthermore, its estimation is sensitive to details in the modeling of intrinsic red noise and of interstellar-medium timing delays in a few pulsars (see the analysis in § 5.2). Notably, in the 12.5-year data set $\gamma_{\mathrm{HD}} = 13/3$ was recovered at $\sim1\sigma$ below the median (NG12gwb); this anomaly is reversed in the newer data set. It is likely that more expansive data sets or more sophisticated chromatic noise models, e.g., next generation Gaussian process models such as those in § 5.1 (Goncharov et al. [59]; Chalumeau et al. [31]; Lamb et al. 2018), will be needed to infer the presence of possible systematic errors in $\gamma_{\mathrm{HD}}$.

Our Bayesian analysis provides evidence that the common-spectrum signal includes Hellings–Downs inter-pulsar correlations. Specifically, the Bayes factor between the $\mathrm{HD}^{\gamma}$ and $\mathrm{CURN}^{\gamma}$ models ranges from 200 (when 14 Fourier frequencies are included in $\Phi_i$) to 1,000 (when 5 frequencies are included, as in NG12gwb). Results are similar for $\mathrm{HD}^{13/3}$ vs. $\mathrm{CURN}^{13/3}$. Figure 2 recapitulates Bayes factors between a variety of models, including some with the alternative spatial-correlation structures discussed in § 5.3. The very peaked $A_{\mathrm{HD}}$ posterior in panel (b) of Figure 1, significantly separated from smaller amplitudes, supports the very large Bayes factor between IRN and $\mathrm{CURN}^{\gamma}$. The 15-year data set favors $\mathrm{HD}^{\gamma}$ over $\mathrm{CURN}^{\gamma}$, and over models with monopolar or dipolar correlations, and it is inconclusive about, i.e., gives roughly even odds for, the presence of spatially correlated signals in addition to $\mathrm{HD}^{\gamma}$.

We can also regard the $\mathrm{HD}^{\gamma}$ vs. $\mathrm{CURN}^{\gamma}$ Bayes factor as a detection statistic in a hypothesis-testing framework, and derive the $p$-value of the observed Bayes factor with respect to its empirical distribution under the $\mathrm{CURN}^{\gamma}$ model. We do so by computing Bayes factors on 5,000 bootstrapped data sets where inter-pulsar spatial correlations are removed by introducing random phase shifts, drawn from a uniform distribution from 0 to $2\pi$, to the common-process Fourier components (Taylor et al. [137]). This procedure alters inter-pulsar correlations to have a mean of zero, while leaving the amplitudes of intrinsic pulsar noise and $\mathrm{CURN}$ unchanged, thus providing a way to test the null hypothesis that no inter-pulsar correlations are present. The resulting background distribution of Bayes factors is shown in the left panel of Figure 3—they exceed the observed value in five of the 5,000 phase shifts ($p = 10^{-3}$). We also performed sky scramble analyses (Cornish & Sampson [36]), which remove the dependence of inter-pulsar spatial correlations on the angular separations between the pulsars by attributing random sky positions to the pulsars. Sky scrambles generate a background distribution for which inter-pulsar correlations are present in the data but they are independent of the pulsars’ angular separations: for this distribution, we find $p = 1.6 \times10^{-3}$. A detailed discussion of sky scrambles and the results of these analyses can be found in App. F.

As in NG12gwb, we also carried out a minimally modeled Bayesian reconstruction of the inter-pulsar correlation pattern, using spline interpolation over seven spline-knot positions. The choice of seven spline-knot positions is based on features of the Hellings–Downs pattern: two correspond to the maximum and minimum angular separations ($0^\circ$ and $180^\circ$, respectively), two are chosen to be at the theoretical zero crossings of the Hellings–Downs pattern ($49.2^\circ$ and $121.8^\circ$), one is at the theoretical minimum ($82.5^\circ$), and the final two are between the end points and zero crossings ($25^\circ$ and $150^\circ$) to allow additional flexibility in the fit. Panel (d) of Figure 1 shows the marginal 1-D posterior densities at these spline-knot positions for a power-law varied-exponent model. The reconstruction is consistent with the over-plotted Hellings–Downs pattern; furthermore, the joint 2-D marginal posterior densities for the knots, not shown in panel (d) of Figure 1, at the HD zero-crossings is consistent with $(0,0)$ within $1\sigma$ credibility.

Optimal statistic analysis

We complement our Bayesian search with a frequentist analysis using the optimal statistic (Anholm et al. [12]; Demorest et al. [39]; Chamberlin et al. [32]), a summary statistic designed to measure correlated excess power in PTA residuals. (Note that there is no accepted definition of “optimal statistic” in modern statistical usage, but the term has become established in the PTA literature to refer to this specific method, so we use it for this reason.) It is enlightening to describe the optimal statistic as a weighted average of the inter-pulsar correlation coefficients

where $\delta t_a^{\mathsf{T}}$ are the residuals of pulsar $a$, and $\mathbf{P}_a = \langle\delta t_a\delta t_a^{\mathsf{T}}\rangle$ is their total auto-covariance matrix. The cross-covariance matrix $\widetilde{\boldsymbol{\Phi}}_{ab}$ encodes the spectrum of the HD-correlated signal, normalized so that $\boldsymbol{\Phi}_{ab} = A^2\Gamma(\xi_{ab})\widetilde{\boldsymbol{\Phi}}_{ab}$ (see Pol et al. [114]), and where elements of $\boldsymbol{\Phi}_{ab}$ are given by Equation (3). Indeed, the $\rho_{ab}$ have expectation value $A^2\Gamma(\xi_{ab})$, but their variance $\sigma_{ab}^2 = (\operatorname{Tr}\mathbf{P}_a^{-1}\widetilde{\boldsymbol{\Phi}}_{ab}\mathbf{P}_b^{-1}\widetilde{\boldsymbol{\Phi}}_{ba})^{-1}+O(A^4)$ is too large to use them directly as estimators. Thus we assemble the optimal statistic as the variance-weighted, $\Gamma$-template-matched average of the $\rho_{ab}$,

This equation represents the optimal estimator of the HD amplitude $A^2$; it can also be interpreted as the best-fit $\hat{A}^2$ obtained by least-squares–fitting the $\rho_{ab}$ to the Hellings–Downs model $\hat{A}^2\Gamma(\xi_{ab})$. Because $\hat{A}^2$ is a function of intrinsic–red-noise and common-process hyperparameters through the $P_a$, we use the results of an initial Bayesian-inference run to refer the statistic to MAP hyperparameters, or to marginalize it over their posteriors. As discussed in Vigeland et al. [147], we obtain more accurate values of the amplitude by this marginalization.

To search for inter-pulsar correlations using the optimal statistic, we evaluate the frequency (the $p$-value) with which an uncorrelated common-spectrum process with parameters estimated from our data set would yield $\hat{A}^2$ greater than we observe. In the absence of a signal, the expectation value of $\hat{A}^2$ is zero, and its distribution is approximately normal. Thus we divide the observed $\hat{A}^2$ by its standard deviation to define a formal signal-to-noise ratio

Figure 4 shows the distribution of this S/N over CURN$^{\gamma}$ and CURN$^{13/3}$ noise-parameter posteriors, with S/Ns of $5 \pm1$ and $4 \pm1$, respectively (means $\pm$ standard deviations across noise-parameter posteriors). We use 14 frequency components to model the signal: the dependence on the number of frequency components is very weak.

Because the distribution of $\hat{A}^2$ is only approximately normal (Hazboun et al. [63]), the S/N of Equation (11) does not map analytically to a $p$-value, and it cannot be interpreted as a “sigma” level. Instead, optimal-statistic $p$-values can be computed empirically by removing inter-pulsar correlations from the 15-year data set with phase shifts (Taylor et al. [137]). We draw random phase offsets from 0 to $2\pi$ for the common-process Fourier components, which is equivalent to making uniform draws from the background distribution of CURN, and ask how often a random choice of phase offsets produces a HD-correlated signal. The right panel of Figure 3 shows the distribution of noise-marginalized S/N over 400,000 phase shifts. There are 19 phase shifts with noise-marginalized S/N greater than observed, with $p = 5 \times10^{-5}$. We compare the phase-shift distribution with backgrounds obtained by simulation (right panel of Figure 3, orange line) and analytic calculation (green line). For the former, we simulate 27,000 CURN$^{\gamma}$ realizations using MAP hyperparameters from the 15-yr data and compute the optimal-statistic S/N for each; for the latter, we evaluate the generalized $\chi^2$ distribution (Hazboun et al. [63]) with median CURN$^{\gamma}$ hyperparameters. Although neither method includes the marginalization over noise-parameter posteriors, we find good agreement with phase shifts, with $p = 1.8 \times10^{-4}$ from simulations, and $p = 1.9 \times10^{-4}$ from the analytic calculation. Finally, we use sky scrambles to compute the $p$-value for the null hypothesis that inter-pulsar correlations are present, but they have no dependence on the angular separation between the pulsars, for which we find $p < 10^{-4}$ (see App. § 14).

Averaging the cross-correlations $\rho_{ab}$ in angular-separation bins with equal numbers of pulsar pairs reveals the Hellings–Downs pattern, as shown in panel (c) of Figure 1 for 15 bins. The $\rho_{ab}$ were evaluated with MAP CURN$^{13/3}$ noise parameters. The black dashed curve traces the expected correlations for an HD-correlated background with the MAP amplitude; the vertical error bars display the expected $1\sigma$ spreads of the binned cross-correlations, accounting for the $\langle\rho_{ab}\rho_{cd}\rangle$ covariances induced by the HD-correlated process (Romano et al. [122]; Allen & Romano [11]). (Neglecting those covariances yields 20–40% smaller spreads. Note that they are not included in $p$-value estimates because those are calculated under the null hypothesis of no spatially correlated process.)

Although each draw from the noise-parameter posterior would generate a slightly different plot, as would different binnings, the quality of the fit seen in Figure 1 provides a visual indication that the excess low-frequency power in the 15-year data set harbors HD correlations. The $\chi^2$ for this 15-bin reconstruction with respect to the Hellings–Downs curve is 8.1, where we account for $\rho_{ab}$ covariance in constructing the bins, and the covariance between bins in constructing the $\chi^2$ (Allen & Romano [11]). This corresponds to a $p$-value of 0.75, calculated using simulations based on the HD$\gamma$ model, or 0.92 if one assumes this value follows a canonical $\chi^2$ with 15 degrees of freedom. These $p$-values are representative of what we find with different binnings: we find $p > 0.3$ when using eight to 20 bins (assuming a canonical $\chi^2$ distribution).

Checks and validation

Prior to analyzing the 15-year data set, we extensively reviewed our data collection and analysis procedures, methods, and tools, in an effort to eliminate contamination from systematic effects and human error. Furthermore, the results presented in § 3 and § 4 are supported by a variety of consistency checks and auxiliary studies. In this section we present those that offer evidence for or against the presence of HD correlations, reveal anomalies, or otherwise highlight features of note in the data: alternative DM modeling (§ 5.1), the spectral content (§ 5.2) and correlation pattern (§ 5.3) of the excess-noise signal, as well as the consistency of our findings across data set “slices,” pulsars, and telescopes (§ 5.4).

Alternative DM models

In this paper and in previous GW searches (e.g., NG12gwb), we model fluctuations in the DM using DMX parameters (a piecewise-constant representation, see NG15). Adopting this DM model as the standard makes it easier to directly compare the results here to those in NG12gwb. An alternative model where DM variations are modeled as a Fourier-domain Gaussian process, DMGP, has been used by Antoniadis et al. [13], Chen et al. [33], and Goncharov et al. [58]. The Fourier coefficients follow a power law similar to those of intrinsic and common-spectrum red noise, but their basis vectors include a $\nu^{-2}$ radio-frequency dependence, and the component frequencies $f_i = i/T$ range through $i = 1$–100. Under the DMGP model we also include a deterministic solar-wind model (Hazboun et al. [67]) and the two chromatic events in PSR J1713+0747 reported in Lam et al. [89] which are modeled as deterministic exponential dips with the chromatic index quantifying the radio-frequency dependence of the dips left as a free parameter. If these chromatic events are not modeled, they raise estimated white noise (Hazboun et al. [66]). A detailed discussion of chromatic noise effects can be found in NG15detchar.

Using the DMGP model in place of DMX has minimal effects on nearly all pulsars in the array. Only PSRs J1713+0747 and J1600−3053 show notable differences in their recovered intrinsic-noise parameters. However, DMGP does affect the parameter estimation of common red noise, as seen in Figure 5, shifting the posterior for

$\gamma$ to higher values that are more consistent with 13/3. Despite this, we still recover HD correlations at the same significance as when we use DMX to model fluctuations in the DM, implying that the evidence reported for the presence of correlations in this work is independent of the choice of DM noise modeling.

Spectral analysis

Adopting power-law spectra for CURN and HD is a useful simplification that reduces the number of fit parameters and yields more informative constraints; furthermore, it is expedient to identify HD$^{13/3}$ with the hypothesis that we are observing the GWB from SMBHBs. Nevertheless, the standard $\gamma= 13/3$ power law for GW inspirals may be altered by astrophysical processes such as stellar and gas friction in nuclei (see, e.g., Merritt & Milosavljević [99] for a review), by appreciable eccentricity in SMBHB orbits (Enoki & Nagashima [48]), and by low-number SMBHB statistics (Sesana et al. [129]). HD$\gamma$ parameter recovery may also be biased if intrinsic pulsar noise is not modeled well by a power law. Indeed, our data show hints of a discrepancy from the idealized HD$^{13/3}$ model: the $\gamma_{\mathrm{HD}}$ posterior in panel (b) of Figure 1 favors slopes much shallower than 13/3, and the HD$\gamma$-to-CURN$\gamma$ Bayes factor drops from 1,000 to 200 when Fourier components at more than five frequencies are included in the model.

We examine the spectral content of the 15-year data set using the CURN$^{\mathrm{free}}$ and HD$^{\mathrm{free}}$ models, which are parametrized by the variances of the Fourier components at each frequency. Their marginal posteriors are shown in the left panel of Figure 6, where bin number $i$ corresponds to $f_i = i/T$, with $T = 16.03$ yr the extent of the data set. For the purpose of illustration, we overlay best-fit power laws that thread the posteriors in a way similar to the factorized PTA likelihood of Taylor et al. [138] and Lamb et al. [90].

We deem excess power, either uncorrelated for CURN${}^{\mathrm{free}}$ or correlated for HD${}^{\mathrm{free}}$, to be observed in a bin when the support of the posterior is concentrated away from the lowest amplitudes. No power of either kind is observed above $f_8$, consistent with the presence of a floor of white measurement noise. Furthermore, no correlated power is observed in bins 6 and 7, where a power-law model would expect a smooth continuation of the trend of bins 1–5 (cf. the dashed fit of Figure 6): this may explain the drop in the Bayes factor. However, correlated power reappears in bin 8, pushing the fit toward shallower slopes. Indeed, repeating the fit by omitting subsets of the bins suggests that the low recovered $\gamma_{\mathrm{HD}}$ is due mostly to bin 8 and to the lower-than-expected correlated power found in bin 1. Obviously, excluding those bins leads to higher $\gamma_{\mathrm{HD}}$ estimates.

To explore deviations from a pure power law that may arise from statistical fluctuations of the astrophysical background or from unmodeled systematics (perhaps related to the timing model), in App. D we relax the normal $c_k$ prior (cf. Equation (3)) to a multivariate Student’s $t$-distribution that is more accepting of mild outliers. The resulting estimate of $\gamma_{\mathrm{CURN}}$ peaks at a higher value and is broader than in CURN$^\gamma$, with posterior medians and 5–95% quantiles of $\gamma_{\mathrm{CURN}}=3.5^{+1.0}_{-1.0}$.

Similarly, spectral turnovers due to interactions between SMBHBs and their environments can result in reduced GWB power at lower frequencies, which might explain the slightly lower correlated power in bin 1. We investigate this hypothesis in App. E using the turnover spectrum of Sampson et al. [124]. For this CURN${}^{\mathrm{turnover}}$ model, the 15-year data favor a spectral bend below 10 nHz (near $f_5$), but the Bayes factor against the standard HD$\gamma$ is inconclusive.

Future data sets with longer time spans and the comparison of our data set with those of other PTAs should help clarify the astrophysical or systematic origin of these possible spectral features.

Alternative correlation patterns

Sources other than GWs can produce inter-pulsar residual correlations with spatial patterns other than HD. For example, errors in the solar-system ephemerides create time-dependent Roemer delays with dipolar correlations (Roebber [119]; Vallisneri et al. [143]), and errors in the correction of telescope time to an inertial timescale (Hobbs et al. [73, 74]) create an identical time-dependent delay for all pulsars (i.e., a delay with monopolar correlations).

Gair et al. [52] showed that, for a pulsar array distributed uniformly across the sky, HD correlations can be decomposed as

where the $P_\ell(\cos\xi_{ab})$ are Legendre polynomials of order $\ell$ evaluated at the pulsar angular separation $\xi_{ab}$. In other words, a HD-correlated signal should have no power at $\ell=0$ or $\ell=1$.

We can perform a frequentist generic correlation search using Legendre polynomials¹¹ with the multiple-component optimal statistic (MCOS; Sardesai & Vigeland [125])—a generalized statistic that allows multiple correlation patterns to be fit simultaneously to the correlation coefficients $\rho_{ab}$. Figure 7 shows the constraints on $A_l^2=A^2g_l$ obtained by fitting the correlations $\rho_{ab}$ to this Legendre series using the MCOS and marginalizing over CURN$^\gamma$ noise-parameter posteriors. The quadrupolar structure of the data is evident, along with a small but significant monopolar contribution.

The same feature from the Legendre decomposition appears if we use the MCOS to search for multiple correlations simultaneously: a multiple regression analysis favors models that contain both significant HD and monopole correlations (see Appendix § 15). From simulations of 15-year-like data sets (see Appendix § 16.§ 1), we find a $p$-value of 0.11 (approx. $2\sigma$) for observing a monopole at this significance or higher with a pure-HD injection of amplitude similar to what we observe. We also perform a model-checking study to assess whether the observed monopole is consistent with the $\mathrm{HD}^{13/3}$ model (see Appendix § 16.§ 2), and we find a $p$-value of 0.11 for producing an apparent monopole when the signal is purely $\mathrm{HD}^{13/3}$. Thus, we conclude that it is possible for a HD-correlated signal to appear to have monopole correlations in an optimal statistic analysis at this significance level.

In contrast, Bayesian searches for additional correlations do not find evidence of additional monopole- or dipole-correlated red noise processes: as shown in Figure 2, the Bayes factors for these processes are $\sim1$. We also perform a general Bayesian search for correlations using a CURN$^{\mathrm{free}}$ + HD$^{\mathrm{free}}$ + MONOPOLE$^{\mathrm{free}}$ + DIPOLE$^{\mathrm{free}}$ model, which allows for independent uncorrelated and correlated components at every frequency bin. We note that this analysis is more flexible than the ones described above, which assume a power-law power spectral density. We find no significant dipole-correlated power at any frequency, and we find monopole-correlated power only in the second frequency bin ($f_2=3.95$ nHz); posteriors of variance for that bin are shown in the right panel of Figure 6.

Motivated by this finding, we perform a search for HD$^\gamma$ + SINUSOID, which includes a deterministic sinusoidal delay (applied to all pulsars alike, as appropriate for a monopole) with free frequency, amplitude, and phase. The sinusoid’s posteriors match the free-spectral analysis in frequency and amplitude; however, the Bayes factor between HD$^\gamma$ + SINUSOID and HD$^\gamma$ calculated using two methods (Hee et al. [69]; Hourihane et al. [75]), is only $\sim1$, so the signal cannot be considered statistically significant. Astrophysically motivated searches for sources that produce sinusoidal or sinusoid-like delays in the residuals, such as an individual SMBHB or perturbations to the local gravitational field induced by fuzzy dark matter (Khmelnitsky & Rubakov [86]), also yield Bayes factors $\sim1$. Thus we conclude that there is some evidence of additional power at 3.95 nHz with monopole correlations; however, the significance in the Bayesian analyses is low, while the optimal-statistic S/N could be produced by a HD-correlated signal. Therefore, we cannot definitively say whether the signal is present, or determine the source. We note that performing an MCOS analysis after subtracting off realizations of a sinusoid using HD$^\gamma$ + SINUSOID posteriors reduces the $(\mathrm{S/N})_{\mathrm{monopole}}$ to $\sim0$, while $(\mathrm{S/N})_{\mathrm{HD}}$ remains unchanged, indicating that this single-frequency monopole-correlated signal is likely causing the nonzero monopole signal observed in the MCOS analysis.

Similar hints of a monopolar signal (though weaker) were found in the NANOGrav 12.5-year data set, unsurprisingly given that it is a subset of the current data set. To exercise due diligence, we audited the correction of telescope time to GPS time at the Arecibo Observatory and at the Green Bank Telescope, and found nothing that could explain our observations. The subsequent steps in the time-correction pipeline rely on very accurate atomic clocks and are unlikely to introduce considerable systematics (Petit [110]). An important test will be whether this signal persists in future data sets. If this monopolar feature is a truly an astrophysical signal, we would expect it to increase in significance as our data set grows. Comparisons with other PTAs and combined IPTA data sets will also provide crucial insight.

Dropout and cross-validation

The GWB is by its nature a signal affecting all of the pulsars in the PTA, although it may appear more significant in some based on their observing time span, noise properties, and on the particular realization of pulsar and Earth contributions ([131]). One way to assess the significance of the GWB in each pulsar is a Bayesian dropout analysis ([6, 19]), which introduces a binary parameter that turns on and off the common signal (or its inter-pulsar correlations) for a single pulsar, leaving all other pulsars unchanged. The Bayes factor associated with this parameter, also referred to as the “dropout factor,” describes how much each pulsar likes to “participate” in the common signal.

Figure 8 plots $\mathrm{CURN}^{\gamma}$ vs. $\mathrm{IRN}$ dropout factors for all 67 pulsars (blue). We find positive dropout factors (i.e., dropout factors $> 2$) for an uncorrelated common process in twenty pulsars, while only one has a dropout factor $< 0.5$. For comparison, in the NANOGrav 12.5-year data set ten pulsars showed positive dropout factors for an uncorrelated common process, while three had negative dropout factors. We also show HD correlations vs. $\mathrm{CURN}^{\gamma}$ dropout factors (orange). For these, the uncorrelated common process is always present in all pulsars, but the cross-correlations for all pulsar pairs involving a given pulsar may be dropped from the likelihood. We find positive factors for HD correlations vs. $\mathrm{CURN}^{\gamma}$ in seven pulsars, while three are negative. We expect more pulsars to have positive dropout factors for $\mathrm{CURN}^{\gamma}$ vs. $\mathrm{IRN}$ than for Hellings–Downs vs. $\mathrm{CURN}^{\gamma}$ because the Bayes factor comparing the first two models is significantly higher than the one comparing the second two models (see Figure 2). Negative dropout factors could be caused by noise fluctuations or they could be an indication that more advanced chromatic noise modeling is necessary ([8]). They could also be caused by the GWB itself, which induces both correlated and uncorrelated noise in the pulsars (the so-called “Earth terms” and “pulsar terms”; [102]).

In addition to Bayes factors, the goodness-of-fit of probabilistic models can be evaluated by assessing their predictive performance ([53]). Specifically, given that the GWB is correlated across pulsars, we can (partially) predict the timing residuals $\boldsymbol{\delta t}_a$ of

pulsar $a$ from the residuals $\boldsymbol{\delta t}_{-a}$ of all other pulsars by way of the “leave-one-out” posterior predictive likelihood (PPL)

where $\boldsymbol{\theta}_a$ are all the parameters and hyperparameters that affect pulsar $a$ in a given model. As discussed in [100], we compare the predictive performance of $\mathrm{CURN}^{13/3}$ and $\mathrm{HD}^{13/3}$ for each pulsar in turn by taking the ratio of the corresponding leave-one-out PPLs. These ratios are closely related to the dropout factors plotted in Figure 8. Multiplying the PPL ratios for all pulsars yields the pseudo Bayes factor (PBF). For the 15-year data set we find $\mathrm{PBF}_{15\mathrm{yr}} = 1{,}400$ in favor of $\mathrm{HD}^{13/3}$ over $\mathrm{CURN}^{13/3}$. The PBF does not have a “betting odds” interpretation, but we obtain a crude estimate of its significance by building its background distribution on 40 $\mathrm{CURN}^{13/3}$ simulations with the MAP $\log_{10} A_{\mathrm{CURN}}$ inferred from the 15-year data set. For all simulations except one, the PBF favors the null hypothesis, and $\log_{10}\mathrm{PBF}_{15\mathrm{yr}}$ is displaced by approx. three standard deviations from the mean $\log_{10}\mathrm{PBF}$.

A different sort of cross-validation relies on evaluating the optimal statistic for temporal subsets of the data set, as in [66]. In the regime where the lowest frequencies of our data are dominated by the GWB, the optimal statistic S/N should grow with the square root of the time span of the data and linearly with the number of pulsars in the array ([130]); in this regime increasing the number of pulsars is the best way to boost PTA sensitivity to the GWB. To verify that this is indeed the case, we analyze “slices” of the data set in six-month increments, starting from a six-year data set. Once a new pulsar accumulates three years of data, we add it to the array. We perform a separate Bayesian CURN$^\gamma$ analysis for each slice and calculate the Hellings–Downs optimal statistic over the noise-parameter posterior. In Figure 9, we plot the S/N distributions against time span and the number of pulsars. As expected, we observe essentially monotonic growth associated with the increase in the number of pulsars.

The signal should also be consistent between timing observations made with Arecibo and GBT. To test this, we analyze the two split-telescope data sets (see App. § 9); both show evidence of common-spectrum excess noise. Figure 10 shows Arecibo (orange) and GBT (green) CURN$^\gamma$ posteriors, which are broadly consistent with each other and with full-data posteriors (blue). Arecibo yields $\log_{10} A = -14.02^{+0.18}_{-0.22}$ and $\gamma= 2.78^{+0.70}_{-0.64}$ (medians with 68% credible intervals), while GBT yields $\log_{10} A = -14.2^{+0.15}_{-0.17}$ and $\gamma= 3.37^{+0.40}_{-0.38}$.

The split-telescope data sets are significantly less sensitive to spatial correlations than the full data set, because they have fewer pulsars and therefore pulsar pairs (see Figure 12 of App. § 9). Nevertheless, we can search them for spatial correlations using the optimal statistic. We find a noise-marginalized Hellings–Downs S/N of 2.9 for Arecibo and 3.3 for GBT, consistent with the split-telescope data sets having about half the number of pulsars as the full data set. The S/Ns for Arecibo and GBT are comparable: while telescope sensitivity, observing cadence, and distribution of pulsars all affect GWB sensitivity, the dominant factor is the number of pulsars because the S/N scales linearly with the number of pulsars but only as $\propto(\sigma\sqrt{c})^{-1/\gamma}$, where $\sigma$ is the residual root-mean-squared, and $c$ is the observing cadence ([130]). We also note that the distributions of angular separations probed by Arecibo and GBT are similar, although GBT observes more pulsar pairs with large angular separations (see Figure 12).

Discussion

In this letter we have reported on a search for an isotropic stochastic GWB in the 15-year NANOGrav data set. A previous analysis of the 12.5-year NANOGrav data set found strong evidence for excess low-frequency noise with common spectral properties across the array, but inconclusive evidence for Hellings–Downs inter-pulsar correlations, which would point to the GW origin of the background. By contrast, the 12.5-year data disfavored purely monopolar (clock-error-like) and dipolar (ephemeris-error-like) correlations. Subsequent independent analyses by the PPTA and EPTA collaborations reported results consistent with ours ([58, 33]), as did the search of a combined data set ([13])—a syzygy of tantalizing discoveries that portend the rise of low-frequency GW astronomy.

We analyzed timing data for 67 pulsars in the 15-year data set (those that span $> 3$ years), with a total time span of 16.03 years, and more than twice the pulsar pairs than in the 12.5-year data set. The common-spectrum stochastic signal gains even greater significance and is detected in a larger number of pulsars. For the first time, we find compelling evidence of Hellings–Downs inter-pulsar correlations, using both Bayesian and frequentist detection statistics (see Figure 1), with false-alarm probabilities of $p = 10^{-3}$ and $p = 5 \times10^{-5}$–$1.9 \times10^{-4}$, respectively (see Figure 3).

The significance of Hellings–Downs correlations increases as we increase the number of frequency components in the analysis up to five, indicating that the correlated signal extends over a range of frequencies. A detailed spectral analysis supports a power-law signal, but at least two frequency bins show deviations that may skew the determination of spectral slope (Figure 6). These discrepancies may arise from astrophysical or systematic effects. Furthermore, slope determination changes significantly using an alternative DM model (Figure 5). The study of spatial correlations with the optimal statistic confirms a Hellings–Downs quasi-quadrupolar pattern (Figure 7 and panel c of Figure 1), with some indications of an additional monopolar signal confined to a narrow frequency range near 4 nHz. However, the Bayesian evidence for this monopolar signal is inconclusive, and we could not ascribe it to any astrophysical or terrestrial source (e.g., an individual SMBHB or errors in the chain of timing corrections).

The GWB is a persistent signal that should increase in significance with number of pulsars and observing time span. This is indeed what we observe by analyzing slices of the data set (see Figure 9). Furthermore, the signal is present in multiple pulsars (Figure 8), and can be found in independent single-telescope data sets (Figure 10). We are preparing a number of other papers searching the 15-year data set for stochastic and deterministic signals, including an all-sky, all-frequency search for GWs from individual circular SMBHBs. This search, together with the same analysis of the 12.5-year data set [21], indicates that the spectrum and correlations we observe cannot be produced by an individual circular SMBHB.

If the Hellings–Downs-correlated signal is indeed an astrophysical GWB, its origin remains indeterminate. Among the many possible sources in the PTA frequency band, numerous studies have focused on the unresolved background from a population of close-separation SMBHBs. The SMBHB population is a direct byproduct of hierarchical structure formation, which is driven by galaxy mergers [e.g. 27]. In a post-merger galaxy, the SMBHs sink to the center of the common merger remnant through dynamical interactions with their astrophysical environment, eventually leading to the formation of a binary [26]. GW emission from a SMBHB at nHz frequencies is quasi-monochromatic because the binaries evolve very slowly. Under the assumption of purely GW-driven binary evolution, the expected characteristic-strain spectrum is $\propto f^{-2/3}$ (or $f^{-13/3}$ for pulsar-timing residuals).

The GWB spectrum may also feature a low-frequency turnover induced by the dynamical interactions of binaries with their astrophysical environment (e.g., with stars or gas, see [15, 128, 99]) or possibly by non-negligible orbital eccentricities persisting to small separations [48]. We find little support for a low-frequency turnover in our data (see App. E).

The GWB amplitude is determined primarily by SMBH masses and by the occurrence rate of close binaries, which in turn depends on the galaxy merger rate, the occupation fraction of SMBHs, and the binary evolution timescale; population models predict amplitudes ranging over more than an order of magnitude [116, 150, 78, 98, 127], under a variety of assumptions. Figure 11 displays a comparison of HD$^\gamma$ parameter posteriors with power-law spectral fits from an observationally constrained semi-analytic model of the SMBHB population constructed with the HOLODECK package [85]. This particular set of SMBHB populations assumes purely GW-driven binary evolution, and uses relatively narrow distributions of model parameters based on literature constraints from galaxy-merger observations (see, e.g., [141]). While the amplitude recovered in our analysis is consistent with models derived directly from our understanding of SMBH and galaxy evolution, it is toward the upper end of predictions implying a combination of relatively high SMBH masses and binary fractions. A detailed discussion of the GWB from SMBHBs in light of our results is given in [3].

In addition to SMBHBs, more exotic cosmological sources such as inflation, cosmic strings, phase transitions, domain walls, and curvature-induced GWs can also produce detectable GWs in the nHz range (see, e.g., [61, 29], and references therein). Similarities in the spectral shapes of cosmological and astrophysical signals make it challenging to determine the origin of the background from its spectral characterization [83]. The question could be settled by the detection of signals from individual loud SMBHBs or by the observation of spatial anisotropies, since the anisotropy expected from SMBHBs are orders of magnitude larger than those produced by most cosmological sources (Caprini & Figueroa 2018; Bartolo et al. 2022. We discuss these models in the context of our results in Afzal et al. [2].

The EPTA and Indian Pulsar Timing Array (InPTA; Joshi et al. [82]), PPTA, and Chinese Pulsar Timing Array (CPTA; Lee [91]) collaborations have also recently searched their most recent data for signatures of a gravitational-wave background (Antoniadis et al. [14]; Reardon et al. [118]; Xu et al. (2023)), and an upcoming IPTA paper will compare the results of these searches.

The IPTA’s forthcoming Data Release 3 will combine the NANOGrav 15-year data set with observations from the EPTA, PPTA, and InPTA collaborations, comprising about 80 pulsars with time spans up to 24 years, and offering significantly greater sensitivity to spatial correlations and spectral characteristics than single-PTA data sets. Future PTA observation campaigns will improve our understanding of this signal and of its astrophysical and cosmological interpretation. Longer data sets will tighten spectral constraints on the GWB, clarifying its origin (Pol et al. [115]). Greater numbers of pulsars will allow us to probe anisotropy in the GWB (Pol et al. [114]) and its polarization structure (see, e.g., Arzoumanian et al. [20], and references therein). The observation of a stochastic signal with spatial correlations in PTA data, suggesting a GWB origin, expands the horizon of GW astronomy with a new Galaxy-scale observatory sensitive to the most massive black-hole systems in the Universe and to exotic cosmological processes.

Acknowledgments

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. G.A., A.A., A.M.A., Z.A., P.T.B., P.R.B., H.T.C., K.C., M.E.D., P.B.D., T.D., E.C.F., W.F., E.F., G.E.F., N.G., P.A.G., J.G., D.C.G., J.S.H., R.J.J., M.L.J., D.L.K., M.K., M.T.L., D.R.L., J.L., R.S.L., A.M., M.A.M., N.M., B.W.M., C.N., D.J.N., T.T.P., B.B.P., N.S.P., H.A.R., S.M.R., P.S.R., A.S., C.S., B.J.S., I.H.S., K.S., A.S., J.K.S., and H.M.W. developed the 15-year data set through a combination of observations, arrival time calculations, data checks and refinements, and timing model development and analysis; additional specific contributions to the data set are summarized in NG15. S.R.T. and S.J.V. led the search and coordinated the writing of this paper. P.T.B., J.S.H., P.M.M., N.S.P., J.S., S.R.T., M.V., and S.J.V. wrote the paper, made the figures, and generated the bibliography. K.P.I., N.L., N.S.P., X.S., J.S., J.P.S., S.R.T., and S.J.V. performed analyses and developed new techniques on a preliminary 14-year data set. P.T.B., B.B., B.D.C., S.C., B.D., G.E.F., K.G., A.D.J., A.R.K., L.Z.K., N.L., W.G.L., N.S.P., S.C.S., L.S., X.S., J.S., J.P.S., S.R.T., C.U., M.V., S.J.V., Q.W., and C.A.W. analyzed preliminary versions of the 15-year data set. J.S.H., J.G., N.L., N.S.P., J.P.S., and J.K.S. performed noise analyses on the data set. P.T.B., B.B., R.B., R.C., M.C., N.J.C., D.D., H.D., G.E.F., K.A.G., S.H., A.D.J., N.L., W.G.L., P.M.M., P.P., N.S.P., S.C.S., X.S., J.S., J.P.S., J.T., S.R.T., M.V., S.J.V., and C.A.W. performed the Bayesian and frequentist analyses presented in this paper. K.G., J.S.H., P.M.M., J.D.R., and S.J.V. computed the background distribution of our detection statistics. P.R.B., P.M.M., N.S.P., and M.V. performed simulations that were used to compute false alarm probabilities. A.M.A., P.T.B., B.B., B.D., S.C., J.A.E., G.E.F., J.S.H., A.D.J., A.R.K., N.L., M.T.L., K.D.O., T.T.P., N.S.P., J.S., J.P.S., J.K.S., S.R.T., M.V., and S.J.V. contributed to the development of our code. P.T.B., N.J.C., J.S.H., A.D.J., T.B.L., P.M.M., J.D.R., and M.V. performed an internal code review. K.C., C.C., J.A.E., J.S.H., D.R.M., M.A.Mc., C.M.F.M., K.D.O., N.S.P., S.R.T., M.V., R.v.H., and S.J.V. were part of a NANOGrav Detection Study Group. M.A.M. and P.N. served on the IPTA Detection Committee. K.C., J.S.H., T.J.W.L., P.M.M., N.S.P., J.D.R., X.S., S.R.T., S.J.V., and C.A.W. responded to the Detection Checklist from the IPTA Detection Committee.

Acknowledgments. The NANOGrav collaboration receives support from National Science Foundation (NSF) Physics Frontiers Center award numbers 1430284 and 2020265, the Gordon and Betty Moore Foundation, NSF AccelNet award number 2114721, an NSERC Discovery Grant, and CIFAR. The Arecibo Observatory is a facility of the NSF operated under cooperative agreement (AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. Méndez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Specifically, it used the Bridges-2 system, which is supported by NSF award number ACI-1928147, at the Pittsburgh Supercomputing Center (PSC). This work was conducted using the Thorny Flat HPC Cluster at West Virginia University (WVU), which is funded in part by National Science Foundation (NSF) Major Research Instrumentation Program (MRI) Award number 1726534, and West Virginia University. This work was also conducted in part using the resources of the Advanced Computing Center for Research and Education (ACCRE) at Vanderbilt University, Nashville, TN. This work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. NANOGrav is part of the International Pulsar Timing Array (IPTA); we would like to thank our IPTA colleagues for their feedback on this paper. We thank members of the IPTA Detection Committee for developing the IPTA Detection Checklist. We thank Bruce Allen for useful feedback. We thank Valentina Di Marco and Eric Thrane for uncovering a bug in the sky scramble code. We thank Jolien Creighton and Leo Stein for helpful conversations about background estimation. L.B. acknowledges support from the National Science Foundation under award AST-1909933 and from the Research Corporation for Science Advancement under Cottrell Scholar Award No. 27553. P.R.B. is supported by the Science and Technology Facilities Council, grant number ST/W000946/1. S.B. gratefully acknowledges the support of a Sloan Fellowship, and the support of NSF under award #1815664. The work of R.B., R.C., D.D., N.La., X.S., J.P.S., and J.T. is partly supported by the George and Hannah Bolinger Memorial Fund in the College of Science at Oregon State University. M.C., P.P., and S.R.T. acknowledge support from NSF AST-2007993. M.C. and N.S.P. were supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. K.Ch., A.D.J., and M.V. acknowledge support from the Caltech and Jet Propulsion Laboratory President’s and Director’s Research and Development Fund. K.Ch. and A.D.J. acknowledge support from the Sloan Foundation. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc./National Radio Astronomy Observatory. Support for H.T.C. is provided by NASA through the NASA Hubble Fellowship Program grant #HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. K.Cr. is supported by a UBC Four Year Fellowship (6456). M.E.D. acknowledges support from the Naval Research Laboratory by NASA under contract S-15633Y. T.D. and M.T.L. are supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. E.C.F. is supported by NASA under award 80GSFC21M0002. G.E.F., S.C.S., and S.J.V. are supported by NSF award PHY-2011772. K.A.G. and S.R.T. acknowledge support from an NSF CAREER award #2146016. The Flatiron Institute is supported by the Simons Foundation. S.H. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1745301. N.La. acknowledges the support from Larry W. Martin and Joyce B. O’Neill Endowed Fellowship in the College of Science at Oregon State University. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). D.R.L. and M.A.Mc. are supported by NSF #1458952. M.A.Mc. is supported by NSF #2009425. C.M.F.M. was supported in part by the National Science Foundation under Grants No. NSF PHY-1748958 and AST-2106552. A.Mi. is supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121 Quantum Universe - 390833306. P.N. acknowledges support from the BHI, funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. K.D.O. was supported in part by NSF Grant No. 2207267. T.T.P. acknowledges support from the Extragalactic Astrophysics Research Group at Eötvös Loránd University, funded by the Eötvös Loránd Research Network (ELKH), which was used during the development of this research. S.M.R. and I.H.S. are CIFAR Fellows. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. J.D.R. also acknowledges support from start-up funds from Texas Tech University. J.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388, and acknowledges previous support by the NSF under award 1847938. C.U. acknowledges support from BGU (Kreitman fellowship), and the Council for Higher Education and Israel Academy of Sciences and Humanities (Excellence fellowship). C.A.W. acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship. O.Y. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-2139292.

Dedication. This paper is dedicated to the memory of Donald Backer: a pioneer in pulsar timing arrays, a term he coined; a discoverer of the first millisecond pulsar; a master developer of pulsar timing instrumentation; a founding member of NANOGrav; and a friend and mentor to many of us.

Facilities: Arecibo, GBT, VLA

Software: acor, astropy ([22]), ceffyl ([90]), chainconsumer ([71]), ENTERPRISE ([47]), enterprise_extensions ([135]), hasasia ((Hazboun et al. 2019)), holodeck ([85]), Jupyter ([87]), libstempo ([142]), matplotlib ([77]), numpy ([62]), PINT ([93]), PTMCMC ([46]), scipy ([148]), Tempo2 ([72])

Appendix

Additional data set details

The observations included in the NANOGrav 15-year data set were performed between July 2004 and August 2020 with the 305-m Arecibo Observatory (Arecibo), the 100-m Green Bank Telescope (GBT), and, since 2015, the 27 25-m antennae of the Very Large Array (VLA). We used Arecibo to observe the 33 pulsars that lie within its declination range ($0^\circ< \delta< +39^\circ$); GBT to observe the pulsars that lie outside of Arecibo’s range, plus J1713+0747 and B1937+21, for a total of 36 pulsars; the VLA to observe the seven pulsars J0437−4715, J1600−3053, J1643−1224, J1713+0747, J1903+0327, J1909−3744, and B1937+21. Six of these were also observed with Arecibo, GBT, or both; J0437−4715 was only visible to the VLA. Figure 12 shows the sky locations of the 67 pulsars used for the GWB search (top) and the distribution of angular separations for the pulsar pairs (bottom).

Initial observations were performed with the ASP (Arecibo) and GASP (GBT) systems, with 64-MHz bandwidth (Demorest [38]). Between 2010 and 2012, we transitioned to the PUPPI (Arecibo) and GUPPI (GBT) systems, with bandwidths up to 800 MHz ([44, 50]). We observe pulsars in two different radio-frequency bands in order to measure pulse dispersion from the interstellar medium: at Arecibo, we use the 1.4 GHz receiver plus either the 430 MHz or 2.1 GHz receiver (and the 327 MHz receiver for early observations of J2317+1439); at GBT, we use the 820 MHz and 1.4 GHz receivers; at the VLA, we use the 1.4 GHz and 3 GHz receivers with the YUPPI system.

In § 5.4 we analyze also two split-telescope data sets: 33 pulsars for Arecibo, and 35 for GBT (excluding J0614−3329, which was observed for less than three years). For the two pulsars timed by both telescopes (J1713+0747 and B1937+21), we partition the timing data between the telescopes and obtain independent timing solutions for each. We do not analyze a VLA-only data set, which would have shorter observation spans and significantly reduced sensitivity.

Bayesian Methods & Diagnostics

The prior probability distributions assumed for all analyses in this paper are listed in Table 1. We use Markov chain Monte Carlo (MCMC) techniques to sample randomly from the joint posterior distribution of our model parameters. Marginal distributions are obtained simply by considering only the parameter of interest

parameter	description	prior	comments
		white noise
$E_k$	EFAC per backend/receiver system	Uniform $[0, 10]$	single-pulsar analysis only
$Q_k$ [s]	EQUAD per backend/receiver system	log-Uniform $[-8.5, -5]$	single-pulsar analysis only
$J_k$ [s]	ECORR per backend/receiver system	log-Uniform $[-8.5, -5]$	single-pulsar analysis only
		intrinsic red noise
$A_{\mathrm{red}}$	red-noise power-law amplitude	log-Uniform $[-20, -11]$	one parameter per pulsar
$\gamma_{\mathrm{red}}$	red-noise power-law spectral index	Uniform $[0, 7]$	one parameter per pulsar
		all common processes, free spectrum
$\rho_i$ [$\mathrm{s}^2$]	power-spectrum coefficients at $f=i/T$	log-Uniform in $\rho_i$ $[-18, -8]$	one parameter per frequency
		all common processes, power-law spectrum
$A$	common process strain amplitude	log-Uniform $[-18, -14]$ ($\gamma=13/3$)	one parameter for PTA
		log-Uniform $[-18, -11]$ ($\gamma$ varied)	one parameter for PTA
$\gamma$	common process power-law spectral index	delta function ($\gamma=13/3$)	fixed
		Uniform $[0, 7]$	one parameter for PTA
		all common processes, broken–power-law spectrum
$A$	broken–power-law amplitude	log-Uniform $[-18, -11]$	one parameter for PTA
$\gamma$	broken–power-law low-freq. spectral index	Uniform $[0, 7]$	one parameter per PTA
$\delta$	broken–power-law high-freq. spectral index	delta function ($\delta=0$)	fixed
$f_{\mathrm{bend}}$ [Hz]	broken–power-law bend frequency	log-Uniform $[-8.7, -7]$	one parameter for PTA
$\ell$	broken–power-law high-freq. transition sharpness	delta function ($\ell=0.1$)	fixed
		all common processes, $t$-process spectrum
$A$	power-law amplitude	log-Uniform $[-18, -11]$	one parameter for PTA
$\gamma$	power-law spectral index	Uniform $[0, 7]$	one parameter per PTA
$\xi_i$	modification factor	Inverse Gamma Distribution	one parameter per frequency
		all common processes, turnover spectrum
$A$	turnover power-law amplitude	log-Uniform $[-18, -11]$	one parameter for PTA
$\gamma$	turnover power-law high-freq. spectral index	Uniform $[0, 7]$	one parameter per PTA
$\kappa$	turnover power-law low-freq. spectral index	Uniform $[0, 7]$	one parameter per PTA
$f_0$ [Hz]	turnover power-law bend frequency	log-Uniform $[-9, -7]$	one parameter for PTA
		all common processes, cross-correlation spline model
$y$	normalized cross-correlation values at spline knots $(10^{-3}, 25, 49.3, 82.5, 121.8, 150, 180)^\circ$	Uniform $[-0.9, 0.9]$	seven parameters for PTA

in each sample. To assess convergence of our MCMC runs beyond visual inspection we use the Gelman–Rubin statistic, requiring $\hat{R} < 1.01$ for all parameters ([56, 146]). We performed most runs discussed in this paper with the PTMCMC sampler ([46]) and postprocessed samples with chainconsumer ([71]).

In NG12gwb we use an analytic approximation for the uncertainty of marginalized-posterior statistics (Wilcox [149]). Here we instead adopt a bootstrap approach: we resample the original MCMC samples (with replacement) to generate new sets that act as independent sampling realizations. We then calculate the distributions of the desired summary statistics (e.g., quantiles, marginalized posterior values) over these sets. From these distributions, we determine central values and uncertainties (either medians and 68% confidence intervals, or means and standard deviations).

We rely on a variety of techniques to perform Bayesian model comparison. The first is thermodynamic integration (e.g., Ogata [107]; Gelman & Meng [54]), which computes Bayesian evidence integrals directly through parallel tempering: we run $N_\beta$ MCMC chains that explore variants of the likelihood function raised to different exponents $\beta$, then approximate the evidence for model $\mathcal{H}$ as

where all likelihoods and posteriors are computed within model $\mathcal{H}$, $\theta$ denotes all of the model’s parameters, and the expectation $\langle\ln p(d|\theta)\rangle_\beta$ is approximated by MCMC with respect to the posterior $p_\beta(\theta|d) \propto p(d|\theta,\mathcal{H})^\beta p(\theta,\mathcal{H})$. The inverse temperatures $\beta$ are spaced geometrically, as is the default in PTMCMC.

To compare nested models, which differ by “freezing” a subset of parameters, we also use the Savage–Dickey density ratio (Dickey [42]): if models $\mathcal{H}$ and $\mathcal{H}_0$ differ by the fact that (say) $\theta_0$ is frozen to 0 in the latter, then $p(d|\mathcal{H}_0)/p(d|\mathcal{H})=p(\theta_0=0|d,\mathcal{H})/p(\theta_0=0|\mathcal{H})$.

When comparing disjoint models with different likelihoods (e.g., HD versus CURN), we use product-space sampling (Carlin & Chib [30]; Godsill [57]). This method treats model comparison as a parameter estimation problem, where we sample the union of the unique parameters of all models, plus a model-indexing parameter that activates the relevant likelihood function and parameter space of one of the sub-models. Bayes factors are then obtained by counting how often the model index falls in each activation region and taking ratios of those counts.

In some situations, it can be difficult to sample a computationally expensive model directly. In these cases, we sample a computationally cheaper approximate distribution and reweight those posterior samples to estimate the posterior for the computationally expensive model (Hourihane et al. [75]). The reweighted posterior can be used in the thermodynamic-integration or Savage–Dickey methods. In addition, the mean of the weights yields the Bayes factor between the expensive and approximate models, which may be of direct interest (e.g., HD can be approximated by CURN). We estimate Bayes-factor uncertainties using bootstrapping and, for product-space sampling, with the Markov-model techniques of [35] and [68].

Broken power-law model

As shown in NG12gwb, the simultaneous Bayesian estimation of white measurement noise and of red-noise processes described by power laws biases the recovery of the spectral index of the latter [88] (Hazboun et al. 2019). Just as in NG12gwb and [13], we impose a high-frequency cutoff on the red-noise processes. To choose the cutoff frequency, we perform inference on our data with a CURN$^\gamma$ model modified so that the common process has power spectral density

then set the cutoff to the MAP $f_{\mathrm{break}}$. Equation (15) is fairly generic, allowing for separate spectral indices at low ($\gamma$) and high ($\delta$) frequencies. The break frequency $f_{\mathrm{break}}$ dictates where the broken power law changes spectral index, while $\ell$ (which we set to 0.1) controls the smoothness of the transition.

The marginal posterior for $f_{\mathrm{break}}$, obtained in the factorized-likelihood approximation using the techniques of Lamb et al. [90], has median and 90% credible region of $3.2^{+5.4}_{-1.2} \times10^{-8}$ Hz, and a MAP value of $2.75 \times10^{-8}$ Hz. The latter is close to $f_{14} = 14/T$ in our frequency basis (with $T$ the total span of the data set), so we use 14 frequencies to model common-spectrum noise processes (see § 2 and NG12gwb).

t-process spectrum model

The free-spectrum analysis of our data (§ 5.2 and Figure 6) shows that the frequency bins at $f_1$, $f_6$, $f_7$, and $f_8$ appear to be in tension with a pure power law, skewing the estimation of $\gamma$ and reducing the HD$^{13/3}$ vs. CURN$^{13/3}$ Bayes factor. Assuming that those frequency components reflect unmodeled systematics or stronger-than-expected statistical fluctuations, we can make our inference more robust to such outliers with a “fuzzy” power-law model that allows the individual $\Phi_i$ to vary more freely around their expected values. To wit, we introduce the $t$-process spectrum (TPS)

where $\Phi_{\mathrm{powerlaw},i}$ follows Equation (6) and $x$ follows the inverse gamma distribution with parameters $\alpha=\beta=1$; the resulting Gaussian mixture yields a Student’s-$t$ distribution for the $\Phi_{\mathrm{TPS},i}$. Figure 13 shows CURN$^\gamma$ power-law posteriors and CURN$^{\mathrm{TPS}}$ modified power-law posteriors, obtained in the factorized-likelihood approximation [138, 90] and compared to CURN$^{\mathrm{free}}$ bin variances. The TPS model is spread more widely and deviates from the perfect power law at bins $f_1$, $f_6$, $f_7$, and $f_8$, as expected. The right panel of Figure 13 shows the joint $\log_{10} A,\gamma$ posteriors for CURN$^\gamma$ and CURN$^{\mathrm{TPS}}$. The latter is more consistent with steeper power laws, and it includes $\gamma=13/3$ at $1\sigma$ credibility.

Turnover model

The final parameterized spectral model that we investigate is motivated by the idea that the dynamics of SMBHBs are influenced by their environments at sub-parsec separations [15, 128, 99]. These interactions affect binary evolution and the resulting spectrum of the GWB. The process of bringing two SMBHs together after galaxy mergers involves a complex chain of interactions: despite significant theoretical work, the lack of observational constraints makes it difficult to draw any conclusions. PTAs, however, provide a unique opportunity to probe the timescale over which two SMBHs evolve from the merger of their galaxies to a bound binary that produces GW signals in the PTA sensitivity band.

When dynamical interactions dominate orbital evolution, binaries will traverse the GW spectrum more quickly, reducing GW emission compared to a GW-driven inspiral. This kind of behavior is captured by the turnover model [124]:

This is qualitatively similar to the broken power law discussed earlier, except that here $f_0$ represents the GW frequency at which typical binary evolution transitions from environmentally dominated (at lower frequencies and wider separations) to GW-dominated (at higher frequencies and smaller separations). The parameter $\kappa$ controls the shape of the spectrum below $f_0$, and depends on the orbital-evolution mechanism. Note that the actual turning point of the spectrum is not at $f_0$ but at $f_{\mathrm{bend}}=f_0\times(3\kappa/4-1)^{1/\kappa}$ (NG9gwb).

Applying this model to our data, we find hints of departures from a pure power law: the transition frequency $f_0$ lies below 10 nHz with 65% credibility, while the bend frequency lies below 10 mHz with 75% credibility.

Nevertheless, Bayesian comparison of this CURN$^{\mathrm{turnover}}$ model with CURN$^{\gamma}$ reports an inconclusive Bayes factor of $1.46 \pm0.02$ in favor of CURN$^{\mathrm{turnover}}$. Furthermore, the estimation of CURN$^{\mathrm{turnover}}$ parameters is sensitive to DM modeling (see § 5.1). While the spectra are broadly consistent whether we use DMX or DMGP to model DM fluctuations, there are differences in the power at certain frequencies that lead to differences in the turnover parameters. This is discussed in greater detail in Agazie et al. [3].

Sky scrambles

In the sky-scramble method (Cornish & Sampson [36]), inter-pulsar correlations are analyzed as if the pulsars occupied random sky positions, with the purpose of creating a background distribution of PTA detection statistics for null-hypothesis testing, in alternative to phase shifts (Taylor et al. [137]; see § 3 and § 4). If a correlated signal is present in the data, phase shifts and sky scrambles actually test different null hypotheses: phase shifts test the hypothesis that no inter-pulsar correlations are present, while sky scrambles assume that inter-pulsar correlations are present at the level measured in the data, but test the hypothesis that these correlations have no dependence on angular separation.

As is the convention in the literature, we require that scrambled overlap reduction functions (ORFs) be independent of each other and of the unscrambled ORF using a match statistic,

where $\Gamma_{ab}$ and $\Gamma'_{ab}$ are two different ORFs. For the sky scrambles used in our analysis, the scrambled ORFs have $\bar{M}<0.1$ with respect to the unscrambled ORF, and $\bar{M}<0.17$ with each other. We generate 10,000 sky scrambles, owing to the difficulty in obtaining large numbers of scrambled ORFs that satisfy the match threshold; because of limitations of computational resources, we obtain our detection statistics for 5,000 of those ORFs. Figure 14 shows the resulting background distributions for the HD$^{\gamma}$-to-CURN$^{\gamma}$ Bayes factor (left panel) and the optimal-statistic S/N (right panel). The Bayes factors exceed the observed value in eight of the 5,000 sky scrambles ($p=1.6\times10^{-3}$), while none of the sky scrambles have noise-marginalized mean S/N greater than observed ($p<10^{-4}$).

We note that the null distribution recovered by the sky scrambles is not very sensitive to the choice of match threshold for $|\bar{M}|\lesssim0.2$. Figure 15 compares the null distributions when the match threshold for all ORFs with each other and with the unscrambled ORF is set to $|\bar{M}|<0.17$ (blue), $|\bar{M}|<0.1$ (orange), and $|\bar{M}|<0.08$ (green). There is very little difference among the distributions; however, imposing a smaller threshold means that fewer sky scrambles can be used (6,043 with $|\bar{M}|<0.1$ and 1,534 with $|\bar{M}|<0.08$, compared to 10,000 with $|\bar{M}|<0.17$), which limits the precision with which the $p$-value can be measured. We find no evidence that the recovered null distribution is biased when including sky scrambles with matches up to 0.17.

Note—All values were computed for the 15-year data set, assuming a power-law power spectral density using the 14 lowest frequency components. Here $\hat{A}^{2}$, S/N, and AIC are marginalized over pulsar noise parameters with fixed $\gamma= 13/3$. The numbers in parentheses represent the mean least-squares errors for the $\hat{A}^{2}_{k}$ coefficients and standard deviations over noise-parameter posteriors for S/Ns. We compute $p(\mathrm{AIC})$ with respect to the model with the lowest mean AIC (i.e., HD + monopole).

Multiple-correlation optimal statistic

The multiple-correlation optimal statistic (MCOS; Sardesai & Vigeland [125]) fits the inter-pulsar correlation coefficients $\rho_{ab}$ with a linear model that includes multiple components with different correlation patterns, but with the same spectral shape. The linear-model coefficients are the squared amplitudes of the components. Within such a model, the significance of each component can be quoted as a S/N given by its best-fit coefficient divided by the fit error. Just as for the noise-marginalized optimal statistic (Vigeland et al. [147]), the posterior distribution of pulsar noise parameters induces a distribution of MCOS statistics.

We fit the 15-year data with models that include HD, monopole, and dipole-correlated components in various combinations. Table 2 lists the noise-marginalized amplitude estimates and S/N for all models. The goodness-of-fit of the models can be compared using the Akaike Information Criterion (AIC; Akaike [7]):

Model	HD Correlations $\hat{A}^{2}_{\mathrm{HD}}$	S/N	Monopole Correlations $\hat{A}^{2}_{\mathrm{mo}}$	S/N	Dipole Correlations Mean $\hat{A}^{2}_{\mathrm{di}}$	S/N	$p(\mathrm{AIC})$
HD only	$6.8(9) \times10^{-30}$	$4(1)$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$3 \times10^{-2}$
Monopole only	$\cdots$	$\cdots$	$1.1(1) \times10^{-30}$	$4(1)$	$\cdots$	$\cdots$	$6 \times10^{-3}$
Dipole only	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$1.5(3) \times10^{-30}$	$4(1)$	$8 \times10^{-4}$
HD + monopole	$5.5(8) \times10^{-30}$	$3.4(8)$	$8(1) \times10^{-31}$	$2.9(8)$	$\cdots$	$\cdots$	$1$
HD + dipole	$5.5(8) \times10^{-30}$	$3.2(7)$	$\cdots$	$\cdots$	$8(2) \times10^{-31}$	$1.7(7)$	$6 \times10^{-2}$
monopole + dipole	$\cdots$	$\cdots$	$8(1) \times10^{-31}$	$2.7(7)$	$9(2) \times10^{-31}$	$1.9(6)$	$1 \times10^{-2}$
HD + monopole + dipole	$5.1(8) \times10^{-30}$	$2.9(6)$	$7(1) \times10^{-31}$	$2.4(6)$	$3(2) \times10^{-31}$	$0.6(4)$	$0.48$

where $k$ is the number of model parameters and $\chi^{2}$ is the fit’s chi-squared, computed without accounting for GW-induced $\rho_{ab}$ correlations. (This can be thought of as a pseudo-Bayes factor, with the factor of $2k$ imposing an Occam penalty.) The relative probability of a model compared to the most-favored model is then given by

where $\mathrm{AIC}_{\min}$ is the minimum AIC across all models. Table 2 lists the AIC probabilities, computed by averaging the AIC of each model over pulsar noise parameters. The HD-correlated model is preferred among the models with a single correlated process. The models with both HD and monopole correlations are preferred among all models: for a model with HD and monopole correlations, we find S/N of $3.4 \pm0.8$ for HD correlations and $2.9 \pm0.8$ for monopolar correlations, while for a model with HD, monopole, and dipole correlations, we find S/N of $2.9 \pm0.6$ for HD correlations, $2.4 \pm0.6$ for monopole correlations, and $0.6 \pm0.4$ for dipole correlations (means $\pm$ standard deviations across noise-parameter posteriors). The statistical significance of these S/Ns can be quantified empirically using simulations of 15-year–like data sets (see App. § 16.§ 1), which report $p$-values $< 10^{-2}$ and $\simeq4 \times10^{-2}$ for the observed mean HD and monopole statistics across data replications with no spatially correlated injections.

As discussed in Sardesai & Vigeland [125], the optimal statistic and the MCOS are metrics of the apparent spatial correlation pattern of the data, but they have a limited ability to identify its actual source. That is because a real HD signal may also excite the monopole optimal statistic and the monopole component of the MCOS; conversely, a real monopolar signal may also excite the HD optimal statistic and the HD component of the MCOS; and so on. The S/Ns quoted in Table 2 quantify how often we would expect to measure the observed value of the optimal statistic if only uncorrelated noise is present, but they do not describe how often one type of correlated noise would produce a given value of the optimal statistic for a different type of correlation.

This effect can be characterized using simulations (see App. § 16.§ 1), which report a $p$-value of $0.11$ for the observed mean monopole statistic when a HD-correlated signal with the MAP 15-year amplitude is included in the simulated data sets. We conclude that there are some indications of a possible monopole-correlated signal in the data with S/N comparable to but smaller than the S/N for HD correlations; however, from simulations we conclude that it is possible for such a signal to appear in an MCOS analysis if only a HD-correlated stochastic process is present.

Multiple-correlation optimal statistic simulations

In this appendix we obtain the distribution of the MCOS over an ensemble of simulated data sets, with the goal of characterizing the probability that the observed S/Ns could have been produced by pulsar noise alone, or by a GWB with HD correlations. Unlike our Bayesian analysis, the MCOS prefers a model that includes both HD and monopolar components. So we are especially interested in asking how frequently we may expect the observed MCOS monopole if the data contain only the GWB. In Apps. H.1 and H.2 we present two different types of simulations: “astrophysical,” where we generate synthetic data with MAP noise parameters inferred from the 15-year data set, both with and without the GWB; and “model checking,” where we create data replications following the HD$^{13/3}$ posteriors for the real data set. Note that neither simulation attempts to account for the monochromatic character of the putative monopolar signal (see § 5.2).

Astrophysical simulations

Following Pol et al. [115], we generate simulated data sets adopting MAP pulsar-noise parameters obtained from the real data independently for each pulsar; these “noise runs” include an additional power-law process to reduce contamination between the putative GWB and the pulsars’ intrinsic red noise (Taylor et al. [138]). We produce 100 realizations each of three different simulations: (i) injecting no spatially correlated power-law GWB or excess uncorrelated common-spectrum noise; (ii) injecting a spatially correlated power-law GWB with amplitude $2.7\times10^{-15}$ and spectral index $13/3$; and (iii) injecting no GWB or common-spectrum noise, but omitting the additional power-law process in the estimation of intrinsic pulsar noise, with the goal of testing how often excess common-spectrum noise is recognized as a spatially correlated GWB. We compute HD + monopole + dipole MCOS S/Ns for all synthetic data sets (see Figure 17). The mean HD S/Ns observed in the real data (see App. § 15) correspond to $p$-values of $<10^{-2}$ for simulations (i) and (iii), and 0.64 for simulation (ii). The mean monopole S/Ns observed in the real data set correspond to $p$-values of $4\times10^{-2}$, $1.1\times10^{-1}$, and $<10^{-2}$ for simulations (i), (ii), and (iii) respectively. We conclude that it is unlikely that we would measure HD correlations at the level observed in real data when no correlated signal is present (simulation (i)) or when only uncorrelated common-spectrum red noise is present (simulation (iii)). In addition, the HD S/Ns obtained from a HD-correlated GWB injection (simulation (ii)) are fully consistent with the S/N observed in real data. By contrast, the observed monopole S/N could have been produced by intrinsic pulsar noise alone, or by a real HD signal.

Model-checking simulations

In App. § 16.§ 1 we have tackled the question of monopole S/N significance using simulations based on real-data MAP estimates $\boldsymbol{\eta}^{\mathrm{MAP}}$ of pulsar-noise and GW parameters. In this appendix we adopt a procedure with a stronger Bayesian flavor, evaluating the MCOS on a population of data replications created using $\mathrm{HD}^{13/3}$ as a generative model with noise hyperparameters $\eta$ drawn from the $\mathrm{HD}^{13/3}$ real-data posterior. This can be seen also as a Bayesian model-checking exercise [55, 53]: if we find that the summary statistic of interest (the monopole MCOS) has a much more extreme value in real data than in data replications, we should suspect that the data model (here $\mathrm{HD}^{13/3}$) is missing something.

We perform the test by drawing 500 parameter vectors $\{\eta^{(k)}\}$ from the $\mathrm{HD}^{13/3}$ real-data posterior; for each $\eta^{(k)}$ we simulate a data set $\delta t^{\mathrm{sim},(k)} \sim p(\delta t \mid\eta^{(k)})$ and compare $\mathrm{MCOS}(\delta t^{\mathrm{sim},(k)};\eta^{(k)})$ with $\mathrm{MCOS}(\delta t;\eta^{(k)})$. Our notation emphasizes the dependence of the MCOS on the pulsar noise parameters through the $\mathbf{P}$ matrices in Equation (9). Figure 18 shows the resulting distribution of monopole S/Ns. The replicated monopole S/N is greater than its observed counterpart for 11% of the draws. Thus, it is plausible that the MCOS could measure the observed monopole S/N in data that contain only a HD-correlated GWB. Conversely, the observed monopole S/N does not by itself suggest that $\mathrm{HD}^{13/3}$ is misspecified.