#  The NANOGrav 15-Year Data Set: Detector Characterization and Noise Budget

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  Bence Bécsy Department of Physics, Oregon State
University, Corvallis, OR 97331, 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
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  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  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
Neil J. Cornish Department of Physics, Montana State University, Bozeman, MT
59717, 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  Brendan Drachler  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  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  Lydia Guertin Department of Physics
and Astronomy, Haverford College, Haverford, PA 19041, 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  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  Luke Zoltan Kelley Department of
Astronomy, University of California, Berkeley, 501 Campbell Hall #3411,
Berkeley, CA 94720, 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 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  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
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 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  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
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 The NANOGrav Collaboration
<EMAIL_ADDRESS>
###### Abstract

Pulsar timing arrays (PTAs) are galactic-scale gravitational wave detectors.
Each individual arm, composed of a millisecond pulsar, a radio telescope, and
a kiloparsecs-long path, differs in its properties but, in aggregate, can be
used to extract low-frequency gravitational wave (GW) signals. We present a
noise and sensitivity analysis to accompany the NANOGrav 15-year data release
and associated papers, along with an in-depth introduction to PTA noise
models. As a first step in our analysis, we characterize each individual
pulsar data set with three types of white noise parameters and two red noise
parameters. These parameters, along with the timing model and, particularly, a
piecewise-constant model for the time-variable dispersion measure, determine
the sensitivity curve over the low-frequency GW band we are searching. We
tabulate information for all of the pulsars in this data release and present
some representative sensitivity curves. We then combine the individual pulsar
sensitivities using a signal-to-noise-ratio statistic to calculate the global
sensitivity of the PTA to a stochastic background of GWs, obtaining a minimum
noise characteristic strain of $7\times 10^{-15}$ at 5 nHz. A power law-
integrated analysis shows rough agreement with the amplitudes recovered in
NANOGrav’s 15-year GW background analysis. While our phenomenological noise
model does not model all known physical effects explicitly, it provides an
accurate characterization of the noise in the data while preserving
sensitivity to multiple classes of GW signals.

gravitational waves — pulsars: general

††facilities: Arecibo, GBT, VLA††software: ENTERPRISE (Ellis et al., 2019),
enterprise_extensions (Taylor et al., 2021), hasasia (Hazboun et al., 2019a),
libstempo (Vallisneri, 2020), matplotlib (Hunter, 2007), PTMCMC (Ellis & van
Haasteren, 2017)

##  1 Introduction

Millisecond pulsars (MSPs) are extremely stable rotators, with long-term
stability comparable to that of accurate atomic clocks (Matsakis et al., 1997;
Hobbs et al., 2012, 2020a). They are therefore uniquely well-suited
astronomical objects to probe the low-frequency window of gravitational waves
(GWs), from nanohertz to microhertz (Sazhin, 1978; Detweiler, 1979; Hellings &
Downs, 1983; Foster & Backer, 1990; Perera et al., 2018).

The most promising sources of low-frequency GWs are inspiralling supermassive
black hole binaries (SMBHBs). The superposition of the GWs produced by the
SMBHB population forms the stochastic GW background (GWB), and the most
massive, closest binaries could be resolved as continuous GWs produced by
individual sources (see Taylor, 2021, for review). The North American
Nanohertz Observatory for Gravitational Waves (NANOGrav; Alam et al., 2021a)
collaboration regularly monitors an array of highly stable MSPs spread across
the sky to achieve the sub-microsecond timing precision required for GW
detection. In a GWB search, a common long-time-correlated (with more power at
low frequency, or “red”) signal should be detected among the pulsars in the
array, combined with a quasi-quadrupolar signature in the angular correlations
between pulsar pairs predicted by General Relativity, the so-called Hellings-
Downs correlations (Hellings & Downs, 1983). NANOGrav reported the detection
of a common red signal in our 12.5-yr pulsar timing data set; however, the
observations were insufficient to exhibit the Hellings-Downs correlations
required to definitively associate this signature with GWs (Arzoumanian et
al., 2020). Similar GW-search results have subsequently been reported by other
PTA collaborations, including the European PTA (EPTA; Chen et al., 2021), the
Parkes PTA (PPTA; Goncharov et al., 2021a), and the International PTA (IPTA;
Antoniadis et al., 2022), a collaboration between NANOGrav, the EPTA, the
PPTA, and the recently formed Indian Pulsar Timing Array (InPTA; Tarafdar et
al., 2022).

In this series of papers presenting the NANOGrav 15-year data release, we
report statistically significant evidence for Hellings-Downs spatial
correlations in the timing residuals for an ensemble of 67 pulsars (see Agazie
et al. 2023b, hereafter NG15gwb). Because of the stochastic nature of the
signal, a careful noise analysis of our experiment is imperative. As we have
done in our previous papers (Arzoumanian et al., 2016, 2018, 2020), we
explicitly model both white (time-uncorrelated) and non-white (time-
correlated) noise components present in the detector and incorporate those
results in additional analysis steps. Unlike in our previous data releases, in
which this topic has been discussed in either the data presentation or GWB
search papers, we devote a separate paper to the noise analysis of our 15-year
data release (see Agazie et al. 2023a, hereafter NG15). This paper describes
and employs a noise analysis technique specifically designed for the data in
this release. Our approach has been developed over our last four data releases
and incorporates empirical expressions of noise sources discussed in Cordes &
Shannon (2010); Cordes (2013); Stinebring (2013); Demorest et al. (2013);
Arzoumanian et al. (2015, 2016); Lam et al. (2016, 2017). Furthermore, we
explicitly outline the structure of the covariance matrix we use since it
incorporates our knowledge of noise source variances and their relation, if
any, to each other.

The structure of the paper is as follows. In §2 we present the known
astrophysical sources of noise in our experiment. We present in §3 a
phenomenological noise model that characterizes the white noise (WN) and its
covariance as well as developing a simple spectral model for red noise (RN)
present in the data, which includes possible contributions from the GWB and
other sources. The covariance matrices that take this knowledge of the noise
and propagate it into later analyses are described in §4. Noise
characterization and its effect on our sensitivity to a GWB are summarized
through transmission functions and sensitivity curves in §5. Finally, we
consider possible future improvements to noise characterization and reduction
in §6. Notation used throughout the paper is listed in Table 3.

##  2 Astrophysical Impacts on Pulse Arrival Times

Pulsar timing at the sub-microsecond level relies on using MSPs with extremely
stable rotation rates as high-precision cosmic clocks. We measure times of
arrival (TOAs) of pulses from a pulsar by producing a summed pulse profile at
every observation epoch and convolving with a pulse profile template using a
matched-filter approach where the observed profile is a scaled and shifted
version of a perfect template shape with additive noise (e.g., Taylor, 1992).
In this way, we are able to track a reference longitude on the pulsar, an
arbitrarily chosen reference meridian. Measured TOAs are then subtracted from
those predicted by a timing model to calculate pulsar timing residuals. Many
known sources of noise both intrinsic and extrinsic to the pulsar impact the
measured TOAs and must be carefully modeled in order to detect the timing
signatures expected due to GWs in the residuals.

In order to keep this paper self-contained, we briefly discuss the numerous
astrophysical perturbations to the TOAs, the fundamental measured quantities
in our experiment. This is a deep, rich, and complicated subject. See Cordes &
Shannon (2010), Cordes (2013), Stinebring (2013), and Verbiest & Shaifullah
(2018) for more detail on all of these effects. Each subsection is followed by
one or a few labels, such as “chromatic,” “achromatic,” etc., describing the
noise. Chromatic (or chromaticity) refers to a dependence on the radio
frequency of the observations. Achromatic features lack this dependence.
Spatially-correlated noise processes show statistical cross-correlations that
may depend on the sky position of the pulsars or possess the same correlation
across all pulsars. These noise processes are particularly important to
understand since PTAs search for the stochastic GWB which is expected to be
spatially correlated according to the Hellings-Downs relationship.

In the following subsections noise sources are discussed roughly in order from
emission mechanisms at the pulsar, following the propagation order, through to
observation of radio pulses at the telescopes. See Figure 1 for a visual
summary of these effects and their physical location relative to the pulsar-
Earth lines of sight.

Figure 1: An illustration of the signal path from the pulsar to the data
products, highlighting a few relevant sources of noise. The pulsar emission
itself is subject to jitter (and in rare cases, glitches), and the source
model includes fits for spin and spindown, astrometric parameters, and
possible binary orbit models. Propagation through the interstellar medium
imposes pulse dispersion and scattering, which are both frequency- and time-
dependent. At the receiver, along with thermal noise, there are effects
related to clock corrections and solar system barycenter parameters, among
others, that need to be modeled in order to produce a data set that includes
pulse times of arrival and time-dependent dispersion measure estimates.

###  2.1 TOA Errors from the Pulsar Emission Process, including Jitter —
Varying Chromaticity

Pheneomenologically, pulse jitter comprises the stochastic phase variation of
single pulses with respect to a mean phase that is locked to the rotation of a
neutron star. It appears to be a ubiquitous feature of radio pulsars and
magnetars but is a much smaller (if not absent) effect in high-energy emission
of pulsars. Pulse jitter limits the fundamental precision of our TOAs at the
emission region itself. While the rotation of the neutron star represents the
clock mechanism for our experiment, the TOAs we measure correspond to radio
emission generated in the pulsar’s magnetosphere. Pulse shapes averaged over
many rotations are generally very stable, allowing for TOA measurements to
precisions $\mathcal{O}(10^{-3})$ the duration of pulse widths. However, the
shapes of individual pulses are not stable. This long-studied phenomenon is
known as pulse jitter, in which single pulse components can vary in phase and
amplitude, with phase variations of order the pulse width (Cordes & Downs,
1985; Lam et al., 2016; Parthasarathy et al., 2021). The net effect is that
any finite sample of pulses will slightly differ in shape compared to another
sample, creating a source of uncertainty in the measured phase of the spot of
longitude that should be uncorrelated between observing epochs and scale with
the number of observed pulses $N_{\rm p}$ as $\propto 1/\sqrt{N_{\rm p}}$.
Jitter is strongly correlated across radio frequencies, though since pulse
profiles vary in shape across frequency, the combined effect on timing due to
jitter can have some frequency dependence, which has been measured for a large
fraction of our MSPs (Lam et al., 2019; Hebel, 2022). However, jitter is
entirely uncorrelated from pulse to pulse, i.e., it acts as WN in the
residuals (Cordes & Downs, 1985).

Giant radio pulses (GRPs) are individual pulses that have been observed to
occur not periodically, show amplitudes many times larger than average and
amplitude distributions that follow a power-law. They are another phenomenon
that arises in a small subset of known pulsars. The millisecond pulsar
B1937+21 was the second emitter of GRPs detected after the Crab pulsar
(Staelin & Reifenstein, 1968; Wolszczan et al., 1984). GRPs occur in specific
ranges of pulse phase and are thus locked to the spin phase of the neutron
star, but their large amplitudes occur stochastically without any distinct
temporal pattern. In the case of B1937+21 its GRPs have pulse widths ranging
from micro- to nanoseconds with intensities that follow a power-law
distribution and show brightness temperatures ranging up to 1039 K for a GRP
of 15 ns pulse width (Sallmen & Backer, 1995; Cognard et al., 1996;
Kinkhabwala & Thorsett, 2000; Soglasnov et al., 2004). GRPs from B1937+21 have
been observed to have little effect on the timing precision of this pulsar.
Extensive studies performed by McKee et al. (2019) show that the GRPs from
B1937+21 provide a timing solution negligibly different from the general
timing solution of the pulsar. Exhaustive studies of lower-level GRP emission
have not been performed for the remaining NANOGrav pulsars though examination
of several individual sources is underway. A frequently discussed pulsar
characteristic for the occurrence of GRPs is a high magnetic field at the
light cylinder (Cognard et al., 1996). Currently, it is not clear whether this
characteristic is connected to the physics of GRPs. Further investigations
need to be carried out before corresponding interpretations can be made.
Another pulsar which also displays a high magnetic field at its light cylinder
and was reported to be the third GRP emitter is PSR B1821-24A (Romani &
Johnston, 2001). Albeit not observed by NANOGrav, it is observed by other
PTAs. No studies of the effect of its GRPs on its timing precision have been
published yet. Although caused by propagation effects in circum-pulsar
environments rather than the emission process itself, highly magnified pulses
were discovered in the MSP B1957+20 (Main et al., 2018) and have been observed
in a variety of eclipsing MSPs (see Lin et al. 2023). Their effect, if any, on
MSP timing precision is largely unexplored.

