# How to Detect an Astrophysical Nanohertz Gravitational-Wave Background

Bence Bécsy Department of Physics, Oregon State University, Corvallis, OR
97331, USA  Neil J. Cornish eXtreme Gravity Institute, Department of Physics,
Montana State University, Bozeman, MT 59717, USA  Patrick M. Meyers Division
of Physics, Mathematics, and Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Luke Zoltan Kelley Department of Astronomy,
University of California at Berkeley, Berkeley, CA 94720, USA  Gabriella
Agazie Center for Gravitation, Cosmology and Astrophysics, Department of
Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Akash Anumarlapudi Center for
Gravitation, Cosmology and Astrophysics, Department of Physics, University of
Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Anne M. Archibald Newcastle
University, NE1 7RU, UK  Zaven Arzoumanian  X-Ray Astrophysics Laboratory,
NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA  Paul T.
Baker Department of Physics and Astronomy, Widener University, One University
Place, Chester, PA 19013, USA  Laura Blecha Physics Department, University of
Florida, Gainesville, FL 32611, USA  Adam Brazier Cornell Center for
Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Cornell Center for Advanced Computing,
Cornell University, Ithaca, NY 14853, USA  Paul R. Brook Institute for
Gravitational Wave Astronomy and School of Physics and Astronomy, University
of Birmingham, Edgbaston, Birmingham B15 2TT, UK  Sarah Burke-Spolaor Sloan
Fellow  Department of Physics and Astronomy, West Virginia University, P.O.
Box 6315, Morgantown, WV 26506, USA  Center for Gravitational Waves and
Cosmology, West Virginia University, Chestnut Ridge Research Building,
Morgantown, WV 26505, USA  J. Andrew Casey-Clyde Department of Physics,
University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046,
USA  Maria Charisi Department of Physics and Astronomy, Vanderbilt University,
2301 Vanderbilt Place, Nashville, TN 37235, USA  Shami Chatterjee Cornell
Center for Astrophysics and Planetary Science and Department of Astronomy,
Cornell University, Ithaca, NY 14853, USA  Katerina Chatziioannou  Division of
Physics, Mathematics, and Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Tyler Cohen Department of Physics, New Mexico
Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA
James M. Cordes Cornell Center for Astrophysics and Planetary Science and
Department of Astronomy, Cornell University, Ithaca, NY 14853, USA  Fronefield
Crawford Department of Physics and Astronomy, Franklin & Marshall College,
P.O. Box 3003, Lancaster, PA 17604, USA  H. Thankful Cromartie NASA Hubble
Fellowship: Einstein Postdoctoral Fellow  Cornell Center for Astrophysics and
Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY
14853, USA  Kathryn Crowter Department of Physics and Astronomy, University of
British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada  Megan
E. DeCesar George Mason University, resident at the Naval Research Laboratory,
Washington, DC 20375, USA  Paul B. Demorest National Radio Astronomy
Observatory, 1003 Lopezville Rd., Socorro, NM 87801, USA  Timothy Dolch
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI
49242, USA  Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA
94602-3017, USA  Elizabeth C. Ferrara Department of Astronomy, University of
Maryland, College Park, MD 20742  Center for Research and Exploration in Space
Science and Technology, NASA/GSFC, Greenbelt, MD 20771  NASA Goddard Space
Flight Center, Greenbelt, MD 20771, USA  William Fiore Department of Physics
and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506,
USA  Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Emmanuel Fonseca
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Gabriel E. Freedman Center for Gravitation, Cosmology and Astrophysics,
Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Nate Garver-Daniels Department of
Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV
26506, USA  Center for Gravitational Waves and Cosmology, West Virginia
University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Peter
A. Gentile Department of Physics and Astronomy, West Virginia University, P.O.
Box 6315, Morgantown, WV 26506, USA  Center for Gravitational Waves and
Cosmology, West Virginia University, Chestnut Ridge Research Building,
Morgantown, WV 26505, USA  Joseph Glaser Department of Physics and Astronomy,
West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Deborah C. Good Department of
Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT
06269-3046, USA  Center for Computational Astrophysics, Flatiron Institute,
162 5th Avenue, New York, NY 10010, USA  Kayhan Gültekin Department of
Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI 48109, USA
Jeffrey S. Hazboun Department of Physics, Oregon State University, Corvallis,
OR 97331, USA  Sophie Hourihane Division of Physics, Mathematics, and
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  Ross
J. Jennings NANOGrav Physics Frontiers Center Postdoctoral Fellow  Department
of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown,
WV 26506, USA  Center for Gravitational Waves and Cosmology, West Virginia
University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Aaron
D. Johnson Center for Gravitation, Cosmology and Astrophysics, Department of
Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Division of Physics, Mathematics, and
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  Megan
L. Jones Center for Gravitation, Cosmology and Astrophysics, Department of
Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Andrew R. Kaiser Department of Physics
and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506,
USA  Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  David L. Kaplan
Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Matthew Kerr Space Science Division,
Naval Research Laboratory, Washington, DC 20375-5352, USA  Joey S. Key
University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA
Nima Laal Department of Physics, Oregon State University, Corvallis, OR 97331,
USA  Michael T. Lam SETI Institute, 339 N Bernardo Ave Suite 200, Mountain
View, CA 94043, USA  School of Physics and Astronomy, Rochester Institute of
Technology, Rochester, NY 14623, USA  Laboratory for Multiwavelength
Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA
William G. Lamb Department of Physics and Astronomy, Vanderbilt University,
2301 Vanderbilt Place, Nashville, TN 37235, USA  T. Joseph W. Lazio  Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA  Natalia Lewandowska Department of Physics,
State University of New York at Oswego, Oswego, NY, 13126, USA  Tyson B.
Littenberg NASA Marshall Space Flight Center, Huntsville, AL 35812, USA
Tingting Liu Department of Physics and Astronomy, West Virginia University,
P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational Waves and
Cosmology, West Virginia University, Chestnut Ridge Research Building,
Morgantown, WV 26505, USA  Duncan R. Lorimer Department of Physics and
Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Jing Luo Deceased
Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George
Street, Toronto, ON M5S 3H4, Canada  Ryan S. Lynch Green Bank Observatory,
P.O. Box 2, Green Bank, WV 24944, USA  Chung-Pei Ma Department of Astronomy,
University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA
94720, USA  Department of Physics, University of California, Berkeley, CA
94720, USA  Dustin R. Madison Department of Physics, University of the
Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA  Alexander McEwen Center
for Gravitation, Cosmology and Astrophysics, Department of Physics, University
of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  James W. McKee E.A. Milne Centre for
Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6
7RX, UK  Centre of Excellence for Data Science, Artificial Intelligence and
Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6
7RX, UK  Maura A. McLaughlin Department of Physics and Astronomy, West
Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Natasha McMann Department of
Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place,
Nashville, TN 37235, USA  Bradley W. Meyers Department of Physics and
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
BC V6T 1Z1, Canada  International Centre for Radio Astronomy Research, Curtin
University, Bentley, WA 6102, Australia  Chiara M. F. Mingarelli Center for
Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY
10010, USA  Department of Physics, University of Connecticut, 196 Auditorium
Road, U-3046, Storrs, CT 06269-3046, USA  Department of Physics, Yale
University, New Haven, CT 06520, USA  Andrea Mitridate Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany  Cherry Ng Dunlap
Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George
St., Toronto, ON M5S 3H4, Canada  David J. Nice Department of Physics,
Lafayette College, Easton, PA 18042, USA  Stella Koch Ocker Cornell Center for
Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Ken D. Olum Institute of Cosmology,
Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA
Timothy T. Pennucci Institute of Physics and Astronomy, Eötvös Loránd
University, Pázmány P. s. 1/A, 1117 Budapest, Hungary  Benetge B. P. Perera
Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA  Nihan S. Pol
Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA  Henri A. Radovan Department of Physics,
University of Puerto Rico, Mayagüez, PR 00681, USA  Scott M. Ransom National
Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
Paul S. Ray Space Science Division, Naval Research Laboratory, Washington, DC
20375-5352, USA  Joseph D. Romano Department of Physics, Texas Tech
University, Box 41051, Lubbock, TX 79409, USA  Shashwat C. Sardesai Center for
Gravitation, Cosmology and Astrophysics, Department of Physics, University of
Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Ann Schmiedekamp Department of
Physics, Penn State Abington, Abington, PA 19001, USA  Carl Schmiedekamp
Department of Physics, Penn State Abington, Abington, PA 19001, USA  Kai
Schmitz Institute for Theoretical Physics, University of Münster, 48149
Münster, Germany  Brent J. Shapiro-Albert Department of Physics and Astronomy,
West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Giant Army, 915A 17th Ave,
Seattle WA 98122  Xavier Siemens Department of Physics, Oregon State
University, Corvallis, OR 97331, USA  Center for Gravitation, Cosmology and
Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Joseph Simon NSF Astronomy and
Astrophysics Postdoctoral Fellow  Department of Astrophysical and Planetary
Sciences, University of Colorado, Boulder, CO 80309, USA  Magdalena S. Siwek
Center for Astrophysics, Harvard University, 60 Garden St, Cambridge, MA 02138
Sophia V. Sosa Fiscella School of Physics and Astronomy, Rochester Institute
of Technology, Rochester, NY 14623, USA  Laboratory for Multiwavelength
Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA
Ingrid H. Stairs Department of Physics and Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada  Daniel R.
Stinebring Department of Physics and Astronomy, Oberlin College, Oberlin, OH
44074, USA  Kevin Stovall National Radio Astronomy Observatory, 1003
Lopezville Rd., Socorro, NM 87801, USA  Abhimanyu Susobhanan Center for
Gravitation, Cosmology and Astrophysics, Department of Physics, University of
Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Joseph K. Swiggum NANOGrav Physics
Frontiers Center Postdoctoral Fellow  Department of Physics, Lafayette
College, Easton, PA 18042, USA  Stephen R. Taylor Department of Physics and
Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235,
USA  Jacob E. Turner Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944,
USA  Caner Unal Department of Physics, Ben-Gurion University of the Negev,
Be’er Sheva 84105, Israel  Feza Gursey Institute, Bogazici University,
Kandilli, 34684, Istanbul, Turkey  Michele Vallisneri Jet Propulsion
Laboratory, California Institute of Technology, 4800 Oak Grove Drive,
Pasadena, CA 91109, USA  Division of Physics, Mathematics, and Astronomy,
California Institute of Technology, Pasadena, CA 91125, USA  Rutger van
Haasteren Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-
Institut), Callinstrasse 38, D-30167, Hannover, Germany  Sarah J. Vigeland
Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Haley M. Wahl Department of Physics
and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506,
USA  Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Caitlin A. Witt
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, Evanston, IL 60208  Adler Planetarium, 1300 S.
DuSable Lake Shore Dr., Chicago, IL 60605, USA  Olivia Young School of Physics
and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA
Laboratory for Multiwavelength Astrophysics, Rochester Institute of
Technology, Rochester, NY 14623, USA

