# The NANOGrav 15 yr Data Set: Chromatic Gaussian Process Noise Models for Six Pulsars

Bjorn Larsen, Chiara M.F Mingarelli, Jeffrey S Hazboun, Aurelien

Chalumeau, Deborah C Good, Joseph Simon, Gabriella Agazie, Akash

Anumarlapudi, Anne M Archibald, Zaven Arzoumanian, et al.

► To cite this version:

Bjorn Larsen, Chiara M.F Mingarelli, Jeffrey S Hazboun, Aurelien Chalumeau, Deborah C Good, et al.. The NANOGrav 15 yr Data Set: Chromatic Gaussian Process Noise Models for Six Pulsars. The Astrophysical Journal, 2024, 972 (1), pp.49. 10.3847/1538-4357/ad5291. hal-04603642

HAL Id: hal-04603642

https://hal.science/hal-04603642v1

Submitted on 13 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The NANOGrav 15 yr Data Set: Chromatic Gaussian Process Noise Models for Six Pulsars

Bjorn Larsen¹, Chiara M. F. Mingarelli¹, Jeffrey S. Hazboun², Aurélien Chalumeau³,⁴, Deborah C. Good⁵, Joseph Simon⁶,⁵⁰, Gabriella Agazie⁷, Akash Anumarlapudi⁷, Anne M. Archibald⁸, Zaven Arzoumanian⁹, Paul T. Baker¹⁰, Paul R. Brook¹¹, H. Thankful Cromartie¹²,⁵¹, Kathryn Crowter¹³, Megan E. DeCesar¹⁴,⁵², Paul B. Demorest¹⁵, Timothy Dolch¹⁶,¹⁷, Elizabeth C. Ferrara¹⁸,¹⁹,²⁰, William Fiore²¹,²², Emmanuel Fonseca²¹,²², Gabriel E. Freedman⁷, Nate Garver-Daniels²¹,²², Peter A. Gentile²¹,²², Joseph Glaser²¹,²², Ross J. Jennings²¹,²²,⁵³, Megan L. Jones⁷, David L. Kaplan⁷, Matthew Kerr²³, Michael T. Lam²⁴,²⁵,²⁶, Duncan R. Lorimer²¹,²², Jing Luo²⁷,⁵⁴, Ryan S. Lynch²⁸, Alexander McEwen⁷, Maura A. McLaughlin²¹,²², Natasha McMann²⁹, Bradley W. Meyers¹³,³⁰, Cherry Ng³¹, David J. Nice³², Timothy T. Pennucci³³, Benetge B. P. Perera³⁴, Nihan S. Pol²⁹, Henri A. Radovan³⁵, Scott M. Ransom³⁶, Paul S. Ray²³, Ann Schmiedekamp³⁷, Carl Schmiedekamp³⁷, Brent J. Shapiro-Albert²¹,²²,³⁸, Ingrid H. Stairs¹³, Kevin Stovall¹⁵, Abhimanyu Susobhanan⁷, Joseph K. Swiggum³²,⁵³, Haley M. Wahl²¹,²², David J. Champion³⁹, Ismaël Cognard⁴⁰,⁴¹, Lucas Guillemot⁴⁰,⁴¹, Huanchen Hu³⁹, Michael J. Keith⁴², Kuo Liu³⁹,⁴³, James W. McKee⁴⁴,⁴⁵, Aditya Parthasarathy³⁹,⁴⁶,⁴⁷, Delphine Perrodin⁴⁸, Andrea Possenti⁴⁸, Golam M. Shaifullah³,⁴,⁴⁸, and Gilles Theureau⁴⁰,⁴¹,⁴⁹

¹ Department of Physics, Yale University, New Haven, CT 06520, USA; bjorn.larsen@yale.edu

² Department of Physics, Oregon State University, Corvallis, OR 97331, USA

³ Dipartimento di Fisica “G. Occhialini,” Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

⁴ INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

⁵ Department of Physics and Astronomy, University of Montana, 32 Campus Drive, Missoula, MT 59812, USA

⁶ Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA

⁷ Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA

⁸ Newcastle University, Newcastle University, NE1 7RU, UK

⁹ X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA

¹⁰ Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA

¹¹ Institute for Gravitational Wave and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK

¹² National Academy of Sciences, Washington, DC 20001, USA

¹³ Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC, V6T 1Z1, Canada

¹⁴ George Mason University, 4400 University Drive, Fairfax, VA 22030, USA

¹⁵ National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA

¹⁶ 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

¹⁸ Department of Astronomy, University of Maryland, College Park, MD 20742, USA

¹⁹ Center for Research and Exploration in Space Science and Technology, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²⁰ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²¹ 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

²³ Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

²⁴ SETI Institute, 339 N Bernardo Avenue Suite 200, Mountain View, CA 94043, USA

²⁵ School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA

²⁶ Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA

²⁷ Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON, M5S 3H4, Canada

²⁸ Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA

²⁹ Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA

³⁰ International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia

³¹ Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada

³² Department of Physics, Lafayette College, Easton, PA 18042, USA

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

³⁴ Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA

³⁵ Department of Physics, University of Puerto Rico, Mayagüez, PR 00681, USA

³⁶ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

³⁷ Department of Physics, Penn State Abington, Abington, PA 19001, USA

³⁸ Giant Army, 915A 17th Avenue, Seattle, WA 98122, USA

³⁹ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany

⁴⁰ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, 45071 Orléans Cedex 02, France

⁴¹ Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

⁴² Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester, M13 9PL, UK

⁴³ Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

⁴⁴ 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

⁴⁶ ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensdijk 4, 7991 PD Dwingeloo, The Netherlands

⁴⁷ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

⁴⁸ INAF-Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius, Italy
<sup>49</sup> Laboratoire Univers et Théories LUTH, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

Received 2024 January 28; revised 2024 May 23; accepted 2024 May 29; published 2024 August 23

**Abstract**

Pulsar timing arrays (PTAs) are designed to detect low-frequency gravitational waves (GWs). GWs induce achromatic signals in PTA data, meaning that the timing delays do not depend on radio frequency. However, pulse arrival times are also affected by radio-frequency-dependent “chromatic” noise from sources such as dispersion measure (DM) and scattering delay variations. Furthermore, the characterization of GW signals may be influenced by the choice of chromatic noise model for each pulsar. To better understand this effect, we assess if and how different chromatic noise models affect the achromatic noise properties in each pulsar. The models we compare include existing DM models used by the North American Nanohertz Observatory for Gravitational waves (NANOGrav) and noise models used for the European PTA Data Release 2 (EPTA DR2). We perform this comparison using a subsample of six pulsars from the NANOGrav 15 yr data set, selecting the same six pulsars as from the EPTA DR2 six-pulsar data set. We find that the choice of chromatic noise model noticeably affects the achromatic noise properties of several pulsars. This is most dramatic for PSR J1713+0747, where the amplitude of its achromatic red noise lowers from $\log_{10}A_{RN} = -14.1_{-0.1}^{+0.1}$ to $-14.7_{-0.5}^{+0.3}$, and the spectral index broadens from $\gamma_{RN} = 2.6_{-0.4}^{+0.5}$ to $\gamma_{RN} = 3.5_{-0.9}^{+1.2}$. We also compare each pulsar’s noise properties with those inferred from the EPTA DR2, using the same models. From the discrepancies, we identify potential areas where the noise models could be improved. These results highlight the potential for custom chromatic noise models to improve PTA sensitivity to GWs.

*Unified Astronomy Thesaurus concepts: Millisecond pulsars (1062); Gravitational wave astronomy (675); Gravitational wave detectors (676); Pulsar timing method (1305); Interstellar medium (847); Astronomy data analysis (1858)*

**1. Introduction**

Pulsar timing arrays (PTAs) are designed to detect low-frequency gravitational waves (GWs). GWs induce small shifts in pulse times-of-arrival (TOAs), which can be measured using a network of millisecond pulsars (Sazhin 1978; Detweiler 1979; Hellings & Downs 1983). Several collaborations around the globe are carrying out searches for GW signals using PTAs, including the North American Nanohertz Observatory for Gravitational waves (NANOGrav; Ransom et al. 2019), the European PTA (EPTA; Desvignes et al. 2016), the Parkes PTA (PPTA; Manchester et al. 2013), the Indian PTA (InPTA; Joshi et al. 2018), the Chinese PTA (CPTA; Lee 2016), and the MeerKAT PTA (MPTA; Miles et al. 2023). Together, NANOGrav, EPTA, PPTA, and InPTA form the International PTA (IPTA; e.g., Verbiest et al. 2016; Perera et al. 2019).

Evidence for a stochastic GW background (GWB) at nanohertz frequencies has recently been presented by NANOGrav (Agazie et al. 2023a, hereafter NG15_GWB), EPTA + InPTA (EPTA Collaboration et al. 2023a), PPTA (Reardon et al. 2023a), and CPTA (Xu et al. 2023), with varying levels of significance but broadly consistent properties across data sets (Agazie et al. 2024a, hereafter IPTA_comp). This evidence is based on the presence of a time-correlated, low-frequency (“red”) noise process common to and spatially correlated between all pulsars across the sky. These spatial correlations

follow the Hellings and Downs (HD) curve, which is the definitive signature of an isotropic GWB (Hellings & Downs 1983).

Among the next major milestones in PTA science is to identify and characterize the source of the GWB, which may be astrophysical, cosmological, or a combination of both (Burke-Spolaor et al. 2019; Mingarelli & Casey-Clyde 2022). A likely source of the GWB is the incoherent superposition of GWs produced by hundreds of thousands of slowly inspiraling supermassive black hole binaries (SMBHBs; e.g., Rajagopal & Romani 1995; Kelley et al. 2017; Agazie et al. 2023b, hereafter NG15_SMBHBs). More exotic sources of the background have also been proposed, such as early-Universe phase transitions, cosmic strings, and relic GWs from inflation (e.g., Caprini & Figueroa 2018; Afzal et al. 2023; Vagnozzi 2023).

Each potential source may be distinguished using the shape of the GWB spectrum inferred from PTA timing residuals (Lasky et al. 2016; Kaiser et al. 2022). For instance, the simplest analytic model of the GWB from SMBHB populations predicts a power-law timing residual spectrum (Equation (6)) with a spectral index $\gamma_{GWB} = 13/3$ (Phinney 2001). However, deviations from this simplified spectrum may result, e.g., from more complicated models of SMBHB evolution (Sesana 2013; Kocsis & Sesana 2011; NG15_SMBHBs), discreteness of the SMBHB population (Sesana et al. 2008; NG15_SMBHBs), or resolvable single sources within the PTA data set (Bécsy et al. 2023).

In order to measure the GWB spectrum as accurately as possible, it is important to account for different sources of noise affecting individual pulsars. Specifically, either overfitting (Hazboun et al. 2020a) or underfitting (Hazboun et al. 2020b; Zic et al. 2022) for pulsar noise processes may bias inferences of a common *uncorrelated* red noise (CURN) process, which encodes the spectrum of the GWB without including interpulsar HD correlations (Arzoumanian et al. 2020; Romano et al. 2021; Taylor et al. 2022). A promising approach to accurately

<figure><img src="image_3.png" /></figure>

<footnote n="50">50 NSF Astronomy and Astrophysics Postdoctoral Fellow.</footnote>

<footnote n="51">51 National Research Council Postdoctoral Associate, USA resident at Naval Research Laboratory, Washington, DC 20375, USA.</footnote>

<footnote n="52">52 Resident at the Naval Research Laboratory, Washington, DC 20375, USA.</footnote>

<footnote n="53">53 NANOGrav Physics Frontiers Center Postdoctoral Fellow.</footnote>

<footnote n="54">54 Deceased.</footnote>
model pulsar noise, first employed by Lentati et al. (2016), is the creation of custom noise models for each pulsar using Bayesian model selection methods. Recently, Reardon et al. (2023b) found use of custom pulsar noise models to significantly influence the recovered spectral characteristics of the CURN in PPTA Data Release 3 (DR3; Zic et al. 2023). Conversely, Chalumeau et al. (2022) found that custom pulsar noise models have a minimal effect on both the spectral characterization of the CURN and the detection statistics for HD correlations using the European PTA Data Release 2 (EPTA DR2) six-pulsar data set (Chen et al. 2021). These differences suggest the importance of custom pulsar noise models for GWB analyses may vary depending on the properties of the data set.

We turn our attention now to the NANOGrav 15 yr data set (Agazie et al. 2023c, hereafter NG15). NG15_GWB and NG15_SMBHBs compared the inferred GWB spectral parameters when changing the dispersion measure (DM) model applied to all pulsars in NG15, where DM variations introduce chromatic (radio-frequency-dependent) timing noise. The choice of DM model in NG15 was found to affect the spectral parameter inferences of the GWB, with a slightly higher $A_{\text{GWB}}$ and lower $\gamma_{\text{GWB}}$ predicted using the standard DM model, “DMX,” than the alternative model, “DMGP.” While the effect is minor (the 68% credible regions of the 2D GWB posteriors overlap using both models; NG15_GWB), the choice of model could still be consequential for astrophysical inferences. It is therefore important to investigate the effect of these different models on a deeper level.

Here we investigate how the choice of chromatic noise model affects the achromatic red noise (ARN) in six pulsars from NG15, noting that a GWB signal contributes to a component of the ARN in each pulsar. We specifically compare the following three models: (i) the standard noise model used by NANOGrav (DMX), (ii) a noise model using Gaussian processes (GPs) for DM variations (DMGP), and (iii) a new noise model using GPs to account for additional chromatic noise alongside DM (CustomGP), similar to the models used in EPTA DR2 (Antoniadis et al. 2023, hereafter EPTA_Noise). We select PSRs J0613−0200, J1012+5307, J1600−3053, J1713+0737, J1744−1134, and J1909−3744 as our focus for this study since they are the same pulsars from the EPTA DR2 six-pulsar data set (Chalumeau et al. 2022). As a consistency check, we compare each pulsar’s noise properties as inferred under the CustomGP model using both NG15 and EPTA DR2. Since the astrophysical noise in each pulsar ought to be consistent regardless of the data set (IPTA_comp), we use the inconsistencies to identify potential improvements to these noise models.

The paper is laid out as follows. In Section 2, we describe the data used here. In Section 3, we present details on relevant noise processes affecting single pulsars. In Section 4, we describe the GP models we use. In Section 5, we present our results, starting with an overview of how (and if) each pulsar’s ARN changes as a function of the noise model, followed by a pulsar-by-pulsar noise breakdown including our comparisons with EPTA DR2. Finally, in Section 6, we discuss our results and provide recommendations for future analyses. Table 1 shows the acronyms and symbols used in this paper.

## 2. Data

### 2.1. The NANOGrav 15 yr Data Set

The NANOGrav 15 yr data set (NG15) contains observations of 68 millisecond pulsars with time spans ranging from 3 to 15 yr. NG15 comprises observations from three radio observatories: the Green Bank Telescope (GBT), the Arecibo Observatory (AO), and the Very Large Array (VLA). All six pulsars studied here are observed by the GBT. PSR J1713+0747 includes additional observations from the AO, while PSRs J1600−3053, J1713+0747, and J1909−3744 also include observations from the VLA. The observations were collected, reduced, and analyzed to produce a best-fit timing model, a set of narrowband and wideband TOAs, and a configuration file for each pulsar (NG15). Here we use the narrowband TOAs, which are derived from many subbands of the radio observing bands. NG15 uses the JPL DE440 solar system ephemeris (Park et al. 2021) and the TT(BIPM2019) timescale in order to correct observatory and terrestrial clocks to an inertial reference frame at the solar system barycenter.