Mode changing is a phenomenon in which a pulsar’s average profile switches
between two or more discrete modes over a variety of sometimes-periodic
timescales, ranging from seconds to days. The most dramatic, and first-
identified, example of mode changing is the abrupt cessation of a pulsar’s
radio emission for several consecutive pulse periods, known as nulling and
first identified by Backer (1970). Nulling has been observed in a large number
of predominantly slow-period pulsars (Huguenin et al., 1970; Page, 1973; Hesse
& Wielebinski, 1974). Nulling has not been observed in any MSP, but is
difficult to detect due to the high sensitivity required to detect MSP single
pulses. However, mode-changing has been observed in three MSPs so far – PSRs
B1957+20 (Mahajan et al., 2018), J0621+1002 (Wang et al., 2021), and
J1909$-$3744 (Miles et al., 2022). Of the three, only PSR J1909$-$3744 is
observed by NANOGrav. The mode changes detected from PSR J1909–3744 are
characterized by a weak and strong mode that are offset in phase and occur on
a single-pulse timescale (Miles et al., 2022). Additionally, PSR J1713+0747
shows evidence of distinct modes of low-amplitude drifting subpulses (Liu et
al., 2016). For both pulsars, these variations should contribute similarly to
WN as pulse-to-pulse jitter. We have shown that noise on timescales of several
hours was shown to integrate down as WN (Dolch et al., 2014; Shapiro-Albert et
al., 2020). Mode-changing would likely be detectable in a larger sample of the
MSPs we observe with sufficient sensitivity. A detailed investigation of the
impact of the nulling/mode changing over the timespans we observe is ongoing.

###  2.2 Rotational Irregularities (Spin Noise) — Achromatic

Spin noise results from rotational irregularities in the pulsars. The
interiors of MSPs are thought to consist of a differentially-rotating
superfluid core surrounded by an iron-rich crust (Langlois et al., 1998).
Rotational irregularities are thought to arise from torques in the pulsar
magnetosphere (Cheng, 1987; Kramer et al., 2006; Cordes & Shannon, 2008; Lyne
et al., 2010; Gao et al., 2016) or coupling between the crust and superfluid
core (Jones, 1990; Melatos & Link, 2014; van Eysden & Link, 2018). The result
is a stochastic variation in the pulse period detectable over timescales of
many years. Fortunately, spin noise is positively correlated with spin
frequency derivative, with root-mean-square contribution of $\sigma_{\rm
spin}\sim\dot{f}_{\rm s}^{1.0}$ over long timescales (Shannon & Cordes, 2010;
Lam et al., 2017; Lower et al., 2020). Since MSPs have very small frequency
derivatives, they typically show negligible spin noise. One counterexample is
the NANOGrav pulsar B1937+21 (Shannon & Cordes, 2010). Beyond rotational
irregularities from neutron-star-specific origins, another proposed cause of
the observed timing variations of this pulsar is an asteroid belt of mass
$\lesssim 0.05M_{\oplus}$ (Shannon et al., 2013).

Glitches  represent a special case of spin noise where the spin of a pulsar
experiences a step-change (see, e.g., Espinoza et al., 2011). Although this
usually results in a very small increase in spin frequency, examples of large
increases (“slow glitches” Shabanova, 2005) have been reported, or even
decreases (so-called “anti-glitches”) in X-ray pulsars and magnetars (Ray et
al., 2019; Archibald et al., 2013). The glitch rate and typical fractional
change of spin frequency show strong inverse correlations with pulsar age
(Fuentes et al., 2017). They are therefore rare in MSPs, although two examples
have been reported: first in the globular cluster pulsar B1821$-$24A (Cognard
& Backer, 2004), and later in the NANOGrav pulsar J0613$-$0200 (McKee et al.,
2016). The fractional changes in spin frequency for these two glitches are
only $8\pm 1\times 10^{-12}$ and $2.5\pm 1\times 10^{-12}$, respectively,
making these the two smallest glitches listed in the Jodrell Bank Glitch
Catalogue111https://www.jb.man.ac.uk/pulsar/glitches.html (Espinoza et al.,
2011). The low glitch amplitudes make glitches very hard to detect in MSPs,
with the PSR J0613$-$0200 glitch occurring in early 1998 but not being
reported until 2016. Note that this glitch occurred prior to the start of our
data set. However, other small and unrecognized glitches from this or other
MSPs could prove problematic for detection of some proposed GW targets of
PTAs, as a GW memory burst (affecting only the pulsar term) is predicted to
induce an apparent spin-frequency step-change in the timing residuals of a
single pulsar that would be indistinguishable from the signature of a small
glitch (Cordes & Jenet, 2012).

###  2.3 Profile Changes in Frequency and Time — Chromatic

Pulse shapes intrinsically vary as a function of frequency and, for the most
distant pulsars, due to multi-path propagation through the interstellar
medium, or scattering; see Section 2.5. We account for this effect in two
different ways. Our traditional narrowband timing method uses a single
template to fit pulses in a number of discrete frequency channels across a
given frequency band and then corrects for constant time offsets between
channels in the timing models with a fitted functional form of the offsets
versus frequency (Agazie et al., 2023a). Our “wideband” method (Pennucci et
al., 2014; Pennucci, 2019) uses a pulse portrait that contains information
about the frequency-dependent pulse shape to fit for the TOA and the time-
dependent dispersion measure, or DM (the integrated column density of free
electrons along the line of sight), delay simultaneously per epoch.

The TOA creation process assumes that the pulse profile of the pulsar is
constant and reproducible when integrating over a suitably long time (at least
thousands of pulses). However, subtle pulse shape variations, likely
attributable to either propagation effects or incorrect polarization
calibration, have been observed for several NANOGrav pulsars over long
timescales (Brook et al., 2018).

In addition, in early 2021 (after the timespan of our data set), the NANOGrav
pulsar J1713+0747 was found to have experienced a drastic change in its pulse
shape on a timescale of less than one day (Lam, 2021; Xu et al., 2021; Meyers
& Chime/Pulsar Collaboration, 2021; Singha et al., 2021a), before relaxing
back towards its original pulse shape over the course of several months
(Jennings et al., 2022a). As the signature in the timing residuals (generated
using templates derived from pre-event data) was similar to that expected from
interstellar medium (ISM) propagation effects (Lin et al., 2021), this shape
change was initially misinterpreted as being accompanied by a step-change in
DM. It was later demonstrated that, as the pulse shape evolution over
frequency was different before and after the event (Singha et al., 2021b), it
is not trivial to measure a change in DM between epochs (Jennings et al.,
2022b), and that the observed chromatic behavior is entirely consistent with a
frequency-dependent pulse shape change unaccompanied by changes in ISM
propagation. This event has necessitated re-examination of the previous
chromatic timing features in PSR J1713+0747 data noted by Lam (2021) as well
as one in PSR J1643$-$1224 (Shannon et al., 2016); it is not yet clear how
much of an impact this phenomenon has on our ability to model pulsar timing
behavior.

###  2.4 Orbital Irregularities — Achromatic

The majority of pulsars observed by NANOGrav are in binary systems, typically
with a white-dwarf companion (Fonseca et al., 2016; Agazie et al., 2023a). A
majority are in orbits well-described by Keplerian parameters, and also post-
Keplerian parameters in many cases. However, several pulsars require different
modeling of their binary systems, which we comment on here. Four (PSRs
J0023+0923, J0636+5128, J1705$-$1903, and J2214+3000) are in low-eccentricity
black-widow systems where tidal and wind effects can cause measurable orbital
variations. Three of these pulsars (excluding PSR J2214+3000) are modeled by
higher-order orbital frequency derivatives. While PSR J1705$-$1903 is newly
added to the NANOGrav PTA and several new higher-order orbital-frequency
derivatives have been measured in (NG15) for the others, the orbital
parameters for the other three measured in common with our 12.5-yr data set
(Alam et al., 2021a) remain stable with additional data (Agazie et al.,
2023a), suggesting excess noise from parameter mis-estimation is low. Only PSR
J1705$-$1903 shows significant RN, with a very shallow spectral index,
suggesting that mis-modeling does not significantly affect our GW analyses.
Note that NANOGrav selects for “well-behaved” systems from a timing
perspective, and so other such black-widow systems may indeed show
significantly more noise due to irregular angular-momentum transfer, though
since the orbital periods are all short compared to our GW signal ($<$1 day
versus years) only a small amount of sensitivity to GWs is lost overall
(Bochenek et al., 2015).

Beyond these systems, PSR J1024$-$0719 is a pulsar in a wide-binary system
(Kaplan et al., 2016; Bassa et al., 2016) such that we cannot feasibly measure
a complete orbit; instead we model the orbital motion as a second derivative
in the pulsar’s spin frequency. In this case, the pulsar shows significant RN
in individual noise and common noise analyses (see Table 2) and so unmodeled
orbital variations may still be contributing to the noise in this pulsar.
Lastly, there are two known pulsars in triple systems (Thorsett et al., 1999;
Ransom et al., 202314), one of which (PSR J0337+1715) is currently being timed
by NANOGrav due to its high-precision TOAs but is not included in NG15. The
timing of this pulsar requires higher-order effects from general relativity to
be included in its timing model and so the modeling procedure involves more
costly likelihood evaluations. Nonetheless, even though a complicated system,
only a small amount of excess noise has been measured beyond the template-
fitting uncertainties (Archibald et al., 2018), with no RN and a simplified DM
fit compared to our procedure (see next subsection), making triple systems
potentially significant contributors to future PTA efforts.

###  2.5 ISM Propagation Effects — Chromatic

Dispersion is the dominant frequency-dependent propagation effect in pulsar
timing resulting from radio pulses traveling through the ionized ISM,
interplanetary medium (IPM; or solar wind), and even Earth’s ionosphere, with
emission being temporally delayed as a function of frequency $\nu$ by an
amount proportional to $\mathrm{DM}/\nu^{2}$. The relative motions of the
Earth, pulsar, and ISM cause the sampled free electron content to vary,
requiring us to estimate a DM for every observing epoch. These dispersive
delays are also covariant with other propagation effects and the frequency-
dependence of the pulse shape itself. While dispersive delays are accounted
for in our timing models, other time-dependent propagation effects are not.

Multi-path propagation through the ISM causes extra time delays that may vary
differently from the $\nu^{-2}$ dispersive delay along with distortions to the
pulse shape that can cause mis-estimation of the TOAs in a variety of ways.
The most visible effect is scattering, manifesting as a broadening of the
pulse shape mathematically described as a convolution with a pulse broadening
function (PBF). The PBF has an approximately exponential tail, with the
broadening scaling approximately as $\nu^{-4.4}$, though the value of the
index is highly dependent on the physics and geometry of the medium. In the
low-scattering regime, the primary impact of scattering is to delay the pulse
with the roughly same frequency scaling (Hemberger & Stinebring, 2008).
However, several of the NANOGrav pulsars show prominent scattering tails,
especially those at high DM and/or those observed at our lowest radio
frequencies.

Another frequency-dependent effect results from the ray paths at different
radio frequencies traversing different sets of electrons, resulting in
slightly different DMs at different frequencies (Cordes et al., 2016). At
sufficiently high levels of precision, this will require a frequency-dependent
characterization of the DM at each epoch, currently outside the scope of our
analysis.

Refraction also results from multi-path propagation and the deflection of the
bulk set of rays through the ISM, causing the observed sky position of the
pulsar to vary. This results in an additional geometric delay to the pulses
proportional to $\nu^{-4}$, with an additional delay caused by the variation
in the angle-of-arrival proportional to $\nu^{-2}$ (Foster & Cordes, 1990),
the latter entirely covariant with the dispersive delay.

Scintillation arises because pulsar images are extremely angularly compact,
and the multiple ray paths interfere coherently with one another. The effect
is most easily visualized observationally as heavy modulation of the dynamic
spectrum, or the intensity of the pulsar as a function of frequency and time.
Normally referred to as diffractive interstellar scintillation (DISS), the
process at a particular epoch can be characterized by a coherence bandwidth,
$\Delta\nu_{\rm d}$, and a coherence timescale, $\Delta t_{\rm d}$, with both
quantities inversely proportional to pulsar distance. DISS changes the PBF
stochastically, resulting in a limit to the fundamental precision of our TOAs
from propagation effects (versus intrinsically at the neutron star for
jitter), with larger uncertainties at lower frequencies (see e.g., in Lam et
al., 2018a). When the total observation time, $T_{\rm obs}$, and total
bandwidth, $B$, are very large with respect to $\Delta t_{\rm d}$ and
$\Delta\nu_{\rm d}$, respectively, the effect of DISS on the timing is
reduced. However, most NANOGrav MSPs are nearby, leading to coherence
timescales of the same order or longer than our approximately 30-minute
observation lengths and coherence bandwidths only factors of a few smaller
than total observing bandwidths. This results in a “finite scintle” effect.
This is independent of pulse signal-to-noise (S/N) to first order (Cordes et
al., 1990; Lam et al., 2016), making the root-mean-square (RMS) timing noise
due to scintillation covariant with jitter in that respect. However, its
strong dependence on frequency allows it to be disentangled. In addition,
measurements of the scintillation parameters can be used as priors on the
analysis, though these are not applied to our current analyses.

###  2.6 Solar Wind Effects — Chromatic (Spatially Correlated)