###### Abstract

Analysis of pulsar timing data have provided evidence for a stochastic
gravitational wave background in the nHz frequency band. The most plausible
source of such a background is the superposition of signals from millions of
supermassive black hole binaries. The standard statistical techniques used to
search for such a background and assess its significance make several
simplifying assumptions, namely: i) Gaussianity; ii) isotropy; and most often
iii) a power-law spectrum. However, a stochastic background from a finite
collection of binaries does not exactly satisfy any of these assumptions. To
understand the effect of these assumptions, we test standard analysis
techniques on a large collection of realistic simulated datasets. The dataset
length, observing schedule, and noise levels were chosen to emulate the
NANOGrav 15-year dataset. Simulated signals from millions of binaries drawn
from models based on the Illustris cosmological hydrodynamical simulation were
added to the data. We find that the standard statistical methods perform
remarkably well on these simulated datasets, despite their fundamental
assumptions not being strictly met. They are able to achieve a confident
detection of the background. However, even for a fixed set of astrophysical
parameters, different realizations of the universe result in a large variance
in the significance and recovered parameters of the background. We also find
that the presence of loud individual binaries can bias the spectral recovery
of the background if we do not account for them.

††software: enterprise (Ellis et al., 2019), enterprise_extensions (Taylor et
al., 2021), QuickCW (Bécsy et al., 2023), corner (Foreman-Mackey, 2016),
libstempo (Vallisneri, 2020), tempo (Nice et al., 2015), tempo2 (Hobbs et al.,
2006), PINT (Luo et al., 2019), matplotlib (Hunter, 2007), astropy (Price-
Whelan et al., 2018; Astropy Collaboration et al., 2013)