### 2.2. EPTA DR2

EPTA DR2 (EPTA Collaboration et al. 2023b, hereafter EPTA DR2) was first made up of six millisecond pulsars (Chen et al. 2021) and later upgraded with more pulsars for a total of 25 millisecond pulsars and a maximum time span of 24 yr. EPTA DR2 comprises observations from six radio telescopes: the Effelsberg Radio Telescope (EFF), the Lovell Telescope (LT), the Mark II Telescope (MK2), the Nançay Radio Telescope (NRT), the Westerbork Synthesis Radio Telescope (WSRT), and the Sardinia Radio Telescope (SRT). All telescopes are also used in tied-array mode to make observations as the Large European Array for Pulsars (LEAP; Bassa et al. 2016). PSR J1909−3744 is only observed by the NRT and SRT due to its lower declination, but the remaining five pulsars are observed by all six telescopes. The EPTA DR2 observations undergo data reduction, combination, outlier and timing analyses to produce a final timing model and set of narrowband TOAs for each pulsar. EPTA DR2 uses the JPL DE440 solar system ephemeris (the same as NG15) and the TT(BIPM2021) timescale.

EPTA DR2 also presents multiple data-set versions, one of which excluded TOAs that had not undergone coherent dispersion and another which included TOAs from InPTA Data Release 1 (DR1; Tarafdar et al. 2022). While there are differences between the noise properties inferred using each data set (see EPTA_Noise), we found the differences were not substantial enough to change the conclusions of our comparison with NG15. As such, we compare our results to EPTA DR2.

### 2.3. Differences between Data Sets

Figure 1 shows each pulsar’s TOAs from NG15 and EPTA DR2, visualized as a function of time and frequency. NG15 has fairly consistent multifrequency coverage over time and a low-frequency floor of 724 MHz for all six pulsars. Since EPTA DR2 has a longer time span and is composed of more telescopes, its level of radio-frequency coverage is more varied. In particular, EPTA DR2 includes TOAs down to 323 MHz from WSRT for PSRs J0613−0200, J1012+5307, J1713+0747, and J1744−1134, although not
<table><caption>**Table 1**
Definitions of Terms and Parameters Commonly Used throughout This Work</caption><thead><tr><th>Category</th><th>Term</th><th>Definition</th></tr></thead><tbody><tr><td>Acronyms (General)</td><td>PTA, TOA, GWB<br>GP, MCMC<br>PSD, NG15</td><td>Pulsar timing array, time of arrival, gravitational wave background<br>Gaussian process, Markov Chain Monte Carlo<br>Power spectral density, NANOGrav 15 yr data set (Agazie et al. 2023c)</td></tr><tr><td>Acronyms (Noise processes)</td><td>DM, ΔDM<br>DMX/DMGP<br>ISM, SW<br>ARN, RN, WN<br>DMX/DMGP/<br>CustomGP</td><td>Dispersion measure, deviation from fiducial DM value<br>Piecewise-constant model/Fourier-basis GP model for DM estimation<br>Interstellar medium, solar wind<br>Achromatic red noise, red noise, white noise<br>Labels full pulsar noise model using DMX/DMGP/<br>DMGP plus additional chromatic terms</td></tr><tr><td>Symbols (General)</td><td>N<sub>ToA</sub>, T<sub>psr</sub><br>Δt, δt<br>δt', ρ<br>b̄/η̄, N<sub>freqs</sub><br>B, N/μ<br>E<sub>1.2</sub></td><td>Number of TOAs, total observation time span of pulsar<br>Time delay (generic), estimated time delay (Appendix D only)<br>Timing residual vector (Section 4 only), timing residual power<br>GP coefficient/hyperparameter vector, number of frequencies in GP Fourier basis<br>Bayes factor, normal/uniform distribution<br>First and second exponential timing events in PSR J1713+0747</td></tr><tr><td>Symbols (Model parameters)</td><td>A/γ/f<br>ν, χ<br>n<sub>Earth</sub><br>A<sub>E</sub>/τ<sub>E</sub>/t<sub>E</sub></td><td>Spectral amplitude/index/frequency<br>Radio frequency, chromatic radio-frequency scaling index<br>Estimated local electron number density<br>Amplitude/timescale/initial time of decaying exponential in PSR J1713+0747</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** TOAs for six pulsars from the NANOGrav 15 yr data set (Agazie et al. 2023c) and EPTA DR2 (EPTA Collaboration et al. 2023b), visualized as a function of time and radio frequency. The radio frequencies used to collect TOAs are important for determining how well each pulsar's chromatic noise can be constrained. Dashed lines are used to mark the times of two known chromatic timing events in PSR J1713+0747 (Lam et al. 2018).</figcaption></figure>

across the entire time span. These differences in radio-frequency coverage will become relevant when accounting for discrepancies in chromatic noise characterization between the two data sets, since chromatic processes induced from the interstellar medium (ISM) introduce larger delays at low radio frequencies.

### 3. Single-pulsar Noise Budget

Here we discuss various sources of noise relevant to millisecond pulsars, which form the basis of NG15. Most pulsars experience rotational irregularities which manifest as spin noise, an ARN process (Verbiest et al. 2009; Cordes & Shannon 2010; Shannon & Cordes 2010). Spin noise is found to be small in millisecond pulsars due to their very small spin frequency derivatives, with notable exceptions such as PSR B1937+21 (Shannon & Cordes 2010; Agazie et al. 2023d, hereafter NG15_Noise). Nevertheless, accounting for spin noise is very important since the GWB also manifests as an ARN process in single pulsars. NG15_Noise finds that 12 pulsars in NG15 still have significant detections of ARN in the presence of a GWB signal. One of these 12, PSR J1012+5307, is among the six we include here. The remaining five pulsars display significant ARN which does not persist in the presence of a GWB signal using NG15, i.e., these pulsars currently show little evidence for both intrinsic spin noise and a GWB signal (NG15_Noise). As such, changes to these five pulsars' ARN properties could impact GWB inferences.

To maximize PTA sensitivity to GWs, it is also important to account for chromatic noise. A major source of chromatic noise comes from DM variations (Rankin & Roberts 1971; You et al. 2007; Jones et al. 2017). DM is defined as the integrated electron column density between the Earth and the pulsar,

$$DM = \int_{0}^{d} n_{e}(l) dl, \qquad (1)$$
where $n_e$ is the free electron density, $l$ defines the Earth-pulsar line of sight, and $d$ is the distance to the pulsar (Lorimer & Kramer 2012). DM may undergo linear, annual, and/or stochastic variations due to the changing line of sight through the ionized ISM (Lam et al. 2016; Jones et al. 2017) and the solar wind (SW; Madison et al. 2019; Tiburzi et al. 2021). DM introduces the following timing delay:

$$ \Delta t_{\text{DM}} = \frac{e^2}{2\pi m_e c} \frac{\text{DM}}{\nu^2} \qquad (2) $$

$$ = (4.15 \ \mu\text{s}) \left( \frac{\text{DM}}{10^{-3} \ \text{pc} \ \text{cm}^{-3}} \right) \left( \frac{1 \ \text{GHz}}{\nu} \right)^2, \qquad (3) $$

where $\nu$ is the radio frequency of the pulse. Millisecond pulsars from NG15 have peak-to-peak DM variations ranging from $\Delta DM \sim 10^{-4}-10^{-2} \ \text{pc cm}^{-3}$ over the full data time span (NG15). As such, the time delays introduced by DM variations often dominate over ARN processes in millisecond pulsars.

The $\nu^{-2}$ dependence in Equation (2) allows DM variations to be decoupled from other processes affecting pulsar timing. However, a number of factors can systematically bias DM measurements, including asynchronous measurements across radio-frequency bands (Lam et al. 2015), finite observing bandwidths (Sosa Fiscella et al. 2024), or the presence of additional chromatic effects (NG15_Noise). Inaccurate DM values may result in ARN with a low spectral index $\gamma_{RN}$ (Cordes & Shannon 2010; Lam et al. 2015; NG15_Noise), which may reduce the pulsar's sensitivity to GW signals.

A secondary source of chromatic noise is interstellar scattering, which results from frequency-dependent refraction of radio pulses through an inhomogeneous ISM (Cordes & Rickett 1998; Hemberger & Stinebring 2008; Lorimer & Kramer 2012). The first-order effect of scattering is to delay the TOA by $\Delta t \propto \nu^{-4}$ for a Gaussian inhomogeneity (Lang 1971) or $\Delta t \propto \nu^{-4.4}$ for a Kolmogorov-turbulent medium (Romani et al. 1986). However, this frequency scaling may vary more substantially depending on the geometry of the scattering medium, with observed scalings ranging from $\nu^{-0.7}$ to $\nu^{-5.6}$ for different pulsar lines of sight (Lewandowski et al. 2015; Turner et al. 2021), and scalings predicted as high as $\nu^{-6.4}$ from simulations (Shannon & Cordes 2017). While noise from time-variable scattering is not expected to affect pulsar timing as strongly as DM variations, many high-DM pulsars have been observed by NANOGrav with large scattering tails (Alam et al. 2021; NG15_Noise). Unmitigated scattering variations may be absorbed as excess white noise (WN), DM noise, ARN, or some combination thereof (Lentati et al. 2016; Shapiro-Albert et al. 2021; NG15_Noise).

Changes to the pulse profile itself also present a potential source of noise in millisecond pulsars. These may result, for example, from polarization calibration errors, scatter broadening of the pulse profile, or intrinsic changes in the pulsar magnetosphere (NG15_Noise). The frequency dependence of pulse profiles is already accounted for in timing analyses using log-polynomial functions in frequency, parameterized by “frequency-dependent” (FD) parameters (NG15). However, these do not have built in time dependence. Brook et al. (2018) observed long-term pulse profile variability in several pulsars from the NANOGrav 11 yr data set (Arzoumanian et al. 2018). Of the pulsars we study here, PSR J1713+0747 was identified to have high profile variability on short timescales. Pulsars may

also experience transient changes to their pulse profiles, with associated delays to their timing residuals (Shannon et al. 2016; Goncharov et al. 2021). A dramatic pulse profile change took place for PSR J1713+0747 in early 2021 (Singha et al. 2021), which introduced chromatic timing delays scaling nonmonotonically with radio frequency (Jennings et al. 2024). While this event is not in NG15, PSR J1713+0747 features two weaker chromatic timing events at earlier times (Lam et al. 2018), with evidence of associated changes to the pulse profile found in one or both events (Goncharov et al. 2021; Lin et al. 2021).

Chromatic noise is not typically mitigated prior to TOA generation in PTA pipelines, but it can be included in the noise model as a GP. Existing alternatives include wideband (Pennucci et al. 2014; Liu et al. 2014) and profile domain timing (Lentati et al. 2017), where chromatic noise mitigation is applied at a different stage. NANOGrav's standard noise model mitigates DM variations using DMX timing model parameters, which fit for the DM value at each observation epoch comprised of multiband observations (Demorest et al. 2013; NANOGrav Collaboration et al. 2015; NG15_Noise). An alternative to DMX is to treat DM variations as a red noise (RN) process, using the formalism of GPs (“DMGP”; e.g., Lentati et al. 2014; van Haasteren & Vallisneri 2014; IPTA_comp), alongside a SW model (e.g., Hazboun et al. 2022). Often a power-law prior is imposed on the DM power spectral density (PSD), constraining how the DM variations may vary over time. This choice is physically motivated by Kolmogorov turbulence in the ISM, which predicts a power-law PSD for DM variations with $\gamma_{DM} = 8/3$ (Keith et al. 2013).

Meanwhile, time-variable scattering is not always modeled explicitly. A power-law GP model with a $\nu^{-4}$ frequency scaling was first introduced by Lam et al. (2018) to mitigate possible scattering delays in PSR J1713+0747. Since then, this GP model has been commonly used as a first-order correction for scattering delays (e.g., Goncharov et al. 2021; Chalumeau et al. 2022; Srivastava et al. 2023). From here on we refer to the $\nu^{-4}$ GP as a “scattering-like” chromatic noise process, since other unmodeled chromatic processes could hypothetically result in preference for this model during the Bayesian analysis. Alternatively, one could attempt to mitigate scattering without assuming a particular frequency scaling, or search for a generic chromatic noise process, by fitting RN processes isolated to single radio-frequency observing bands (i.e., “band” noise; Lentati et al. 2016; Goncharov et al. 2021; Chalumeau et al. 2022.)

Additional chromatic noise processes not considered here include decorrelation of pulse jitter over radio frequency (Lam et al. 2019), frequency-dependent DM due to multipath propagation effects (Cordes et al. 2016), and low-level radio-frequency interference. While these processes may introduce additional timing errors, they are also difficult to measure and expected to primarily affect the WN level in most pulsars (NG15_Noise). The ionosphere may also become a significant source of noise at very low radio frequencies ($\nu < 100$ MHz; de Gasperin et al. 2018).

## 4. Gaussian Process Models

We contextualize our models throughout this work in the framework of GPs. GPs are flexible mathematical models which represent a series of values $\bar{y}$ (e.g., pulsar timing
residuals) as samples from a multivariate Gaussian distribution:

$$p(\vec{y}) = \mathcal{N}(\vec{m}, \mathbf{C}), \quad (4)$$

where $\vec{m}$ is the mean vector and $\mathbf{C}$ is the covariance matrix (Rasmussen & Williams 2006). GPs are particularly useful for modeling stochastic signals in astronomical time series, as otherwise unknown information about the functional form of the stochastic signal is represented by the off-diagonal elements of $\mathbf{C}$ (Aigrain & Foreman-Mackey 2023).

We summarize the implementation of GPs used in PTA analyses (e.g., van Haasteren & Vallisneri 2014; Taylor 2021; NG15_Noise). Our set of timing residuals $\delta\vec{t}$ are represented as a combination of deterministic terms (comprising the mean vector $\vec{m}$), WN terms, and rank-reduced terms (which are themselves GPs). The WN and rank-reduced GPs are used to build the covariance matrix as

$$\mathbf{C} = \mathbf{N} + \mathbf{TBT}^T. \quad (5)$$

Here $\mathbf{N}$ is a block-diagonal WN matrix. $\mathbf{T}$ is a $(N_{\text{TOA}} \times N_{\text{b}})$ design matrix containing a series of $N_{\text{b}}$ basis functions. $\mathbf{B} = \langle \vec{b} \vec{b}^T \rangle$ is a diagonal matrix encoding the variance of GP coefficients $\vec{b}$, which are given Gaussian (hyper)priors $p(\vec{b}|\vec{\eta}) = \mathcal{N}(0, \mathbf{B}(\vec{\eta}))$ with hyperparameters $\vec{\eta}$. We first obtain hyperparameter posteriors $p(\vec{\eta}|\delta\vec{t})$ using Markov Chain Monte Carlo (MCMC) sampling, while marginalizing over $p(\vec{b}|\vec{\eta})$. We subsequently draw samples of our GP coefficients from the conditional probability distribution $p(\vec{b}|\vec{\eta}, \delta\vec{t})$ (see, e.g., Laal et al. 2023; Meyers et al. 2023).