The solar wind is a stream of charged particles which escape from the solar
corona due to their high ($\sim$1 keV) energies (Marsch, 2006). This leads to
an ionized IPM with an electron number density that decays with increasing
distance from the Sun. Throughout the Earth’s orbit the line of sight to a
given pulsar cuts through different parts of the IPM, resulting in an annual
contribution to the pulsar’s DM. The DM variations caused by the solar wind
are subsumed in the generic DM variation model, DMX, used by NANOGrav;
however, this excess DM can be modeled by assuming a spherically-symmetric
electron number density as a function of Sun-pulsar separation angle $\theta$
(Edwards et al., 2006)

$\text{DM}_{\text{sw}}(\theta)=(4.85\times
10^{-6}\textrm{pc})\,n_{0}\frac{\theta}{\sin\theta},$ (1)

where $n_{0}$ is the nominal electron number density at a distance of 1 au,
typically assumed to be several particles per cm3 in pulsar timing codes and
measured to be in the same range (e.g., Madison et al., 2019; Hazboun et al.,
2022). The maximum amplitude (i.e., at the smallest angular separation from
the Sun) of this periodic contribution to the DM depends on the pulsar’s angle
from the ecliptic plane, with pulsars very close to the ecliptic (latitude
$|\beta|\lesssim 10^{\circ}$) showing very sharp peaks around the minimum
solar separation (Jones et al., 2017; Donner et al., 2020). The solar wind
therefore represents both a chromatic and spatially-correlated signal (Tiburzi
et al., 2016) among pulsars in a timing array.

However, the solar wind is neither spherically symmetric nor static over the
course of the 11-year solar cycle (e.g., Issautier et al., 2001), making the
simple spherical model in Equation 1 inadequate to describe the solar wind
contribution to DM for pulsars at our DM precisions (Madison et al., 2019;
Tiburzi et al., 2019; Hazboun et al., 2022). The model also does not
differentiate between the densities of the ‘fast’ and ‘slow’ solar wind
(Tiburzi et al., 2021), leading to large annual changes in amplitude. As shown
in Hazboun et al. (2022) our current practice of measuring and removing the DM
delay on a per-observation basis allows us to account for the changing DMs
along each line of sight except in extreme circumstances where the Sun-pulsar
separation angle is extremely small and dual-receiver observations are not
available. NANOGrav removes data for which the DM variation is likely to be
too rapid to reliably represent it with DMX segments (see §4.1 in NG15 for
more details), but retains the raw data from observing pulsars near the cusp
of closest line of sight approach to the Sun for solar wind studies.

###  2.7 Solar System Ephemeris — Achromatic (Spatially Correlated)

TOAs are measured at an observatory and then must be transformed to the quasi-
inertial reference frame of the Solar System barycenter (SSB) (see e.g.,
Lorimer & Kramer, 2005 for full details). This is done through the use of a
planetary ephemeris, notably the Development Ephemeris (DE; Standish, 1982)
maintained by NASA’s Jet Propulsion Laboratory, most recently DE441 (Park et
al., 2021), and the Observatoire de Paris-maintained INPOP ephemerides (Fienga
et al., 2009), most recently INPOP19a (Fienga et al., 2019). The dominant
source of uncertainty in these ephemerides arises from inaccurate measurements
of planetary masses, particularly the outer planets, each of which contributes
a sinusoid to the timing residuals with period equal to the planetary orbit.
This signature will be spatially correlated among pulsars, roughly as a
dipolar signal, in a PTA (Tiburzi et al., 2016). The orbital periods of the
giant planets range from 11.9 – 164.8 yr, corresponding to frequencies of
$\sim 0.4$ – 2.7 nHz, similar to the GWB frequency range that is probed by
PTAs. The impact on PTA analyses is modeled by perturbing the orbital
parameters in the ephemeris used to derive PTA noise limits and comparing with
other models, as detailed in, e.g., Vallisneri et al. (2020) and Chen et al.
(2021). The dependence of the GW statistics in NANOGrav’s 11-yr data set on
the version of the ephemeris used was the catalyst for the development of
these methods (Arzoumanian et al., 2018). These models allow the analyses to
be bridged between different versions of the solar system ephemeris to obtain
equivalent results. Arzoumanian et al. (2020) and Vallisneri et al. (2020)
showed a diminishing dependency on the choice of ephemeris for PTA GWB results
and all tests of different ephemeris versions on NG15, including using a
Bayesian ephemeris model (Vallisneri et al., 2020), have shown insignificant
differences in parameter recovery. In addition, PTA data sets have been used
to identify non-GW common noise signals and have enabled limits to be placed
on the masses of the outer planets, large asteroids, and undiscovered planets
in the outer Solar System (Champion et al., 2010; Caballero et al., 2018; Guo
et al., 2019).

###  2.8 Clock Errors — Achromatic (Spatially Correlated)

Pulsar data are referred to an observatory time standard when they are
produced, but are transformed into Barycentric Dynamical Time (TDB) for timing
calculations and comparison across data sets (Luo et al., 2021; Hobbs &
Edwards, 2012). Observatory clocks are not perfect and so corrections are
applied at various stages of this transformation. See Luo et al. (2021) for a
detailed discussion of clock corrections and transformations used in pulsar
timing. Any errors in these transformations, or drifts between time standards,
can manifest in pulsar data as the same offset in all data sets with the clock
error (Hobbs et al., 2012; Miles et al., 2023). These offsets are spatially
correlated in the sense that all affected pulsars will be impacted by the same
shift at the same time. The spatial correlation function is therefore
monopolar, i.e., the same for all pulsar sky positions. See Tiburzi et al.
(2016) for a complete discussion of clock errors and other spatially
correlated noise processes in pulsar timing data.

###  2.9 Measurements at the Observatories — Varying Chromaticity

Radiometer noise from the observing systems used at our telescopes dominates
our measurement uncertainty for our lowest S/N pulsars but becomes a less
significant component at higher S/N; we are radiometer-noise-dominated for
most of our pulsars at most epochs. Random fluctuations due to radiometer
noise should be time and frequency independent. However, the system
temperature, and thus the S/N, is chromatic as it depends on frequency-
dependent contributions from the receiver bandpass, Galactic background, etc.
We have shown that the uncertainties from our template matching procedure
follow the expectations from radiometer noise for lower S/N (see Appendix B in
Arzoumanian et al. 2015 and note that we correct for this effect in our more
recent data sets) and deviate due to jitter and diffractive interstellar
scintillation at the higher S/Ns (Lam et al., 2016), regardless of the
frequency. Therefore, even though the system temperature will be frequency-
dependent, our per-channel modeling of the template-fitting error should
contain no known biases.

Incorrect polarization calibration could be a source of error in our data. For
every observation we correct for differential gain and phase variations in the
two chiralities of polarization of the received radio waves as well as the
changing parallactic angle. The latter effect is easily calculated based on
the source’s apparent position whereas the former is corrected for based on
measurements of a noise diode taken prior to each observation and an
unpolarized quasar roughly once per month by which we correct for variations
in the noise source power. Mis-calibration causes alterations to the
polarization-summed profiles used in timing (Stinebring et al., 1984) which
would then be a chromatic source of noise given the frequency dependence of
the (polarization) profiles. If polarization calibration is incomplete, for
example when a calibrator source is not observed at a particular epoch, we
would expect some contribution to chromatic WN.

Strong radio-frequency interference (RFI) is removed in several stages of our
pipeline as discussed in NG15. Low-level RFI will still remain and perturb the
TOA estimates in the template-matching procedure. RFI can take many forms –
remaining narrowband RFI will perturb individual frequency channels but
broadband RFI can perturb all of the TOAs across the band. RFI can also be
impulsive or periodic. Both can be removed from individual sub-integrations if
strong enough to identify. If not, in either case RFI will be folded at the
pulse period and diminished, but still affect the TOAs as a source of
possibly-chromatic WN.

##  3 Phenomenological Noise Model

Independently fitting for each and every source of noise identified in Section
2 across all pulsars in the PTA would be extremely challenging. In order to
efficiently characterize the noise of our detector, and allow noise modeling
from unknown sources, we therefore use a phenomenological model that takes
these effects into account while reducing the number of fit parameters
dramatically. The model has two main parts, distinguished by the timescale and
spectral characteristics of the noise being modeled. The first part models the
white, or uncorrelated in time, noise that has a “flat” contribution to the
power spectra of timing residuals. We model this noise using three parameters
which increase the uncertainties on the TOAs by accounting for WN unaccounted
for in the template fitting process (Arzoumanian et al., 2015). The second
part uses Gaussian process regression (Williams & Rasmussen, 2006) to model
low-frequency, time-correlated, or red, stochastic processes in the data. In
theory, both of these models could be included in the covariance matrix for
the TOA data, however, it is more numerically expeditious to separate out the
Gaussian process formalism for real data analyses. See Section 4 for details
about how this is accomplished.

###  3.1 White Noise Model

The uncertainties that are initially associated with pulsar TOAs $\sigma_{\rm
S/N}$ are due to the finite S/N of the matched-filtering process used to
calculate them (Taylor, 1992). A template of the pulse shape, built from many
observations, is convolved in the Fourier domain with the summed pulse profile
from a single observing session, under the assumption that the data comprise a
scaled and shifted version of the template added to WN (Lommen & Demorest,
2013). However, there are other sources of WN (e.g., jitter, scattering, RFI,
etc.; see Section 2) that do not adhere to the assumptions of matched
filtering (i.e., that the pulse shape is a copy of the template). These
various noise effects cumulatively induce variance in the TOAs, and it is
nontrivial to disentangle the distinct contributions of each noise source to
the total TOA uncertainties. We therefore employ an empirical WN model for
pulsar timing data that inflates the $\sigma_{\rm S/N}$ measured from the
pulse template-matching process using three parameters.

Three WN parameters are used to adjust the TOA uncertainties in order to
accurately reflect WN present in the data. This process yields a
reduced-$\chi^{2}$ near unity for the fit to the timing residuals, if the
timing model is complete and accurate. Various differences between pulsar
timing backends and radio observatory receivers make it necessary to give
different values of these WN parameters to each receiver/backend combination.
These parameters thereby encode the trustworthiness of TOAs from each
receiver/backend combination, down-weighting the TOAs from combinations where
effects in addition to template matching reduce the reliability of the data.
These three WN terms – EFAC ($\mathcal{F}$), EQUAD ($\mathcal{Q}$) and ECORR
($\mathcal{J}$) – come together with receiver/backend combination $re/be$
dependence as

$C_{ij}=\mathcal{F}^{2}(re/be)[\sigma_{{\rm
S/N},i}^{2}+\mathcal{Q}^{2}(re/be)]\;\delta_{ij}+\mathcal{J}^{2}(re/be)\;\mathcal{U}_{ij}$
(2)

where the $i,j$ denote TOA indices across all observing epochs, $\delta_{ij}$
is the Kronecker delta and we omit the dependence on receiver and backend,
$re/be$, from here on for simplicity. While EFAC and EQUAD only add to the
diagonal of $C$, where $C_{ij}$ are the elements of the covariance matrix to
be discussed in Section 4, the ECORR terms are block diagonal for single
observing epochs. ECORR is modeled using a block diagonal matrix,
$\mathcal{U}$, with values of $1$ for TOAs from the same observation and zeros
for all other entries.

Historically, the EFAC (error factor, $\mathcal{F}$) parameter was the first
WN parameter added to the pulsar timing covariance matrix. EFAC is a scale
factor on the $\sigma_{\rm S/N}$. The increase in uncertainties from EFAC
attempts to account for underestimated template matching errors from low S/N
ratio observing epochs and template mismatches due to pulse profile
variability. This parameter has tended towards $1$ as pulsar backends have
improved, more high dynamic range (i.e., 8 or more bit sampling) digital
systems have been implemented, and the treatment of low S/N TOAs has improved
(Alam et al., 2021a, b).

To include additional WN, the EQUAD ($\mathcal{Q}$) parameter222Note that
there are two distinct definitions of EQUAD in the literature, depending on
whether EFAC multiplies only $\sigma_{\rm S/N}$, referred to as the TempoNest
convention for the paper (Lentati et al., 2014) where it was first used, or
whether it multiplies the sum, in quadrature of $\sigma_{\rm S/N}$ and EQUAD.
All three of the main pulsar timing packages and the Enterprise software stack
use the latter convention, laid out in Equation 2. (first used in Nice &
Taylor 1995) is added in quadrature following the usual rules of noise
propagation. EQUAD encompasses additional sources of WN not accounted for in
the TOA uncertainties, and not modeled by EFAC.

The most recent WN parameter added to the pulsar covariance matrix is ECORR
(error correlated in _radio_ frequency, $\mathcal{J}$), which is also added in
quadrature and specifically models noise correlated across radio frequencies
within a single observing epoch333Here an observing epoch is effectively one
observation with a single receiver. Multiple observations within a single MJD
are not treated as correlated.. ECORR was first used in Arzoumanian et al.
(2015) to mitigate noise correlated across the narrowband TOAs obtained for a
single observing epoch. In part this parameter became necessary due to
NANOGrav’s data acquisition strategy of using multiple TOAs across the full
observing band for chromatic mitigation. While ECORR is not strictly “white
noise” from the perspective of being only diagonal in the covariance matrix,
its spectral contribution is effectively white, i.e., constant in the
frequency domain, for the frequencies important to GW searches. The
correlation timescale modeled is only across pulses emitted from the pulsar
within a relatively short ($\sim 30\;{\rm min}$) observation. ECORR partially
accounts for intra-band correlations caused by pulse jitter, but also includes
other short-timescale correlations across radio frequency, which can be
produced by, e.g., RFI and time-variable scattering (Lam et al., 2019;
Shapiro-Albert et al., 2021).