##  1 Introduction

Pulsar timing arrays (PTAs) monitor millisecond pulsars over a timespan of
decades. The times-of-arrival (TOAs) of these radio pulses are sensitive to
various physical processes, including variations in the proper distance
between the pulsars and the observer due to passing gravitational waves (GWs).
PTAs recently reported evidence for a nanohertz stochastic GW background (GWB)
in their most recent datasets (Agazie et al., 2023a; Antoniadis et al., 2023a;
Reardon et al., 2023; Xu et al., 2023). While such a background can be
explained by exotic phenomena (Afzal et al., 2023), their most plausible
source is a large collection of supermassive black hole binaries (SMBHBs,
Agazie et al. 2023b).

Several methods have been developed to search for a GWB in PTA datasets. While
the GWB is expected to first appear as a common uncorrelated red noise process
(CURN), the definitive signature of a GWB is considered to be the
characteristic Hellings-Downs (HD, Hellings & Downs 1983) angular correlation
pattern between pulsars. The most notable frequentist method is the so called
optimal statistic (OS, Chamberlin et al. 2015), which is an unbiased estimator
of the GWB amplitude and its variance based on the cross-correlations between
pulsars and the theoretical correlation pattern. It also allows one to define
a signal-to-noise ratio (S/N), which is used as a detection statistic.
Bayesian methods, on the other hand, model the background using Fourier basis
Gaussian processes (Lentati et al., 2013; van Haasteren & Vallisneri, 2014).
These can be constructed either with (HD model) or without HD correlations
between pulsars (CURN model). One can calculate the Bayes factor between these
two models, which is the most important quantity in quantifying the level of
Bayesian evidence for a GWB.

As with any statistical data analysis, one needs to make various assumptions
when calculating either the OS S/N or the HD vs. CURN Bayes factor. Most
commonly one assumes that the GWB is: i) Gaussian; ii) isotropic; iii)
described by a power-law spectrum111This also means we implicitly assume that
the signal is stationary.. However, if the GWB is produced by a finite number
of SMBHBs, none of these assumptions hold exactly (see e.g. Bécsy et al.
2022). Methods have been developed that relax some of these assumptions.
Assumption iii) is the one most often relaxed either by modeling the spectrum
as a two-component so called broken power-law model or allowing the amplitudes
of different frequency components to vary freely (see e.g. Sampson et al.
2015; Taylor et al. 2017; Agazie et al. 2023a, c; Meyers et al. 2023; Lamb et
al. 2023; Gersbach et al. 2023). There are also specific searches looking for
an anisotropic GWB (Mingarelli et al., 2013; Taylor & Gair, 2013; Romano et
al., 2015; Gair et al., 2015; Taylor et al., 2020; Ali-Haïmoud et al., 2020;
Pol et al., 2022; Agazie et al., 2023d), thus relaxing assumption ii). One way
of relaxing the assumption of Gaussianity is by modeling the background as a
t-process, where the GWB power at each frequency follows a Student’s
t-distribution instead of a Gaussian distribution, thus allowing more
flexibility (see Appendix D in Agazie et al. 2023a). While such generalized
search algorithms are available, the flagship results in most GWB searches
continue to rely upon these assumptions (see e.g. Agazie et al. 2023a;
Antoniadis et al. 2023a; Reardon et al. 2023). This is in part due to their
simplicity, which makes them easier to compute. In addition, it is also widely
expected that their assumptions are good enough to make them strong and robust
methods.

The expectation that models employing these assumptions are capable of
detecting GWBs, even when these assumptions are broken is based largely on
Cornish & Sesana (2013) and Cornish & Sampson (2016). In Cornish & Sesana
(2013) it was shown that even a single SMBHB produces HD correlations in an
array of isotropically distributed pulsars - showing that we can trade
isotropy of the background for isotropy of the pulsar array. In Cornish &
Sampson (2016) the authors studied both realistic GWBs based on SMBHB
population models and also idealized Gaussian, isotropic GWBs with a power-law
spectrum. They found that the significance with which one can detect these is
similar, except for cases with an unrealistically low number of binaries. In
this paper we carry out a similar analysis with a number of improvements:

  1. 1.

We use a more up-to-date SMBHB population model;

  2. 2.

Instead of pulsars with evenly sampled data, we use the actual observation
times from the NANOGrav 15-year dataset (Agazie et al., 2023e);

  3. 3.

In addition to white noise, we include pulsar red noise as well, and we set
the level of those based on the noise properties found in the NANOGrav 15yr
dataset (Agazie et al., 2023c);

  4. 4.

Instead of simple frequency-domain analysis, we use the actual software used
in Agazie et al. (2023a), that works in the time-domain, and takes into
account covariances with the timing model.

Our analysis can also be considered as an extension of the consistency checks
carried out in Johnson et al. (2023). There, these analysis pipelines were
extensively stress tested to make sure they perform properly under the
assumption that the data analyzed conforms to the models used. We extend the
scope to see how these pipelines perform if their assumptions are not met
perfectly, as it is expected to be the case when analyzing real data.

The rest of the paper is organized as follows. In Section 2, we describe our
realistic population-based simulated datasets and give more details about the
statistical methods employed on them. In Section 3, we present our results,
and in Section 4 we provide a brief conclusion and discuss future work.

##  2 Methods

###  2.1 Realistic simulated datasets

Our simulated datasets are based on the pulsars and their measured noise
properties in the NANOGrav 15-year dataset (Agazie et al., 2023e), following
the simulation framework of Pol et al. (2021). We use the libstempo software
package, and we simulate datasets with the actual observing times from the
real dataset, but to reduce the data volume we only keep one observation per
epoch. We simulate data with the white and red noise properties fixed to the
maximum a posteriori values inferred from the NANOGrav 15-year dataset (Agazie
et al., 2023c).

To simulate the contribution of a realistic GWB, we use an SMBHB population
model implemented in the holodeck software package (Kelley et al., 2023),
based on the Illustris cosmological hydrodynamical simulation (Vogelsberger et
al., 2014). This model applies a post processing step to account for small-
scale physics not captured by Illustris, uses a phenomenological model to
describe binary evolution (for details see Agazie et al. 2023b), and produces
a simulated list of binaries in the Universe (Kelley et al., 2017a, b, 2018).
To add the contribution of these binaries, we model the brightest 1000
binaries in each frequency bin individually and add the rest as a stochastic
background. This was shown to produce equivalent results to a full simulation,
but with drastically reduced runtimes (Bécsy et al., 2022). We produce
hundreds of realizations of these simulations by drawing new binary
populations to account for cosmic variance, and randomizing their extrinsic
parameters (sky location, inclination, polarization angle, phases). We then
analyze them with both Bayesian and frequentist statistical methods, which are
described in the next sections.