Next we define the components of our noise model, which we apply to all six pulsars. Additional deterministic signals applied only to PSR J1713+0747 are presented in Section 4.5.

### 4.1. Red Noise

Chromatic and achromatic RN processes are defined as rank-reduced GPs using a sine-cosine Fourier basis (Lentati et al. 2013). The PSD of the Fourier coefficients $\vec{b}$ are parameterized by a power-law prior with hyperparameters $\vec{\eta} = A, \gamma$ (van Haasteren & Levin 2013):

$$S_{\text{RN}}(f_i) = \frac{A_{\text{RN}}^2}{12\pi^2} \left( \frac{f_i}{\text{yr}^{-1}} \right)^{-\gamma_{\text{RN}}} \text{yr}^3. \quad (6)$$

Here we use a log-uniform prior $\log_{10}U(10^{-18}, 10^{-10})$ to sample $A$, and a uniform prior $U(0, 7)$ to sample $\gamma$. We set the frequencies of the Fourier basis as integer multiples of the reciprocal of the pulsar's observation time $f_i = i/T_{\text{psr}}$, where $i = 1, 2, 3, \dots, N_{\text{freqs}}$, where $2N_{\text{freqs}}$ is the size of the Fourier basis, and $T_{\text{psr}}$ is the pulsar's observation time span.

We also model ARN using a “free-spectral” PSD, where the power at each frequency $f_i$ is a separate parameter $\rho_i$ (Lentati et al. 2013). To match EPTA_Noise, we use a log-uniform prior $\log_{10}U(10^{-10}, 10^{-4})$ to sample each $\rho_i$ (in units of seconds). This method is useful to gauge the presence of noise across the spectrum, without making any assumptions about the form of the PSD.

To specify chromatic noise, we scale the Fourier basis by a frequency-dependent factor:

$$\mathbf{T}_{\text{RN}} \rightarrow \mathbf{T}_{\text{RN}}\left(\frac{\nu}{1400 \text{ MHz}}\right)^{-\chi}, \quad (7)$$

where $\nu$ is the radio frequency and $\chi$ is the chromatic index (Goncharov et al. 2021). We select larger values of $N_{\text{freqs}}$ for chromatic processes, which generally have shallow spectra and are more easily decoupled from WN at higher frequencies than ARN. We use the following values of these parameters by default: $\chi = 0, N_{\text{freqs}} = 30$ for ARN; $\chi = 2, N_{\text{freqs}} = 100$ for DM noise; and $\chi = 4, N_{\text{freqs}} = 150$ for scattering-like chromatic noise. The base values for $N_{\text{freqs}}$ are chosen for consistency with other PTA literature, where $N_{\text{freqs}} = 30$ (100) for ARN (DM noise) is a common, albeit arbitrary, choice, and $N_{\text{freqs}} = 150$ was the only value favored for scattering-like noise in Chalumeau et al. (2022).

We also calculate Bayes factors (Kass & Raftery 1995) comparing a model with RN versus a model without RN. We use these RN Bayes factors, $\mathcal{B}^{\text{RN}}$, to quantify the statistical evidence for each RN process under the given modeling assumptions. We calculate each $\mathcal{B}^{\text{RN}}$ using the Savage-Dickey density ratio (Dickey 1971), approximated as the prior-to-posterior ratio at the lower bound of the prior distribution ($\log_{10}A_{\text{RN}} = -18$). In many cases, $\mathcal{B}^{\text{RN}}$ cannot be calculated using the Savage-Dickey approximation due to lack of MCMC samples consistent with $\log_{10}A_{\text{RN}} = -18$. In these cases we place a lower limit of $\log_{10} \mathcal{B}^{\text{RN}} > 3$, as they correspond to a statistically significant detection of RN.

### 4.2. Solar Wind

The SW may contribute substantially to DM variations as the Earth-pulsar line of sight cuts through different regions of the heliosphere over the course of each year (Lommen et al. 2006). For a time-independent $1/r^2$ SW density profile, the SW's contribution to DM variations goes as

$$\text{DM}_{\text{SW}} = n_{\text{Earth}}(1 \text{ au}) \frac{\pi - \theta_i}{\sin \theta_i}, \quad (8)$$

where $n_{\text{Earth}}$ is the SW free electron density at 1 au, and $\theta_i$ is the angle between the Earth-Sun line of sight and the Earth-pulsar line of sight (Splaver et al. 2005). For pulsars close to the ecliptic, $\sin \theta_i$ may become very small at the nearest conjunction of the Sun and the pulsar, leading to larger and more peaked annual spikes in the DM time series (Madison et al. 2019; Hazboun et al. 2022).

We include this SW model as a deterministic signal in our Bayesian analysis, with $n_{\text{Earth}}$ fit independently for each pulsar using a uniform prior $U(0, 30)$ in units of cubic centimeters. Modifications to this model can be made to account for time-dependent or nonspherical SW density profiles (You et al. 2007; Hazboun et al. 2022). As such, we emphasize this model acts only as a first-order correction for SW effects. Similarly as for RN, we use the Savage-Dickey density ratio to calculate Bayes factors, gauging how sensitive each pulsar is to detecting the SW (see Appendix A).

### 4.3. Timing Model Perturbations

To account for covariances between noise model parameters and timing model parameters, we vary a linear approximation of the timing model (van Haasteren & Levin 2013). The linearized timing response from perturbations to each of the best-fit timing model parameters make up the timing model design matrix (van Haasteren 2013; Taylor 2021). The coefficients $\vec{b}$ corresponding to the amplitude of these
perturbations are given Gaussian priors of effectively infinite variance to mimic improper uniform priors.

When using DMX timing model parameters, these priors ensure the DM estimated at each epoch is constrained only by the fit to the data. When instead modeling DM variations as a RN process, we remove the DMX parameters and replace them with the DM1 and DM2 timing model parameters, which parameterize a linear and quadratic trend in DM over time (Lentati et al. 2014). These are needed to account for long-timescale DM variations below the fundamental frequency $f_1 = 1/T_{psr}$ of the DM RN process.

### 4.4. White Noise

TOA uncertainties are initially estimated based on radio-meter noise, which affects the pulse profile signal-to-noise ratio (Lorimer & Kramer 2012). We model WN by modifying these estimated uncertainties using three parameters: EFAC ($\mathcal{F}$), EQUAD ($\mathcal{Q}$), and ECORR ($\mathcal{J}$). These parameters are designed to model errors in estimates of template-matching uncertainties, independent measurement noise, and pulse jitter (Cordes & Downs 1985; Cordes & Shannon 2010), respectively. An independent set of these parameters is fit for each unique receiver/ backprop pair (NG15\_Noise). Mathematically, these effects are represented in the following elements of $N$:

$$ \langle n_{i,\mu} n_{j,\nu} \rangle = \mathcal{F}_{\mu}^2 (\sigma_i^2 \delta_{i,j} \delta_{\mu,\nu} + \mathcal{Q}_{\mu}^2 \delta_{i,j} \delta_{\mu,\nu}) + \mathcal{J}_{\mu}^2 \delta_{e(i),e(j)} \delta_{\mu,\nu}, \quad (9) $$

where $i, j$ label each TOA, $\mu, \nu$ label each receiver/ backend pair, $e(i), e(j)$ label all TOAs within the same observation epoch, and $\sigma_i$ are the original TOA errors. The $\mathcal{J}_{\mu}$ terms make $N$ a block-diagonal matrix, whose inverse we calculate using the Sherman-Morrison formula. To match EPTA\_Noise, we use a uniform prior $\mathcal{U}(0.1, 5)$ for all $\mathcal{F}_{\mu}$ and a log-uniform prior $\log_{10} \mathcal{U}(10^{-9}, 10^{-5})$ in units of seconds for all $\mathcal{Q}_{\mu}$ and $\mathcal{J}_{\mu}$.

### 4.5. PSR J1713+0747 Chromatic Events

PSR J1713+0747 exhibited unusual timing events near MJDs 54750 and 57510 (Keith et al. 2013; Demorest et al. 2013; Lam et al. 2018). These events each manifest as a sudden “dip” in the apparent DM value which gradually returns back to a previous level. GP analyses of PSR J1713+0747 (e.g., Lam et al. 2018; Hazboun et al. 2020b; Goncharov et al. 2021) have modeled these noise transients using decaying exponential functions:

$$ \Delta t = -A_E \Theta(t_E) \exp\left(-\frac{t}{\tau_E}\right) \left(\frac{\nu}{1400 \text{ MHz}}\right)^{-\chi_E}, \quad (10) $$

where $A_E$ is the amplitude, $\Theta(t_E)$ is a Heaviside step function centered at the initial time of the event, $\tau_E$ is the decay timescale, $\chi_E$ is the chromatic scaling index, and $\nu$ is the radio frequency. We further refer to each dip as $E_1$ and $E_2$, respectively. We use log-uniform priors $\log_{10} \mathcal{U}(10^{-10}, 10^{-2})$ for $A_{E_{1,2}}$ and $\log_{10} \mathcal{U}(10^0, 10^{3.5})$ for $\chi_{E_{1,2}}$. We use uniform priors $\mathcal{U}(54650, 54850)$ for $t_{E_1}$ and $\mathcal{U}(57490, 57530)$ for $t_{E_2}$. We treat the chromatic indices $\chi_{E_{1,2}}$ in two ways: (i) we hold them at fixed $\chi_{E_{1,2}} = 2$ to model $E_{1,2}$ as DM events (Lam et al. 2018), and (ii) we sample $\chi_{E_{1,2}}$ with a uniform prior $\mathcal{U}(0, 7)$ (see, e.g., Goncharov et al. 2021; Chalumeau et al. 2022).

<table><caption>**Table 2**
The Three Noise Models We Use for Each Pulsar from NG15</caption><tr><td>DMX</td><td>RN (achromatic); TM (DMX); WN</td></tr><tr><td>DMGP</td><td>RN (achromatic, DM); SW; TM (DM1, DM2); WN; DM exp. dips (J1713+0747 only)</td></tr><tr><td>CustomGP</td><td>RN (achromatic, DM, scattering-like); SW; TM (DM1, DM2); WN; Custom exp. dips (J1713+0747 only)</td></tr></table>

**Notes.** Bolded terms indicate key changes from one noise model to the next. Model components are detailed throughout Section 4.

### 4.6. Composite Models

Table 2 summarizes the set of three composite noise models we apply to these six pulsars in NG15, labeled DMX, DMGP, and CustomGP. DMX labels the standard NANOGray noise model, which includes WN, ARN, and DMX parameters. DMGP labels a model in which DMX parameters are removed and replaced with the following components: a DM RN GP, the DM1 and DM2 parameters, and the deterministic SW model, as well as deterministic exponential dips scaling as $\Delta t \sim \nu^{-2}$ for PSR J1713+0747. We use the comparison of DMX and DMGP to assess if DMX parameters produce similar results as time-correlated DM models for these pulsars.

CustomGP extends the DMGP model by including additional nondispersive chromatic noise processes used by EPTA\_Noise for pulsars from EPTA DR2. Namely, CustomGP includes the addition of a $\chi = 4$ scattering-like chromatic RN process for all pulsars. It also uses $N_{freqs} = 150$ for ARN in PSR J1012+5307 and for DM noise in PSR J1909−3744 as these processes favored a large number of Fourier modes in EPTA\_Noise. Furthermore, the chromatic indices $\chi_{E_{1,2}}$ of the deterministic dips in PSR J1713+0747 are allowed to vary as free parameters instead of being fixed to $\chi_{E_{1,2}} = 2$.