###  3.2 Red Noise Model

Various sources of noise in Section 2, such as variations in the electron
density of the ISM, instabilities in the spin of the pulsars, corrections to
Earth-based clock systems, corrections to the solar system ephemeris model,
and, of course, TOA shifts due to the stochastic GWB are time-correlated
across long timescales. The characterization of these effects in the frequency
domain, e.g., a power spectral density, shows that the effect has more power
at lower (redder) frequencies, hence RN. In most cases the theoretical models
for these sources of noise have a power spectral density $P(f)$ that follows a
power law with frequency $f$ Cordes & Shannon (2010), $P(f)\sim
A^{2}f^{-\gamma}$, with amplitude $A$ and where the “redness” is dictated by
the explicit minus sign in the exponent and $\gamma>0$. More complex models
have been investigated, such as turnover models, but the straight power law is
still considered an accurate and effective model for MSPs (Goncharov et al.,
2020).

There are two methods by which RN can be included in a pulsar timing analysis.
One way is to use the Wiener-Khinchin theorem to write the correlations
between different TOAs, $(t_{i},t_{j})$,

$C^{\rm RN}_{ij}=\int_{0}^{f_{\rm Nyq}}{\rm d}\\!f\cos[2\pi
f(t_{i}-t_{j})]P(f),$ (3)

where $C^{\rm RN}_{ij}$ are the elements of the covariance matrix as noted
above and $f_{\rm Nyq}$ is the Nyquist frequency. Note that, strictly
speaking, the integral in Equation 3 does not converge for all forms of
$P(f)$, in which case a low-frequency cutoff must be used, see van Haasteren
et al. (2009) for an exhaustive discussion. The details of the RN model are
discussed more in Section 4.

Alternatively, the perturbations of the TOAs due to RN can be modeled directly
with a Fourier basis and a set of coefficients, $F\mathbf{a}$, where $F$ is an
$(N_{\rm TOA}\times 2N_{\rm freq})$ Fourier design matrix, and $\mathbf{a}$ is
a $2N_{\rm freq}$ vector with 2 coefficients for each of $N_{\rm freq}$
frequencies included. As will be discussed in Section 4, the form of $P(f)$
can be dictated by imposing various functional forms on the $\mathbf{a}$
coefficients. This method is at the heart of PTA searches for the GWB (Lentati
et al., 2016, 2013; van Haasteren & Vallisneri, 2014) and achromatic (in radio
frequency) RN models used for individual pulsars. In practice, the methods of
Gaussian process regression are used in the analyses to accurately include the
stochastic nature of these signals, and the explicit Fourier basis-modeling
will eventually be incorporated into the covariance matrix, but in a different
form, see Section 4. It is these analyses that are used to fit for RN
parameters used in the covariance matrix for pulsar timing packages. While
NANOGrav uses the DMX model (NG15) in this analysis, these same Gaussian
methods are used by other PTA collaborations to model the variations in DM by
including a dependence on the radio frequency (Lentati et al., 2016).

##  4 The NANOGrav Pulsar Covariance Matrix

PTA data analysis is done using the timing residuals $\delta\mathbf{t}$,
produced by subtracting times of arrival predicted by a timing model
$\mathbf{t}_{M}$ from the TOAs $\mathbf{t}$:

$\delta\mathbf{t}=\mathbf{t}-\mathbf{t}_{M}=\mathbf{t}_{D}-\mathbf{t}_{M}+F\mathbf{a}+\mathbf{n},$
(4)

where the bolded symbols are column vectors. The difference between the actual
underlying, deterministic, non-GW delays and the timing model,
$\mathbf{t}_{D}-\mathbf{t}_{M}\equiv M\bm{\epsilon}$, is represented by a
linear-order Taylor series in the (assumed small) perturbations to the model
parameters, $\delta\mathbf{p}\equiv\bm{\epsilon}$, and the design matrix, $M$,
where

$M_{ij}\equiv\left(\frac{\partial t_{M,i}}{\partial
p_{j}}\right)_{\mathbf{p}=\mathbf{p}_{0}}$ (5)

is evaluated (usually analytically, numerically otherwise) at the best-fit
parameters, $\mathbf{p}_{0}$, from a linear least squares analysis. The use of
the linearized timing models was introduced in Ellis et al. (2013) and is
especially expeditious for full PTA analyses. However, full timing model fits
can also be done as part of the Bayesian searches (Lentati et al., 2014;
Vigeland & Vallisneri, 2014, Kaiser et al. in prep). The $F\mathbf{a}$ term
encodes the long-timescale correlated (red) noise, both for individual pulsar
data sets and the GWB, see Section 3. Finally, $\mathbf{n}$ is the WN
remaining in the residuals, assumed to have a multivariate Gaussian
distribution,

$\mathcal{P}(\mathbf{n})=\frac{\exp\left(-\frac{1}{2}\mathbf{n}^{\mathrm{T}}C^{-1}\mathbf{n}\right)}{\sqrt{\det\left(2\pi
C\right)}},$ (6)

with covariance
$C(\mathbf{E})\equiv\langle\mathbf{n},\mathbf{n}^{\mathrm{T}}\rangle$, which
is a function of the set of WN parameters, e.g.,
$\mathbf{E}\equiv(\mathcal{F},\mathcal{Q},\mathcal{J})$. We construct the
likelihood function for a PTA by combining the timing residuals, $C$, and
Equation 6:

$\mathcal{P}(\delta\mathbf{t}|\bm{\epsilon},\mathbf{a},\mathbf{E})=\frac{\exp\left(-\frac{1}{2}(\delta\mathbf{t}-M\bm{\epsilon}-F\mathbf{a})^{\mathrm{T}}C^{-1}(\delta\mathbf{t}-M\bm{\epsilon}-F\mathbf{a})\right)}{\sqrt{\det\left(2\pi
C\right)}}\;.$ (7)

When conducting a search for GWs in pulsar timing data we are not concerned
with the values of the timing parameter perturbations included in
$\bm{\epsilon}$. Using a Gaussian distribution for the prior probability
distribution of the parameter perturbations,

$\mathcal{P}(\bm{\epsilon})=\frac{\exp\left(-\frac{1}{2}\bm{\epsilon}^{\mathrm{T}}X^{-1}\bm{\epsilon}\right)}{\sqrt{\det\left(2\pi
X\right)}},$ (8)

and defining $X\equiv\langle\bm{\epsilon},\bm{\epsilon}^{\mathrm{T}}\rangle$,
one can marginalize over the timing model parameter perturbations (see
Appendix A.1)

$\displaystyle\mathcal{P}(\delta\mathbf{t}|\mathbf{a},\mathbf{E})$
$\displaystyle=$
$\displaystyle\int\mathcal{P}(\delta\mathbf{t}|\bm{\epsilon},\mathbf{a},\mathbf{E})\mathcal{P}(\bm{\epsilon}){\rm
d}(\bm{\epsilon})$ (9) $\displaystyle=$
$\displaystyle\frac{\exp\left(-\frac{1}{2}(\delta\mathbf{t}-F\mathbf{a})^{T}D^{-1}(\delta\mathbf{t}-F\mathbf{a})\right)}{\sqrt{\det\left(2\pi
D\right)}},$ (10)

where $D\equiv C+MXM^{\mathrm{T}}$ (Eq. A3).

In a similar fashion the Fourier basis coefficients are also marginalized over
assuming that the RN, whether intrinsic to the pulsar or a GWB, is a Gaussian
process,

$\mathcal{P}(\mathbf{a})=\frac{\exp\left(-\frac{1}{2}\mathbf{a}^{\mathrm{T}}\varphi^{-1}\mathbf{a}\right)}{\sqrt{\det\left(2\pi\varphi\right)}},$
(11)

where $\varphi\equiv\langle\mathbf{a},\mathbf{a}^{\mathrm{T}}\rangle$ is a
$(2N_{\rm freq}\times 2N_{\rm freq})$ diagonal matrix with dimensions dictated
by the number of frequencies, $N_{\rm freq}$, used in the analysis and
describing the variance of the RN coefficients. In principle these
coefficients could be used as free parameters, however they are usually
parameterized as a power law,

$\varphi_{i}=\frac{A^{2}}{12\pi^{2}}\left(\frac{f_{i}}{f_{\rm
yr}}\right)^{-\gamma}\;{\rm yr}^{3}{\rm
d}\\!f\sim\frac{A^{2}}{12\pi^{2}\;T_{\rm span}}\left(\frac{f_{i}}{f_{\rm
yr}}\right)^{-\gamma}\;{\rm yr}^{3},$ (12)

or some other functional form for the power spectral density. Note that
$f_{\rm yr}=1~{}\mathrm{yr}^{-1}$, and the $\sim$ assumes the usual constant
$\mathrm{d}f\sim 1/T_{\rm span}$, where $T_{\rm span}$ is the total timespan
(years) of the data set. Marginalizing over the $\mathbf{a}$ coefficients
gives the final form for a single pulsar likelihood,

$\displaystyle\mathcal{P}(\delta\mathbf{t}|\mathbf{E})$ $\displaystyle=$
$\displaystyle\int\mathcal{P}(\delta\mathbf{t}|\mathbf{a},\mathbf{E})\mathcal{P}(\mathbf{a}){\rm
d}(\mathbf{a})$ (13) $\displaystyle=$
$\displaystyle\frac{\exp\left(-\frac{1}{2}\delta\mathbf{t}^{\mathrm{T}}N^{-1}\delta\mathbf{t}\right)}{\sqrt{\det\left(2\pi
N\right)}},$

where $N\equiv D+F\varphi F^{\mathrm{T}}$, again see Appendix A.1 (Eq. A7).

In a single pulsar noise analysis using this likelihood, one searches over all
of the parameters $(A,\gamma,\mathbf{E})$; however, the WN parameters are
fixed to their maximum likelihood values during full PTA GW searches. The WN
parameters are independent of the parameters varied in the GW search, and the
necessary matrix inversion444For further mathematical discussion of this
inversion, see Appendix A.1. is prohibitively expensive from a computational
standpoint. The marginalizations above were developed to reduce this
computational expense to a single inversion of a $(2N_{\rm freq}\times 2N_{\rm
freq})$ matrix (van Haasteren & Vallisneri, 2014).

###  4.1 Phenomenological Covariance Matrix

As noted above, the phenomenological WN model is a function of the parameters
$\mathbf{E}\equiv(\mathcal{F},\mathcal{Q},\mathcal{J})$, which together
account for the noise processes the NANOGrav-observed pulsars; see Table 1.
Collectively, the WN in each of our pulsar timing models can be expressed as a
covariance matrix where each of $N_{\rm TOA}$ TOAs is specified with both its
recorded time $t$ and its radio frequency $\nu$. This covariance matrix is
then fit globally over all epochs and over all a pulsar’s frontend-backend
combinations. For a set of TOAs on a given epoch, we write the noise model in
the form