###  2.2 Bayesian methods

To carry out a Bayesian statistical analysis of our datasets, we use the
following likelihood function (for more details see, e.g., Eq. (7.35) in
Taylor 2021):

$\log L=-\frac{1}{2}\left[\delta t^{T}C^{-1}\delta t+\log\det(2\pi C)\right],$
(1)

where $\delta t$ are the timing residuals, and $C$ is the noise covariance
matrix that includes the effects of both white and red noise and the
linearized timing model. The matrix $C$ depends on parameters describing the
noise. In our analysis we fix the white noise parameters to their true values.
We model the spectrum of intrinsic pulsar red noise with a power-law model, so
that its power spectral density (PSD) in the $i$th pulsar is:

$S_{i}(f)=\frac{A_{i}^{2}}{12\pi^{2}}\left(\frac{f}{f_{\rm
yr}}\right)^{-\gamma_{i}}\ f_{\rm yr}^{-3},$ (2)

where $A_{i}$ and $\gamma_{i}$ are the amplitude and the spectral slope of the
RN in the $i$th pulsar, respectively, and $f_{\rm yr}=({\rm yr})^{-1}$ is the
frequency corresponding to a period of a year. We marginalize over $A_{i}$ and
$\gamma_{i}$ in our analysis. Similarly, we model the timing residuals induced
by the GWB as a red noise process with a power-law spectrum, such that its
cross-power spectral density between the $i$th and $j$th pulsar is:

$S_{ij}(f)=\Gamma_{ij}\frac{A_{\rm GWB}^{2}}{12\pi^{2}}\left(\frac{f}{f_{\rm
yr}}\right)^{-\gamma_{\rm GWB}}\ f_{\rm yr}^{-3},$ (3)

where $\Gamma_{ij}$ is the overlap reduction function (ORF) between the two
pulsars. For a GWB, the ORF is given by the HD curve (Hellings & Downs, 1983),
which introduces correlations between pulsars depending on their angular
separation on the sky. Another useful model is a so called common uncorrelated
red noise (CURN), which assumes $\Gamma_{ij}=\delta_{ij}$, so no correlation
between pulsars. This model is much more efficient to calculate, because $C$
is block-diagonal. In practice, the posteriors on $A_{\rm GWB}$ and
$\gamma_{\rm GWB}$ are similar between the CURN and HD models. This allows one
to sample the simpler CURN model and reweight the posterior samples to get
fair draws from the HD posterior (for details see Hourihane et al. 2022). The
mean of the weights also serves as a measurement of the Bayes factor between
the HD and CURN models, which is the primary quantity one needs to measure to
claim a detection of the GWB.

###  2.3 Frequentist methods

While GWB searches usually focus on Bayesian methods, frequentist methods can
also be useful, particularly given their significantly faster computational
speed. The most common frequentist method used in searches for a GWB is the so
called optimal statistic (OS, Chamberlin et al. 2015), which is an unbiased
estimator of the GWB amplitude. It is constructed by matching the expected ORF
to the measured cross-correlation values under the assumption of a (typically
power-law) PSD model. This results in the estimated GWB amplitude ($\hat{A}$)
and its standard deviation ($\sigma_{A}$). The ratio of these gives the
signal-to-noise ratio, S/N, which can be used as a detection statistic for the
presence of a GWB.

One disadvantage of the OS is that one needs to know the PSD describing the
GWB and the intrinsic red noise of each pulsar, which is not known a priori in
a real dataset. One can measure the PSD using Bayesian methods, and use a
point estimate (usually the maximum a posteriori) to calculate the OS.
However, this does not take into account the uncertainty of the measured
spectrum. To do so, Vigeland et al. (2018) developed the so called noise-
marginalized OS, which calculates the OS for random draws of the spectrum from
a Bayesian posterior, thus producing a distribution of OS values.
Subsequently, Vallisneri et al. (2023) have also proposed a related method for
interpreting the OS in the context of posterior predictive checking.

The OS can also be used to calculate the S/N of processes described by
different ORFs. Two commonly calculated ORFs are the monopole and dipole,
which could arise e.g. from systematic clock errors and ephemeris errors,
respectively. One limitation of the standard OS, is that since these ORFs are
usually not orthogonal, processes with less preferred ORFs can still produce
high S/N values. To alleviate this problem, Sardesai & Vigeland (2023)
developed the so called multi-component OS, which can allow multiple processes
with different ORFs to fit the data simultaneously. As a result, different
processes can be better decoupled. In this study we calculate S/N values with
the noise-marginalized multi-component OS method to assess the frequentist
significance of the GWB, along with monopole or dipole correlated processes.

##  3 Results

To produce robust results, we analyzed 773 different realizations of our
simulated 15-yr dataset. The median (thick lines) and mean (thin lines) GWB
spectra over these realizations is shown in Figure 1, both for all the
binaries (red solid lines) and for all except the brightest binaries in each
bin (green dashed lines). The shaded region represents amplitudes between the
5th and 95th percentile in each bin. We also show a spectrum with
$\gamma=13/3$ as a reference. Note that this has a slope between the median
and mean lines for all sources, which give $\gamma=4.6$ and $\gamma=4.2$ in a
linear fit to the first five frequency bins, respectively. We can also see
that the median spectrum is similar to the mean spectrum once the brightest
binary is removed. However, for the total spectrum, the mean spectrum is
shallower and follows the expected 13/3 spectral slope better. This is due to
the fact that the median is insensitive to outliers, and thus we lose the
information about the bright sources.

Figure 1: Mean and median GWB spectra over 773 realizations. Also shown is the
spectrum after removing the brightest binary’s contribution in each frequency
bin. The shaded regions show the range of spectra between the 5th and 95th
percentile. We also show a power-law spectrum with a canonical 13/3 slope as a
reference.

###  3.1 Stochastic background