Accounting for all of these effects, our CustomGP model tests for the same extent of noise processes as included in EPTA\_Noise. However, CustomGP is not quite equivalent to the models from EPTA\_Noise. This is partly due to intrinsic differences between NG15 and EPTA DR2. For instance, we do not use the same values of $N_{freqs}$ for RN processes as were selected in EPTA\_Noise. Each data set features a different observation time span and cadence for each pulsar, therefore the values $N_{freqs}$ favored for pulsars from EPTA DR2 are unlikely to be optimal for the same pulsars from NG15. Instead, we ensure the ARN and chromatic noise spectra extend to at least the same high-frequency cutoff $f_{max} = N_{freqs}/T_{psr}$ as the favored models from EPTA\_Noise. One additional difference is that where EPTA\_Noise fixes $n_{Earth} = 7.9 \text{ cm}^{-3}$ (Madison et al. 2019), we allow it to vary as a free parameter for each pulsar separately. Furthermore, EPTA\_Noise only includes scattering-like chromatic noise for PSR J1600−3053, and does not include ARN for PSR J1600−3053, while we have both ARN and scattering-like chromatic noise processes in the six NG15 pulsars using CustomGP. To account for these differences, we modified the favored models from EPTA\_Noise to include ARN, include scattering-like chromatic noise, and fit $n_{Earth}$ as a free parameter for all six pulsars using EPTA DR2. These modifications were made to ensure fair comparison with NG15 using CustomGP, but this
<table><caption>**Table 3**
Estimated Noise Parameters and Bayes Factors for Six Pulsars in NG15 under All Three Modeling Assumptions</caption><thead><tr><th rowspan="2">Signal/Parameter</th><th rowspan="2">Model</th><th colspan="6">Pulsar</th></tr><tr><th>J0613-0200</th><th>J1012+5307</th><th>J1600-3053</th><th>J1713+0747</th><th>J1744-1134</th><th>J1909-3744</th></tr></thead><tbody><tr><td rowspan="12">Achromatic red noise</td><td>log<sub>10</sub>A</td><td>DMX<br>DMGP<br>CustomGP</td><td>-13.8<sup>+0.3</sup><sub>-0.4</sub><br>-13.8<sup>+0.2</sup><sub>-0.3</sub><br>-14.2<sup>+0.4</sup><sub>-0.6</sub></td><td>-12.64<sup>+0.06</sup><sub>-0.06</sub><br>-12.81<sup>+0.07</sup><sub>-0.07</sub><br>-12.85<sup>+0.05</sup><sub>-0.05</sub></td><td>-13.5<sup>+0.2</sup><sub>-0.6</sub><br>-14.5<sup>+0.8</sup><sub>-0.9</sub><br>-14.2<sup>+0.5</sup><sub>-0.7</sub></td><td>-14.1<sup>+0.1</sup><sub>-0.1</sub><br>-14.1<sup>+0.1</sup><sub>-0.2</sub><br>-14.7<sup>+0.3</sup><sub>-0.5</sub></td><td>-14.1<sup>+0.4</sup><sub>-0.6</sub><br>-14.3<sup>+0.4</sup><sub>-0.6</sub><br>-15.2<sup>+1.0</sup><sub>-1.8</sub></td><td>-14.5<sup>+0.3</sup><sub>-0.4</sub><br>-14.6<sup>+0.3</sup><sub>-0.4</sub><br>-14.7<sup>+0.3</sup><sub>-0.4</sub></td></tr><tr><td>γ</td><td>DMX<br>DMGP<br>CustomGP</td><td>3.1<sup>+0.9</sup><sub>-0.7</sub><br>3.1<sup>+0.8</sup><sub>-0.6</sub><br>4.0<sup>+1.3</sup><sub>-1.0</sub></td><td>0.8<sup>+0.3</sup><sub>-0.3</sub><br>1.1<sup>+0.3</sup><sub>-0.3</sub><br>1.2<sup>+0.2</sup><sub>-0.2</sub></td><td>1.7<sup>+1.6</sup><sub>-0.8</sub><br>3.8<sup>+2.0</sup><sub>-1.9</sub><br>3.9<sup>+1.7</sup><sub>-1.4</sub></td><td>2.6<sup>+0.5</sup><sub>-0.4</sub><br>2.6<sup>+0.5</sup><sub>-0.4</sub><br>3.5<sup>+1.2</sup><sub>-0.9</sub></td><td>3.6<sup>+1.4</sup><sub>-1.2</sub><br>3.8<sup>+1.4</sup><sub>-1.2</sub><br>3.5<sup>+1.9</sup><sub>-2.1</sub></td><td>4.1<sup>+1.0</sup><sub>-0.9</sub><br>4.1<sup>+1.0</sup><sub>-0.8</sub><br>4.4<sup>+1.0</sup><sub>-0.8</sub></td></tr><tr><td>log<sub>10</sub>B</td><td>DMX<br>DMGP<br>CustomGP</td><td>>3<br>>3<br>>3</td><td>>3<br>>3<br>>3</td><td>2.4<br>>3<br>1.7</td><td>>3<br>>3<br>>3</td><td>2.1<br>1.6<br>-0.1</td><td>>3<br>>3<br>>3</td></tr><tr><td rowspan="6">DM noise</td><td>log<sub>10</sub>A</td><td>DMGP<br>CustomGP</td><td>-13.38<sup>+0.04</sup><sub>-0.04</sub><br>-13.7<sup>+0.1</sup><sub>-0.3</sub></td><td>-13.17<sup>+0.05</sup><sub>-0.04</sub><br>-14.1<sup>+0.6</sup><sub>-2.3</sub></td><td>-13.12<sup>+0.05</sup><sub>-0.05</sub><br>-13.20<sup>+0.08</sup><sub>-0.11</sub></td><td>-13.82<sup>+0.05</sup><sub>-0.05</sub><br>-13.80<sup>+0.05</sup><sub>-0.05</sub></td><td>-13.46<sup>+0.04</sup><sub>-0.04</sub><br>-13.55<sup>+0.07</sup><sub>-0.11</sub></td><td>-13.63<sup>+0.04</sup><sub>-0.04</sub><br>-13.67<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>γ</td><td>DMGP<br>CustomGP</td><td>1.9<sup>+0.2</sup><sub>-0.2</sub><br>2.8<sup>+0.7</sup><sub>-0.4</sub></td><td>1.4<sup>+0.2</sup><sub>-0.2</sub><br>2.6<sup>+1.9</sup><sub>-0.9</sub></td><td>2.2<sup>+0.2</sup><sub>-0.2</sub><br>2.6<sup>+0.3</sup><sub>-0.2</sub></td><td>1.8<sup>+0.2</sup><sub>-0.2</sub><br>1.9<sup>+0.2</sup><sub>-0.2</sub></td><td>1.5<sup>+0.2</sup><sub>-0.2</sub><br>1.7<sup>+0.4</sup><sub>-0.2</sub></td><td>1.5<sup>+0.1</sup><sub>-0.1</sub><br>1.6<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>log<sub>10</sub>B</td><td>DMGP<br>CustomGP</td><td>>3<br>>3</td><td>>3<br>0.1</td><td>>3<br>>3</td><td>>3<br>>3</td><td>>3<br>>3</td><td>>3<br>>3</td></tr><tr><td rowspan="3">Scattering-like chromatic noise</td><td>log<sub>10</sub>A</td><td>CustomGP</td><td>-14.01<sup>+0.05</sup><sub>-0.06</sub></td><td>-13.75<sup>+0.05</sup><sub>-0.06</sub></td><td>-13.59<sup>+0.05</sup><sub>-0.05</sub></td><td>-14.22<sup>+0.06</sup><sub>-0.05</sub></td><td>-14.23<sup>+0.09</sup><sub>-0.13</sub></td><td>-14.67<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>γ</td><td>CustomGP</td><td>1.5<sup>+0.2</sup><sub>-0.2</sub></td><td>1.5<sup>+0.2</sup><sub>-0.2</sub></td><td>1.6<sup>+0.2</sup><sub>-0.2</sub></td><td>1.3<sup>+0.2</sup><sub>-0.2</sub></td><td>1.3<sup>+0.4</sup><sub>-0.3</sub></td><td>0.5<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>log<sub>10</sub>B</td><td>CustomGP</td><td>>3</td><td>>3</td><td>>3</td><td>>3</td><td>3.0</td><td>2.6</td></tr></tbody></table>

**Notes.** Noise parameters are presented using the median and 68.3% Bayesian credible intervals (referenced here as $1\sigma$ regions), and Bayes factors indicating statistical detection significance of the given signal are calculated from our posterior distributions using the Savage-Dickey approximation. If a parameter is bolded, that means the parameter's $1\sigma$ region estimated under the current model is inconsistent with the $1\sigma$ region estimated under the previous model (from one row above). If a Bayes factor is bolded, that means the Bayes factor estimated under the current model is at least an order of magnitude different from the Bayes factor estimated under the previous model.

does not noticeably alter the ARN parameter estimation results from EPTA_Noise.

## 5. Results

Table 3 presents the medians and 68.3% Bayesian credible intervals ($1\sigma$ regions) of inferred noise parameters and Bayes factors for each noise process. These are tabulated for our six pulsars in NG15 under the three different noise models: DMX, DMGP, and CustomGP. Bolded parameter values indicate cases where going from one model to the next results in a discrepancy (significant at a >1$\sigma$ level) between noise parameters. Bolded Bayes factors indicate cases where the Bayes factor changed by over an order of magnitude, i.e., if the detection significance of a noise process has substantially dropped. The inferred SW electron density from each pulsar is reported in Appendix A. Interestingly, using NG15 we find the presence of $\chi=4$ scattering-like chromatic noise is supported by a Bayes factor $log_{10}B^{chrom} \ge 2.6$ for all six pulsars using model CustomGP, despite a significant scattering-like chromatic noise detection in only one pulsar from EPTA DR2 (PSR J1600-3053; Chalumeau et al. 2022; EPTA_Noise). For several pulsars, these chromatic noise amplitudes are substantially higher using NG15 than the upper limits set using EPTA DR2. Including scattering-like chromatic noise in CustomGP changes the estimated DM noise parameters by >1$\sigma$ for PSRs J0613-0200 and J1012+5307. A deeper investigation comparing with measurements of pulse

broadening or scintillation would be needed to confirm the origin of these variations.

To complement Table 3, Figure 2 shows the 2D posterior distributions for $log_{10}A_{RN}$ and $\gamma_{RN}$ for our six-pulsar sample of NG15, under all three modeling assumptions. All contours enclose 68.3% ($1\sigma$) and 95.4% ($2\sigma$) 2D credible intervals. NG15_Noise showed that out of these six pulsars, only PSR J1012+5307 shows evidence for additional ARN on top of a GWB signal. To highlight this, the GWB parameters, inferred from NG15_GWB using the DMX model, are also shown in Figure 2. With the exception of PSR J1012+5307, the ARN and GWB parameters are similar for every pulsar, indicating the GWB makes up a substantial portion of these pulsars' ARN budget.

We first consider how the recovered ARN parameters are affected by switching from DMX to DMGP. Each model results in notably different ARN posteriors for PSRs J1012+5307 and J1600-3053, but no major differences for the remaining pulsars (Figure 2). The change to $log_{10}A_{RN}$ for PSR J1012+5307 is significant at a >1$\sigma$ level (Table 3). These findings coincide with much larger differences in DM recovery for PSRs J1012+5307 and J1600-3053 than the remaining pulsars (Appendix B).

Transitioning from DMX to CustomGP yields further changes to the single-pulsar ARN parameters. In general, switching to CustomGP results in lower $log_{10}A_{RN}$ and higher $\gamma_{RN}$ (Figure 2). As a counterexample, PSR J1600-3053 favors a slightly higher $log_{10}A_{RN}$ going from DMGP to CustomGP. Overall, the ARN properties of PSR J1909-3744 remain the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Accounting for nondispersive chromatic noise using the *CustomGP* model noticeably affects ARN posteriors for multiple pulsars within NG15. Posterior distributions for the six pulsars' ARN parameters using the DMX (solid purple), DMGP (dashed green), and *CustomGP* (solid blue) models. Since the GWB makes up a portion of all pulsars' ARN budget, the GWB parameters inferred from NG15_GWB using the DMX model for all pulsars is included in gray. The inferred ARN posteriors change the most substantially for PSRs J1713+0747 and J1012+5307 using *CustomGP*. PSRs J0613-0200, J1600-3053, and J1744-1134 also feature noticeable changes to their ARN parameters at a less significant level. PSR J1909-3744's ARN parameters are the least affected by the choice of noise model. The apparent shift of PSR J1713+0747's ARN parameters away from the inferred GWB parameters indicates its choice of noise model is likely to affect GWB spectral characterization.</figcaption></figure>

most similar under all three models. Meanwhile, PSR J1713+0747's change is the most dramatic, as it experiences a significant ($>1\sigma$) decrease to $\log_{10} A_{\text{RN}}$ and supports a much broader range of $\gamma_{\text{RN}}$ values. Furthermore, PSR J1713+0747's ARN parameters using DMX and DMGP are highly constrained near the GWB parameters measured in NG15_GWB, but become less consistent with the measured GWB parameters when using *CustomGP*, favoring instead a lower amplitude and higher spectral index. Notably, the alternative noise model used in NG15_GWB, which included the *CustomGP* model for PSR J1713+0747 alongside DMGP for the remaining pulsars, also resulted in a shift toward lower amplitude and higher spectral index of the common noise. As such, our results signpost PSR J1713+0747's noise model as a strong contributor to this change. However, the GWB parameter inference is dependent on information from 61 additional pulsars not studied here. As such, a direct quantification of these impacts on GWB characterization will require a more careful analysis using the full PTA, which is the subject of an upcoming work.

To assess the performance of the models, we next compare the following cases on a pulsar-by-pulsar basis: (i) DMX applied to NG15, (ii) *CustomGP* applied to NG15, and (iii) *CustomGP* applied to EPTA DR2. We do not include the model DMGP in this comparison as it is intermediary to the more disparate models DMX and *CustomGP* (although we do compare DM estimates using DMX and DMGP in Appendix B). We present a separate figure for each pulsar (starting from

Figure 3), displaying parameter posterior parameters, spectra, and time-domain GP realizations for ARN, DM variations, and $\chi^2 = 4$ scattering-like chromatic noise. For model DMX, we display only ARN and the time series of DMX parameters. Time-domain realizations of DM variations using *CustomGP* include all stochastic and deterministic contributions to DM to allow fair comparison with the DMX time series. Spectra and time-domain realizations are visualized using medians and 68% (1$\sigma$) Bayesian credible intervals, each computed from 100 GP realizations. In cases where the 68% regions were difficult to make out by eye, we additionally added the 95% regions, distinguished from the 68% regions with a different plot style. Furthermore, since EPTA DR2 has a longer time span than NG15 for several pulsars (Figure 1), this can result in differing ARN properties between EPTA DR2 and NG15 for the same pulsar. To provide a more useful comparison to NG15, we generated time-domain realizations from EPTA DR2 using the full data time spans for each pulsar, but excised the portion of the realizations before the start of the NG15 time spans. We then fit out a quadratic in the shortened ARN realizations from EPTA DR2 to account for the covariance between ARN and pulsar spindown which could not have been resolved using NG15. All chromatic noise spectra are referenced to a delay at 1400 MHz.

### 5.1. PSR J0613-0200

Applying *CustomGP* to PSR J0613-0200 results in a steeper ARN spectrum than using DMX (Figure 3). In the time-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** PSR J0613−0200: a variation in the ARN, uniquely observed using NG15, is mitigated using *CustomGP*. Noise properties are displayed by column: posterior parameter distributions (left), spectra (middle), and time-domain GP realizations (right). Noise types are displayed by row: ARN (top), DM variations (middle), and scattering-like chromatic noise with $\chi = 4$ (bottom). Data/model in use are displayed by color: NG15 using DMX (purple), NG15 using *CustomGP* (blue), and EPTA DR2 (yellow). All chromatic noise spectra are referenced to a delay at 1400 MHz. A variation present in the time-domain ARN realizations near MJD 56000 using DMX is now absorbed by the chromatic noise model using *CustomGP*. This ARN variation is not present in EPTA DR2. The ARN spectrum is steeper and more consistent with EPTA DR2 when applying *CustomGP*. An annual trend in the time series of DMX parameters is now absorbed by the scattering-like chromatic model using *CustomGP*. However, the scattering-like chromatic noise amplitude from NG15 is above the upper limit set by EPTA DR2.</figcaption></figure>

domain ARN realizations, this change corresponds to a variation just before MJD 56000 becoming flat when using *CustomGP*. When applying *CustomGP*, this variation is classified as chromatic, as variations of similar width appear in the DM and scattering-like chromatic variations, with opposite sign. It is plausible that this variation was falsely characterized as achromatic using DMX, as no such variation is present in EPTA DR2 ARN (Chalumeau et al. 2022; *IPTA_comp*), nor does it appear to be present in ARN from PPTA DR2 (Goncharov et al. 2021) or PPTA DR3 (Reardon et al. 2023b; *IPTA_comp*). Switching to *CustomGP* also reduces the Bayes factor for excess power in PSR J0613−0200’s ARN free spectrum just above a frequency of $f = 1 \text{ yr}^{-1}$ (Figure 11, Appendix C). *EPTA_Noise* notably does not favor the inclusion of any Fourier modes at or above $f = 1 \text{ yr}^{-1}$.

The DMX time series for PSR J0613−0200 shows evidence of annual sinusoidal variations, which could result from a steep DM gradient along the line of sight as the Earth orbits around the Sun (Keith et al. 2013; Jones et al. 2017). This annual DM trend manifests using both DMX and DMGP (Appendix B). However, when applying model *CustomGP*, the annual DM trend disappears and instead manifests in the scattering-like chromatic noise. This is evidenced by a peak at $f = 1 \text{ yr}^{-1}$ in the posterior scattering-like chromatic noise spectrum

(Figure 3). An annual scattering trend is supported by measurements of annual scintillation arc variability in this pulsar (Main et al. 2020, 2023; Liu et al. 2023). However, no annual chromatic noise in J0613−0200 appears to be present using EPTA DR2 (Chalumeau et al. 2022). Furthermore, the scattering-like chromatic noise parameters estimated using NG15 lie above the upper limit set by EPTA DR2, as evidenced by the lack of overlap in their scattering-like chromatic noise parameters and spectra. This discrepancy warrants further investigation, as the inclusion of $\chi = 4$ chromatic noise is entirely responsible for the changes to PSR J0613−0200’s ARN properties.

## 5.2. PSR J1012+5307