$C^{\mathrm{epoch}}=\begin{bmatrix}\mathcal{F}^{2}[\sigma_{\rm
S/N}^{2}(\nu_{1})+\mathcal{Q}^{2}]+\mathcal{J}^{2}&\mathcal{J}^{2}&\cdots&\mathcal{J}^{2}\\\
\mathcal{J}^{2}&\mathcal{F}^{2}[\sigma_{\rm
S/N}^{2}(\nu_{2})+\mathcal{Q}^{2}]+\mathcal{J}^{2}&\cdots&\mathcal{J}^{2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{J}^{2}&\mathcal{J}^{2}&\cdots&\mathcal{F}^{2}[\sigma_{\rm
S/N}^{2}(\nu_{N})+\mathcal{Q}^{2}]+\mathcal{J}^{2}\end{bmatrix},$ (14)

where $\mathcal{F}$, $\mathcal{Q}$, and $\mathcal{J}$ are the EFAC, EQUAD, and
ECORR described in Section 3. This matrix is block diagonal, with the entries
of $\mathcal{J}^{2}$ connecting frequency channels from the same observations
in blocks. The WN covariance matrix can then be written down as

$C=\begin{bmatrix}C^{\mathrm{epoch}}(t_{1})&0&\cdots&0\\\
0&C^{\mathrm{epoch}}(t_{2})&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&C^{\mathrm{epoch}}(t_{N})\end{bmatrix},$ (15)

where the noise is independent between epochs, but the noise parameters are
fit globally over all TOAs from a single pulsar. If one wants to include the
RN in the covariance matrix the total noise would then be $C\rightarrow
C+C^{\rm RN}$. This approach is used in pulsar timing software when the best-
fit amplitude and spectral index for RN in a given pulsar have already been
fit for in an earlier analysis (Coles et al., 2011).

It is tempting to try and map the phenomenological noise model to a signal
covariance matrix made up of all of the effects from Section 2. However, in
practice, there are a few reasons for not doing so in a GW analysis. First of
all, the complexity of these models makes this task unwieldy for GW searches,
especially given that some of these effects are not directly measurable.
Second, the noise measured through the phenomenological model is known to be
larger than the sum of the various physical effects, so a signal based
covariance matrix might not include all of the noise in our detector. In
Section A.2 we sketch a signal covariance matrix to show how some of the terms
from Section 2 would enter into the analysis.

##  5 PTA Noise Budget

###  5.1 Comparison of Physical and Phenomenological Models

The noise in Section 2 enters into the phenomenological noise model through a
diverse set of physical processes and any particular physical effect may enter
more than one part of the phenomenological model. From the standpoint of
detecting GWs with PTAs the most salient aspect of this complex relationship
is that the phenomenological model includes the known physical sources of
noise and is pliable enough to include unknown sources as well.

That pliability is demonstrated in Table 1 where we provide a list of the
various noise sources from Section 2, and relationships between the
phenomenological parameters in which they are included, along with a number of
other descriptors that can be used to sort the sources of noise.

Table 1: Connecting Noise Terms: Astrophysical noise sources, along with a summary of their attributes and connections to the phenomenological noise model, are listed. Short timescales are roughly the order of the observation length or less, while long timescales refer to periods between observation epochs. Noise Type | Symbol | Origin | Phenomenological Model Component | Timescale | Spatially Correlated
---|---|---|---|---|---
Radiometer | $\sigma_{\rm S/N}$ | Telescope | $\sigma_{\rm S/N}$ | Short |
Jitter | $\sigma_{\rm J}$ | Pulsar | $\mathcal{J}$ | Short |
Diffractive Interstellar Scintillation | $\sigma_{\rm DISS}$ | ISM |  $\mathcal{Q}$, $\mathcal{J}$ | Short |
DM Mis-estimation | $\delta\mathrm{DM}$ | ISM | RN/DMX | Long |
Solar Wind DM Mis-estimation | $\delta\mathrm{DM}_{\mathrm{sw}}$ | IPM | RN/DMX | Long | $\checkmark$
Frequency-dependent DM | $\mathrm{DM}(\nu)$ | ISM | RN/DMX | Long |
RFI |  | Telescope | $\mathcal{J}$ | Short |
Polarization Mis-calibration |  | Telescope | $\mathcal{J}$ | Both |
Scattering |  | ISM | $\mathcal{F},\mathcal{J}$ | Short |
Solar System Barycenter Mis-modeling |  | Solar System | RN | Long | $\checkmark$
Clock Errors |  | Telescope | RN | Long | $\checkmark$

####  5.1.1 DM Mis-estimation due to Asynchronous Measurements

In order to estimate DM on each epoch, we make TOA measurements across a wide
frequency range. We used two or more receiver bands at one or more telescopes
to cover the total range. Data were recorded with each receiver independently,
with observations occurring within 30 minutes of each other at Arecibo but
within several days at Green Bank due to efficiency constraints in changing
the receivers. DMs can change rapidly (i.e., on timescale of a day or less)
due to the passage of the line of sight through the inhomogeneous ISM, changes
in the solar wind, etc. In the presence of this stochastic variation, Lam et
al. (2015) estimated the effect of asynchronous measurement, finding an
induced perturbation on the infinite-frequency arrival time with a shallow
($\gamma=2/3$ for a Kolmogorov medium) RN spectral index. The RMS error is
given as (see Eq. 14 in Lam et al. 2015)

$\sigma_{\delta t_{\infty}}\approx
6.5~{}\mathrm{ns}\left(\frac{\nu}{\mathrm{GHz}}\right)^{-2}\left(\frac{r^{2}}{r^{2}-1}\right)\left(\frac{t_{\rm
sep}}{\Delta t_{\rm d,1~{}GHz}/1000~{}\mathrm{s}}\right)^{5/6}.$ (16)

This expression assumes a simplified model where pulsar observations are made
at only two spot frequencies $\nu$ and $\nu^{\prime}$, where
$r\equiv\nu/\nu^{\prime}$ with $\nu>\nu^{\prime}$, $t_{\rm sep}$ is the time
spacing between the observations, and $\Delta t_{\rm d,1~{}GHz}$ is the
scintillation timescale at a radio frequency of 1 GHz. For pulsars observed at
Green Bank with $r=1.5~{}\mathrm{GHz}/0.8~{}\mathrm{GHz}=1.875$, a time
between observations (in days) of $t_{\rm sep}\approx$ 2, and a typical
scintillation timescale of approximately 1000 s at 1 GHz, we see that the RMS
error from this effect is approximately 7 ns. These errors are therefore
generally smaller than the intrinsic RN we measure for most pulsars. For PSR
J1713+0747, observed with all three telescopes, reducing the impact of
rapidly-varying DM requires the window over which we fit a single DM to be
shorter than typical (e.g., $\lesssim$1 day compared to six days) to avoid a
larger DM mis-estimation error.

####  5.1.2 DM Mis-estimation due to Additional Chromatic Effects

There are a host of unmitigated chromatic effects in our timing model, though
many or all are expected to be small except for specific lines of sight.
Scattering is expected to be the next dominant effect. It is highly frequency
dependent, broadening the width of the PBF of the interstellar medium by an
amount that is $\propto\nu^{-4.4}$ for a Kolmogorov medium. When this
scattering delay, $t_{\rm scatt}\ll W$, where $W$ is the pulse width, the
first-order effect is to delay the TOA by $t_{\rm scatt}$. However, when
$t_{\rm scatt}$ is comparable to $W$, there is distortion of the pulse shape
and the time delay becomes a non-linear function of the scattering delay
(Hemberger & Stinebring, 2008). We are likely only in this regime for one
NANOGrav pulsar, J1903+0327, with a spin period of 2.2 ms and scattering delay
of $\sim$150 $\mu$s at 1.5 GHz (Geiger & Lam, 2022).

As discussed in §2.5, DM itself is a frequency-dependent quantity (Cordes et
al., 2016). Again considering a simplified model where pulsar observations are
made at only two spot frequencies $\nu$ and $\nu^{\prime}$, we can write the
TOA perturbation as

$\delta
t_{\infty}=\frac{-r^{2}t_{C,\nu}-r^{2}\epsilon_{\nu}+t_{C,\nu^{\prime}}+\epsilon_{\nu^{\prime}}}{r^{2}-1},$
(17)

where $t_{C,\nu}$ and $t_{C,\nu^{\prime}}$ are non-dispersive chromatic delays
at frequencies $\nu$ and $\nu^{\prime}$, and $\epsilon_{\nu}$ and
$\epsilon_{\nu^{\prime}}$ are additive, frequency-dependent errors at
frequencies $\nu$ and $\nu^{\prime}$. The additional error is uncorrelated
between different frequencies for radiometer noise but is highly correlated
for intrinsic pulse jitter, at least over modest frequency separations (see
below). The RMS timing error is found by calculating the square root of the
variance. Even when there is no WN, i.e.,
$\epsilon_{\nu}=\epsilon_{\nu^{\prime}}=0$, there is a non-zero perturbation
that affects our timing. Frequency-dependent DM affects us at the 10s of ns
level (Lam et al., 2018a) for most pulsars at the frequencies we observe,
though for a few at the highest DMs, that number can be substantially higher.

####  5.1.3 Accounting for Jitter

As noted in Section 3, the ECORR parameter was added to the pulsar signal
covariance matrix to account for frequency-correlated WN. When examined over
NANOGrav’s 12.5-yr data set pulsars, we find that a global fiducial estimate
for the RMS uncertainty from jitter is (Lam et al., 2019)

$\sigma_{\rm J}\approx 60~{}\mathrm{ns}\left(\frac{P_{\rm
s}}{5~{}\mathrm{ms}}\right)^{3/2}\left(\frac{T_{\rm
obs}}{30~{}\mathrm{min}}\right)^{-1/2}$ (18)

for a pulsar with a $P_{\rm s}=5$-ms spin period and a $T_{\rm obs}=30$ minute
observing time, if we take the 1500 MHz estimate as a representative value.
This uncertainty level makes up a non-negligible fraction of our known WN
budget. In its general use, for narrowband timing, the ECORR estimation is
higher than the physically-derived constraints in Lam et al. (2019) (see also
the comparison between ECORR and jitter in Lam et al., 2016), which prevents
potential underestimation of the frequency-correlated noise . Likewise, the
in-depth study of jitter in PTA pulsars revealed that jitter typically
decorrelates given a sufficiently broad frequency range (Shannon et al.,
2014). This decorrelation bandwidth is currently not incorporated into the
ECORR parameter because current receivers are narrow enough to avoid
decorrelation; however, with the implementation of wideband receivers in the
near future (see Section 6), this decorrelation bandwidth may fall within the
band of these new receivers.

####  5.1.4 The Impact of Polarization Mis-calibration

For any given narrowband profile, polarization gain mis-calibration will cause
a TOA uncertainty dependent on the fractional gain error $\varepsilon$ and the
degree of circular polarization $\pi_{\rm V}$ for a pulsar with pulse width
$W$ as (Cordes et al., 2004)

$\sigma_{\rm pol}\approx
1~{}\mu\mathrm{s}~{}\left(\frac{\varepsilon}{0.1}\right)\left(\frac{\pi_{\rm
V}}{0.1}\right)\left(\frac{W}{100~{}\mu\mathrm{s}}\right).$ (19)

The fiducial values given are representative of the NANOGrav pulsars and data
(Lam et al., 2018a) though the $\pi_{\rm V}$ and $W$ values can still vary
quite widely (e.g., van Straten, 2006; Gentile et al., 2018). Cross coupling
induces a false circular polarization that adds to the measured circular
polarization. While a seemingly large value, the uncertainty is on the
narrowband TOAs and not the infinite-frequency arrival time (e.g., the epoch-
averaged TOAs), and so the value should be reduced by the square root of the
number of TOAs, nearly an order of magnitude reduction. Polarization mis-
calibration uncertainties in the NANOGrav data set have not been fully
explored but are expected to add a of order 100 ns uncertainty to our excess
noise budget.

###  5.2 Red Noise Analysis

Since the strongest GW signal expected in PTAs manifests as a long-time
correlated stochastic signal, understanding RN in pulsar timing data sets is
of the utmost importance for mitigating false positives and/or parameter mis-
estimation. The use of power-law RN models in pulsar timing analyses is
ubiquitous and predates modern PTA data analysis techniques (Blandford et al.,
1984; Cordes & Downs, 1985). There are a number of more complex models for RN
in PTA data sets (Hazboun et al., 2020b; Lam et al., 2018b; Chalumeau et al.,
2022; Goncharov et al., 2021b), but we will demonstrate shortly that the
power-law model is a reasonable and effective choice by carrying out an
additional noise analysis.

####  5.2.1 Accounting for Red Noise

Searches for GWs in PTA data sets start with single-pulsar noise analyses and
then move on to full analyses that include all of the pulsar data sets. In the
case of the GWB this involves a correlated search across all pulsar data sets.
The main reason for the individual noise analyses is to search for the WN
parameters described in section Section 3. They involve full $N_{\rm
TOA}\times N_{\rm TOA}$ ($\mathcal{O}(10^{4})$ narrowband TOAs per pulsar, and
a factor of $\sim$30 reduction on average for the wideband data) matrix
inversions and upwards of 21 parameters per pulsar. Hence they are too costly
to do across all pulsars at once in a full Bayesian search. It was shown in
early, smaller data sets (Arzoumanian et al., 2015; Demorest et al., 2013)
that the WN parameters do not change substantially when the analysis is done
across other pulsars, as long as there is a RN model included in both
analyses, since this model can substantially change the amount of WN in some
pulsars. The presence of two RN analyses, one for individual pulsars and one
in the full GWB search, allows us to do an accounting of the RN in each pulsar
and keep track for which pulsars the RN power moves from the individual noise
channel into the common signal.

In Table 2 we show the pulsars that have significant RN detections in their
individual noise analyses. Our criterion is the same as for NG15; we consider
the detection to be significant when the pulsar’s RN has a Bayes factor
$\mathcal{B}>100$, using the Savage-Dickey (Dickey, 1971) approximation, when
possible. The italicized entries in Table 2 denote the pulsars that no longer
have significant detections of RN in the full PTA analysis, i.e., the power in
the single pulsar RN model has moved into the common channel555Depending on
the search this can either be a spatially-correlated process, as in searches
for the GWB, or just a “common” process, described by only the power spectral
density and not the spatial correlations. See NG15gwb for more details.. The
bolded entries denote the pulsars that continue to have significant RN
detections.

Table 2: Individual RN model parameter values and 68% credible intervals for pulsars with significant detections of RN. Pulsars with italicized entries only have significant detections in individual pulsar noise analyses while bolded entries show RN even when a common process is included. Note that $A_{\rm RN}$ is in strain amplitude to match the GW analyses, hence it is unit-less. Pulsar | $\log_{10}A_{\rm RN}$ | $\gamma_{\rm RN}$
---|---|---
B1855+09 | $\mathit{-14.0_{-0.4}^{+0.3}}$ | $\mathit{3.9_{-0.8}^{+1.0}}$
B1937+21 | $\mathbf{-13.6_{-0.1}^{+0.1}}$ | $\mathbf{4.0_{-0.3}^{+0.4}}$
B1953+29 | $\mathbf{-12.8_{-0.3}^{+0.2}}$ | $\mathbf{1.8_{-0.7}^{+1.1}}$
J0030+0451 | $\mathit{-14.4_{-0.5}^{+0.4}}$ | $\mathit{4.6_{-0.9}^{+1.1}}$
J0437$-$4715 | $\mathit{-13.4_{-0.2}^{+0.2}}$ | $\mathit{0.5_{-0.4}^{+0.6}}$
J0610$-$2100 | $\mathbf{-12.9_{-0.5}^{+0.3}}$ | $\mathbf{4.1_{-1.9}^{+2.0}}$
J0613$-$0200 | $\mathit{-13.8_{-0.3}^{+0.3}}$ | $\mathit{3.1_{-0.7}^{+0.9}}$
J1012+5307 | $\mathbf{-12.6_{-0.1}^{+0.1}}$ | $\mathbf{0.8_{-0.3}^{+0.3}}$
J1600$-$3053 | $\mathit{-13.5_{-0.6}^{+0.2}}$ | $\mathit{1.6_{-0.7}^{+1.5}}$
J1614$-$2230 | $\mathit{-14.9_{-0.8}^{+1.0}}$ | $\mathit{4.7_{-2.0}^{+1.6}}$
J1643$-$1224 | $\mathbf{-12.3_{-0.1}^{+0.1}}$ | $\mathbf{0.9_{-0.4}^{+0.4}}$
J1705$-$1903 | $\mathbf{-12.6_{-0.1}^{+0.1}}$ | $\mathbf{0.5_{-0.3}^{+0.4}}$
J1713+0747 | $\mathit{-14.1_{-0.1}^{+0.1}}$ | $\mathit{2.6_{-0.4}^{+0.5}}$
J1738+0333 | $\mathit{-14.6_{-0.6}^{+0.8}}$ | $\mathit{5.2_{-1.8}^{+1.3}}$
J1744$-$1134 | $\mathit{-14.1_{-0.6}^{+0.4}}$ | $\mathit{3.6_{-1.2}^{+1.4}}$
J1745+1017 | $\mathbf{-11.9_{-0.1}^{+0.1}}$ | $\mathbf{2.4_{-0.5}^{+0.6}}$
J1747$-$4036 | $\mathbf{-12.6_{-0.2}^{+0.1}}$ | $\mathbf{2.4_{-0.7}^{+1.0}}$
J1802$-$2124 | $\mathbf{-12.2_{-0.2}^{+0.2}}$ | $\mathbf{1.8_{-0.6}^{+0.7}}$
J1853+1303 | $\mathit{-13.5_{-0.4}^{+0.2}}$ | $\mathit{2.3_{-0.7}^{+1.1}}$
J1903+0327 | $\mathbf{-12.2_{-0.1}^{+0.1}}$ | $\mathbf{1.5_{-0.4}^{+0.4}}$
J1909$-$3744 | $\mathit{-14.5_{-0.4}^{+0.3}}$ | $\mathit{4.1_{-0.9}^{+1.0}}$
J1918$-$0642 | $\mathit{-13.8_{-0.7}^{+0.4}}$ | $\mathit{2.7_{-1.0}^{+1.5}}$
J1946+3417 | $\mathbf{-12.5_{-0.1}^{+0.1}}$ | $\mathbf{1.4_{-0.4}^{+0.5}}$
J2145$-$0750 | $\mathbf{-12.9_{-0.1}^{+0.1}}$ | $\mathbf{0.6_{-0.4}^{+0.5}}$
J2234+0611 | $\mathit{-13.9_{-0.9}^{+0.5}}$ | $\mathit{3.2_{-1.9}^{+2.7}}$

Figure 2 shows the RN power law posteriors for the same pulsars included in
Table 2. Those italicized in Table 2 are grayed in Figure 2, along with the
overlaid posterior for a common process across all of the pulsars. Note that
the gray contours follow a well known covariance trend between the spectral
index and RN amplitude. See NG15gwb for more details, along with a way to
ameliorate the covariance in parameter recoveries.

Figure 2: Two-dimensional parameter posterior contours for the power-law RN
models used in the pulsar noise models of the NG15 data set. The power
spectral density is parameterized by a spectral index and amplitude. The black
contours show the RN parameters of pulsars that have significant individual RN
in both their single pulsar analyses and the full PTA analysis. The gray
contours show RN parameters of pulsars that have significant individual RN
detections only in the single pulsar analyses, i.e., the power in the
individual-pulsar model has moved into the common model. The blue concentric
contours show the parameter recovery from a full PTA common process analysis
(NG15gwb). Posteriors for pulsars without significant individual RN detections
are not included for clarity.

The results are broadly as one would expect. The pulsars that possess
significant RN detections in their individual noise analysis, that lie along a
similar slope of covariance to the common process recovery, lose their RN
power to the common process. This is corroborated by a “leave one out”
analysis in NG15gwb, where individual pulsars are dropped out of a common
uncorrelated red process (curn) analysis and an odds ratio is calculated
between the analysis with $N$ pulsars versus $N-1$ pulsars. All of the grayed-
out pulsars have positive dropout factors, i.e. the analysis with that pulsar
is favored over the analysis without that pulsar, except for PSR J0437$-$4715
and PSR J1713+0747. These two pulsars are known to have challenging noise
properties that are discussed below.

It is important to note that while there is obviously another intriguing
cluster of posteriors that are a few orders of magnitude larger than the
recovered curn, the RN in those pulsars does not seem to be common and would
be excluded by the pulsars in the gray cluster that do not see this RN.
Various Bayesian analyses find no curn in that part of parameter space
(NG15gwb). Given the shallow spectral index this red noise is likely due to
mis-modeled chromatic noise (Cordes & Shannon, 2010). While only the
individual noise analysis results are shown in Table 2 there is no change in
the median RN values at the precision reported, and only slight changes in the
credible intervals.

The individual RN models are working as expected. For pulsars where the RN
seems to be a part of the curn the RN models are a stand-in for the common
process in the individual noise analyses, but do not hold onto the RN during
the full PTA analysis666The choice of priors is very important to ensure that
this is the case (Hazboun et al., 2020a). where the common process instead
models the RN. For pulsars where the RN seems to be truly intrinsic to only a
single pulsar data set the RN models are successfully modeling that noise,
which is _not_ able to contaminate the common signal. One interesting
characteristic of the grayed-out posteriors is that they cover a large span in
the spectral index, but all follow a roughly common slope. Pulsars with short
time spans, e.g., PSR J0437$-$4715 with only $4.8$ yrs of data and the gray
posterior against the left side of Figure 2, may have shallow spectral index
recoveries because they lack enough sensitivity at low frequencies where the
GWB will rise above their WN floor.

Next we present three separate noise analyses for four separate representative
pulsars777Figures and analyses are included for all pulsars at the end of this
manuscript. in Figure 3 through Figure 6. Each figure contains all of the
components of the power-law model, a free spectral analysis, and what we refer
to as an “excess noise analysis,” discussed below. The plots show RMS
fluctuations of the timing residuals, i.e., the amplitude of the noise in the
residuals across the frequencies shown, in units of $\log_{10}(\mathrm{s})$.
We convert the parameter posterior probabilities from all noise models into
amplitude for ease of comparison against the precision of the measurements.
The most generic noise model uses a separate set of $\varphi$ coefficients
(Equation 12), allowed to vary freely, without reference to any particular
power spectral density model. This so-called free spectral model (Lentati et
al., 2016; Hazboun et al., 2020b) is a Bayesian spectrogram of the pulsar
data. Here the spectra have been recovered from each pulsar’s data across 30
frequencies, ranging from $1/T_{\rm span}$ to $30/T_{\rm span}$, where $T_{\rm
span}$ is the full time span of the PTA888The choice of these 30 frequencies
is based on the usual Nyquist frequency sampling considerations and extends up
to frequencies higher than where a GWB would be detectable due to the WN floor
(Lentati et al., 2013). Unfortunately, the use of these models, with 30
frequency _parameters_ for each pulsar, is currently not feasible for full PTA
analyses due to computational limitations.. The solid orange “violin” plots
show the posteriors, frequency by frequency, for the free spectral analysis.
Significant posteriors are those where the violin plot is separated completely
from the frequency axis, as in the second to lowest frequency free spectral
(solid orange) posterior of Figure 3. Insignificant posteriors have broad
posteriors all the way down to the minimum values, while slightly significant
detections have thin posteriors extending to the minimum.

The free spectral model posteriors represent a “raw” spectrum, which does not
assume anything about the power in a pulsar data set. However, when we search
for common uncorrelated processes across the pulsars a simpler power law model
is used for both the RN and the curn. The recovered curn from NG15gwb, plotted
the same for all pulsars, is shown as the dashed red line. The horizontal gray
dashed line shows the WN power spectral density of the residuals (subscript
R), converted into these units as ${\rm WN}=\sqrt{P^{\rm(WN)}_{R}/T_{\rm
span}}$. The individual RN is shown with a dotted green line. This line may
not be visible for pulsars where the power law amplitude is small or not very
significant. Since an individual RN model is used for every pulsar, regardless
of significance, we have included the maximum likelihood values of the RN for
all pulsars. We do not include the RN from individual pulsar noise analyses in
these figures for clarity. See Table 2 below for a list of pulsars significant
detections of RN. The solid black line shows the total of the common process
RN, the WN, and any individual RN model and represents the total noise power
(in these units) for the pulsar using these models.

Sources of RN were studied in Cordes & Shannon (2010). Large RN with a
shallower spectral index is thought to be due to modeling errors of time-
correlated chromatic effects, while steeper spectral indices are thought to
originate from achromatic processes either intrinsic to the pulsar system or
due to a stochastic GWB. The pulsars can broadly be separated into four
categories depending on what type of noise they are dominated by in their
lowest frequencies.

> Common Process Dominated: Figure 3 shows the noise budget for PSR
> J1909$-$3744, one of the most sensitive pulsars in the NANOGrav array. This
> pulsar has very low WN power spectral density, has a significant detection
> of power-law RN in its individual noise analysis, and shows significant
> recovery of power at a few low frequencies of the free spectral noise
> analysis. As denoted by the italicized text in Table 2, most of the RN power
> moves into the curn in a full PTA analysis, hence it is dominated by RN that
> appears to be a part of the common process. See Section 5.2.1 for a full
> accounting of the RN across the PTA.
>
> Looking at the free spectral posteriors, one might ask whether a power-law
> model is sufficient for such a pulsar, but in fact the excess noise analysis
> shows no significant detections of additional noise. The standard model
> accounts for power in the 30 frequencies considered here. This is true
> across all of the pulsars in the array.  Figure 3: The excess timing
> residual delay as a function of frequency for PSR J1909$-$3744\. The total
> noise (solid black line) includes WN (gray dashed line), common RN (red
> dashed line), and individual RN (green dashed line). The free spectral model
> (solid orange posteriors) does not dictate any relationship between the
> amplitude of the power at different frequencies. Holding the total noise
> model (solid black line) parameters constant and searching for additional
> noise using a free spectral model results in the excess noise shown in blue.
> The vertical dashed line denotes a frequency of 1 yr-1. Bayes factors for
> these parameters, shown in the bottom panel, are fairly insignificant across
> all frequencies.
>
> White Noise Dominated: Figure 4 shows a representative of a second category
> of pulsar – one dominated by WN at the lowest frequencies. PSR J0509+0856
> only has $3.6$ years of data, just making the cut for a pulsar included in
> NANOGrav gravitational analyses. One does not expect sensitivity to RN from
> this pulsar at the lowest frequencies considered here, since they are based
> on the longest time span covered by the full PTA data set. One can see that
> there are no significant detections of power at individual frequencies in
> either the free spectral model or the excess noise model.  Figure 4: The
> excess timing residual delay as a function of frequency for PSR J0509+0856.
> See Figure 3 for details. Bayes factors for these parameters, shown in the
> bottom panel, are insignificant across all frequencies. WN dominates this
> pulsar’s short-timespan data set.
>
> Shallow Red Noise Dominated: Figure 5 shows a third type of pulsar, PSR
> J1903+0327, also dominated by RN at the lowest frequencies, but which does
> not move into the common channel during a full PTA analysis. PSR J1903+0327
> is our highest DM pulsar and so significant unmodeled chromatic propagation
> effects are expected to impact the timing (see e.g., Geiger & Lam, 2022).
> Here again, the power-law RN model is effective at modeling most of the
> noise picked up by the free spectral model, as seen by the insignificant
> Bayes factors.  Figure 5: The excess timing residual delay as a function of
> frequency for PSR J1903+0327. See Figure 3 for details. Bayes factors for
> these parameters, shown in the bottom panel, are fairly insignificant across
> all frequencies. The lowest frequencies are dominated by RN intrinsic to the
> pulsar.
>
> Steep Red Noise Dominated: Figure 6 shows a fourth type of pulsar, also
> dominated by RN at the lowest frequencies, that also does not move into the
> common channel during a full PTA analysis. PSR B1937+21 is a well-known
> example of high amplitude, steep-spectral index RN which appears to be
> intrinsic to the pulsar system. Here again, the power-law RN model is
> effective at modeling most of the noise picked up by the free spectral
> model, as seen by the insignificant Bayes factors.  Figure 6: The excess
> timing residual delay as a function of frequency for PSR B1937+21. See
> Figure 3 for details. Bayes factors for these parameters, shown in the
> bottom panel, are fairly insignificant across all frequencies. The lowest
> frequencies are dominated by RN intrinsic to the pulsar.

Noise budget figures for all of the pulsars used in Agazie et al. (2023b) are
included as Figures 14–76.

####  5.2.2 Comparison of Time-Correlated Noise Models

In order to understand the quantities of noise that might be unaccounted for
by the power-law model in our data set, we run an additional noise analysis on
each pulsar. This analysis fixes the WN, power-law RN, and GWB parameters,
i.e., the solid black line in Figure 3 through Figure 4, and adds in an
additional analysis, using a separate set of $\varphi$ coefficients, allowed
to vary freely. This additional model is effectively the same as the free
spectral analysis mentioned above, but this model now accounts for any excess
noise unaccounted for in the usual noise model. As we will see, it is mostly
superfluous when the power-law model is used. The WN parameters are from the
individual pulsar noise analyses discussed in Section 3.1. In addition we use
the 2-dimensional posterior power-law RN maximum likelihood values from a full
PTA search for the curn along with the median value for the curn amplitude
from that same analysis. The objective is to understand the noise unmodeled by
the WN, RN and curn stochastic processes. As with the free spectral model the
analyses are done over 30 frequencies ranging from $1/T_{\rm span}$ to
$30/T_{\rm span}$. The new free spectral coefficients, referred to as “excess
noise,” quantify how much power at various frequencies in each pulsar’s data
set is not modeled by the WN + RN + curn model.

The results from this excess noise analysis are compared to the results of the
other noise analyses for the four representative pulsars in Figure 3 to Figure
6. The excess noise posteriors are shown as the hollow blue violins. The
bottom panel shows the Bayes factor, calculated using a Savage-Dickey
approximation (Dickey, 1971), for the excess noise parameters.

Figure 7 summarizes the excess noise results from all of the 67 pulsars in the
NANOGrav PTA. We show the individual Bayes factors for each pulsar across the
30 frequencies considered in our GWB search. Again, the solid orange dots show
the Bayes factor for the free spectral parameters while the blue circles show
the Bayes factors for the excess noise parameters. The orange dots with arrows
at the top represent posteriors where the detection is too significant to use
the Savage-Dickey approximation. The number of such detections is noted with a
digit in the circle. The main takeaway from Figure 7 is that there is no case
where a significant free spectral detection corresponds to a significant
excess noise detection, revealing that the power-law model used across the PTA
is a sufficient model, given the sensitivity of the data, for mitigating RN in
individual pulsars and detecting RN from a common process.

Figure 7: Free spectral (orange) and excess noise (blue) model Bayes factors
for all 67 pulsars used in the 15-yr GWB search. The orange dots with arrows
represent posteriors where the detection is too significant to use the Savage-
Dickey approximation. The blue and orange circles represent the same
frequencies, but the blue circles are offset slightly to the left of the
orange for ease of viewing. The vertical dashed black line demarcates the 14th
frequency. Only the lowest 14 frequencies were used in the GWB search.

###  5.3 Sensitivity Curves

Detection sensitivity curves are commonly used in the GW community to
summarize all aspects of detector characterization into a single “figure of
merit”. They are often used by the broader astrophysics community to assess
the detectability of various GW sources. Here we use the Python package
hasasia (Hazboun et al., 2019a), based on the formalism developed in Hazboun
et al. (2019b), to calculate sensitivity curves for the 67 pulsars in the
NANOGrav PTA and then combine these into a sensitivity curve for the
stochastic GWB. This combination is based on the signal-to-noise ratio of the
GWB optimal statistic developed in Anholm et al. (2009); Chamberlin et al.
(2015); Rosado et al. (2015).

The effects of the timing model on searches for various astrophysically
interesting signals, especially noise correlated on long timescales, has been
known for a long time (Blandford et al., 1984). These effects can be
represented by a transmission function, $\mathcal{T}(f)$, that encodes the
fraction of power transmitted through the timing model. In this way the
transmission function acts as the transfer function for a pulsar (Blandford et
al., 1984; Cordes & Downs, 1985). A slightly different, but equivalent, method
(van Haasteren & Levin, 2013) than described in Section 4 is used to
marginalize over the timing model parameters and build the transmission
function. Here the $G$-matrix, derived from the timing model design matrix,
$M$, encodes the information about the timing model fit and acts to project
the data and covariance matrix into a basis orthogonal to the timing model
(see Hazboun et al. (2019b) for details). The transmission function can be
written in terms of $G$,

$\mathcal{T}(f)\equiv\frac{1}{N_{\rm
TOA}}\sum_{i,j}\left(GG^{\mathcal{T}}\right)_{ij}e^{i2\pi f(t_{i}-t_{j})}$
(20)

and calculated for any pulsar in the array. In Figure 11$a$ the transmission
function for PSR J1909$-$3744 is shown. The transmission function has a few
interesting features. First, the fit for the rotation frequency and its
derivative (the spin and spindown) of the pulsar pulls a quadratic polynomial
out of the data, which acts as a high pass filter, and limits sensitivity at
the lowest frequencies. As the time span of the pulsar data set increases, the
spindown parameters are fit more precisely and the frequency at which the
transmission “turns down” moves to lower and lower frequencies. Second, the
single frequency fits for the sky position/proper motion of the pulsar
($f=1~{}{\rm yr}^{-1}$) and parallax ($f=2~{}{\rm yr}^{-1}$) remove power in a
narrow band around those frequencies. If the pulsar has a binary period within
the frequency range, there will be another dip in the transmission function
for that fit as well. The width of these dips is proportional to $1/T_{\rm
span}$, and therefore narrows the longer a pulsar is timed. Lastly, the DM
variation model used removes power across the entire GW frequency band. The
DMX model constitutes $>200$ parameters for a few of the pulsars timed and
diminishes the power by a factor of $\sim 3$ across the frequencies searched
for GWs. Therefore the transmission function does not asymptote to
$\mathcal{T}(f)\sim 1$ at the highest frequencies. (Note that any DM variation
model with the same frequency resolution would remove a similar amount of
power.)

Figure 8: (a) Figure 9: (b) Figure 10: (c)

Figure 11: Various spectra for PSR J1909$-$3744\. The transmission function is
shown in (a) as a dashed orange line, while perfect transmission
($\mathcal{T}(f)=1$) is shown as the dotted black line. The various frequency-
dependent features are discussed in the text. Various power spectral densities
of the residuals, $P_{R}(f)$, are shown in (b). The black dotted line shows
the WN power spectral density $P^{\rm(WN)}_{R}$. The orange dashed line shows
$P^{\rm(WN)}_{R}$ with the effects of the timing model (TM) included. The gray
dashed line shows the GWB reproduced using the median values of the curnγ
model recovered in NG15gwb. The blue solid line shows the full power spectral
density of the pulsar and includes the effects of the timing model, WN,
individual pulsar RN, if significant, and the GWB. In plot (c) $P_{R}(f)$ has
been converted into units of characteristic strain, $h_{c}(f)$, for easy
comparison with other GW detectors.

The next step in calculating a pulsar sensitivity curve is to calculate
$\mathcal{N}^{-1}$, the noise-weighted transmission function.

$\mathcal{N}^{-1}(f)\equiv\frac{1}{2T_{\rm
span}}\sum_{i,j}\left[G(G^{T}CG)^{-1}G^{T}\right]_{ij}e^{i2\pi
f(t_{i}-t_{j})},$ (21)

which can also be thought of as the timing model-marginalized power spectral
density, i.e., the Fourier domain equivalent of the timing model-marginalized
covariance matrix. In fact, $\mathcal{N}=P_{R}(f)$ is the power spectral
density of the pulsar’s residuals, and hence the power output of a single arm
of our Galactic-scale GW detector. The power spectral density can be
transformed to units of GW strain by taking into account the response function
of TOAs to GWs.

Figure 11$b$ shows three separate curves in residual power. The dotted line
shows the WN power spectral density, $P^{\rm(WN)}_{R}$, where
$P^{\rm(WN)}_{R}=\frac{1}{2T_{\rm
span}}\sum_{i,j}C^{\rm(WN)}_{ij}e^{i2\pi(t_{i}-t_{j})}$. This only includes
the WN described in Section 3.1 and shows what the WN power would look like
without the effects of the transmission function. The orange dashed line shows
the effects of the transmission function on the WN power. It is calculated by
replacing $C$ in Equation 21 with $C^{\rm(WN)}$. The solid blue line includes
all of the modeled noise, modeling the same power as shown in the black solid
line of Figure 3.

The power spectral density in terms of pulsar residuals, $P_{R}$, can be
converted into strain power spectral density,

$S(f)=\frac{P_{R}(f)}{\mathcal{R}}=12\pi^{2}f^{2}P_{R}(f),$ (22)

by taking into account the pulsar response function to GWs. The sky-averaged
response function useful for GWB characterization is simply

$\mathcal{R}(f)=\frac{1}{12\pi^{2}f^{2}},$ (23)

where a factor of three comes from the sky averaging and the $(2\pi f)^{2}$
stems from the fact that we are working with timing residuals instead of
Doppler shifts in the pulse frequency, where $h_{\rm GWB}=\Delta f_{\rm
s}/f_{\rm s}$. See Hazboun et al. (2019b); Sesana et al. (2004) for more
details about the pulsar residual response function.

The strain power spectral density can then, in turn, be converted into units
of characteristic strain, the units in which sensitivity curves are often
plotted, via $h_{c}(f)\equiv\sqrt{fS(f)}=\pi
f^{\frac{3}{2}}\sqrt{12P_{R}(f)}$. Figure 11$c$ shows the strain power
spectral density for J1909$-$3744 converted into units of characteristic
strain. In these units there is a distinctive positive slope in the curve at
higher frequencies that goes like $\sim f^{\frac{3}{2}}$. It is also evident
that the GWB acts as a noise floor for an individual pulsar at some
frequencies. It is important to note that the correlated Earth-term power in
the GWB as measured at the Earth (“Earth-term”) is the signal that we are
searching for, but the power at the pulsar (“pulsar-term”) is uncorrelated
across the pulsars (Cornish & Sesana, 2013), and hence needs to be included in
the noise budget and calculations of the sensitivity.

The individual pulsar strain power spectral densities can be combined using
the optimal statistic signal-to-noise ratio in the frequency domain to
construct a detection sensitivity curve for the NANOGrav PTA999These curves
are constructed using a statistic specific to a stochastic GWB. Curves for
single, resolvable sources can also be constructed using a matched filter
statistic (Hazboun et al., 2019b). (Hazboun et al., 2019b), shown in Figure
12. The blue curve is the fiducial curve to compare the stochastic
gravitational signals against since it contains the self noise due to the
background itself. The orange curve, which only includes WN, shows the
importance of including the GWB and RN in calculating the sensitivity.

Figure 12: Stochastic GWB sensitivity curves. The solid blue curve shows the
characteristic strain sensitivity as a function of GW frequency for the
NANOGrav 15-year data set and includes contributions from the WN, RN, and GWB
self noise. The gray line shows the recovered GWB (NG15gwb), while the dashed
orange curve is the sensitivity curve only including the contributions from
WN. The importance of including time-correlated noise can be seen in the
overly optimistic orange dashed curve. The vertical lines show frequencies at
one and two times the inverse of the data set time span.

The sensitivity curve in Figure 12 is calculated for any possible spectral
shape of the stochastic GWB. In order to understand the detectability of that
background one needs to integrate the GWB power spectral density against the
full “effective” power spectral density in the frequency domain. Following
Thrane & Romano (2013) we determine the power law-integrated sensitivity curve
by calculating the amplitude of a power-law GWB we would detect with an
$\mathrm{S/N}=5$, as reported by the optimal statistic in NG15gwb. In Figure
13 we show the detection sensitivity curve in blue and the power law-
integrated curve in black101010In order to make the full curve aesthetically
complete a wide range of spectral indices are used.. Two particular spectral
index values are highlighted: $\alpha=-2/3$ ($\gamma=13/3$), the value
theoretically expected for a population of binaries where the loss of energy
to GWs dominates their evolution and $\alpha=-0.18$ ($\gamma=3.35$), the
median value for $\alpha$ obtained from the GWB search in NG15gwb.

Figure 13: Power-law integrated sensitivity curve. The full stochastic
sensitivity curve for PTAs combines the individual pulsar sensitivities using
the denominator of the optimal statistic Chamberlin et al. (2015); Hazboun et
al. (2019b). As in Thrane & Romano (2013) one can integrate the sensitivity
curve against a set of power law stochastic backgrounds, here shown as the set
of gray lines with $\alpha=-\frac{7}{4}$ to $\alpha=\frac{5}{4}$, and solve
for the amplitude needed, $\hat{A}^{\alpha}_{\mathrm{GWB}}$, to meet a given
threshold in S/N. The power-law integrated sensitivity curve, solid black, is
the maximum of the set of those lines at each frequency. The blue line shows
the GWB sensitivity, as described in Figure 12, and the orange and green lines
show the power-law GWBs detectable with $\mathrm{S/N}=5$ for two discrete
values of $\alpha$, described in the text.

The power law-integrated calculation gives us another important accounting of
the noise budget for the NANOGrav detector. This analysis demonstrates that by
folding together our knowledge of the noise characterization of the individual
pulsars and the parameters for the common process, we can estimate the
recovered amplitude, $\hat{A}^{\alpha}_{\mathrm{GWB}}$, for the two values of
the GWB spectral index used in NG15gwb.

##  6 Improvements and Future Outlook

###  6.1 Near-Future Wideband Noise Mitigation Strategies

Over the last several decades, the bandwidths of the radio telescope receivers
used for pulsar observations have steadily increased, with many receivers now
capable of observing over more than one or even several octaves of bandwidth.
Some examples of receivers used for pulsar timing array science are the
Parkes-Murriyang Ultra-Wideband Low Receiver (UWL; 0.7-‐4.2 GHz, Hobbs et al.
2020b), the Effelsberg Ultra-Broadband Receiver (UBB; 1.3–6
GHz111111https://www.mpifr-bonn.mpg.de/research/fundamental/ubb), and, for
NANOGrav, a 0.7‐-4.0 GHz Ultra Wideband Receiver (UWBR; Bulatek & White 2020)
being commissioned at the GBT, as well as the octave capability of CHIME
(400–800 MHz, CHIME/Pulsar Collaboration et al. 2021). Processing the data
streams remains a “big data” problem in which raw voltage data is not
desirable.

For pulsar timing applications, data is first coherently dedispersed to remove
DM delays in real time. Traditionally, TOAs were then calculated over a set of
discrete frequency channels. Comparing TOAs over wide frequency ranges then
yielded deviations of DM from the fiducial value used in the dedispersion
kernel. However, this traditional procedure ignored the pulse profile
evolution apparent across wide frequency ranges which is covariant with the
DM. This has long been accounted for in NANOGrav timing models through fitting
frequency-dependent (FD) pulse profile parameters to the channelized TOAs. The
move toward even wider and continuously sampled bandwidths, however, favors a
newer method of timing TOAs, by representing the profile as a 2D function that
takes into account a fiducial DM, pulse profile frequency evolution, and
scattering in a procedure known as “pulse portraiture” (Pennucci et al. 2014,
Liu et al. 2014) which reduces covariances. Starting with the 12.5-year data
release, NANOGrav has been creating both narrowband and wideband data sets
(Alam et al. 2021a, Alam et al. 2021b). With each new data release, we expect
the wideband data set to become more and more important, as data from the
newer wideband receivers will become a larger fraction of each existing data
set. Additionally, the portraiture framework will be critical for probing
frequency-dependent DM behavior, expected for some pulsars across very wide
bandwidths (Cordes et al., 2016).

Not all near-future noise mitigation strategies depend on new hardware. We
highlight a signal processing technique known as cyclic spectroscopy (CS).
Demorest (2011) and Walker et al. (2013) used CS as a deconvolution method to
analyze data from PSR B1937+21, both resolving an intrinsic profile as well as
revealing scintillation features with better resolution than that derived from
filterbank techniques. Additional work since has shown how we can improve
pulsar scintillation and TOA measurements (Archibald et al., 2014; Palliyaguru
et al., 2015; Dolch et al., 2021; Turner et al., 2023), where CS is most
effective for pulsars with high profile S/N and/or high scattering timescales,
as is the case for PSR B1937+21. While PSR B1937+21 is more highly scattered
than most current PTA pulsars, future telescopes and upcoming receivers on
existing telescopes will likely discover new, distant pulsars for which CS
deconvolution is applicable. In addition, some currently known pulsars may
become PTA-quality.

We expect significant S/N improvement both from the highly resolved pulse
scintillation structure enabled by CS and from improved RFI mitigation due to
frequency binning. This will be especially important for the heavy RFI-
contamination expected in the UWBR context. The computational limit, however,
is that CS requires either saved baseband data or cyclic spectra written in
real time, much like pulse profiles are currently written in real time after
the application of coherent dedispersion. For the second (preferable) case, a
real-time CS-backend on the GBT’s UWBR is currently under development.

###  6.2 Noise Mitigation for Different GW Sources

PTA collaborations continually strive to understand and mitigate the noise in
pulsar data sets. A deep understanding of the noise allows PTA collaborations
to forecast the sensitivity of our detector to various GW signals and even
tune it towards specific sources in some cases (Liu et al., 2023; Speri et
al., 2023). Figure 12 clearly demonstrates the strong effect of the stochastic
GWB on the sensitivity of our detector. As PTAs enter the post-GWB detection
era, the GWB will in fact be a foreground source of noise, relative to the
many individual GW sources searched for by PTAs. While continued resolution of
the background (Pol et al., 2021) will allow us to mitigate up to half of the
power in the GWB, half of the power in any given pulsar is uncorrelated across
the pulsars, i.e., the pulsar term. Flagship GW analyses, e.g. for resolvable
supermassive binary black holes (Arzoumanian et al., 2023), already
necessitate joint searches for a GWB and the deterministic signal to mitigate
false positives. Searches for steeper spectral index stochastic backgrounds
will require continued long time span observations (Kaiser et al., 2022),
while searches for resolvable binary signals will be served by better
sensitivity at higher frequencies (Liu et al., 2023).

In recent years a number of new strategies for noise modeling of pulsar timing
data sets have been developed (Lam et al., 2018b; Chalumeau et al., 2022;
Goncharov et al., 2021b) to mitigate against higher-order chromatic
propagation effects like scattering and telescope issues (i.e., “band noise”
and “system noise”). The latter noise types are characterized using power law
RN models exclusively in a particular frequency band or for a specific pulsar
backend. Largely these advances come in the form of Gaussian process models,
built into the Bayesian noise analyses and also used in the full PTA GW
searches. In Falxa et al. (2023) it was shown that including a large number of
frequencies in the DM modeling was necessary to mitigate against spurious
high-frequency single-source detections. These efforts reveal the need for and
potential of these noise mitigation strategies for increasing sensitivity to
high-frequency, deterministic GW sources.

###  6.3 The Future of PTA Noise Mitigation

New noise mitigation strategies for PTAs offer a two-fold benefit for the low-
frequency search for GWs. Mitigation techniques that involve new data will
allow for better and more sensitive searches of GWs from sources in the higher
end of our frequency band. Strategies that allow for mitigation of noise
across the $15+$ year data set allow us to extract more sensitivity out of the
data that we already possess. These techniques, as well as the in-depth
characterization of PTAs as Galactic-scale detectors, will allow PTAs to reach
their full potential as sensitive nanohertz GW observatories.

##  7 Conclusion

The North American Nanohertz Observatory for Gravitational Waves is a
Galactic-scale GW detector made up of some of the most sensitive radio
telescopes in the world and nature’s most stable clocks. The astrophysical
makeup of our detector requires careful noise analysis, starting from the
underlying physical processes adding uncertainty to our measurements. It also
requires that we use a well-motivated phenomenological model in order to
incorporate all sources of noise, whether explicitly accounted for or not, in
order to be as conservative as possible in our characterization.

We have carried out a full characterization of the noise in our nanohertz GW
detector, including a full accounting of RN in the pulsars and new analyses to
corroborate the use of power-law RN models. These noise analyses allow us to
construct detection sensitivity curves for the stochastic GWB which are in
agreement with the results of NG15gwb. We find that our phenomenological model
faithfully encapsulates the noise in our detector.

The NANOGrav collaboration receives support from National Science Foundation
(NSF) Physics Frontiers Center award numbers 1430284 and 2020265, the Gordon
and Betty Moore Foundation, 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. NANOGrav is part of the International Pulsar Timing Array
(IPTA); we would like to thank our IPTA colleagues for their help with this
paper. Particular thanks to Joris Verbiest for extensive conversations and
comments concerning early drafts of this manuscript. This work was 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. 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. 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.C. 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. A.D.J. and M.V. acknowledge
support from the Caltech and Jet Propulsion Laboratory President’s and
Director’s Research and Development Fund. A.D.J. acknowledges support from the
Sloan Foundation. 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.

## 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.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-yr 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. J.S.H.
coordinated the writing of this paper. S.C., T.D., L.G., J.S.H., M.T.L.,
N.Le., J.W.M., M.A.M., S.K.O, B.B.B.P., and D.R.S. wrote the paper and
collected the bibliography. J.S.H., N.La., J.S., and X.S carried out the noise
analyses. J.S.H carried out the sensitivity analysis.

## References

  * Agazie et al. (2023a) Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023a, ApJ, doi: 10.3847/2041-8213/acda9a
  * Agazie et al. (2023b) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023b, ApJ, doi: 10.3847/2041-8213/acdac6
  * Alam et al. (2021a) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021a, ApJS, 252, 4, doi: 10.3847/1538-4365/abc6a0
  * Alam et al. (2021b) —. 2021b, ApJS, 252, 5, doi: 10.3847/1538-4365/abc6a1
  * Anholm et al. (2009) Anholm, M., Ballmer, S., Creighton, J. D. E., Price, L. R., & Siemens, X. 2009, Phys. Rev. D, 79, 084030, doi: 10.1103/PhysRevD.79.084030
  * Antoniadis et al. (2022) Antoniadis, J., Arzoumanian, Z., Babak, S., et al. 2022, MNRAS, 510, 4873, doi: 10.1093/mnras/stab3418
  * Archibald et al. (2014) Archibald, A. M., Kondratiev, V. I., Hessels, J. W. T., & Stinebring, D. R. 2014, ApJ, 790, L22, doi: 10.1088/2041-8205/790/2/L22
  * Archibald et al. (2018) Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73, doi: 10.1038/s41586-018-0265-1
  * Archibald et al. (2013) Archibald, R. F., Kaspi, V. M., Ng, C.-Y., et al. 2013, Nature, 497, 591, doi: 10.1038/nature12159
  * Arzoumanian et al. (2015) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2015, ApJ, 813, 65, doi: 10.1088/0004-637X/813/1/65
  * Arzoumanian et al. (2016) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2016, ApJ, 821, 13, doi: 10.3847/0004-637X/821/1/13
  * Arzoumanian et al. (2018) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47, doi: 10.3847/1538-4357/aabd3b
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34, doi: 10.3847/2041-8213/abd401
  * Arzoumanian et al. (2023) Arzoumanian, Z., Baker, P. T., Blecha, L., et al. 2023, arXiv e-prints, arXiv:2301.03608, doi: 10.48550/arXiv.2301.03608
  * Backer (1970) Backer, D. C. 1970, Nature, 228, 42, doi: 10.1038/228042a0
  * Bassa et al. (2016) Bassa, C. G., Janssen, G. H., Stappers, B. W., et al. 2016, MNRAS, 460, 2207, doi: 10.1093/mnras/stw1134
  * Blandford et al. (1984) Blandford, R., Narayan, R., & Romani, R. W. 1984, Journal of Astrophysics and Astronomy, 5, 369, doi: 10.1007/BF02714466
  * Bochenek et al. (2015) Bochenek, C., Ransom, S., & Demorest, P. 2015, ApJ, 813, L4, doi: 10.1088/2041-8205/813/1/L4
  * Brook et al. (2018) Brook, P. R., Karastergiou, A., McLaughlin, M. A., et al. 2018, ApJ, 868, 122, doi: 10.3847/1538-4357/aae9e3
  * Bulatek & White (2020) Bulatek, A., & White, S. 2020, in American Astronomical Society Meeting Abstracts, Vol. 235, American Astronomical Society Meeting Abstracts #235, 175.17
  * Caballero et al. (2018) Caballero, R. N., Guo, Y. J., Lee, K. J., et al. 2018, MNRAS, 481, 5501, doi: 10.1093/mnras/sty2632
  * Chalumeau et al. (2022) Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538, doi: 10.1093/mnras/stab3283
  * 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
  * Champion et al. (2010) Champion, D. J., Hobbs, G. B., Manchester, R. N., et al. 2010, ApJ, 720, L201, doi: 10.1088/2041-8205/720/2/L201
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970, doi: 10.1093/mnras/stab2833
  * Cheng (1987) Cheng, K. S. 1987, ApJ, 321, 805, doi: 10.1086/165673