Figure 2: Median recovered $\log_{10}A_{\rm GWB}$ (referenced at $(10\ {\rm
yr})^{-1}$) and $\gamma_{\rm GWB}$ values from 14-frequency CURN analysis of
773 realizations (blue). Also shown are the values we get from a simple least-
squares fit to the first 5 (orange) and first 14 (green) frequency bins of the
simulated spectra, and the parameters reported in Agazie et al. (2023a) based
on the NANOGrav 15-year dataset (red). Contours represent 99.7% levels
(3$\sigma$), and dots show the medians. 1-dimensional marginal distributions
are also shown, with 95.4% intervals (2$\sigma$) highlighted for some of the
distributions. While the NANOGrav results show a $\gamma_{\rm GWB}$ value
lower than the $13/3$ value expected from circular GW-driven GWB (gray), it is
well within the spread of $\gamma_{\rm GWB}$ values found in our simulated
datasets.

We carry out a Bayesian search for a 14-frequency CURN process in each
realization. Figure 2 shows the distribution of recovered $A_{\rm GWB}$ and
$\gamma_{\rm GWB}$ values (medians of posterior samples, blue contours) along
with the ones estimated by a simple least squares fit to the theoretical
spectrum using either the first 5 (orange) or the first 14 (green) frequency
bins. Interestingly, fitting to the first 5 frequencies gives very similar
values to the actual Bayesian CURN analysis, even though that uses the first
14 frequencies. This is due to the fact that the Bayesian analysis is
dominated by the first few frequencies, where the signal is the strongest.
Conversely, a 14-frequency fit to the spectra provides unrealistically
accurate measurements, as it does not take into account the effect of red and
white noise, the latter of which dominates the GWB signal at higher
frequencies. Note that all of these distributions center roughly around the
theoretically expected $\gamma_{\rm GWB}=13/3$ value, and the amplitude to
which this population was calibrated. However, they show a significant spread
indicating that the lower $\gamma_{\rm GWB}$ value and higher amplitude found
in the NANOGrav 15-year GWB search (red contours) is compatible with the
binary population in our simulations. Note that this is consistent with the
findings of Agazie et al. (2023a), where they found that the NANOGrav results
are compatible with an SMBHB interpretation based on a comparison between the
recovered $A_{\rm GWB}$ and $\gamma_{\rm GWB}$ values and fits to SMBHB
population models.222Note that Agazie et al. (2023b) carried out much more
rigorous analyses, which also showed that the GWB is consistent with an SMBHB
interpretation. As we can see, Bayesian spectral analysis of population-based
simulated datasets yield similar conclusions. This is not surprising given
that we see a simple fit to the spectrum gives very similar results to a full
analysis.

Figure 3: Distribution of HD vs. CURN Bayes factors for all 773 realizations
(blue) and for 14 realizations with a recovered amplitude consistent with the
amplitude recovered in the NANOGrav 15-year analysis at the 90% confidence
level (purple). Also shown are the means of these distributions (dashed
vertical lines), and the Bayes factor reported in Agazie et al. (2023a) based
on the NANOGrav 15-year dataset (red). While this particular astrophysical
model predicts amplitudes and Bayes factors that tend to be lower than the one
found in the NANOGrav 15-year dataset, if we select a subset of these that
produce a consistent amplitude, they also show good consistency in terms of
Bayes factor.

For each realization, we also produce HD posteriors using the reweighting
method (Hourihane et al., 2022). We calculate the Bayes factor between the HD
and CURN models which are shown as the blue histogram in Figure 3. Note the
large variance of the Bayes factor, even though we fixed the astrophysical
parameters of the simulation (more than eight orders of magnitude between the
lowest and highest Bayes factor). To understand the source of this variance,
it is worth looking at the recovered Bayes factors as a function of the
amplitude of the GWB. Figure 4 shows the Bayes factors as a function of
$\log_{10}A_{0.1yr^{-1}}$, i.e. the amplitude of the GWB referenced at a
frequency of $0.1yr^{-1}$ instead of the usual $1yr^{-1}$. We use this
reference point, since there the amplitude is less covariant with $\gamma_{\rm
GWB}$. We can see that while the Bayes factor is correlated with the
amplitude, even at a given fixed amplitude we see significant scatter in the
Bayes factor. Thus even for universes described by the same astrophysics and a
GWB with the same amplitude, the recovered Bayes factor shows a significant
variance. This is consistent with previous findings on both astrophysical and
simple power-law backgrounds (Cornish & Sampson, 2016; Hourihane et al.,
2022). Thus this is a generic property of any stochastic background, not
specifically of astrophysical backgrounds.

Figure 4: HD vs. CURN Bayes factors as a function of the recovered GWB
amplitude referenced at $f=0.1{\rm yr}^{-1}$ (blue dots). While there is
correlation between the amplitude and the Bayes factor, the Bayes factors show
significant variance even at a given amplitude. Also shown is the Bayes factor
and amplitude reported in Agazie et al. (2023a) based on the NANOGrav 15-year
dataset (red).

Figure 3 also shows the histogram of Bayes factors for realizations that show
an amplitude consistent at the 90% level with that reported in the NANOGrav
15-year GWB search (Agazie et al., 2023a), along with the actual measured
Bayes factor in that search. We can see that the simulations show about four
orders of magnitude scatter even in this narrow amplitude range, but they
predict Bayes factors that are roughly consistent with the real data. Note
that this is not a rigorous consistency check of the 15-year results with the
astrophysical predictions, since our astrophysical models tend to produce a
lower amplitude than the one found in the 15-year dataset. Thus, the mean of
our simulations produces a lower amplitude and lower Bayes factor than found
in the 15-year data. However, once we filter out only those realizations that
produce an amplitude close to the one found in the real data, we find that the
Bayes factors are consistent with the one measured on real data. This stands
to show that our simplified analysis methods are capable of detecting the GWB
with a significance reported in (Agazie et al., 2023a) even if the GWB is
coming from a finite number of binaries.

Figure 5: Distribution of mean OS S/Ns, HD vs. CURN Bayes factors, and
recovered GWB amplitude referenced at $f=0.1{\rm yr}^{-1}$. Grey lines
represent zero S/Ns and unit Bayes factors, and red lines represent the values
reported in Agazie et al. (2023a) based on the NANOGrav 15-year dataset. We
also indicate the Pearson correlation coefficient for every pair where the
associated p-value is below $10^{-3}$.