PSR J1012+5307 uniquely features a shallow ARN spectrum, present at a much higher amplitude than the GWB across the spectrum (NG15_Noise). Increasing the number of ARN Fourier modes $N_{freqs} \to 150$ in model *CustomGP* (Chalumeau et al. 2022) results in a slightly steeper, more constrained ARN spectrum, since power is now spread across more frequencies (Figure 4). The posteriors for $\gamma_{RN}$ are almost an exact match for NG15 and EPTA DR2. In the time-domain ARN realizations there are many features common to both data
<figure><img src="image_7.png" /><figcaption>**Figure 4.** PSR J1012+5307: each noise model fails to mitigate anticorrelations between chromatic noise and a shallow-spectrum ARN process in NG15. J1012+5307 features a peculiar shallow-spectrum ARN process, modeled in *CustomGP* using $N_{freqs} = 150$, following Chalumeau et al. (2022). Using *CustomGP*, NG15 and EPTA DR2 have the same ARN spectral index and share many common features in the time domain, including sharp spikes nears MJDs 56100 and 58350. However, the NG15 ARN spectrum features excess power near 50 nHz, which is not seen in EPTA DR2. ARN variations near MJDs 55000 and 56000 are anticorrelated with chromatic variations using NG15, and the scattering-like chromatic noise detection made by NG15 is also well above the upper limit set by EPTA DR2. See Figure 3 caption for panel descriptions.</figcaption></figure>

sets, including sharp and sudden delays near MJDs 56100 and 58350. However, $log_{10}A_{RN}$ is higher for NG15 than EPTA DR2. Furthermore, NG15_Noise showed that PSR J1012+5307 features excess noise at $f \sim 50$ nHz. Here this is produced as a bump in the NG15 ARN spectra using both DMX and *CustomGP*, but this bump is not present in the EPTA DR2 ARN spectrum.

Using both DMX and *CustomGP*, we notice simultaneous, anticorrelated variations in both the time-domain ARN and chromatic noise realizations, particularly near MJD 55000. Where DMX features simultaneous ARN and DM variations, *CustomGP* shows simultaneous ARN and $\chi=4$ scattering-like chromatic variations, with little support for power-law DM noise (Table 3). The anticorrelated nature of these chromatic and achromatic variations is highly unusual, and signifies a high chance of chromatic mismodeling (see Appendix D). To diagnose the issue using *CustomGP*, we highlight a sudden achromatic delay near MJD 56000, which corresponds to a sudden advance in the scattering-like chromatic noise at the same time. At $\nu = 800$ MHz, the median chromatic advance near MJD 56000 corresponds to $\Delta t_{chrom} \sim -3.36$ µs, which is nearly the same amplitude as the achromatic delay near MJD 56000 ($\Delta t_{ARN} \sim 3.71$ µs). As such, at the time of the event, both noise processes together ($\Delta t_{chrom} + \Delta t_{ARN}$) effectively cancel out in the lowest-frequency band. The achromatic delay remains in the higher-radio-frequency bands, while the inferred

chromatic delay decays down to $|\Delta t| < 360$ ns at $\nu > 1400$ MHz. Interestingly, no evidence for these anticorrelated variations is observed using EPTA DR2. Additionally, the scattering-like chromatic noise parameters are in major tension, as $log_{10}A_{chrom}$ from NG15 is over an order of magnitude larger than the upper limits set by EPTA DR2.

### 5.3. PSR J1600-3053

The choice of chromatic noise model has a noticeable effect on the inferred ARN parameters (Figure 2). DMX detects a shallow-spectrum ARN process, while the spectrum is much steeper using *CustomGP*. In the time domain, short-timescale ARN fluctuations with DMX are replaced by short-timescale scattering-like chromatic variations (with opposite sign) using *CustomGP*. This is especially clear near MJDs 55000, 56500, 57800, and 58000. The sign change may result from the DM model overcompensating for an unmodeled scattering delay at low radio frequencies (see Appendix D for an example of this). Interestingly, the ARN spectrum is not only steeper using *CustomGP*, but *amplified* at the lowest Fourier mode (Figure 5, top panel). This begins to raise the power-law ARN posteriors inferred using *CustomGP* above the upper limits from EPTA DR2.

Both NG15 and EPTA DR2 agree on the presence of $\chi=4$ scattering-like chromatic noise in PSR J1600-3053 with
<figure><img src="image_8.png" /><figcaption>**Figure 5.** PSR J1600−3053: ARN is highly sensitive to the chromatic noise parameters, which are inconsistent between NG15 and EPTA DR2. CustomGP mitigates short-timescale chromatic noise, resulting in a steeper ARN spectrum than DMX. NG15 and EPTA DR2 both agree on the presence of strong DM and scattering-like variations, but the scattering-like chromatic noise amplitude and DM noise spectral index are both lower using NG15. In the time domain, several chromatic features are present in both data sets, but neither data set agrees on whether the DM or scattering-like chromatic term should mitigate them. See Figure 3 caption for panel descriptions.</figcaption></figure>

$\log_{10} \beta^{\text{chrom}} > 2.2$. However, Figure 5 shows the chromatic noise parameters estimated by NG15 and EPTA DR2 using CustomGP are not consistent. EPTA DR2 favors a scattering-like chromatic noise spectrum with a similar spectral index but a higher amplitude than NG15. Meanwhile, NG15’s DM noise spectrum deviates from a pure power law, as it is similar to EPTA DR2 below 10 nHz but becomes more shallow past 10 nHz. In the time domain, several events (a spike near MJD 57500, a large bump between MJDs 55200 and 56500, and a dip near MJD 56700) are characterized as scattering-like chromatic events by EPTA DR2 but as DM events by NG15, or vice versa. These discrepancies suggest PSR J1600−3053 could benefit from a modified chromatic noise model, e.g., a model with a varied radio-frequency dependence.

## 5.4. PSR J1713+0747

Using CustomGP, PSR J1713+0747’s ARN parameters change significantly ($\log_{10} A_{\text{RN}}$ by >1σ; see Table 3) over the use of DMX and DMGP. We find that allowing $\chi_{E_2}$ to vary as a free parameter in CustomGP is directly responsible for this change. Figure 6 shows at the time of $E_2$ there is a noticeable dip in both the DMX time series and the time-domain ARN realizations, i.e., the ARN and DM variations are coupled at the time of $E_2$. Allowing $\chi_{E_2}$ to vary successfully decouples the event from the ARN and DM variations, as evidenced by the lack of excess noise in the time-domain ARN and DM

realizations at the time of $E_2$ using CustomGP. Furthermore, the power in the eighth frequency bin of PSR J1713+0747’s ARN free spectrum experiences a major drop in detection significance when using CustomGP (Appendix C, Figure 11). This may be relevant for GWB characterization, as the eighth frequency bin of the GWB free spectrum has been identified as a driver in pushing the fit for $\gamma_{\text{GWB}}$ to lower values (NG15_GWB; Agazie et al. 2024b).

Despite using the same models for NG15 and EPTA DR2 with CustomGP, there are several inconsistencies between the two data sets. For one, the recovered scattering-like chromatic noise parameters using NG15 are above the upper limit set by EPTA DR2. Furthermore, the recovered ARN amplitude is lower in NG15 than EPTA DR2. IPTA_comp reports a similar discrepancy between NG15 and a version of EPTA DR2 without legacy data at the 1.4σ level. The comparison of ARN realizations shows that the ARN is less consistent between the two data sets prior to MJD 57000, while it is more consistent after MJD 57000. Understanding the nature of these differences may be useful to improve PSR J1713+0747’s noise modeling in the future.

Figure 7 shows posterior distributions for the exponential dip model parameters, inferred using both NG15 and EPTA DR2. For $E_1$, EPTA DR2 favors a high amplitude, short decay timescale, and a chromatic index $\chi_{E_1} = 4.2_{-1.0}^{+1.1}$, whereas NG15 favors a lower amplitude, longer recovery timescale, and a
<figure><img src="image_9.png" /><figcaption>**Figure 6.** PSR J1713+0737: simultaneous ARN and DM variations near the second transient timing event are decoupled using *CustomGP*. The times of events $E_1$ and $E_2$ are indicated by the dashed lines. Using *DMX*, $E_2$ manifests as a dip in both the NG15 DMX time series and time-domain ARN realizations. Modeling $E_2$ in *CustomGP* using a free-chromatic index successfully decouples this event from both the ARN and DM noise. This lowers the value of $\log_{10}A_{RN}$ and widens the posterior on $\gamma_{RN}$. However, the ARN and scattering-like chromatic noise parameters in NG15 are in tension with EPTA DR2. See Figure 3 caption for panel descriptions.</figcaption></figure>

chromatic index $\chi_{E1} = 3.0_{-0.3}^{+0.3}$. The NG15 $E_1$ posteriors also have long tails and are covariant, requiring a larger amplitude and a smaller chromatic index as the decay timescale becomes smaller. The inconsistent characterization of $E_1$ between data sets is explained by the uneven properties of each data set near $t_{E1}$. EPTA DR2 features multiple TOAs from the NRT, WSRT, and EFF with an average cadence of ~3.5 days, but only at 1400 MHz and above. NG15 features lower-frequency (800 MHz) TOAs at MJD 54765 from the GBT, but these are not followed up with higher-frequency (1400, 2300 MHz) TOAs from the GBT and AO until MJD 54819.

The NG15 and EPTA DR2 posteriors for $E_2$ are more consistent with each other, only featuring differences in their variances (Figure 7, right panel). NG15 includes many subbanded TOAs and thus a high-radio-frequency resolution near $t_{E2}$, which explains why its posterior on $\chi_{E2}$ is more constrained than EPTA DR2. The chromatic index of $E_2$ is $\chi_{E2} = 1.1_{-0.5}^{+0.5}$ using EPTA DR2 and $\chi_{E2} = 1.37_{-0.09}^{+0.09}$ using NG15. These values are also consistent with the $\chi = 1.15_{-0.19}^{+0.18}$ estimated by Goncharov et al. (2021) using PPTA DR2, supporting the proposal therein that these events originate from the pulsar's magnetosphere. Furthermore, the observation that $0 < \chi_{E2} < 2$ critically explains why this event manifests as both excess ARN and DM noise using the DMX and DMGP models.

### 5.5. PSR J1744-1134

We begin by noticing some differences in the ARN properties of EPTA DR2 and NG15 when using DMX. Specifically, in the top-middle panel of Figure 8, the ARN spectrum using NG15 with DMX is at a higher amplitude than the spectrum using EPTA DR2 across the region of frequency space where the spectra overlap. In the time domain, there are peaks and troughs in NG15's ARN realizations near MJDs 54000 and 58000 that are not present using EPTA DR2. Note that PSR J1744-1134's ARN is dominated in EPTA DR2 by the lowest-frequency bin, which would manifest in the time domain as a single sinusoidal trend if using the EPTA DR2 data from before MJD 53000 (Chalumeau et al. 2022).

Using *CustomGP* reduces the ARN Bayes factor and the median amplitude of ARN over both the DMX and DMGP models (Table 3). In particular, the dip in the ARN realizations near MJD 58000 using DMX is now absorbed by the DM variations using *CustomGP*, with a corresponding bump in the scattering-like chromatic noise to counterbalance. As a result, the ARN properties are also now more consistent with EPTA DR2 (Figure 8). That said, the new ARN posteriors still overlap entirely with the ARN posteriors inferred under DMX and DMGP (as well as the GWB posteriors), and as such the ARN process detected using DMX is not ruled out.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Posterior distributions for PSR J1713+0747's exponential dip parameters inferred using model *CustomGP* with both *NG15* and *EPTA DR2*, where the dashed line indicates $\chi = 2$ as expected for DM variations. Both data sets provide consistent estimates on the second dip ($E_2$) model parameters but differ in characterization of the first dip ($E_1$).</figcaption></figure>

Information from additional frequency bands in future data sets will help us better decouple and understand these signals.

The $\chi = 4$ scattering-like chromatic noise detected in NG15 is also above the upper limits set by EPTA DR2. This discrepancy appears to be less significant in PSR J1744−1134 than the other pulsars, as the scattering-like chromatic noise spectra of NG15 and EPTA DR2 still overlap at the lowest Fourier modes (Figure 8). However, the time-domain scattering-like chromatic noise realizations allowed by EPTA DR2 do not line up with the NG15 scattering-like chromatic noise realizations.

## 5.6. PSR J1909−3744

Unlike the preceding pulsars, there is very little change to PSR J1909−3744's ARN parameters when switching from DMGP to *CustomGP*, aside from a very minute increase in $\gamma_{RN}$ (Table 3). There is also no change to PSR J1909−3744's free spectrum below $f = 1 \text{ yr}^{-1}$ (Appendix C, Figure 11). This indicates the inferred ARN signal in PSR J1909−3744 is robust to the choice of chromatic noise model. The ARN is also extremely consistent with EPTA DR2 (Figure 9).

Using *CustomGP*, the NG15 DM noise spectrum features an excursion from the power-law prior at $f = 1 \text{ yr}^{-1}$, indicating support for an annual DM process. We also detect $\chi = 4$ scattering-like chromatic variations in PSR J1909−3744, characterized by a nearly flat spectrum ($\gamma_{chrom} = 0.5_{-0.3}^{+0.3}$). This was also found by Srivastava et al. (2023) using InPTA DR1 (Tarafdar et al. 2022). While EPTA DR2 does not detect a significant scattering-like noise process, this introduces no tension with NG15, since the NG15 scattering-like chromatic noise posteriors are below the EPTA DR2 upper limits. The NG15 time-domain scattering-like chromatic noise realizations feature three sharp spikes that all align with annual DM cusps from the deterministic SW model. This suggests either that the scattering-like chromatic model is capturing

scattering variations from within the heliosphere, or that it is compensating for excess DM from a time-variable SW electron density.

# 6. Summary and Discussion

We compare three different chromatic noise models on a subsample of six pulsars from NG15. Since five out of these six pulsars' ARN processes are spectrally similar to the GWB signal (NG15_Noise), we pay special attention to the model dependence of the ARN. Out of these five, PSR J1713+0747 is the only pulsar whose ARN parameters change significantly (at a $>1\sigma$ level) using the *CustomGP* model over the DMX model. The change is directly linked to the modeling of its nondispersive chromatic timing event near MJD 57510 (Lam et al. 2018; Hazboun et al. 2020b). Since the GWB spectrum inferred in NG15_GWB is similar to the ARN in just PSR J1713+0747 using DMX, this change is likely to impact GWB spectral characterization.

Alongside PSR J1713+0747, several of these pulsars favor steeper ARN spectra under the *CustomGP* models. This result is unsurprising, since we detect excess nondispersive chromatic noise in these pulsars, and unmitigated chromatic noise is expected to result in low-spectral-index ARN (Cordes & Shannon 2010; NG15_Noise). However, both NG15_GWB and Reardon et al. (2023b) found preference for steeper GWB spectra using alternative pulsar noise models similar to the DMGP and *CustomGP* models used here. Since the GWB spectrum is used to inform astrophysical interpretations of the signal (Phinney 2001; NG15_SMBHBs; Agazie et al. 2024b), we plan to continue exploring the effects of chromatic noise models on the full PTA in an upcoming work. Furthermore, it will be worth investigating if mitigation of excess chromatic noise in the data set could improve the measurements of HD cross-correlations between pulsar pairs, as suggested by Di Marco et al. (2024).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** PSR J1744−1134: use of *CustomGP* reduces the discrepancy in ARN characterization between NG15 and EPTA DR2. When applying *CustomGP* to NG15, the detection of ARN is reduced to an upper limit. This improves the consistency in ARN characterization between data sets, as indicated by their ARN spectra. However, the detection of scattering-like chromatic noise in NG15 is above the upper limit set by EPTA DR2. See Figure 3 caption for panel descriptions.</figcaption></figure>

It is also interesting that the changes to the ARN properties under the *CustomGP* model are much smaller in some pulsars (e.g., PSR J1909−3744) than in others. As such, it is possible that applying the *CustomGP* model on some pulsars could have little to no impact on GWB analyses, even if *CustomGP* provides a better fit to the pulsar’s TOAs. In a future work, the impact of the choice of noise model for individual pulsars on GWB characterization could be assessed, e.g., from factorized likelihood methods (Taylor et al. 2022; Lamb et al. 2023). This could then be used to expand upon existing model selection methods (e.g., Lentati et al. 2016; Hazboun et al. 2020a; Goncharov et al. 2021) by only using custom noise models for pulsars that measurably impact GWB spectral characterization (or other results of interest). Since the *CustomGP* model is more computationally expensive than the DMX model, this could help reduce the computational burden of cross-correlation analyses of the GWB (which require simultaneous analysis of all pulsars) with custom noise models.

Custom pulsar noise models are also useful tools to study ISM processes. For instance, the distribution of $\gamma_{DM}$ values can be used to assess if the ISM is consistent with the model of Kolmogorov turbulence, which predicts $\gamma_{DM} = 8/3$ (Keith et al. 2013; Lentati et al. 2016; Goncharov et al. 2021). Here we find $\gamma_{DM} < 8/3$ when using DMGP for these six pulsars in NG15. However, for PSRs J0613−0200 and J1600−3053, consistency with $\gamma_{DM} = 8/3$ improves using *CustomGP* (Table 3). Separating DM variations into a stochastic GP and

a SW component is also very useful for estimating the SW electron density, $n_{\text{Earth}}$ (Appendix A; Madison et al. 2019; Hazboun et al. 2022; Kumar et al. 2022; Reardon et al. 2023b). Since the model we use here assumes $n_{\text{Earth}}$ is constant, a time-variable SW will still induce excess noise in the DM noise spectrum at harmonics of $f = 1 \text{ yr}^{-1}$ (Hazboun et al. 2022). We observe this here for PSR J1909−3744, and expect this effect to be much more pronounced for pulsars close to the ecliptic. Introducing an additional GP to vary $n_{\text{Earth}}$ over time is a promising method to mitigate this effect (Hazboun et al. 2022; Nitu et al. 2024).

To validate our results, we compare the inferred ARN properties from both NG15 and EPTA DR2 using the *CustomGP* model, which is based on the favored models from *EPTA_Noise*. This analysis complements *IPTA_comp*, who have also compared pulsar noise properties from recent PTA data sets using the same noise models. For most pulsars, we find the inferred ARN properties are consistent using *CustomGP*. Furthermore, applying *CustomGP* instead of DMX to three pulsars in NG15 (PSRs J0613−0200, J1012+5307, and J1744−1134) alleviated some discrepancies in ARN characterization between the two data sets. This strengthens our confidence that the *CustomGP* model is improving ARN estimation for these pulsars. However, PSR J1713+0747 is an exception: The ARN detected using NG15 is lower in amplitude than the ARN detected using EPTA DR2. This discrepancy was already found by *IPTA_comp*, but here
<figure><img src="image_12.png" /><figcaption>**Figure 9.** PSR J1909−3744: ARN is robust to the choice of chromatic noise model. PSR J1909−3744’s ARN in NG15 experiences little to no change going from DMX to CustomGP, indicating any errors in chromatic noise mitigation are small and decoupled from ARN. ARN characterization is also remarkably consistent between NG15 and EPTA DR2. Spikes in the NG15 time series for scattering-like chromatic noise are also observed near annual DM cusps, indicating they likely originate from the heliosphere. See Figure 3 caption for panel descriptions.</figcaption></figure>

we confirm its existence using the *CustomGP* model. *EPTA_Noise* also recover similarly high ARN amplitudes for PSR J1713+0747 using both alternative versions of *EPTA DR2*. Analysis of a future combined data set (IPTA Data Release 3, DR3) may help resolve this inconsistency and assess any potential impacts on GWB inferences. It would also be useful to extend the comparison to other data sets, e.g., by applying the PPTA DR3 noise models (Reardon et al. 2023b) to the NG15 data and comparing the results.

This comparison of data sets reveals another major discrepancy in the scattering-like chromatic noise: Four out of six pulsars’ scattering-like noise amplitudes estimated using NG15 are above the upper limits set by EPTA DR2. Throughout this work, we have used $\Delta t \propto \nu^{-\chi}$ with $\chi = 4$ to describe scattering-like chromatic noise. However, an index $\chi < 4$ may resolve the above discrepancy. We suggest this because the scattering-like chromatic noise detected using NG15 is only above the upper limit set by EPTA DR2 for pulsars where EPTA DR2 contains low-radio-frequency ($\sim$300 MHz) TOAs from the WSRT (Figure 1). If these low-frequency TOAs are responsible for ruling out a $\chi = 4$ scattering-like noise process, a lower chromatic index would likely reduce the delay at the lowest radio frequencies and therefore raise the upper limits in EPTA DR2. Supporting evidence for scattering scaling indices $\chi < 4$ has been found in several past studies (Lewandowski et al. 2015; Levin et al. 2016; Turner et al. 2021), and the scaling indices may also vary

over time (Bansal et al. 2019; Liu et al. 2022). Other chromatic processes, such as low-level profile variations, also need not have $\chi = 4$ dependence or even a power-law dependence on radio frequency. If the scattering-like chromatic noise GP is absorbing some additional chromatic process, this is also a viable explanation for the discrepancy. Ultimately, the $\chi = 4$ scattering-like chromatic noise GP is still a valuable phenomenological component of the noise model, as it is still favored with a large Bayes factor for all six pulsars.

In particular, we suggest that cases of chromatic noise processes displaying time correlations with one another, or “competing” with each other, may highlight cases of chromatic mismodeling. Assuming a two-radio-frequency measurement of the TOA, modeled as an ARN and DM process, an excess achromatic delay will be introduced with the opposite sign of any unmodeled scattering delay (Appendix D). This occurs since the unmodeled scattering delay is primarily absorbed by the DM model, rather than the ARN. We suggest this is potentially taking place for PSRs J0613−0200 and J1600−3053 in NG15, as based on our comparisons of the ARN and chromatic noise realizations using the DMX and *CustomGP* models. We similarly find that incomplete modeling of the second chromatic event ($E_2$) in PSR J1713+0747 (which is found to have a chromatic index of $\chi = 1.37_{-0.09}^{+0.09}$) results in excess ARN and DM noise at the exact time of the event. Furthermore, we find using both DMX and *CustomGP* that the ARN and scattering-like chromatic noise in PSR J1012+5307
displays anticorrelations over time. The behavior is consistent with the presence of a $\chi < 0$ chromatic process as explored at the end of Appendix D. It is unknown to us what type of physical process this may correspond to. However, this idea is supported by EPTA_Noise, who found $\chi = -0.65_{-0.41}^{+0.46}$ after a free-chromatic analysis of PSR J1012+5307 from the joint EPTA DR2 and InPTA DR1 data set (Tarafdar et al. 2022; EPTA DR2).

The GP models used here could benefit from further advances, several of which we will test on more pulsars from NG15 in an upcoming work. As discussed, it will be important to assess the evidence of time-variable scattering at chromatic indices other than $\chi = 4$. This could be assessed using a GP with $\chi$ as a free variable (e.g., Goncharov et al. 2021; Srivastava et al. 2023; EPTA_Noise). Furthermore, for pulsars close to the ecliptic, it will be important to implement more sophisticated GP models for time variability of the SW density (e.g., Hazboun et al. 2022; Nitu et al. 2024). GP models could also help mitigate other effects such as frequency-dependent jitter (Lam et al. 2019; Kulkarni et al. 2024), frequency-dependent DM variations due to multipath propagation (Cordes et al. 2016), or nonstationary noise (Ellis & Cornish 2016). Performing WN model selection (e.g., Srivastava et al. 2023; Miles et al. 2023) or including additional ARN processes at high fluctuation frequencies (Reardon et al. 2023b) may also reduce the noise floor at higher GW frequencies. Additionally, PSR J1713+0747 displayed a dramatic pulse profile change in early 2021 (Singha et al. 2021), which introduced timing delays that scale nonmonotonically with radio frequency (Jennings et al. 2024). Mitigating this event at the level of timing residuals would require a more sophisticated chromatic model than those previously used for PSR J1713+0747.

We highlight that the six pulsars we investigate here are only a small subsample of the full NG15 data set, and are not necessarily representative of the whole. Many pulsars (e.g., PSRs B1937+21 and J1903+0327) have very different noise properties due either to factors intrinsic to the pulsar or to its location in the ISM (NG15_Noise). Pulsars closer to the ecliptic will be impacted more strongly by the choice of SW model (Tiburzi et al. 2021). Many pulsars in NG15 also have shorter time spans than the six we study here. Additionally, all six pulsars here have been observed by the GBT, while many pulsars in NG15 have been primarily observed by the AO. As such, investigating more pulsars should reveal new discoveries about the implications of chromatic noise modeling choices.

Finally, upcoming data-set improvements are very promising for chromatic noise mitigation prospects. The GBT has been upgraded with the VERSatile GBT Astronomical Spectrometer (Bussa & VEGAS Development Team 2012) and an ultra-wide-bandwidth receiver capable of observations up to 3.8 GHz. Work is underway to install a cyclic spectroscopy backend at the GBT, with the goal to instantly remove scattering effects before any further timing analysis has taken place (Dolch et al. 2021; Turner et al. 2023). Even if scattering cannot be removed in some pulsars, the high-frequency resolution enabled by cyclic spectroscopy can allow more accurate regular measurement of a pulsar's scintillation bandwidth, and thus the scattering delay (Dolch et al. 2021). CHIME/Pulsar will also provide observations in a 400–800 MHz bandwidth with daily cadence (CHIME/Pulsar Collaboration et al. 2021). Taken together, these developments

should allow for highly precise modeling of DM and other chromatic processes in future NANOGrav data sets. Finally, a future IPTA data set (IPTA DR3) will combine data from all PTAs together to maximize the data cadence, timing baselines, sky coverage, and effective radio-frequency coverage achievable using current data sets. This will allow excellent mitigation of chromatic noise processes and further improve PTA sensitivity to GWs.

Acknowledgments

*Author contributions.* This paper is the result of the work of many people and uses data from over a decade of pulsar timing observations. B.L. performed the analysis, created all figures and tables, and wrote a majority of the text. C.M.F.M. conceived of the project, supervised the analysis, helped develop the manuscript, and provided advice on figures and interpretation. C.M.F.M., J.S.H., and D.C.G. provided continuous advice and feedback on results throughout the project. A.C. assisted with the analysis of EPTA DR2 data and interpretation of the noise analysis results. J.S.H. and J.S. provided additional expertise on pulsar noise modeling as well as extensive resources and codes which were adapted to perform the analysis. 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.D., D.C.G., P.A.G., J.G., 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.N., B.B.P.P., N.S.P., H.A.R., S.M.R., P.S.R., A.S., C.S., B.J.S.A., I.H.S., K.S., A.S., J.K.S., and H.M.W. all ran observations and developed timing models for the NANOGrav 15 yr data set. A.C., D.J.C., I.C., L.G., H.H., M.J.K., K.L., J.W.M., A.Pa., D.P., A.Po., G.M.S., and G.T. contributed observations, timing, and noise models for EPTA DR2.

We thank the anonymous reviewer for several useful comments which helped us improve the quality of this paper. We additionally thank Jacob Turner, Sarah Vigeland, Jeremy Baier, Priyamvada Natarajan, Joris Verbiest, James Andrew Casey-Clyde, London Willson, Jonathan Nay, Małgorzata Curyło, Deven Bhakta, and NANOGrav's Detection and Noise Budget Working Groups for useful input. B.L. was supported in part by NASA CT Space Grant PTE Federal Award No. 80NSSC20M0129. C.M.F.M. was supported in part by the National Science Foundation under grant Nos. NSF PHY-1748958 and AST-2106552. A.C. acknowledges support from the Paris Île-de-France Region and financial support provided under the European Union's H2020 ERC Consolidator Grant "Binary Massive black hole Astrophysics" (B. Massive, grant agreement No. 818691). 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. P.R.B. is supported by the Science and Technology Facilities Council, grant No. ST/W000946/1. Support for H.T.C. is provided by NASA through the NASA Hubble Fellowship Program grant No. HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. Pulsar research at UBC is supported by an NSERC Discovery Grant and by CIFAR. K.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 No. 80GSFC21M0002. G.E.F. is supported by NSF award PHY-2011772. D.R.L. and M.A.M. are supported by NSF No. 1458952. M.A.M. is supported by NSF No. 2009425. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. 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. N.S.P. was supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. 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. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. I.C., L.G., and G.T. acknowledge financial support from “Programme National de Cosmologie and Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA, and CNES, France. I.C., L.G., and G.T. acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. Pulsar research at Jodrell Bank Centre for Astrophysics is supported by an STFC Consolidated Grant (grant No. ST/T000414/1).

Software: ENTERPRISE (Ellis et al. 2020), enterprise_extensions (Taylor et al. 2021), PTMCMC (Ellis & van Haasteren 2017), la_forge (Hazboun 2020), PINT (Luo et al. 2021), numpy (Harris et al. 2020), matplotlib (Hunter 2007), corner (Foreman-Mackey 2016), seaborn (Waskom 2021).

## Appendix A
### Single-pulsar Solar Wind Electron Density Estimates

We estimate the SW electron density at 1 au ($n_{Earth}$) independently for each pulsar using our time-independent, $1/r^2$ density profile SW model. PTA estimates of $n_{Earth}$ have been performed more comprehensively elsewhere (e.g., Madison et al. 2019; Hazboun et al. 2022; Reardon et al. 2023b), but we report our own estimates of $n_{Earth}$ to detail the similarities and differences among the models and data sets considered here. Table 4 compares $n_{Earth}$ medians and 68.3% credible intervals (1σ regions) as well as log-scaled SW Bayes factors ($log_{10} B^{SW}$) estimated using NG15 (under models DMGP and CustomGP) and EPTA DR2 for each pulsar. For reference, we also show the ecliptic latitude (ELAT) value for each pulsar, since the SW introduces a larger DM correction for pulsars closer to the ecliptic. PSRs J0613-0200 and J1012+5307 are relatively far from the ecliptic, and low values of $log_{10} B^{SW}$ indicate the SW is not well detected using these pulsars. Meanwhile, PSR J1713+0747 is more precisely timed and has more TOAs, so it can still constrain the SW despite having ELAT = 30.7°.