We also calculate the noise marginalized multi-component optimal statistic for
all realizations. Figure 5 shows the distribution of the mean monopole,
dipole, and HD S/Ns, along with the HD vs. CURN Bayes factors and
$\log_{10}A_{0.1yr^{-1}}$. We also indicate the values found in the NANOGrav
15-year dataset on Figure 5 (red markers). Note that these are all roughly
consistent with our simulations. Black lines represent zero S/Ns and Bayes
factors. Note that as expected both monopole and dipole S/N values are
centered around zero, while the HD S/N distribution is offset towards positive
values.

We can see some interesting correlations between values shown on Figure 5. To
quantify these, we calculate the Pearson correlation coefficient ($r$) for
each parameter pair and show them on the figure where the correlation is
significant (p-value $<10^{-3}$). The strongest correlation of $r=0.82$ can be
observed between HD S/N and HD vs. CURN BF. This is reassuring, since these
are both meant to quantify the significance of the presence of HD cross-
correlations in the data. In fact, Vallisneri et al. (2023) gives an
approximate mapping between these two quantities. Both the HD S/N and HD vs.
CURN BF values are also correlated with $\log_{10}A_{0.1yr^{-1}}$, with
$r=0.39$ and $r=0.43$, respectively. These relatively low correlation values
are in line with our finding above, that a large variation remains in
detection statistics even at fixed amplitude. The fact that the amplitude is
more correlated with the Bayes factor than with the HD S/N can potentially be
explained by the the fact that the Bayesian analysis uses uses both auto-
correlation and cross-correlation information, while the OS relies solely on
cross-correlations. Additionally, both monopole and HD S/N is anticorrelated
with the dipole S/N. However, there seems to be no correlation between the HD
and monopole S/N.

###  3.2 Individual binaries

Figure 6: Distribution of GWB spectrum parameters ($\gamma_{\rm GWB}$,
$\log_{10}A_{0.1yr^{-1}}$) based on CURN analysis, and S/N and $\log_{10}f$ of
the highest-S/N CW source in each realization. Red lines represent the values
reported in Agazie et al. (2023a) based on the NANOGrav 15-year dataset. We
also indicate the Pearson correlation coefficient for every pair where the
associated p-value is below $10^{-3}$.

We investigate detection prospects of individual binaries and how they relate
to background properties by following Bécsy et al. (2022): we calculate the
expected S/N of binaries defined as $\mathrm{S/N}=\sqrt{(s|s)}$, where $s$ is
the simulated CW signal, i.e. we calculate the inner product of the waveform
with itself. We do so for each binary with randomly assigned extrinsic
parameters (sky location, inclination angle, polarization angle, phases), and
we find the binary with the highest S/N in each realization. Figure 6 shows
the S/N and GW frequency of that best binary in each realization, along with
the $\log_{10}A_{0.1yr^{-1}}$ and $\gamma_{\rm GWB}$ values we find in our
CURN analysis. Similarly to Bécsy et al. (2022) we find that the the
highest-S/N sources tend to be found at moderate frequencies around 2-10 nHz.
We find a higher mean S/N, which is not surprising given the inclusion of
additional pulsars and updated noise models of NANOGrav 15-year pulsars.

We calculate the Pearson correlation coefficient ($r$) for each parameter pair
and show it on Figure 6 where the correlation is significant (p-value
$<10^{-3}$). The most significant correlation is shown between the CW S/N and
$\log_{10}A_{0.1yr^{-1}}$, with $r=0.51$. This can be attributed to the fact
that the CURN analysis does not include a CW model, so if a CW is present, the
recovered amplitude is expected to be biased high. We also find significant
anti-correlation between the CW S/N and $\gamma_{\rm GWB}$ ($r=-0.26$). This
is probably due to the fact that $\gamma_{\rm GWB}$ can be biased low by a
significant individual binary at moderate frequencies as we will see below.
Finally, the only parameter that correlates with the frequency of the loudest
binary is $\gamma_{\rm GWB}$, showing a negative correlation with $r=-0.34$.
This is understandable by the fact that the higher the binary’s frequency is,
the more it can bias the GWB slope low.

To further investigate the covariance between individual binaries and the GWB,
we analyzed a single realization with the model of an individual binary and
CURN (CURN+CW). For this analysis we used QuickCW (Bécsy et al., 2022, 2023),
a software package that builds on the enterprise (Ellis et al., 2019) and
enterprise_extensions (Taylor et al., 2021) libraries, but uses a custom
likelihood calculation and a Markov chain Monte Carlo (MCMC) sampler tailored
to the search for individual binaries. We selected this realization randomly
under the constraints that it has a CW S/N$>$5 and a monopole S/N$>$1.5. The
selected realization has a monopole S/N of 1.8, a dipole S/N of $-1.33$, an HD
S/N of 2.29, an HD vs CURN BF of 77, a $\log_{10}A_{0.1yr^{-1}}=-14.4$, a CW
SNR of 6.2, and a CW frequency of 9 nHz.

Analyzing this dataset with the CURN+CW model we found that the CW signal was
clearly detected, with a CURN+CW vs CURN BF of $\sim$25\. Moreover, the
parameter recovery of the CURN process was significantly affected by the
inclusion of the CW model. Figure 7 shows the distribution of $\gamma_{\rm
GWB}$ and $\log_{10}A$ under the CURN (red), HD (orange), CURN+CW (black), and
HD+CW (blue) models. Posteriors with models including HD were produced by
reweighting (Hourihane et al., 2022). We can see that in the case of the CURN-
only and HD-only models, the unmodeled binary biases the $\gamma_{\rm GWB}$
recovery low compared to the canonical 13/3 value, while the CURN+CW and HD+CW
models produce $\gamma_{\rm GWB}$ posteriors centered on values close to 13/3.
This is understandable, as the CURN model can partially model the CW signal by
introducing more power at higher frequencies where the CW is present. This is
also consistent with the correlation between $\gamma_{\rm GWB}$ and CW
properties we found on Figure 6. While we cannot draw definitive conclusions
from a single realization, this result suggests that modeling individual
binaries along with the GWB may be crucial for unbiased spectral
characterization of the GWB. We leave a detailed investigation of this problem
to a future study.