Switching from DMGP to CustomGP results in slight changes to the estimated $n_{Earth}$ parameters. These changes are most significant for PSRs J1600-3053 and J1713+0747, as the median value of $n_{Earth}$ estimated under CustomGP lies just outside the 68.3% credible interval under DMGP. The $n_{Earth}$ estimates are fairly consistent using NG15 and EPTA DR2, although for PSRs J1600-3053, J1713+0747, J1744-1134, and J1909-3744, $n_{Earth}$ estimates are slightly higher using NG15. Estimates of $n_{Earth}$ from NG15 and EPTA DR2 are the most different for PSR J1600-3053. However, this difference in PSR J1600-3053 could be related to the different chromatic noise properties estimated using NG15 and EPTA DR2 (Section 5.3).

<caption>Table 4</caption>

Comparisons of Solar Wind Electron Density Parameters and Bayes Factors from Each Pulsar and Model/Data-set Combination: NG15 with DMGP, NG15 with CustomGP, and EPTA DR2

<table><thead><tr><th rowspan="2">PSR</th><th rowspan="2">ELAT<br>(deg)</th><th colspan="2">NG15 DMGP</th><th colspan="2">NG15 CustomGP</th><th colspan="2">EPTA DR2</th></tr><tr><th>n<sub>Earth</sub><br>(cm<sup>-3</sup>)</th><th>log<sub>10</sub> B<sup>SW</sup></th><th>n<sub>Earth</sub><br>(cm<sup>-3</sup>)</th><th>log<sub>10</sub> B<sup>SW</sup></th><th>n<sub>Earth</sub><br>(cm<sup>-3</sup>)</th><th>log<sub>10</sub> B<sup>SW</sup></th></tr></thead><tbody><tr><td>J0613-0200</td><td>-25.4</td><td>2.7<sup>+1.3</sup><sub>-1.3</sub></td><td>-0.2</td><td>1.8<sup>+1.3</sup><sub>-1.1</sub></td><td>-0.8</td><td>3.3<sup>+2.6</sup><sub>-2.1</sub></td><td>-0.5</td></tr><tr><td>J1012+5307</td><td>38.8</td><td>8.6<sup>+4.9</sup><sub>-4.5</sub></td><td>0.2</td><td>4.7<sup>+4.3</sup><sub>-3.2</sub></td><td>-0.4</td><td>6.4<sup>+1.9</sup><sub>-1.9</sub></td><td>1.3</td></tr><tr><td>J1600-3053</td><td>-10.1</td><td>6.2<sup>+0.8</sup><sub>-0.8</sub></td><td>&gt;3</td><td>5.2<sup>+0.9</sup><sub>-0.9</sub></td><td>&gt;3</td><td>2.9<sup>+0.7</sup><sub>-0.7</sub></td><td>2.6</td></tr><tr><td>J1713+0747</td><td>30.7</td><td>4.9<sup>+0.6</sup><sub>-0.6</sub></td><td>&gt;3</td><td>5.9<sup>+0.8</sup><sub>-0.8</sub></td><td>&gt;3</td><td>4.0<sup>+1.0</sup><sub>-1.1</sub></td><td>1.7</td></tr><tr><td>J1744-1134</td><td>11.8</td><td>4.1<sup>+0.5</sup><sub>-0.5</sub></td><td>&gt;3</td><td>4.2<sup>+0.5</sup><sub>-0.5</sub></td><td>&gt;3</td><td>3.2<sup>+0.7</sup><sub>-0.7</sub></td><td>2.5</td></tr><tr><td>J1909-3744</td><td>-15.2</td><td>3.7<sup>+0.4</sup><sub>-0.4</sub></td><td>&gt;3</td><td>3.8<sup>+0.4</sup><sub>-0.4</sub></td><td>&gt;3</td><td>3.3<sup>+0.4</sup><sub>-0.4</sub></td><td>&gt;3</td></tr></tbody></table>

Notes. Parameters are presented using the median and 68.3% Bayesian credible intervals (referenced here as 1σ regions), and Bayes factors are calculated from our posterior distributions using the Savage-Dickey approximation.
## Appendix B
### Comparison of Dispersion Measure Time Series

In Figure 10, we compare DM(t) values as estimated under DMX and DMGP for the six pulsars in NG15. The top panels show the estimated deviations to each pulsar's fiducial DM value over time, in units of $10^{-3}$ pc cm$^{-3}$. These are represented by the time series of DMX parameters using the DMX model, and 100 GP realizations of DM variations using the DMGP model. Qualitatively, both DM models produce similar trends for all pulsars, especially once transient events in PSR J1713+0747 are accounted for using DMGP (Lam et al. 2018; Wang et al. 2019; Hazboun et al. 2020b).

The bottom panels of Figure 10 show the difference between the medians of each estimated DM value over time, ΔDM.

These are calculated by subtracting the medians of the time-domain DM realizations from the DMX time series. We report the rms ΔDM for each pulsar. We find ΔDM RMS > $10^{-4}$ pc cm$^{-3}$ for PSRs J1012+5307 and J1600–3053, while ΔDM RMS < $10^{-4}$ pc cm$^{-3}$ for the remaining pulsars. DM estimation errors are known to induce deviations to TOAs at infinite radio frequency (Cordes & Shannon 2010; Lam et al. 2015), therefore the larger ΔDM rms values in PSRs J1012+5307 and J1600–3053 help explain why only only these two pulsars have significantly different ARN parameters using DMX and DMGP (Table 3, Figure 2). There are also some subtle time correlations in each pulsar's ΔDM, which may be useful to study in future analyses.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of DM variations (ΔDM) recovery using the two models, DMX and DMGP, for each pulsar. Both models recover qualitatively similar trends in DM for all pulsars, however the difference in estimated DM values is largest for PSRs J1012+5307 and J1600–3053. Top panels: time series of DMX parameters (black) superimposed with 100 DM GP realizations (turquoise). Bottom panels: difference in estimated DM over time and rms difference between ΔDM values.</figcaption></figure>
## Appendix C
### Achromatic Red Noise Free Spectra

We generate Bayesian power spectra for each pulsar using a free-spectral PSD, which allows the power in each frequency bin to vary as a free parameter (Lentati et al. 2013). In Figure 11, we compare the power spectra for all six pulsars using DMX and CustomGP. The top panels show the posteriors for the log-scaled rms timing residual power $log_{10} \rho_i$ at each frequency $f_i$, alongside a power law using the maximum-likelihood values of $log_{10} A_{RN}$ and $\gamma_{RN}$ from each model. The bottom panels show the log-scaled Bayes factors $log_{10} B_i$ for the presence of excess power in each frequency bin, measured using the Savage-Dickey density ratio (Dickey 1971). These may be interpreted as a measure of the consistency of each $\rho_i$

with zero excess power. In cases where there are a lack of samples consistent with $log_{10} \rho_i = -10$, we set a lower limit of $log_{10} B_i > 1.7$.

The choice of chromatic noise model influences the free spectra of all six pulsars. These changes are most interesting for PSRs J1713+0747 and J1909-3744 since they each have the lowest residual rms power overall. For PSR J1713+0747, the changes to $log_{10} B_i$ indicate the change to $\rho_8$ (the eighth frequency bin) is most consequential for the changes to this pulsar's ARN properties. For PSR J1909-3744, the free spectra inferred using both models are nearly identical at low frequencies below $f = 1/\text{yr}$; above $f = 1/\text{yr}$, CustomGP reduces the rms power at a few frequencies.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Power spectra and Bayes factors using a free-spectral PSD for each pulsar under the two models, DMX and CustomGP. Top panels: the posteriors for the log-scaled residual power $log_{10} \rho_i$ at each frequency, in units of excess timing delay. Bottom panels: log-scaled Bayes factors $log_{10} B_i$, indicating the statistical significance of excess power in each frequency bin. For each pulsar, changes to $\rho_i$ are isolated to a few key frequencies, but still influence the inferred power-law noise parameters (Figure 2).</figcaption></figure>
## Appendix D

### Systematic Errors from Unmodeled Chromatic Effects

We demonstrate how unmodeled chromatic effects may bias estimates of ARN and DM variations in a simplified analytic case (see also Lentati et al. 2016; Lam et al. 2020; Sosa Fiscella et al. 2024). For simplicity, we will assume a TOA is measured only at two radio frequencies, $\nu_0$ and $\nu_1$, where $\nu_1 > \nu_0$. Let us define the true frequency-dependent timing delay $\Delta t(\nu)$ from ARN, DM, and scattering as

$$ \Delta t(\nu) = \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}^{-2} + \Delta t_{chrom} \tilde{\nu}^{-\chi}, $$

where $\tilde{\nu}$ is a dimensionless frequency scaled to some reference frequency $\nu_{ref}$, $\Delta t_{ARN}$ is the delay from achromatic processes (i.e., spin noise or GWs), $\Delta t_{DM}$ is the delay due to DM at the reference frequency $\nu_{ref}$, $\Delta t_{chrom}$ is the delay due to scattering at the reference frequency $\nu_{ref}$, and $\chi$ is the scattering scaling index. We may assume that $\chi = 4$, however this calculation works for any $\chi > 2$. We will assume our two frequencies $\nu_0$ and $\nu_1$ are widely separated, such that $\Delta t_{ARN}$ is small (in comparison to chromatic errors) at $\nu_0$,

$$ \begin{aligned} \Delta t(\nu_0) &= \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}_0^{-2} + \Delta t_{chrom} \tilde{\nu}_0^{-\chi} \\ &\cong \Delta t_{DM} \tilde{\nu}_0^{-2} + \Delta t_{chrom} \tilde{\nu}_0^{-\chi}, \end{aligned} $$

and $\Delta t_{chrom}$ is small at $\nu_1$,

$$ \begin{aligned} \Delta t(\nu_1) &= \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}_1^{-2} + \Delta t_{chrom} \tilde{\nu}_1^{-\chi} \\ &\cong \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}_1^{-2}. \end{aligned} $$

Now we estimate the time delay by modeling it as the sum of only a DM and an achromatic process. We will define the total estimated delay as

$$ \delta t(\nu) = \delta t_{ARN} + \delta t_{DM} \tilde{\nu}^{-2}, $$

where $\delta t_{ARN}$ and $\delta t_{DM}$ are the estimated ARN and DM delays. Again, we will write this down at our two frequencies $\nu_0$ and $\nu_1$, assuming $\delta t_{ARN}$ is small at $\nu_0$,

$$ \begin{aligned} \delta t(\nu_0) &= \delta t_{ARN} + \delta t_{DM} \tilde{\nu}_0^{-2} \\ &\cong \delta t_{DM} \tilde{\nu}_0^{-2}, \end{aligned} $$

while at $\nu_1$ we have exactly

$$ \delta t(\nu_1) = \delta t_{ARN} + \delta t_{DM} \tilde{\nu}_1^{-2}. $$

If we have only measured the TOA at $\nu_0$ and $\nu_1$, then our two model parameters $\delta t_{DM}$ and $\delta t_{RN}$ can perfectly fit our data such that $\delta t(\nu_0) = \Delta t(\nu_0)$ and $\delta t(\nu_1) = \Delta t(\nu_1)$. We can then use the measurement at $\nu_0$ to determine how the estimated DM delay $\delta t_{DM}$ relates to the true delays:

$$ \begin{aligned} \delta t(\nu_0) &= \Delta t(\nu_0), \\ \delta t_{DM} \tilde{\nu}_0^{-2} &= \Delta t_{DM} \tilde{\nu}_0^{-2} + \Delta t_{chrom} \tilde{\nu}_0^{-\chi}, \\ \delta t_{DM} &= \Delta t_{DM} + \Delta t_{chrom} \tilde{\nu}_0^{-(\chi-2)}. \end{aligned} $$

Meanwhile, we can use the measurement at $\nu_1$ to determine how the estimated achromatic delay relates to the true delays:

$$ \begin{aligned} \delta t(\nu_1) &= \Delta t(\nu_1), \\ \delta t_{ARN} + \delta t_{DM} \tilde{\nu}_1^{-2} &= \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}_1^{-2}, \\ \delta t_{ARN} &= \Delta t_{ARN} + \Delta t_{DM} \tilde{\nu}_1^{-2} \\ &\quad - (\Delta t_{DM} + \Delta t_{chrom} \tilde{\nu}_0^{-(\chi-2)}) \tilde{\nu}_1^{-2} \\ &= \Delta t_{ARN} - \Delta t_{chrom} \tilde{\nu}_0^{-(\chi-2)} \tilde{\nu}_1^{-2}. \end{aligned} $$

This shows that if there is an unmodeled scattering delay with $\chi > 2$, the DM will be shifted by $\Delta t_{chrom} \tilde{\nu}_0^{-(\chi-2)}$, while the achromatic delay will be shifted by $-\Delta t_{chrom} \tilde{\nu}_0^{-(\chi-2)} \tilde{\nu}_1^{-2}$. In other words, unmodeled scattering variations may manifest as an excess ARN process proportional to the true scattering-induced delay with an opposite sign. Additionally, the excess achromatic delay will tend to zero as $\tilde{\nu}_1$ becomes very large. While this result is based on a simplified model of the TOA with measurements at only two radio frequencies, it is possible this effect may still arise in real data. A more rigorous quantification will be left for future work.

One can also perform a similar calculation in the case of an unmodeled chromatic process with $0 < \chi < 2$. The result is that the DM and achromatic delays will each be overestimated as

$$ \begin{aligned} \delta t_{DM} &= \Delta t_{DM} + \Delta t_{chrom} \tilde{\nu}_0^{2-\chi}, \\ \delta t_{ARN} &= \Delta t_{ARN} + \Delta t_{chrom} \tilde{\nu}_1^{-2} (\tilde{\nu}_1^{2-\chi} - \tilde{\nu}_0^{2-\chi}). \end{aligned} $$

The difference $\tilde{\nu}_1^{2-\chi} - \tilde{\nu}_0^{2-\chi}$ will always be positive for $\chi < 2$. For a hypothetical chromatic process with $\chi < 0$:

$$ \begin{aligned} \delta t_{DM} &= \Delta t_{DM} - \Delta t_{chrom} \tilde{\nu}_1^{-\chi} \tilde{\nu}_0^2, \\ \delta t_{ARN} &= \Delta t_{ARN} + \Delta t_{chrom} \tilde{\nu}_1^{-\chi}. \end{aligned} $$

This shows that an unmodeled chromatic process with $0 < \chi < 2$ may manifest as excess DM and achromatic delays with the same sign, whereas if $\chi < 0$ then the excess DM and achromatic delays will again have opposite sign.

## ORCID iDs

Bjorn Larsen https://orcid.org/0000-0001-6436-8216
Chiara M. F. Mingarelli https://orcid.org/0000-0002-4307-1322
Jeffrey S. Hazboun https://orcid.org/0000-0003-2742-3321
Aurélien Chalumeau https://orcid.org/0000-0003-2111-1001
Deborah C. Good https://orcid.org/0000-0003-1884-348X
Joseph Simon https://orcid.org/0000-0003-1407-6607
Gabriella Agazie https://orcid.org/0000-0001-5134-3925
Akash Anumarlapudi https://orcid.org/0000-0002-8935-9882
Anne M. Archibald https://orcid.org/0000-0003-0638-3340
Paul T. Baker https://orcid.org/0000-0003-2745-753X
Paul R. Brook https://orcid.org/0000-0003-3053-6538
H. Thankful Cromartie https://orcid.org/0000-0002-6039-692X
Kathryn Crowter https://orcid.org/0000-0002-1529-5169
Megan E. DeCesar https://orcid.org/0000-0002-2185-1790
Paul B. Demorest https://orcid.org/0000-0002-6664-965X
Timothy Dolch https://orcid.org/0000-0001-8885-6388
Elizabeth C. Ferrara https://orcid.org/0000-0001-7828-7708
William Fiore https://orcid.org/0000-0001-5645-5336

Emmanuel Fonseca https://orcid.org/0000-0001-8384-5049

Gabriel E. Freedman https://orcid.org/0000-0001-7624-4616

Nate Garver-Daniels https://orcid.org/0000-0001-6166-9646

Peter A. Gentile https://orcid.org/0000-0001-8158-683X

Joseph Glaser https://orcid.org/0000-0003-4090-9780

Ross J. Jennings https://orcid.org/0000-0003-1082-2342

Megan L. Jones https://orcid.org/0000-0001-6607-3710

David L. Kaplan https://orcid.org/0000-0001-6295-2881

Matthew Kerr https://orcid.org/0000-0002-0893-4073

Michael T. Lam https://orcid.org/0000-0003-0721-651X

Duncan R. Lorimer https://orcid.org/0000-0003-1301-966X

Jing Luo https://orcid.org/0000-0001-5373-5914

Ryan S. Lynch https://orcid.org/0000-0001-5229-7430

Alexander McEwen https://orcid.org/0000-0001-5481-7559

Maura A. McLaughlin https://orcid.org/0000-0001-7697-7422

Natasha McMann https://orcid.org/0000-0002-4642-1260

Bradley W. Meyers https://orcid.org/0000-0001-8845-1225

Cherry Ng https://orcid.org/0000-0002-3616-5160

David J. Nice https://orcid.org/0000-0002-6709-2566

Timothy T. Pennucci https://orcid.org/0000-0001-5465-2889

Benetge B. P. Perera https://orcid.org/0000-0002-8509-5947

Nihan S. Pol https://orcid.org/0000-0002-8826-1285

Henri A. Radovan https://orcid.org/0000-0002-2074-4360

Scott M. Ransom https://orcid.org/0000-0001-5799-9714

Paul S. Ray https://orcid.org/0000-0002-5297-5278

Ann Schmiedekamp https://orcid.org/0000-0003-4391-936X

Carl Schmiedekamp https://orcid.org/0000-0002-1283-2184

Brent J. Shapiro-Albert https://orcid.org/0000-0002-7283-1124

Ingrid H. Stairs https://orcid.org/0000-0001-9784-8670

Kevin Stovall https://orcid.org/0000-0002-7261-594X

Abhimanyu Susobhanan https://orcid.org/0000-0002-2820-0931

Joseph K. Swiggum https://orcid.org/0000-0002-1075-3837

Haley M. Wahl https://orcid.org/0000-0001-9678-0299

David J. Champion https://orcid.org/0000-0003-1361-7723

Ismaël Cognard https://orcid.org/0000-0002-1775-9692

Lucas Guillemot https://orcid.org/0000-0002-9049-8716

Huanchen Hu https://orcid.org/0000-0002-3407-8071

Michael J. Keith https://orcid.org/0000-0001-5567-5492

Kuo Liu https://orcid.org/0000-0002-2953-7376

James W. McKee https://orcid.org/0000-0002-2885-8485

Aditya Parthasarathy https://orcid.org/0000-0002-4140-5616

Delphine Perrodin https://orcid.org/0000-0002-1806-2483

Andrea Possenti https://orcid.org/0000-0001-5902-3731

Golam M. Shaifullah https://orcid.org/0000-0002-8452-4834

Gilles Theureau https://orcid.org/0000-0002-3649-276X

## References

Afzal, A., Agazie, G., Anumarlapudi, A., et al. 2023, ApJL, 951, L11

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJL, 951, L8

Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023c, ApJL, 951, L9

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023b, ApJL, 952, L37

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023d, ApJL, 951, L10

Agazie, G., Antoniadis, J., Anumarlapudi, A., et al. 2024a, ApJ, 966, 105

Agazie, G., Baker, P. T., Bécsy, B., et al. 2024b, arXiv:2404.07020

Aigrain, S., & Foreman-Mackey, D. 2023, ARA&A, 61, 329

Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021, ApJS, 252, 5

EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023, A&A, 678, A49

Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJL, 905, L34

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37

Bansal, K., Taylor, G. B., Stovall, K., & Dowell, J. 2019, ApJ, 875, 146

Bassa, C. G., Janssen, G. H., Karuppusamy, R., et al. 2016, MNRAS, 456, 2196

Bécsy, B., Cornish, N. J., Meyers, P. M., et al. 2023, ApJ, 959, 9

Brook, P. R., Karastergiou, A., McLaughlin, M. A., et al. 2018, ApJ, 868, 122

Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019, A&ARV, 27, 5

Bussa, S. & VEGAS Development Team 2012, AAS Meeting, 219, 446.10

Caprini, C., & Figueroa, D. G. 2018, CQGra, 35, 163001

Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538

Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970

CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5

Cordes, J. M., & Downs, G. S. 1985, ApJS, 59, 343

Cordes, J. M., & Rickett, B. J. 1998, ApJ, 507, 846

Cordes, J. M., & Shannon, R. M. 2010, arXiv:1010.3785

Cordes, J. M., Shannon, R. M., & Stinebring, D. R. 2016, ApJ, 817, 16

de Gasperin, F., Mevius, M., Rafferty, D. A., Intema, H. T., & Fallows, R. A. 2018, A&A, 615, A179

Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94

Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341

Detweiler, S. 1979, ApJ, 234, 1100

Di Marco, V., Zic, A., Shannon, R. M., & Thrane, E. 2024, MNRAS, 532, 4026

Dickey, J. M. 1971, Ann. Math. Stat., 41, 204

Dolch, T., Stinebring, D. R., Jones, G., et al. 2021, ApJ, 913, 98

Ellis, J., & van Haasteren, R. 2017, jellis18/PTMCMCSampler: Official Release v1.0.0, Zenodo, doi:10.5281/zenodo.1037579

Ellis, J. A., & Cornish, N. J. 2016, PhRvD, 93, 084048

Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2020, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE v3.0.0, Zenodo, doi:10.5281/zenodo.4059815

EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023a, A&A, 678, A50

EPTA Collaboration, Antoniadis, J., Babak, S., et al. 2023b, A&A, 678, A48

Foreman-Mackey, D. 2016, JOSS, 1, 24

Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021, MNRAS, 502, 478

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hazboun, J. S. 2020, La Forge, v0.3.0, Zenodo, doi:10.5281/zenodo.4152550

Hazboun, J. S., Simon, J., Madison, D. R., et al. 2022, ApJ, 929, 39

Hazboun, J. S., Simon, J., Siemens, X., & Romano, J. D. 2020a, ApJL, 905, L6

Hazboun, J. S., Simon, J., Taylor, S. R., et al. 2020b, ApJ, 890, 108

Hellings, R. W., & Downs, G. S. 1983, ApJL, 265, L39

Hemberger, D. A., & Stinebring, D. R. 2008, ApJL, 674, L37

Hunter, J. D. 2007, CSE, 9, 90

Jennings, R. J., Cordes, J. M., Chatterjee, S., et al. 2024, ApJ, 964, 179

Jones, M. L., McLaughlin, M. A., Lam, M. T., et al. 2017, ApJ, 841, 125

Joshi, B. C., Arumugasamy, P., Bagchi, M., et al. 2018, JApA, 39, 51

Kaiser, A. R., Pol, N. S., McLaughlin, M. A., et al. 2022, ApJ, 938, 115

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Keith, M. J., Coles, W., Shannon, R. M., et al. 2013, MNRAS, 429, 2161

Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017, MNRAS, 471, 4508

Kocsis, B., & Sesana, A. 2011, MNRAS, 411, 1467

Kulkarni, A. D., Shannon, R. M., Reardon, D. J., et al. 2024, MNRAS, 528, 3658

Kumar, P., White, S. M., Stovall, K., Dowell, J., & Taylor, G. B. 2022, MNRAS, 511, 3937

Laal, N., Lamb, W. G., Romano, J. D., et al. 2023, PhRvD, 108, 063008

Lam, M. T., Cordes, J. M., Chatterjee, S., & Dolch, T. 2015, ApJ, 801, 130

Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, ApJ, 821, 66

Lam, M. T., Ellis, J. A., Grillo, G., et al. 2018, ApJ, 861, 132

Lam, M. T., Lazio, T. J. W., Dolch, T., et al. 2020, ApJ, 892, 89

Lam, M. T., McLaughlin, M. A., Arzoumanian, Z., et al. 2019, ApJ, 872, 193

Lamb, W. G., Taylor, S. R., & van Haasteren, R. 2023, PhRvD, 108, 103019

Lang, K. R. 1971, ApJ, 164, 249
Lasky, P. D., Mingarelli, C. M. F., Smith, T. L., et al. 2016, PhRvX, 6, 011035

Lee, K. J. 2016, ASP Conf. Ser. 502, Frontiers in Radio Astronomy and FAST
Early Sciences Symposium 2015, ed. L. Qain & D. Li (San Francisco, CA:
ASP), 19

Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004

Lentati, L., Alexander, P., Hobson, M. P., et al. 2013, PhRvD, 87, 104021

Lentati, L., Kerr, M., Dai, S., et al. 2017, MNRAS, 466, 3706

Lentati, L., Shannon, R. M., Coles, W. A., et al. 2016, MNRAS, 458, 2161

Levin, L., McLaughlin, M. A., Jones, G., et al. 2016, ApJ, 818, 166

Lewandowski, W., Kowalińska, M., Kijak, J., et al. 2015, MNRAS, 449, 1570

Lin, F. X., Lin, H.-H., Luo, J., et al. 2021, MNRAS, 508, 1115

Liu, K., Desvignes, G., Cognard, I., et al. 2014, MNRAS, 443, 3752

Liu, Y., Main, R. A., Verbiest, J. P. W., et al. 2023, SCPMA, 66, 119512

Liu, Y., Verbiest, J. P. W., Main, R. A., et al. 2022, A&A, 664, A116

Lommen, A. N., Kipphorn, R. A., Nice, D. J., et al. 2006, ApJ, 642, 1012

Lorimer, D. R., & Kramer, M. 2012, Handbook of Pulsar Astronomy
(Cambridge: Cambridge Univ. Press)

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

Madison, D. R., Cordes, J. M., Arzoumanian, Z., et al. 2019, ApJ, 872, 150

Main, R. A., Antoniadis, J., Chen, S., et al. 2023, MNRAS, 525, 1079

Main, R. A., Sanidas, S. A., Antoniadis, J., et al. 2020, MNRAS, 499, 1468

Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017

Meyers, P. M., Chatzioioannou, K., Vallisneri, M., & Chua, A. J. K. 2023,
PhRvD, 108, 123008

Miles, M. T., Shannon, R. M., Bailes, M., et al. 2023, MNRAS, 519, 3976

Mingarelli, C. M. F., & Casey-Clyde, J. A. 2022, Sci, 378, 592

NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ,
813, 65

Nițu, I. C., Keith, M. J., Tiburzi, C., et al. 2024, MNRAS, 528, 3304

Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ,
161, 105

Pennucci, T. T., Demorest, P. B., & Ransom, S. M. 2014, ApJ, 790, 93

Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS,
490, 4666

Phinney, E. S. 2001, arXiv:astro-ph/0108028

Rajagopal, M., & Romani, R. W. 1995, ApJ, 446, 543

Rankin, J. M., & Roberts, J. A. 1971, in IAU Symp. 46, The Crab Nebula, ed.
R. D. Davies & F. Graham-Smith (Cambridge: Cambridge Univ. Press), 114

Ransom, S., Brazier, A., Chatterjee, S., et al. 2019, BAAS, 51, 195

Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine
Learning (Cambridge, MA: MIT Press)

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023a, ApJL, 951, L6

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023b, ApJL, 951, L7

Romani, R. W., Narayan, R., & Blandford, R. 1986, MNRAS, 220, 19

Romano, J. D., Hazboun, J. S., Siemens, X., & Archibald, A. M. 2021, PhRvD,
103, 063027

Sazhin, M. V. 1978, SvA, 22, 36

Sesana, A. 2013, CQGra, 30, 224014

Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192

Shannon, R. M., & Cordes, J. M. 2010, ApJ, 725, 1607

Shannon, R. M., & Cordes, J. M. 2017, MNRAS, 464, 2075

Shannon, R. M., Lentati, L. T., Kerr, M., et al. 2016, ApJL, 828, L1

Shapiro-Albert, B. J., Hazboun, J. S., McLaughlin, M. A., & Lam, M. T. 2021,
ApJ, 909, 219

Singha, J., Surnis, M. P., Joshi, B. C., et al. 2021, MNRAS, 507, L57

Sosa Fiscella, S. V., Lam, M. T., Arzoumanian, Z., et al. 2024, ApJ, 966, 95

Splaver, E. M., Nice, D. J., Stairs, I. H., Lommen, A. N., & Backer, D. 2005,
ApJ, 620, 405

Srivastava, A., Desai, S., Kolhe, N., et al. 2023, PhRvD, 108, 023008

Tarafdar, P., Nobleson, K., Rana, P., et al. 2022, PASA, 39, e053

Taylor, S. R. 2021, arXiv:2105.13270

Taylor, S. R., Baker, P. T., Hazboun, J. S., Simon, J., & Vigeland, S. J. 2021,
enterprise_extensions v2.4.3, GitHub, https://github.com/nanograv/
enterprise_extensions

Taylor, S. R., Simon, J., Schult, L., Pol, N., & Lamb, W. G. 2022, PhRvD, 105,
084049

Tiburzi, C., Shaifullah, G. M., Bassa, C. G., et al. 2021, A&A, 647, A84

Turner, J. E., McLaughlin, M. A., Cordes, J. M., et al. 2021, ApJ, 917, 10

Turner, J. E., Stinebring, D. R., McLaughlin, M. A., et al. 2023, ApJ, 944, 191
Vagnozzi, S. 2023, JHEAp, 39, 81

van Haasteren, R. 2013, MNRAS, 429, 55

van Haasteren, R., & Levin, Y. 2013, MNRAS, 428, 1147

van Haasteren, R., & Vallisneri, M. 2014, PhRvD, 90, 104012

Verbiest, J. P. W., Bailes, M., Coles, W. A., et al. 2009, MNRAS, 400, 951

Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458
1267

Wang, H., Taylor, S. R., & Vallisneri, M. 2019, MNRAS, 487, 3644

Waskom, M. 2021, JOSS, 6, 3021

Xu, H., Chen, S., Guo, Y., et al. 2023, RAA, 23, 075024

You, X. P., Hobbs, G., Coles, W. A., et al. 2007, MNRAS, 378, 493

Zic, A., Hobbs, G., Shannon, R. M., et al. 2022, MNRAS, 516, 410

Zic, A., Reardon, D. J., Kapur, A., et al. 2023, PASA, 40, e049