Figure 7: Distribution of GWB parameters under the CURN (red), HD (orange),
CURN+CW (black), and HD+CW (blue) models for a particular realization with a
clearly detectable binary ($\mathrm{S/N}=6.2$). Not modeling the CW signal
biases the spectral slope low, while once we include it we recover the
expected $\gamma=13/3$ value. Note that in this particular example, the CURN
and HD models show very similar posteriors, while the spectral recovery is
slightly different between the CURN+CW and HD+CW models.

It is important to note that while the NANOGrav 15-year analysis recovers a
$\gamma_{\rm GWB}$ value lower than 13/3, that does not significantly change
with the inclusion of an individual binary model (Agazie et al., 2023f), and
the $\gamma_{\rm GWB}$ recovery is instead more affected by the particular
noise models employed (Agazie et al., 2023a). This is unlike the EPTA DR2
analysis, where the inclusion of an individual binary completely changes the
GWB recovery (Antoniadis et al., 2023b).

Note that in this particular realization, the recovered spectra under the CURN
and HD models are practically identical (red and orange on Figure 7). This is
not particularly surprising, as this forms the basis of the reweighting
technique (Hourihane et al., 2022). In the presence of an individual binary,
the inclusion of correlations does have a small effect on the recovered
spectrum (compare black and blue on Figure 7), but the difference is small
enough that the reweighting technique can be effective.

In addition, the significance of the individual binary does not change
significantly after reweighting. We show the distribution of some individual
binary model parameters on Figure 8. Black shows posteriors from the CURN+CW
run, while blue shows those reweighted to HD+CW. Red indicates the true
parameters of the signal. In both cases, the binary is clearly detected. The
reweighted posterior distributions are similar to the original ones, but are
less constrained. We can also see that while the CURN+CW posteriors fail to
recover some of the model parameters, the HD+CW posteriors are generally
consistent with those. While this needs to be studied in more detail in a
systematic injection study, these results are promising for being able to
extend the reweighting technique to models including individual binaries.

Figure 8: Distribution of CW parameters for a particular realization from a
CURN+CW analysis (black) and those reweighted to HD+CW (blue). Red lines
represent the true parameter values.

##  4 Conclusion and future work

We used state-of-the art statistical data analysis techniques to search for a
stochastic GW background in realistic simulated datasets. We employed both
Bayesian and frequentist techniques to estimate the parameters of the
background and assess its significance. Our simulated datasets were produced
based on an astrophysical population of supermassive black hole binaries, and
as such the resulting stochastic backgrounds were anisotropic, non-Gaussian,
and had non-power-law spectra. We found that standard analysis techniques were
able to detect the characteristic Hellings-Downs correlations in these
backgrounds, despite the fact that they assume isotropic, Gaussian, power-law
backgrounds.

These analyses were also able to correctly characterize the spectrum of these
backgrounds most of the time. However, we found the presence of a strong
individual binary can bias the recovered amplitude high, and the recovered
spectral slope low. We calculated the signal-to-noise ratio of the loudest
individual binary in each realization and examined how that relates to the
background properties.

The fact that standard statistical analysis techniques used by the PTA
community can indeed detect the signals from a realistic SMBHB population is
reassuring. On the other hand, the possibility that individual binaries can
bias spectral recovery is cause for concern, and suggests that joint analyses
of individual binaries and the background might become necessary to avoid
biased inference. Future studies investigating the interplay between these two
source types will be crucial in making sure that any astrophysical
interpretation is robust. We also plan to investigate the expected level of
anisotropy in the stochastic background based on these realistic simulations
using methods described e.g. in Agazie et al. (2023d) and Gardiner & Kelley
(2023). While the results presented in this paper are based on a fixed
astrophysical model, investigating the effect of astrophysical parameters on
the resulting background and individual binary prospects is also being
investigated (Gardiner & Kelley, 2023).

##  5 Acknowledgements

We appreciate the support of the NSF Physics Frontiers Center Award
PFC-1430284 and the NSF Physics Frontiers Center Award PFC-2020265.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. M.C. 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. Pulsar research at UBC is supported by an NSERC Discovery
Grant and by CIFAR. 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 number 80GSFC21M0002. G.E.F., S.C.S., and S.J.V. are supported by
NSF award PHY-2011772. 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. The work of N.La. and X.S. is
partly supported by the George and Hannah Bolinger Memorial Fund in the
College of Science at Oregon State University. 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.M. are supported by NSF #1458952. M.A.M. 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. 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. S.R.T.
acknowledges support from an NSF CAREER award #2146016. 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.

## References

  * Afzal et al. (2023) Afzal, A., Agazie, G., Anumarlapudi, A., et al. 2023, ApJ, 951, L11, doi: 10.3847/2041-8213/acdc91
  * Agazie et al. (2023a) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.16220, doi: 10.48550/arXiv.2306.16220
  * Agazie et al. (2023c) Agazie, G., et al. 2023c, ApJ
  * Agazie et al. (2023d) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023d, arXiv e-prints, arXiv:2306.16221, doi: 10.48550/arXiv.2306.16221
  * Agazie et al. (2023e) Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023e, ApJ, 951, L9, doi: 10.3847/2041-8213/acda9a
  * Agazie et al. (2023f) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023f, ApJ, 951, L50, doi: 10.3847/2041-8213/ace18a
  * Ali-Haïmoud et al. (2020) Ali-Haïmoud, Y., Smith, T. L., & Mingarelli, C. M. F. 2020, Phys. Rev. D, 102, 122005, doi: 10.1103/PhysRevD.102.122005
  * Antoniadis et al. (2023a) Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023a, arXiv e-prints, arXiv:2306.16214, doi: 10.48550/arXiv.2306.16214
  * Antoniadis et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.16226, doi: 10.48550/arXiv.2306.16226
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bécsy et al. (2022) Bécsy, B., Cornish, N. J., & Digman, M. C. 2022, Phys. Rev. D, 105, 122003, doi: 10.1103/PhysRevD.105.122003
  * Bécsy et al. (2022) Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022, Astrophys. J., 941, 119, doi: 10.3847/1538-4357/aca1b2
  * Bécsy et al. (2023) Bécsy, B., Digman, M., & Cornish, N. J. 2023, QuickCW (v1.0.1).  https://github.com/nanograv/QuickCW/tree/v1.0.1
  * Chamberlin et al. (2015) Chamberlin, S. J., Creighton, J. D. E., Siemens, X., et al. 2015, Phys. Rev. D, 91, 044048, doi: 10.1103/PhysRevD.91.044048
  * Cornish & Sampson (2016) Cornish, N. J., & Sampson, L. 2016, Phys. Rev. D, 93, 104047, doi: 10.1103/PhysRevD.93.104047
  * Cornish & Sesana (2013) Cornish, N. J., & Sesana, A. 2013, Class. Quant. Grav., 30, 224005, doi: 10.1088/0264-9381/30/22/224005
  * Ellis et al. (2019) Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE.  http://ascl.net/1912.015
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Gair et al. (2015) Gair, J. R., Romano, J. D., & Taylor, S. R. 2015, Phys. Rev. D, 92, 102003, doi: 10.1103/PhysRevD.92.102003
  * Gardiner & Kelley (2023) Gardiner, E., & Kelley, L. Z. 2023, in preparation
  * Gersbach et al. (2023) Gersbach, K., Meyers, P., Taylor, S., Joseph, R., et al. 2023, in preparation
  * Hellings & Downs (1983) Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39, doi: 10.1086/183954
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655, doi: 10.1111/j.1365-2966.2006.10302.x
  * Hourihane et al. (2022) Hourihane, S., Meyers, P., Johnson, A., Chatziioannou, K., & Vallisneri, M. 2022, arXiv e-prints, arXiv:2212.06276, doi: 10.48550/arXiv.2212.06276
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Johnson et al. (2023) Johnson, A. D., Meyers, P. M., Baker, P. T., et al. 2023, arXiv e-prints, arXiv:2306.16223, doi: 10.48550/arXiv.2306.16223
  * Kelley et al. (2017a) Kelley, L. Z., Blecha, L., & Hernquist, L. 2017a, MNRAS, 464, 3131, doi: 10.1093/mnras/stw2452
  * Kelley et al. (2017b) Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017b, MNRAS, 471, 4508, doi: 10.1093/mnras/stx1638
  * Kelley et al. (2018) —. 2018, MNRAS, 477, 964, doi: 10.1093/mnras/sty689
  * Kelley et al. (2023) Kelley, L. Z., et al. 2023, in preparation
  * Lamb et al. (2023) Lamb, W. G., Taylor, S. R., & van Haasteren, R. 2023, arXiv e-prints, arXiv:2303.15442, doi: 10.48550/arXiv.2303.15442
  * Lentati et al. (2013) Lentati, L., Alexander, P., Hobson, M. P., et al. 2013, Phys. Rev. D, 87, 104021, doi: 10.1103/PhysRevD.87.104021
  * Luo et al. (2019) Luo, J., Ransom, S., Demorest, P., et al. 2019, PINT: High-precision pulsar timing analysis package, Astrophysics Source Code Library, record ascl:1902.007.  http://ascl.net/1902.007
  * Meyers et al. (2023) Meyers, P. M., Chatziioannou, K., Vallisneri, M., & Chua, A. J. K. 2023, arXiv e-prints, arXiv:2306.05559, doi: 10.48550/arXiv.2306.05559
  * Mingarelli et al. (2013) Mingarelli, C. M. F., Sidery, T., Mandel, I., & Vecchio, A. 2013, Phys. Rev. D, 88, 062005, doi: 10.1103/PhysRevD.88.062005
  * Nice et al. (2015) Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo: Pulsar timing data analysis, Astrophysics Source Code Library, record ascl:1509.002.  http://ascl.net/1509.002
  * Pol et al. (2022) Pol, N., Taylor, S. R., & Romano, J. D. 2022, ApJ, 940, 173, doi: 10.3847/1538-4357/ac9836
  * Pol et al. (2021) Pol, N. S., Taylor, S. R., Kelley, L. Z., et al. 2021, ApJ, 911, L34, doi: 10.3847/2041-8213/abf2c9
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Reardon et al. (2023) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6, doi: 10.3847/2041-8213/acdd02
  * Romano et al. (2015) Romano, J. D., Taylor, S. R., Cornish, N. J., et al. 2015, Phys. Rev. D, 92, 042003, doi: 10.1103/PhysRevD.92.042003
  * Sampson et al. (2015) Sampson, L., Cornish, N. J., & McWilliams, S. T. 2015, Phys. Rev. D, 91, 084055, doi: 10.1103/PhysRevD.91.084055
  * Sardesai & Vigeland (2023) Sardesai, S. C., & Vigeland, S. J. 2023, arXiv e-prints, arXiv:2303.09615, doi: 10.48550/arXiv.2303.09615
  * Taylor (2021) Taylor, S. R. 2021, arXiv e-prints, arXiv:2105.13270, doi: 10.48550/arXiv.2105.13270
  * Taylor et al. (2021) Taylor, S. R., Baker, P. T., Hazboun, J. S., Simon, J., & Vigeland, S. J. 2021, enterprise_extensions.  https://github.com/nanograv/enterprise_extensions
  * Taylor & Gair (2013) Taylor, S. R., & Gair, J. R. 2013, Phys. Rev. D, 88, 084001, doi: 10.1103/PhysRevD.88.084001
  * Taylor et al. (2017) Taylor, S. R., Simon, J., & Sampson, L. 2017, Phys. Rev. Lett., 118, 181102, doi: 10.1103/PhysRevLett.118.181102
  * Taylor et al. (2020) Taylor, S. R., van Haasteren, R., & Sesana, A. 2020, Phys. Rev. D, 102, 084039, doi: 10.1103/PhysRevD.102.084039
  * Vallisneri (2020) Vallisneri, M. 2020, libstempo: Python wrapper for Tempo2, Astrophysics Source Code Library, record ascl:2002.017.  http://ascl.net/2002.017
  * Vallisneri et al. (2023) Vallisneri, M., Meyers, P. M., Chatziioannou, K., & Chua, A. J. K. 2023, arXiv e-prints, arXiv:2306.05558, doi: 10.48550/arXiv.2306.05558
  * van Haasteren & Vallisneri (2014) van Haasteren, R., & Vallisneri, M. 2014, Phys. Rev. D, 90, 104012, doi: 10.1103/PhysRevD.90.104012
  * Vigeland et al. (2018) Vigeland, S. J., Islo, K., Taylor, S. R., & Ellis, J. A. 2018, Phys. Rev. D, 98, 044003, doi: 10.1103/PhysRevD.98.044003
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Xu et al. (2023) Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024, doi: 10.1088/1674-4527/acdfa5