Refined Mass and Geometric Measurements of the
High-mass PSR J0740+6620

E. Fonseca, H.T. Cromartie, T.T. Pennucci, P.S. Ray, A.Yu. Kirichenko, S.M.

Ransom, P.B. Demorest, I.H. Stairs, Z. Arzoumanian, Lucas Guillemot, et al.

► To cite this version:

E. Fonseca, H.T. Cromartie, T.T. Pennucci, P.S. Ray, A.Yu. Kirichenko, et al. Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620. The Astrophysical Journal Letters, 2021, 915 (1), pp.L12. 10.3847/2041-8213/ac03b8. hal-03210324

HAL Id: hal-03210324

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

Submitted on 15 Jul 2025

**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
Refined Mass and Geometric Measurements of the High-Mass PSR J0740+6620

E. FONSECA,¹,²,³,⁴ H. T. CROMARTIE,⁵,⁶ T. T. PENNUCCI,⁷,⁸ P. S. RAY,⁹ A. YU. KIRICHENKO,¹⁰,¹¹ S. M. RANSOM,⁷

P. B. DEMOREST,¹² I. H. STAIRS,¹³ Z. ARZOUMANIAN,¹⁴ L. GUILLEMOT,¹⁵,¹⁶ A. PARTHASARATHY,¹⁷ M. KERR,⁹

I. COGNARD,¹⁵,¹⁶ P. T. BAKER,¹⁸ H. BLUMER,³,⁴ P. R. BROOK,³,⁴ M. DECESAR,¹⁹ T. DOLCH,²⁰,²¹ F. A. DONG,¹³

E. C. FERRARA,²²,²³,²⁴ W. FIORE,³,⁴ N. GARVER-DANIELS,³,⁴ D. C. GOOD,¹³ R. JENNINGS,²⁵ M. L. JONES,²⁶

V. M. KASPI,¹,² M. T. LAM,²⁷,²⁸ D. R. LORIMER,³,⁴ J. LUO,²⁹ A. MCEWEN,²⁶ J. W. MCKEE,²⁹ M. A. MCLAUGHLIN,³,⁴

N. MCMANN,³⁰ B. W. MEYERS,¹³ A. NAIDU,³¹ C. NG,³² D. J. NICE,³³ N. POL,³⁰ H. A. RADOVAN,³⁴

B. SHAPIRO-ALBERT,³,⁴ C. M. TAN,¹,² S. P. TENDULKAR,³⁵,³⁶ J. K. SWIGGUM,³³ H. M. WAHL,³,⁴ AND W. W. ZHU³⁷

¹Department of Physics, McGill University, 3600 rue University, Montréal, QC H3A 2T8, Canada

²McGill Space Institute, McGill University, 3550 rue University, Montréal, QC H3A 2A7, Canada

³Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506-6315, USA

⁴Center for Gravitational Waves and Cosmology, Chestnut Ridge Research Building, Morgantown, WV 26505, USA

⁵Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

⁶NASA Hubble Fellowship Program Einstein Postdoctoral Fellow

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

⁸Institute of Physics, Eötvös Lorán University, Pázmány P. s. 1/A, 1117 Budapest, Hungary

⁹U.S. Naval Research Laboratory, Washington, DC 20375, USA

¹⁰Universidad Nacional Autónoma de México, Instituto de Astronomía, AP 106, Ensenada 22800, BC, México

¹¹Ioffe Institute, Politeknichnicheskaya 26, St. Petersburg, 194021, Russia

¹²National Radio Astronomy Observatory, P.O. Box O, Socorro, NM 87801, USA

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

¹⁴Astrophysics Science Division, Code 662, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹⁵Laboratoire de Physique et Chimie de l'Environnement et de l'Espace – Université d'Orléans / CNRS, F-45071 Orléans Cedex 02, France

¹⁶Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, F-18330 Nançay, France

¹⁷Max Planck Institute for Radio Astronomy, Auf dem H "u gel 69, D-53121 Bonn, Germany

¹⁸Department of Physics and Astronomy, Widener University, Chester, PA 19013, USA

¹⁹Department of Physics and Astronomy, George Mason University, Fairfax, VA 22030, resident at U.S. Naval Research Laboratory, Washington, D.C. 20375, USA

²⁰Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA

²¹Eureka Scientific, Inc. 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA

²²Department of Astronomy, University of Maryland, College Park, MD 20742, USA

²³Center for Exploration and Space Studies (CRESST)

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

²⁵Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

²⁶Center for Gravitation, Cosmology, and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, 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

²⁹Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada

³⁰Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA

³¹University of Oxford, Sub-Department of Astrophysics, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH, United Kingdom

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

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

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

³⁵Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Mumbai, 400005, India

³⁶National Centre for Radio Astrophysics, Post Bag 3, Ganeshkhind, Pune, 411007, India

³⁷National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China

(Received XXX; Revised YYY; Accepted ZZZ)

Submitted to ApJ Letters
ABSTRACT

We report results from continued timing observations of PSR J0740+6620, a high-mass, 2.8-ms radio pulsar in orbit with a likely ultra-cool white dwarf companion. Our data set consists of combined pulse arrival-time measurements made with the 100-m Green Bank Telescope and the Canadian Hydrogen Intensity Mapping Experiment telescope. We explore the significance of timing-based phenomena arising from general-relativistic dynamics and variations in pulse dispersion. When using various statistical methods, we find that combining ~1.5 years of additional, high-cadence timing data with previous measurements confirms and improves upon previous estimates of relativistic effects within the PSR J0740+6620 system, with the pulsar mass $m_p = 2.08_{-0.07}^{+0.07} M_\odot$ (68.3% credibility) determined by the relativistic Shapiro time delay. For the first time, we measure secular variation in the orbital period and argue that this effect arises from apparent acceleration due to significant transverse motion. After incorporating contributions from Galactic differential rotation and off-plane acceleration in the Galactic potential, we obtain a model-dependent distance of $d = 1.14_{-0.15}^{+0.17}$ kpc (68.3% credibility). This improved distance confirms the ultra-cool nature of the white dwarf companion determined from recent optical observations. We discuss the prospects for future observations with next-generation facilities, which will likely improve the precision on $m_p$ for J0740+6620 by an order of magnitude within the next few years.

Keywords: pulsars: general – neutron stars: masses – pulsars: timing – etc.

1. INTRODUCTION

The masses and radii of neutron stars are a central focus in many observational high-energy experiments. Members of the neutron-star population with the largest masses are particularly important. A key application of these measurements is in testing proposed neutron-star equations of state (EoSs), of which there are many that incorporate baryonic and/or “exotic” compositions not accessible through terrestrial experiments (Özel & Freire 2016). For many decades, the primary means of measuring precise neutron star masses was through timing of radio pulsars in orbital systems that exhibit relativistic dynamical processes (e.g., Stairs 2003). Observatories sensitive to X-ray and gravitational radiation have begun to yield similar estimates of masses and radii, compelling statistical inference of fundamental physics and EoS constraints from neutron stars into a new and exciting era (e.g., Abbott et al. 2018; Miller et al. 2019; Riley et al. 2019).

A recent example of important observational constraints arises from the timing of PSR J0740+6620, a 2.8-ms millisecond pulsar (MSP) in a near-circular, 4.7-day orbit discovered by the Green Bank North Celestial Cap Survey (GBNCC; Lynch et al. 2018). Measurement of the relativistic Shapiro time delay (Shapiro 1964) in the PSR J0740+6620 system allowed for an estimate of the pulsar mass, $m_p = 2.14_{-0.09}^{+0.10} M_\odot$ (68.3% credibility; Cromartie et al. 2020). An optical search for the white-dwarf companion to PSR J0740+6620 yielded a detection of a likely ultra-cool counterpart (Beronya et al. 2019), though a firm association is dependent on a reliable measure of distance that was not obtainable at the

time. These radio and optical-photometric data were shown to be consistent with binary-stellar-evolution simulations of a progenitor system undergoing conservative mass transfer within a Hubble time, regardless of whether irradiation feedback was prominent in system evolution (Echeveste et al. 2020).

Along with PSRs J0348+0423 (Antoniadis et al. 2013) and J1614−2230 (e.g., Demorest et al. 2010), PSR J0740+6620 resides at the high-mass extremum of the known Galactic pulsar population. In addition, and while subject to greater modeling uncertainty, observations of “redback” pulsars and their irradiated companions clearly suggest a substantial population of high-mass neutron stars (e.g., van Kerkwijk et al. 2011; Linares et al. 2018; Romani et al. 2021). Gravitational-wave observations of GW190814 indicate a high-mass compact object that could be the heaviest neutron star known (Abbott et al. 2020). A small yet growing population, high-mass neutron stars serve as the most important constraints on the maximum allowed masses and, by extension, the EoS models that predict these measurements (e.g., Hu et al. 2020; Raaijmakers et al. 2020; Tsokaros et al. 2020). Further study of the high-mass pulsar population will also allow for an ensemble estimation of the maximum mass, as well as the study of distribution moments that reflect underlying formation processes, before or during the spin-up to millisecond rotational periods, that produce high-mass neutron stars (e.g., Antoniadis et al. 2016).

Recently, the Neutron Star Interior Composition Explorer (NICER) measured the mass and radius of the isolated PSR J0030+0451 through modeling of X-ray ther-
mal emission from hotspots on the neutron-star surface (Bogdanov et al. 2019; Miller et al. 2019; Raaijmakers et al. 2019; Riley et al. 2019). Constraints on neutron-star radii from NICER observations are strengthened with a priori knowledge of their masses. For sufficiently bright pulsars in binary systems, typically with system inclination $i \gtrsim 60^\circ$, mass and geometric measurements can be achieved with high-precision pulsar timing. In late 2019, NICER confirmed X-ray pulsations from PSR J0740+6620 (Wolff et al. 2021) and began a dedicated observing program to infer the stellar compactness via modeling of its pulsed X-ray emission. We endeavor to improve upon the mass measurement of Cromartie et al. (2020) and constrain other physical parameters in support of NICER observations and their constraints on the radius.

In this work, we present observations, analysis and results from ongoing timing observations of PSR J0740+6620 with two different telescopes. In Section 2, we overview the instruments and processing methods used to acquire pulse arrival-time data for PSR J0740+6620. In Section 3, we describe the procedure used to simultaneously model arrival times obtained with different instruments and receivers. In Section 4, we outline the theoretical framework used in order to obtain robust constraints on the Shapiro delay and distance to PSR J0740+6620. In Section 5, we discuss various aspects of our analysis and results, including dispersion measure (DM) modeling and impacts our measurements have on recent optical observations. Finally, in Section 6, we summarize our results and speculate about future improvement of our measurements.

## 2. OBSERVATIONS & REDUCTION

Continued observations of PSR J0740+6620 were conducted using the 100-m Green Bank Telescope (GBT) and the Canadian Hydrogen Intensity Mapping Experiment (CHIME) telescope. A summary of observing parameters and statistics is presented in Table 1.

### 2.1. Telescopes

Nearly all data from the single-dish GBT were acquired as part of the ongoing observing program of the North American Nanohertz Observatory for Gravitational Waves (NANOGrav<fnref n="1" />), which monitors an array of nearly 100 MSPs for signatures of gravitational radiation at nanohertz frequencies. Details of the NANOGrav observing program are given by Alam et al. (2021a,b), and we provide a summary here.

At the GBT, raw telescope voltage data were coherently dedispersed and processed in real time through the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al. 2008). During its operation, GUPPI produced archives of time-integrated (or “folded”) pulse profiles, and for all elements of the Stokes polarization vector, over 10-s time bins. Folded profiles were collected using two different radio receivers at the GBT, with the receivers centred at radio frequencies of 820 MHz and 1400 MHz; these receivers possessed full bandwidths of 200 MHz and 800 MHz, respectively. Moreover, the spectrum of each folded profile was resolved over 128 and 512 channels when using the 820-MHz and 1400-MHz receivers, respectively.

The GBT data set analyzed in this work also contains supplemental observations acquired by Cromartie et al. (2020) near and during superior conjunction in the PSR J0740+6620 binary system. Targeted observations have been shown to increase the measurement significance of the Shapiro time delay, which is at its maximum amplitude during superior conjunction (orbital phase of 0.25, relative to the ascending node; e.g., Demorest et al. 2010; Pennucci 2015). The targeted observations of PSR J0740+6620 occurred over three sessions: one 6-hour session occurred on MJD 58368 using the 820-MHz receiver during orbital phase of 0.25; and two 5-hour sessions were undertaken on MJDs 58448 and 58449 using the 1400-MHz receiver, corresponding to orbital phases of 0.15 and 0.25, respectively.

The CHIME telescope is a static radio interferometer comprised of four half-cylinder reflectors and a total of 1,024 dual-polarization antennas sensitive to the 400–800 MHz range. Digitized telescope voltage data from all CHIME feeds are coherently averaged using time-dependent phase delays corresponding to 10 different sources, instantaneously yielding 10 “tracking” voltage timeseries. The CHIME/Pulsar backend (CHIME/Pulsar Collaboration et al. 2021) is a ten-node computing system that receives these 10 independent, complex voltage streams – at resolutions of 2.56-μs and 1,024 frequency channels – and uses dspsr (van Straten & Bailes 2011) to perform real-time coherent dedispersion and folding. The resultant products are similar in structure to those acquired with GUPPI, though they contain 1,024 channels evaluated across the 400 MHz bandwidth of CHIME for each recorded profile. At the declination of PSR J0740+6620, acquisitions typically occur with ~ 33-minute duration in order to record a full transit at the CHIME telescope.

### 2.2. Observing Cadences

<footnote n="1">¹ nanograv.org</footnote>
<table><caption>**Table 1.** A summary of observing parameters and data sets of PSR J0740+6620 analyzed in this work.</caption><thead><tr><th>Parameter</th><th>GBT/1400-MHz</th><th>GBT/820-MHz</th><th>CHIME</th></tr></thead><tbody><tr><td>Time Range (MJD)</td><td>56675-58945</td><td>56640-58944</td><td>58517-58975</td></tr><tr><td>Frequency Range (MHz)</td><td>1151-1885</td><td>722-919</td><td>400-800</td></tr><tr><td>Number of Channels, raw</td><td>512</td><td>128</td><td>1024</td></tr><tr><td>Number of Channels, downsampled</td><td>64</td><td>64</td><td>32</td></tr><tr><td>Number of Observing Epochs</td><td>187</td><td>142</td><td>263</td></tr><tr><td>Number of TOAs, narrowband</td><td>8,208</td><td>5,191</td><td>4,862</td></tr><tr><td>Number of TOAs, wideband</td><td>209</td><td>154</td><td>263</td></tr><tr><td>Min./Median/Max. σ<sub>TOA</sub>, narrowband (μs)</td><td>0.24/2.65/32.7</td><td>0.16/3.16/24.9</td><td>0.49/5.21/62.8</td></tr><tr><td>Min./Median/Max. σ<sub>TOA</sub>, wideband (μs)</td><td>0.05/0.39/2.04</td><td>0.04/0.51/3.10</td><td>0.14/0.97/3.28</td></tr></tbody></table>

Under the nominal NANOGgrav program, multi-frequency observations of PSR J0740+6620 with the GBT typically lasted for ~25 minutes per receiver and occurred with a monthly cadence; observations with each receiver were acquired within ~3 days of one another. After 2018 September 24 (i.e., MJD 58385), observations of PSR J0740+6620 were prioritized to occur during every scheduled NANOGgrav session, which led to 3–5 multi-frequency observations within a one-week period of time during each month. See Alam et al. (2021a,b) for additional discussion on the NANOGgrav observing program.

Folded data were acquired with the CHIME/Pulsar coherent-dedispersion backend on a near-daily basis from 2019 February 3 to 2020 May 6 (i.e., MJD range 58517–58975). Exceptions to the near-daily cadence largely correspond to ~one-week periods where the CHIME correlator ceased operation for software upgrades and hardware maintenance. Such activities only occurred 3 – 4 times per year during the time span of CHIME data presented in this work.

## 2.3. Offline Processing, Calibration and Downsampling

All data acquired with the GBT and CHIME/Pulsar were processed offline using the psrchive<fnref n="2" /> suite of analysis utilities, through the use of separate psrchive-based procedures tailored to the NANOGgrav and CHIME/Pulsar data sets. For NANOGgrav data on PSR J0740+6620, we used the nanopipe pipeline<fnref n="3" /> – which is employed in ongoing NANOGgrav analysis of all its sources – for the following manipulations: mask-based excision of radio frequency interference (RFI); calibration of flux density and polarization informa-

tion using on/off-source observations of the compact radio source J1445+0958 (B1442+101) modulated with a pulsed noise diode; downsampling to frequency resolution as low as 3.125 MHz per channel; and full integration of profiles across the duration of each scan.

CHIME/Pulsar data were processed using a similar procedure, but with three modifications: (1) the Stokes profiles were not calibrated due to ongoing work in developing calibration methods for CHIME/Pulsar data; (2) spectra were downsampled to 32 channels with 12.5 MHz bandwidth, in order to retain adequate signal to noise (S/N) in as many channels as possible; and (3) ~1 minute of profile data acquired before and after the transit of PSR J0740+6620 at CHIME were discarded prior to integration of profile data over time. The third modification was performed in order to minimize overweighting of CHIME/Pulsar profiles recorded at low frequencies, where the primary beam of CHIME is wider and is thus sensitive to transiting pulsars for longer periods of time.

## 2.4. Computation of Arrival Times

We generated two distinct sets of times-of-arrival (TOAs) and TOA uncertainties ($\sigma_{TOA}$) from processed GBT and CHIME/Pulsar data: a set of frequency-resolved TOAs, which we refer to as “narrowband” TOAs; and a set of “wideband” TOAs, which extract a single arrival time measurement from each observation.

The narrowband TOAs were measured using the Fourier phase-gradient technique described in Taylor (1992); the observed total-intensity (i.e., Stokes-I) profile from each frequency channel is cross-correlated with a de-noised template to measure a best-fit phase offset. The template profile is constructed for each telescope and receiver combination by averaging all available per-receiver timing data, and then wavelet-smoothing the result (e.g., NANOGgrav Collaboration et al. 2015). These

<footnote n="2">² http://psrchive.sourceforge.net/</footnote>

<footnote n="3">³ https://github.com/demorest/nanopipe</footnote>
techniques of template and TOA generation form a stan-
dard protocol in pulsar astronomy, and were employed
by Cromartie et al. (2020) in their analysis of the PSR
J0740+6620 system.

In this work, we obtained a maximum of 64 narrow-
band TOAs for each receiver bandwidth and all pro-
cessed, RFI-cleaned observations. Narrowband TOAs
derived from GBT data were further cleaned by remov-
ing outlier timing data using the NANOGrav analysis
procedure outlined by Alam et al. (2021a), where TOAs
were discarded from analysis based on: low profile fi-
delity (i.e., with a pulse profile S/N < 8); corrupted
calibration data; a probabilistic classification of out-
lier TOAs based on prior timing models (Vallisneri &
van Haasteren 2017); and manual TOA excision. Nar-
rowband TOAs derived from CHIME/Pulsar data were
cleaned using a similar procedure, though only con-
sisted of rejecting low-S/N profiles and performing small
amounts of manual TOA excision. Work is underway to
fully integrate CHIME/Pulsar data into the NANOGrav
analysis infrastructure and will be presented in future
analyses.

Wideband TOAs are a more recent innovation (Liu
et al. 2014; Pennucci et al. 2014; Pennucci 2019), but
have been incorporated in the recently released 12.5-yr
NANOGrav data set (Alam et al. 2021b). The wideband
TOAs used for our analysis were generated using the
procedure described by Pennucci et al. (2014)<fnref n="4" />, which
we briefly summarize as follows. A single, wideband
TOA and instantaneous DM were estimated for each ob-
servation using a generalized Fourier phase-gradient al-
gorithm that constrains the frequency-dependent phase
offsets. In order to derive a “portrait” template for each
receiver band, we first integrated all available data and
computed a high-S/N template; the three mean por-
traits are conjoined and shown in Figure 1. The de-
noised portraits, computed for each receiver band from
the averaged templates, were then derived through a
combination of principal component analysis, wavelet
smoothing, and spline interpolation. The only TOAs
excised from the wideband data set were those from low-
S/N observations (with S/N < 25; see Alam et al. 2021b,
for further explanation).

3. TIMING METHODS

We constructed separate models of timing variations
in PSR J0740+6620 using the narrowband and wide-
band data sets described in Section 2.4. For each
TOA set, we combined NANOGrav and CHIME/Pulsar

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

**Figure 1.** A composite spectrum of mean pulse profiles for PSR J0740+6620 across the three observed CHIME and GBT bands (Table 1). For each receiver, we integrated all data across observing time and down-sampled in frequency and phase to match the CHIME/Pulsar resolution (12.5 MHz channels, 1024 phase bins). The amplitudes of each profile are scaled such that all profiles have unit integrated flux density. The profiles from each distinct band have been aligned for display against an arbitrary template of fixed Gaussian components, and the frequency-averaged profile from each band is displayed in the top panel. The horizontal line and bar in the top panel both denote timescale, with the line spanning 1 ms and the bar spanning 100 µs. The large gap around 1000 MHz is the gap between the GBT bands, and the narrow gap around 720 MHz exists because the upper 20% of the CHIME band was excised for this plot. Other small gaps are places where there was persistent RFI or band roll-off.

TOAs and used the tempo pulsar-timing package<fnref n="5" />
for conversion and generalized least-squares modeling
(Coles et al. 2011) in terms of parameters that are spe-
cific to relevant physical processes as well as instrumen-
tal effects. A summary of the timing procedure and
model is provided below.

<footnote n="5">⁵ https://ascl.net/1509.002</footnote>

<footnote n="4">⁴ https://github.com/pennucci/PulsePortraiture</footnote>
### 3.1. Pre-Modeling Conversion of TOAs

For the CHIME/Pulsar data set, a small portion of the earliest TOAs were impacted by timing offsets that occurred upon restarts of the CHIME correlator. These offsets arose due to improper packaging of timing data required for downstream, real-time processing and determination of UTC timestamps. The relevant correlator software was amended to resolve this error, and all CHIME/Pulsar data acquired after MJD 58600 no longer contain achromatic timing offsets between observing cycles. We used configuration metadata recorded by the CHIME software infrastructure to determine the exact time steps introduced by two correlator offsets, and used those values as corrections to data acquired between MJDs 58517–58543.<fnref n="6" />

### 3.2. Construction of the Timing Model

The barycentric TOAs were then modeled against a superposition of analytic time delays describing various effects of physical and instrumental origin (e.g., Backer & Hellings 1986). Parameters of these time delays consist of terms describing: pulsar spin and its spin-down evolution; astrometric effects, which yield measures of position, proper motion, and apparent variations of position due to timing parallax ($\bar{\omega}$); piecewise-constant estimates of DM across the time span of the combined data set, referred to below as “DMX”; and orbital motion.

We used the ELL1 binary timing model (Lange et al. 2001) to model the near-circular Keplerian motion of PSR J0740+6620 in its orbit with the companion white dwarf. The five Keplerian elements in our model consist of: the orbital period ($P_b$); the projected semi-major axis along the line of sight ($x$); the epoch of passage through the ascending node of the binary system ($T_{asc}$); and the “Laplace-Lagrange” eccentricity parameters, $\epsilon_1$ and $\epsilon_2$, that describe departures from circular motion. The ELL1 parameters are related to the traditional orbital elements – eccentricity ($e$), argument of periastron ($\bar{\omega}$), and epoch of passage through periastron ($T_0$) – in the following manner:

$$e = \sqrt{\epsilon_1^2 + \epsilon_2^2} \qquad (1)$$

$$\omega = \arctan(\epsilon_1/\epsilon_2) \qquad (2)$$

$$T_0 = T_{asc} + P_b \times \omega/(2\pi). \qquad (3)$$

Following Cromartie et al. (2020), we also included the Shapiro delay in our modeling of the extended data set,

described by the “range” ($r$) and “shape” ($s$) parameters. While $s = \sin i$, where $i$ is the inclination of the binary system relative to the plane of the sky, the exact expression of $r$ depends on the assumed theory of gravitation (Damour & Taylor 1992); we assumed that general relativity is valid at the level of timing precision achieved in this work, which leads to $r = T_\odot m_c$, where $T_\odot = GM_\odot/c^3 = 4.925490947 \ \mu s$ and $m_c$ is the mass of the companion star. With estimates of $m_c$ and $\sin i$, the pulsar mass ($m_p$) is determined by the Keplerian mass function,

$$f_m = \frac{4\pi^2 x^3}{T_\odot P_b^2} = \frac{(m_c \sin i)^3}{(m_p + m_c)^2}. \qquad (4)$$

The significance of all model parameters was determined using the F-test criterion, as described in Alam et al. (2021a) for standard NANOGrav analysis, when comparing best-fit models that included or ignored each relevant set of parameters. We chose a $p$-value threshold of 0.0027 to reject the null hypothesis that adding or removing a parameter results in a variation of the residuals consistent with noise. Under Gaussian interpretation, a threshold of 0.0027 corresponds to a $3\sigma$ deviation. Parameters were included in the model if their inclusion resulted in rejection of the null hypothesis based on this test. We explored fitting for secular variations in the orbital elements using the F-test criterion, and found that the time derivative of the orbital period ($\dot{P}_b$) possessed sufficient statistical significance for inclusion as a degree of freedom. We interpret and discuss the implications of the $\dot{P}_b$ measurement in several sections below.

### 3.3. Frequency Evolution of Pulse Profiles

One of the key differences between the commonly-used narrowband methods of TOA estimation and wideband TOA estimation is their treatment of intrinsic profile evolution across receiver bandpasses. A single, achromatic template profile is used to measure our narrowband TOAs. Systematic time delays arise in the TOAs as a function of frequency in the presence of significant variation of the profile shape, which can be seen in Figure 1. We modeled these frequency-dependent (FD) variations in our narrowband TOA data set using the heuristic model developed by NANOGrav Collaboration et al. (2015), where the associated time delay $\Delta t_{FD} = \sum_{i=1}^n c_i \ln(f/1\text{GHz})^i$. The coefficients $c_i$ are free parameters in tempo, and we included three FD coefficients – with the number of coefficients determined using the F-test procedure described above – in all timing models derived from narrowband TOAs.

For the wideband TOAs, the use of a high-fidelity frequency-dependent template ameliorates the timing

<footnote n="6">⁶ The values of these offsets are noted in the “-to” metadata flag in the TOA data made publicly available in this work.</footnote>
biases introduced by profile evolution. The FD parameters were not significant when checked by the F-test. However, it is necessary to include three timing model parameters that quantify an offset in the average DM measured in each receiver band due to template profile misalignment; these DM parameters are analogous to the common phase “JUMP” parameters in tempo (see the “DMJUMP” parameters described in Alam et al. 2021b).

### 3.4. Analysis of Noise Properties

A common procedure in pulsar timing is the analysis of noise properties in TOA residuals (the differences between measured TOAs and their values predicted by the best-fit timing model). This type of analysis involves using a number of heuristic parameters such that the uncertainty-normalized residuals are normally distributed with unit variance. One set of these parameters adjusts the TOA uncertainties by a multiplicative factor ($F_k$), and a second set of these parameters is added to the TOA uncertainties in quadrature ($Q_k$). The subscript $k$ denotes the unique data subset for a set of frontend/backend combinations, and corresponds to the three subsets presented in Table 1 for this work. The factor $F_k$ accounts for incorrect estimation of $\sigma_{\text{TOA}}$ due to mismatch between the data and template profiles, while $Q_k$ characterizes underlying additive white noise (i.e., statistical fluctuations with constant power spectral density across all times and frequencies). A third set of parameters ($C_k$) encapsulates underlying noise processes that are uncorrelated in time but instead are fully correlated across observing frequency, and thus applies to simultaneously observed, multi-frequency TOAs. One source for $C_k$ is the stochastic scatter introduced from pulse “jitter” (e.g., NANOGgrav Collaboration et al. 2015; Lam et al. 2017).

We used the same methodology employed by Alam et al. (2021a,b) for determining the optimal values of {$F_k, Q_k, C_k$} for the three narrowband data subsets listed in Table 1. For all timing solutions, we used the enterprise Bayesian pulsar timing suite (Ellis et al. 2020) for sampling the {$F_k, Q_k, C_k$} terms while marginalizing over all other free parameters in the timing model. Following Alam et al. (2021b), a slightly different noise model is used in the wideband analysis. $C_k$ cannot be modeled for wideband TOAs because there are no simultaneously observed multi-frequency TOAs. However, in addition to {$F_k, Q_k$} for the TOA uncertainties, analogous error-scaling values of $F_k$ for the DM uncertainties were also modeled. In all cases, the timing solution was refined until convergence using tempo by

applying the noise model from the enterprise analysis.

We also explored the significance of temporal correlations in TOA residuals for PSR J0740+6620. Such “red” noise has been seen in other NANOGgrav MSPs and is understood to reflect irregularities in pulsar-spin rotation (Shannon & Cordes 2010). However, we determined that red noise is not prominent in the timing of PSR J0740+6620, as its estimation with enterprise yielded a Bayes significance factor of ~1. We therefore did not include terms that quantify red noise in the noise model for PSR J0740+6620, as NANOGgrav sets a Bayes factor threshold of 100 for including red-noise parameters into timing models (e.g., NANOGgrav Collaboration et al. 2015).

## 4. ANALYSIS

We explored several methods for statistical analysis described below to obtain robust credible intervals on key parameters of the PSR J0740+6620 system, using the modeling procedure described in Section 3. The following analysis methods are agnostic to the use of narrowband or wideband TOAs. However, we chose to use the wideband TOAs in the work presented below given the consistency in parameters from narrowband and wideband TOAs in Table 2. Moreover, the following methods benefit from relaxed computational requirements afforded by the considerably smaller size (~1/20) of the wideband TOA data set. The timing residuals and DM time series from the wideband data set are shown in Figure 2; a zoom-in of the overlap between the GBT and CHIME/Pulsar data sets is provided in Figure 3.

### 4.1. DMX Models

Table 2 presents two tempo models, derived from independent fits to the narrowband and wideband TOA data sets, when using the same DMX method employed by Cromartie et al. (2020). However, the different observing cadences within each TOA subset can lead to an uneven weighting of DMX estimates that impact the measurement of low-amplitude timing effects. In order to assess these impacts, we generated three sets of timing models that used slightly different DMX bins across the combined data set. Two of these DMX models used a single bin width of 3 and 6.5 days, respectively. In the third model, we used a hybrid scheme where the modeling of data acquired between MJDs 56640 and 58400 used a DMX model with bin size of 6.5 days, and all data acquired after MJD 58400 (when heightened-cadence observations began) were modeled with a DMX model using a bin size of 3 days.

These three DMX binning choices form a small subset of possible DM variation models. Alternative methods
<table><caption>**Table 2.** A summary of timing parameters for PSR J0740+6620 when using DMX to estimate DM values<sup>a</sup> in 6.5-day bins.</caption><thead><tr><th>Parameter<sup>b</sup></th><th>Narrowband TOAs</th><th>Wideband TOAs</th><th>Difference (σ<sub>max</sub>)<sup>c</sup></th></tr></thead><tbody><tr><td colspan="4"><strong>Astrometry</strong></td></tr><tr><td>Ecliptic longitude, λ (degrees)</td><td>103.759135333(12)</td><td>103.759135338(13)</td><td>0.4</td></tr><tr><td>Ecliptic latitude, β (degrees)</td><td>44.102478368(13)</td><td>44.102478361(13)</td><td>0.5</td></tr><tr><td>Proper motion in λ, μ<sub>λ</sub> (mas yr<sup>-1</sup>)</td><td>-2.735(14)</td><td>-2.737(15)</td><td>0.1</td></tr><tr><td>Proper motion in β, μ<sub>β</sub> (mas yr<sup>-1</sup>)</td><td>-32.48(2)</td><td>-32.48(2)</td><td>0.0</td></tr><tr><td>Timing parallax, ϖ (mas)</td><td>1.04(18)</td><td>0.87(19)</td><td>0.9</td></tr><tr><td colspan="4"><strong>Spin</strong></td></tr><tr><td>Spin frequency, ν<sub>s</sub> (s<sup>-1</sup>)</td><td>346.5319964608338(3)</td><td>346.5319964608337(3)</td><td>0.0</td></tr><tr><td>Time rate of change in frequency, ṽ<sub>s</sub> (10<sup>-15</sup> s<sup>-2</sup>)</td><td>-1.463874(11)</td><td>-1.463870(11)</td><td>0.4</td></tr><tr><td colspan="4"><strong>Binary Motion</strong></td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>4.76694461933(8)</td><td>4.76694461936(8)</td><td>0.4</td></tr><tr><td>Projected semi-major axis, x (lt-s)</td><td>3.97755608(10)</td><td>3.97755607(11)</td><td>0.0</td></tr><tr><td>First Laplace-Lagrange eccentricity parameter, ε<sub>1</sub> (10<sup>-6</sup>)</td><td>-5.68(3)</td><td>-5.70(3)</td><td>0.7</td></tr><tr><td>Second Laplace-Lagrange eccentricity parameter, ε<sub>2</sub> (10<sup>-6</sup>)</td><td>-1.833(18)</td><td>-1.840(19)</td><td>0.3</td></tr><tr><td>Epoch of ascending-node passage, T<sub>asc</sub></td><td>57804.731308893(17)</td><td>57804.731308895(18)</td><td>0.1</td></tr><tr><td>Time rate of change in period, ṡ<sub>b</sub> (10<sup>-12</sup> s s<sup>-1</sup>)</td><td>1.2(2)</td><td>1.2(2)</td><td>0.2</td></tr><tr><td>Companion mass, m<sub>c</sub> (M<sub>☉</sub>)</td><td>0.251(5)</td><td>0.253(6)</td><td>0.3</td></tr><tr><td>Sine of inclination angle, sin i</td><td>0.99909(12)</td><td>0.99908(13)</td><td>0.1</td></tr><tr><td colspan="4"><strong>Profile Evolution</strong></td></tr><tr><td>First coefficient of FD expansion, c<sub>1</sub> (10<sup>-5</sup>)</td><td>-3.2(2)</td><td>n/a</td><td>n/a</td></tr><tr><td>Second coefficient of FD expansion, c<sub>2</sub> (10<sup>-5</sup>)</td><td>-2.0(3)</td><td>n/a</td><td>n/a</td></tr><tr><td>Third coefficient of FD expansion, c<sub>3</sub> (10<sup>-5</sup>)</td><td>-1.0(2)</td><td>n/a</td><td>n/a</td></tr><tr><td colspan="4"><strong>Fit Configuration and Statistics</strong></td></tr><tr><td>Reference epoch for spin, astrometry (MJD)</td><td>57807</td><td>57807</td><td>n/a</td></tr><tr><td>Terrestrial clock standard</td><td>TT(BIPM2019)</td><td>TT(BIPM2019)</td><td>n/a</td></tr><tr><td>Solar system ephemeris</td><td>DE438</td><td>DE438</td><td>n/a</td></tr><tr><td>Barycentric time scale</td><td>TDB</td><td>TDB</td><td>n/a</td></tr><tr><td>Degrees of freedom</td><td>18,041</td><td>1,099</td><td>n/a</td></tr><tr><td>Goodness of fit, χ<sup>2</sup></td><td>18,183</td><td>1,072</td><td>n/a</td></tr><tr><td>Daily-averaged, weighted RMS (μs)</td><td>0.27</td><td>0.28</td><td>n/a</td></tr></tbody></table>

<sup>a</sup>We do not list DM values in this table due to the large number of DMX bins used to model DM variations. We refer the reader to Figures 2 and 3 for graphical representations of the best-fit DMX data.

<sup>b</sup>All parameter uncertainties, listed in parentheses, denote the 68.3% (i.e., 1 σ) confidence intervals in the preceding digit(s) as obtained from tempo.

<sup>c</sup>Differences between best-fit parameter estimates derived from narrowband and wideband TOAs are listed in units of the larger of the two statistical uncertainties (σ<sub>max</sub>).

for modeling DM variation include the Fourier decomposition of frequency-resolved TOAs as implemented in the tempoNest Bayesian analysis suite (Lentati et al. 2014). We nonetheless chose to use DMX as this method is able to resolve discrete, stochastic variations in DM that cannot otherwise be adequately modeled with, for example, a Taylor expansion of the DM.

## 4.2. *Secular Variation in the Orbital Period*

All fitted timing models exhibit a significant and consistent $\dot{P}_b$. Our analysis represents the first time such variations have been detected in the PSR J0740+6620 system. Several mechanisms can yield apparent or intrinsic variations in $P_b$, such as those due to: energy loss from quadrupole-order gravitational radiation ($\langle\dot{P}_b\rangle_{GR}$; e.g., Damour & Taylor 1992); differential rotation in the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Timing residuals from the best-fit timing model for wideband TOAs (top), and DM time series (bottom). The colors indicate the receiver and telescope of the observation in both panels. The best-fit DMX model parameters (grey crosses) are determined by both the wideband DM measurements and the wideband TOAs. The plotted ΔDM values are offset from the nominal zero value of 14.9631 pc cm⁻³. These wideband residual data and DMX model parameters correspond to the timing model that used the “hybrid” DMX-binning scheme discussed in Section 4.1. A zoom-in of the CHIME/Pulsar data is provided in Figure 3.</figcaption></figure>

Galaxy ((Ṗb)DR; e.g., Nice & Taylor 1995); off-plane acceleration in the Galactic gravitational potential ((Ṗb)z; Kuijken & Gilmore 1989); and apparent acceleration due to transverse motion ((Ṗb)μ; Shklovskii 1970). We assumed that the total observed variation is therefore (Ṗb)obs = (Ṗb)GR + (Ṗb)DR + (Ṗb)z + (Ṗb)μ, where

$$ (\dot{P}_b)_{\text{GR}} = -\frac{192\pi}{5} \frac{(n_b T_{\odot})^{5/3}}{(1 - e^2)^{7/2}} \frac{m_p m_c}{m_{\text{tot}}^{1/3}} \left(1 + \frac{73}{24}e^2 + \frac{37}{96}e^4\right), \quad (5) $$

$$ (\dot{P}_b)_{\text{DR}} = -P_b \cos b \left(\frac{\Theta_0^2}{cR_0}\right) \left(\cos l + \frac{\kappa^2}{\sin^2 l + \kappa^2}\right), \quad (6) $$

$$ (\dot{P}_b)_z = -1.08 \times 10^{-19} \frac{P_b}{c} \left( \frac{1.25z}{[z^2 + 0.0324]^{1/2}} + 0.58z \right) \sin b, \quad (7) $$

$$ (\dot{P}_b)_\mu = \frac{\mu^2 d}{c} P_b, \quad (8) $$

and where: $n_b = 2\pi/P_b$ is the orbital frequency; $m_{\text{tot}} = m_p + m_c$; l and b are the Galactic longitude and latitude, respectively; $\Theta_0$ and $R_0$ are the Galactocentric circular-speed and distance parameters for the Solar-System barycenter, respectively (e.g., Reid et al. 2014); $\kappa = (d/R_0) \cos b - \cos l$; $z = d \sin b$ is the projected

vertical distance of the pulsar-binary system from the Galactic plane; and $\mu = \sqrt{\mu_{\lambda}^2 + \mu_{\beta}^2}$ is the magnitude of proper motion. In this work, we used $\Theta_0 = 236.9(4.2)$ km s⁻¹ and $R_0 = 8.178(26)^7$ kpc as determined by Gravity Collaboration et al. (2019).

In Equation 7 we chose a model for quantifying the Galactic gravitational potential that is traditionally used in pulsar-timing studies (Kuijken & Gilmore 1989). While other models have been developed and studied in recent pulsar literature, the ~20% uncertainty in (Ṗb)obs is too large for resolving statistically meaningful differences between model predictions (Pathak & Bagchi 2018).

While expected to be negligible, we included the (Ṗb)GR term in our analysis for completeness, as it has been observed in other pulsar-binary systems. Hu et al. (2020) presented equations for additional sources of Ṗb corrections that can eventually be observed through pulsar timing, arising from mechanisms such as mass loss

<footnote n="7">⁷ The uncertainty we report for $R_0$ is the composite value determined after adding the statistical and systematic uncertainties reported by Gravity Collaboration et al. (2019) in quadrature.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 3.** A zoomed-in view of the ~1 yr overlap of GBT and CHIME/Pulsar data from Figure 2, sharing the same legend. Worthy of note are: the relative cadences of observation; the smooth, shallow, and linear trend in DM over time; and the statistical weight of CHIME/Pulsar DM measurements as compared to its TOAs, relative to the GBT measurements.

Equations 5-8 demonstrate that $(\dot{P}_b)_{\text{obs}}$ is ultimately a function of $m_p$, $m_c$, and $d$.<fnref n="8" /> To check for self-consistency in our best-fit timing models, we performed a Monte Carlo analysis of the relation $(\dot{P}_b)_{\text{obs}} = (\dot{P}_b)_{\text{GR}} + (\dot{P}_b)_{\text{DR}} + (\dot{P}_b)_z + (\dot{P}_b)_\mu$ by sampling a set of values for ${\{ (\dot{P}_b)_{\text{obs}}, m_p, m_c, R_0, \Theta_0, \mu \}}$, based on uncertainties obtained in this work or as reported in past literature, and solving for $d$ using a Newton-Raphson method. The distribution for $d$ obtained from a Monte Carlo simulation is shown in Figure 4 along with the PDF for $d$ obtained from measurement of $\bar{\omega} = 1/d$. The overlapping PDFs demonstrate that the two independent timing effects yield statistically consistent estimates of $d$. The width of the PDF obtained from $\dot{P}_b$ is dominated by its measurement uncertainty, and width contributions from sampling the ${\{R_0, \Theta_0, \mu\}}$ terms are negligible. We therefore held the ${\{R_0, \Theta_0, \mu\}}$ terms fixed to their best-fit values in subsequent analysis.

The PDF derived from $\bar{\omega}$ was not corrected for the Lutz-Kelker bias in pulsar timing-parallax measurements, which arises in low-precision measurements of $\bar{\omega}$ and can lead to statistical underestimates of $d$ (Verbiest

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

**Figure 4.** Independent estimates of distance from the extended timing data set for PSR J0740+6620. Shown above are distributions of distance to the PSR J0740+6620 binary system determined from the measured timing parallax (line) and Monte Carlo calculations of $(\dot{P}_b)_{\text{obs}}$.

and higher-order corrections predicted by GR. However, these additional terms are at least five orders of magnitude smaller than the current uncertainty in $\dot{P}_b$, and we therefore chose to ignore those terms in subsequent calculations.

<footnote n="8">⁸ $(\dot{P}_b)_{\text{obs}}$ is also a function of several well-measured timing parameters, such as the orbital elements and astrometry. However, we held such values fixed to their best-fit estimates due to their high statistical significance.</footnote>
et al. 2010). However, the consistency between the two distributions of d for PSR J0740+6620 – derived from independent timing effects – indicates that the Lutz-Kelker bias is insignificant in our estimate of $\bar{\omega}$, and we therefore ignore this bias in subsequent analysis.

### 4.3. Likelihood Analysis of the Shapiro Delay and System Distance

For all timing models presented here, we used the $\chi^2$-grid method developed by Splaver et al. (2002) for determining posterior probability density functions (PDFs) and robust credible intervals for {$m_p, m_c, i, \text{and } d$}. The procedure we used for PDF and credible-interval estimation is described as follows:

* selected a set of values for $m_c, \cos i = \sqrt{1 - \sin^2 i}$, and $d$, each from uniform distributions;

* computed values of $\sin i, m_p, \bar{\omega} = 1/d$, and $(\dot{P}_b)_{\text{tot}}$ using Equations 4-8, based on the current location in the $(m_c, \cos i, d)$ phase space;

* held the corresponding values of {$\bar{\omega}, m_c, \sin i, (\dot{P}_b)_{\text{tot}}$} fixed in the timing model;

* re-fitted the timing model, allowing all timing-model parameters not defined on the grid to remain unconstrained during the fit;

* recorded the best-fit $\chi^2$ value, and repeated the above steps for different values of $(m_c, \cos i, d)$.

We ultimately obtained a three-dimensional grid of $\chi^2$ values computed over uniform steps in $(m_c, \cos i, d)$, which were then mapped to a likelihood function $p(\text{data}, \text{DMX}_j | m_c, \cos i, d) \propto \exp(-\Delta\chi^2/2)$ where $\Delta\chi^2 = \chi^2 - \min(\chi^2)$. The notation “DMX_j” refers to one of the three DMX models, labeled with subscript j, that was generated for this work. We then used Bayes’ theorem to compute the posterior distribution $p(m_c, \cos i, d|\text{data}, \text{DMX}_j)$ for uniform priors on the physical parameters and choice of DMX model.

Systematic uncertainties in the $(m_p, m_c, i, d)$ parameters may arise due to choices in modeling DM variations. We used Bayesian model averaging (e.g., Hoeting et al. 1999) to obtain PDFs and credible intervals that better reflect the model-independent distributions, in order to address DMX-model uncertainty. The Bayesian model averaging method defines the model-independent PDF as a weighted summation of model-dependent PDFs; in the case of the three-dimensional posterior distribution described above, the model-averaged PDF is

<figure><img src="image_6.png" /><figcaption>Figure 5. A comparison of normalized posterior PDFs of $m_p$ estimated from the $\chi^2$-grid method applied to timing models of different DMX binning choices (dashed lines). For reference, the resultant posterior PDF obtained using the Bayesian model averaging method described in Section 4.3 is shown as the red-solid line, and largely overlaps the orange-dashed curve.</figcaption></figure>

$$p(m_c, \cos i, d|\text{data}) = \sum_j p(m_c, \cos i, d|\text{data}, \text{DMX}_j) \times p(\text{DMX}_j|\text{data}) \quad (9)$$

where $p(\text{DMX}_j|\text{data})$ is the conditional posterior PDF for model DMX_j. In the absence of any preference in DMX modeling, all DMX models are equally likely and Equation 9 reduces to a straightforward averaging of the three normalized PDFs. We used Equation 9 to compute model-averaged posterior PDFs and their corresponding credible intervals for all four gridded parameters. The DMX_j and model-averaged posterior PDFs for $m_p$ are shown in Figure 5.

Posterior PDFs and credible intervals for $m_p$ were derived from the original three-dimensional posterior PDF by noting that

$$p(m_p, \cos i, d|\text{data}) = p(m_c, \cos i, d|\text{data}) \left| \frac{\partial m_c}{\partial m_p} \right|, \quad (10)$$

$$\text{where } \frac{\partial m_c}{\partial m_p} = \frac{2f_m m_{\text{tot}}}{3m_c^2 \sin^3 i - 2f_m m_{\text{tot}}}, \quad (11)$$

and where Equation 11 is determined by the Keplerian mass function ($f_m$). Using these original and translated PDFs, we computed credible intervals by marginalizing over the relevant subset of parameters, and also calculated two-dimensional posterior PDFs for all possible parameter pairs. We ignored the uncertainty in $f_m$ for all PDF calculations due to its negligible contribution
to PDF widths (relative uncertainty ~10⁻⁷; Cromartie et al. 2020). Table 3 lists the credible intervals obtained using the methods described above.

## 5. DISCUSSION

The combined timing data set for PSR J0740+6620 contains new information on temporal and DM variations, as well as improved constraints of effects previously reported by Cromartie et al. (2020). We discuss various aspects of these features in detail below, using data products generated with the techniques outlined in Section 4. These data products are publicly available and include an initial set of posterior-PDF estimates used by NICER to inform their modeling of surface thermal emission with our improved estimates for the PSR J0740+6620 system parameters (Fonseca et al. 2021). Our updated constraints on mass, geometric and distance parameters were used by NICER to constrain the radius of PSR J0740+6620; their application is described elsewhere (Miller et al. 2021; Riley et al. 2021; Raaijmakers et al. 2021).

### 5.1. *Timing Properties in the CHIME/Pulsar Era*

All best-fit timing models developed for this work yielded consistent weighted root mean square (RMS) residuals: ~0.3 µs for the GBT/1400-MHz data; ~0.4 µs for GBT/820-MHz; and ~1 µs for CHIME/Pulsar data. The factor of ~2 difference in RMS residual between CHIME/Pulsar and GBT data is consistent with the expected timing precision based on observed sensitivity of the CHIME/Pulsar system; a recent timing analysis of PSR J0645+5158, another MSP observed with the GBT and CHIME/Pulsar, yielded similar RMS statistics (CHIME/Pulsar Collaboration et al. 2021).

The white-noise model of the combined TOA set, determined using the methods outlined in Section 3.4, is consistent with the model developed by Cromartie et al. (2020). Moreover, the CHIME/Pulsar data set yields TOA-uncertainty scale and quadrature values comparable to those obtained for the NANOGgrav data set. Therefore, while slightly less precise than its NANOGgrav counterpart, the CHIME/Pulsar instrument is producing TOAs with noise properties that are consistent with behavior observed when using other observatories.

The GBT timing data currently possess greater statistical weight on parameter constraints due to their higher timing precision and larger number of TOAs. It is worth noting that a timing analysis of the CHIME/Pulsar data set on its own yields robust (albeit weaker) measurement of the Shapiro delay, with $m_c = 0.28(2)$ M☉ and $\sin i = 0.9989(8)$ that correspond to $m_p = 2.4(3)$

M☉. These CHIME/Pulsar estimates are consistent with the combined-set values obtained using the methods discussed in Section 4. The high-cadence nature of CHIME/Pulsar observations, along with the 4.8-day orbit of PSR J0740+6620, has led to a statistically significant constraint on the Shapiro delay with only ~1 year of timing data, considerably faster than the rate that was achieved with the GBT. Nonetheless, a combination of high cadence and high sensitivity is the only way to meaningfully improve the mass and geometric estimates of the PSR J0740+6620 system. We discuss these prospects further below.

### 5.2. *DM Variations Toward PSR J0740+6620*

As shown in Figure 2, the DM of PSR J0740+6620 varies smoothly and slowly across the full data set, with few outlying points and no sign of annual variations due to interactions with free electrons of the Solar wind. These features are consistent with those reported by Donner et al. (2020), who analyzed a 4.8-yr data set of PSR J0740+6620 acquired using the German LOng Wavelength (GLOW) consortium of telescope stations built for the LOw Frequency ARray (LOFAR), which were all sensitive around ~150 MHz. The lack of periodic solar-wind variations is expected given the high ecliptic latitude of the PSR J0740+6620 system ($\beta \approx 44^{\circ}$, see Table 2), which leads to negligible traversal of the pulsar signal through the circumsolar medium.

The differences in frequency coverage and timing precision between the GBT and CHIME/Pulsar data sets lead to observable differences in their DM measurements. In particular, while the GBT data yield superior RMS timing residuals in both receiver bands than those obtained with CHIME/Pulsar, the CHIME/Pulsar data nonetheless yield better precision in DM measurements; the median uncertainty in GBT DMs determined with *PulsePortraiture* is ~8 × 10⁻⁴ pc cm⁻³, while the same measure in the CHIME/Pulsar data set is ~2 × 10⁻⁴ pc cm⁻³. The median DMX measurement uncertainty in the CHIME/Pulsar era has comparable precision to that obtained with GLOW by Donner et al. (2020), ~4 × 10⁻⁵ pc cm⁻³. Further analysis of the DM timeseries for PSR J0740+6620 will be the subject of future works.

### 5.3. *Updated Estimates of the Shapiro Delay Parameters*

We performed the χ²-grid and PDF analyses outlined in Section 4.3 for all three DMX timing solutions derived from the combined NANOGgrav+CHIME/Pulsar data set for PSR J0740+6620. A summary of the credible intervals derived from all posterior PDFs is presented
<table><caption>**Table 3.** Credible intervals for the Shapiro-delay and distance parameters in the J0740+6620 system.</caption><thead><tr><th>Parameter</th><th>Cromartie et al. (2020)<sup>a</sup></th><th>DMX=3.0</th><th>DMX=6.5</th><th>Hybrid-DMX</th><th>Model-Averaged</th></tr></thead><tbody><tr><td>Pulsar mass, m<sub>p</sub> (M<sub>☉</sub>)</td><td>2.14<sup>+0.10</sup><sub>-0.09</sub></td><td>2.08<sup>+0.07</sup><sub>-0.07</sub></td><td>2.07<sup>+0.07</sup><sub>-0.06</sub></td><td>2.10<sup>+0.07</sup><sub>-0.07</sub></td><td>2.08<sup>+0.07</sup><sub>-0.07</sub></td></tr><tr><td>Companion mass, m<sub>c</sub> (M<sub>☉</sub>)</td><td>0.260<sup>+0.008</sup><sub>-0.007</sub></td><td>0.253<sup>+0.006</sup><sub>-0.005</sub></td><td>0.252<sup>+0.005</sup><sub>-0.005</sub></td><td>0.254<sup>+0.006</sup><sub>-0.006</sub></td><td>0.253<sup>+0.006</sup><sub>-0.005</sub></td></tr><tr><td>System inclination, i (degrees)</td><td>87.38<sup>+0.2</sup><sub>-0.2</sub></td><td>87.53<sup>+0.17</sup><sub>-0.18</sub></td><td>87.61<sup>+0.16</sup><sub>-0.16</sub></td><td>87.53<sup>+0.17</sup><sub>-0.18</sub></td><td>87.56<sup>+0.17</sup><sub>-0.18</sub></td></tr><tr><td>Distance, d (kpc)</td><td>&lt; 4.4<sup>b</sup></td><td>1.16<sup>+0.17</sup><sub>-0.15</sub></td><td>1.11<sup>+0.17</sup><sub>-0.15</sub></td><td>1.15<sup>+0.17</sup><sub>-0.15</sub></td><td>1.14<sup>+0.17</sup><sub>-0.15</sub></td></tr></tbody></table>

<sup>a</sup>Values initially determined by Cromartie et al. (2020) are listed for comparison. The $\chi^2$ gridding procedure used by Cromartie et al. (2020) only considered one DMX model, with maximum bin size of 6.5 days.

<sup>b</sup>Neither $\dot{P}_b$ nor timing parallax – and thus the distance – were significantly constrained with the data set presented by Cromartie et al. (2020). We instead quote here the $2\sigma$ upper limit obtained from their reported constraint.

in Table 3. The model-averaged PDFs and credible intervals are shown in Figure 6.

For all four parameters defined on each $\chi^2$ grid, the three posterior PDFs are consistent with each other at 68.3% credibility as well as with the estimates made by Cromartie et al. (2020). Moreover, the extent of each credible interval is largely unaffected by the choice of DMX bin width. This consistency indicates that the choice of DMX model with different bin widths makes no statistical difference to the significance of our estimates in the PSR J0740+6620 system. However, it is also likely that this consistency is specific to PSR J0740+6620 due to its slowly-varying DM evolution being adequately modeled with coarser piecewise-constant DMX models.

Our 68.3% credible interval of the model-averaged $m_p = 2.08^{+0.07}_{-0.07}$ M$_{☉}$ is reduced by ~30% in comparison to the estimate made by Cromartie et al. (2020). The model-averaged $m_p$ remains the largest of all other precisely measured pulsar masses determined with the Shapiro delay to date. The 95.4% lower bound on the model-averaged $m_p$, 1.95 M$_{☉}$, is similar to the 1.96 M$_{☉}$ lower bound initially determined by Cromartie et al. (2020). The lower limit on the maximum mass of neutron stars from the PSR J0740+6620 system therefore remains unchanged in our analysis.

The model-averaged estimate of the companion mass, $m_c = 0.253^{+0.006}_{-0.005}$ M$_{☉}$, remains largely consistent with the prediction from expected correlations between companion masses and orbital sizes that arise due to extended periods of mass transfer (Tauris & Savonije 1999). The model-averaged credible interval on $m_c$ is in improved agreement with correlation parameters that define the $m_c - P_b$ relation for low-metallicity progenitors (i.e., metallic mass fraction $Z \sim 10^{-3}$ or lower), which was noted by Cromartie et al. (2020) to yield an expected $m_c \sim 0.25$ M$_{☉}$. Echeveste et al. (2020) demonstrated through numerical calculations that only

metal-poor progenitors with helium interiors can donate matter and evolve to yield the current masses of the PSR J0740+6620 system, and our improved measurement of $m_c$ further supports this conclusion.

## 5.4. Distance to the PSR J0740+6620 System

We obtained a constrained, model-averaged distance of $d = 1.14^{+0.17}_{-0.15}$ kpc from the combined NANOGrav and CHIME/Pulsar TOA data set when using the $\chi^2$-grid method described in Section 4.3. This updated distance estimate is consistent with the distance of $d \approx 0.9$ kpc derived from the observed mean DM, placement within the Milky Way Galaxy, and the model of Galactic electron number density developed by Yao et al. (2017); the distance estimated by the number-density model of Cordes & Lazio (2002), $d \approx 0.6$ kpc, is less consistent, though underlying systematic uncertainties for both electron-density models correspond to ~30% on $d$ and thus reduce the tension.

Our model-averaged estimate of $d$ is statistically consistent at the $2\sigma$ level with the marginal estimate of $d$ derived from $\bar{\omega} = 0.5(3)$ mas made by Cromartie et al. (2020). However, these two estimates are in tension with the first estimate of $d = 0.4^{+0.2}_{-0.1}$ kpc made from initial NANOGrav timing of PSR J0740+6620 (derived from $\bar{\omega} = 2.3(6)$; Arzoumanian et al. 2018). An important distinction in these estimates is that the timing solution presented by Arzoumanian et al. (2018) yielded only a marginal detection of the Shapiro delay parameters. Arzoumanian et al. (2018) also employed an “approximate” orthometric parameterization of the Shapiro delay to model the relativistic effect as a Fourier expansion about the orbital period using a finite number of harmonic terms. As noted by Freire & Wex (2010), prominent Shapiro delays from highly inclined binary systems – as was first established for PSR J0740+6620 by Cromartie et al. (2020) – are best modeled using the analytically exact expression predicted by general relativity,
<figure><img src="image_7.png" /></figure>

**Figure 6.** A triangle plot of two-dimensional, model-averaged posterior PDFs (off-diagonal) and marginalized, one-dimensional PDFs (on-diagonal) of parameters that characterize the Shapiro delay and distance to the PSR J0740+6620 system, derived using the methods outlined in Section 4.3. Vertical red lines in each diagonal panel represent the median value (solid) and 68.3% credible interval (dashed). Inner- to outer-most contours denote the 68.3%, 95.4% and 99.7% regions of credibility, respectively.
<figure><img src="image_8.png" /></figure>

**Figure 7.** Color–magnitude diagram with white dwarf evolutionary sequences, evaluated at the *r'* and *i'* photometric bands. Black and green curves demonstrate the cooling predictions for white dwarfs with hydrogen and helium atmospheres, respectively (Holberg & Bergeron 2006; Kowalski & Saumon 2006; Tremblay et al. 2011; Bergeron et al. 2011). For each colored set of curves, different tracks represent masses 0.2, 0.3, 0.4, 0.6, 0.8 and 1.0 $M_{\odot}$, with the mass increasing from upper to lower curves. Cooling ages along each track are marked by different symbols. The position of the companion to PSR J0740+6620 is indicated by the crosses: the dashed lower and middle crosses correspond to the latest DM distance estimate *d* = 0.93 kpc and the initial timing parallax distance *d* = 0.4$^{+0.2}_{-0.1}$ kpc from Arzoumanian et al. (2018), and the red cross demonstrates the companion location based on the new *d* = 1.14$^{+0.17}_{-0.15}$ kpc. The blue-shaded region shows the initial magnitude-distance uncertainty range as presented by Beronya et al. (2019), and the green-shaded region indicates the new uncertainty range corresponding to *d* = 1.14$^{+0.17}_{-0.15}$ kpc derived in this work. The two regions slightly overlap, and the green region covers the blue one.

instead of a finite-term Fourier expansion that is more
appropriate for low-inclination systems. The combina-
tion of a sparse, low-cadence data set and sub-optimal
modeling likely led to an inaccurate distance estimate
determined for PSR J0740+6620 by Arzoumanian et al.
(2018).

Finally, we reassessed the optical properties of the
white-dwarf companion to PSR J0740+6620 using our
direct measurement of *d*. Beronya et al. (2019) con-
cluded, based on their derived magnitudes and colors of
the optical counterpart and the range of contemporane-
ous timing parallaxes and DM-based distances, that the

pulsar companion is an ultra-cool, helium-atmosphere
white dwarf with an effective temperature < 3,500 K
and cooling age >5 Gyr. However, Beronya et al. (2019)
were not able to constrain the mass of the companion
due to the significant distance ambiguity. Observations
with the Gaia satellite, designed for astrometric mea-
surements, are currently not possible as the *r'* and *i'*
apparent magnitudes of the white-dwarf companion to
PSR J0740+6620 are ~26 (Beronya et al. 2019). These
observed magnitudes are well above the limiting *g'* mag-
nitude of ~19 for Gaia Data Release 2.<fnref n="9" />

Figure 7 shows a comparison of the absolute magni-
tude and color of the optical source with cooling predic-
tions for white dwarfs with hydrogen and helium atmo-
spheres, and over a range of masses<fnref n="10" /> (Holberg & Berg-
eron 2006; Kowalski & Saumon 2006; Tremblay et al.
2011; Bergeron et al. 2011). The range of prior con-
straints on *d* led to a large range in absolute magnitude,
shown as the blue-shaded region in Figure 7, that encom-
passes cooling tracks of white dwarfs with 0.2 − 1 $M_{\odot}$.
Our updated measurement of *d* restricts the absolute
magnitude to the green-shaded region in Figure 7 where
only cooling curves for helium-atmosphere white dwarfs
with masses ~0.2 − 0.3 $M_{\odot}$ remain consistent with ob-
servations. This result is in agreement with the initial
estimate of $m_c$ made by Cromartie et al. (2020), and
with our updated $m_c = 0.253_{-0.005}^{+0.006} M_{\odot}$. The improved
synergy in optical- and radio-based estimates of intrin-
sic parameters for the PSR J0740+6620 system favors
our new measurement of *d*, and further strengthens the
counterpart association made by Beronya et al. (2019).

Future optical and infrared studies with next-
generation telescopes, combined with improved distance
uncertainties from radio-timing observations, will be
useful for further binary evolution modelling (Echeveste
et al. 2020). The increasing timespan of the PSR
J0740+6620 data set will further improve the measure-
ment of $\dot{P}_b$ and thus yield a precise kinematic distance,
as was recently obtained for the PSR J1909−3744 binary
system (Liu et al. 2020). Cooling predictions for white
dwarfs with different composition, as well as stronger
constraints on the distance and optical magnitudes, will
allow for constraints on other properties of the white-
dwarf companion, such as its hydrogen abundance.

6. CONCLUSIONS

In this work, we extended the timing data set on the
high-mass PSR J0740+6620 to incorporate additional
GBT data acquired as part of the NANOGrav project,

<footnote n="9">⁹ https://www.cosmos.esa.int/web/gaia/dr2</footnote>

<footnote n="10">¹⁰ http://www.astro.umontreal.ca/~bergeron/CoolingModels</footnote>
and to include the initial data set being generated by CHIME/Pulsar. We generated timing solutions for PSR J0740+6620 that were derived from both narrowband and wideband arrival times collected with the GBT and CHIME/Pulsar. For each data set, we also explored using three different piecewise-constant models of DM variations to assess sensitivity of the improved mass and geometric measurements on such modeling choices.

We found that all solutions were statistically consistent with one another at 68.3% credibility, with minimal (i.e., ~1%) variation in credible intervals of the Shapiro delay parameters when comparing between several DM-evolution models. We obtained a model-averaged estimate of $m_p = 2.08_{-0.07}^{+0.07}$ M$_\odot$, $m_c = 0.253_{-0.005}^{+0.006}$ M$_\odot$ and $i = 87.56_{-0.18}^{+0.17^\circ}$, consistent with the estimates made by Cromartie et al. (2020). Our constraint on $m_p$ is improved by ~30%, though the lower limit we derived for the maximum mass of neutron stars remains unchanged from that determined by Cromartie et al. (2020). For the first time, our analysis also yielded a significant measurement of $\dot{P}_b$ which we argued is due to apparent acceleration from transverse motion. We found that the model-averaged $d = 1.14_{-0.15}^{+0.17}$ kpc when using the observed $\dot{P}_b$ as a constraint on the system distance, using the framework described in Section 4.

We performed an independent timing analysis of TOA data recorded with the Nançay Radio Telescope (NRT) over a 7-yr timespan. All timing parameters presented in our work were statistically consistent with those estimated from NRT arrival times. The NRT data, along with forthcoming observations with the GBT and CHIME/Pulsar, will be the subject of future studies.

There are four ways in which our current estimates of the masses, geometry and distance can be significantly improved: continued, high-cadence timing with existing facilities; observations with “ultra-wideband” radio receivers, which will become possible at the GBT in the next few years; observations with forthcoming telescopes of greater sensitivity and broadband coverage; and measurement of additional orbital variations that can be related to one or more of the fundamental parameters considered in this work. The most promising avenue for meaningful improvement on current constraints is the use of planned next-generation facilities. Several future radio observatories – such as the Canadian Hydrogen Observatory and Radio-transient Detector (CHORD; Vanderlinde et al. 2019), 2000-dish Deep Synoptic Array (DSA-2000; Hallinan et al. 2019), and the next-generation Very Large Array (ngVLA; e.g., Chatterjee 2018) – will provide key opportunities in pulsar science with their heightened sensitivity and increased receiver bandwidths.

For pulsar timing, the radiometer equation dictates the TOA uncertainty scales as $\sigma_{\text{TOA}} \propto S_{\text{sys}}/\sqrt{t_{\text{obs}}\Delta f}$, where: $S_{\text{sys}}$ is the system-equivalent flux density of the observatory; $t_{\text{obs}}$ is the observation timespan, and $\Delta f$ is the receiver bandwidth.<fnref n="11" /> The expected improvement of TOA precision between CHIME/Pulsar (with $S_{\text{sys}} \approx 50$ Jy; CHIME/Pulsar Collaboration et al. 2021) and CHORD observations can be estimated based on current projections of CHORD design specification, with $S_{\text{sys}} = 9$ Jy and $\Delta f = 1200$ MHz (Vanderlinde et al. 2019). For observations of the same $t_{\text{obs}}$, the band/time-averaged TOA uncertainty obtained with CHORD is $\sigma_{\text{TOA}} \sim 0.1$ µs when using CHIME/Pulsar values listed in Table 1. The DSA-2000 and ngVLA have similar planned bandwidths, but a factor of ~5 smaller $S_{\text{sys}}$ than that of CHORD. With CHORD projected to be operational by 2025, it is likely that daily TOA uncertainties of 0.1 µs and lower will be regularly achieved with next-generation observatories as soon as several years from now.

A reduction in $\sigma_{\text{TOA}}$ by a factor of 10 or greater will lead to a similarly large improvement in the parameters of the Shapiro time delay. We therefore expect the credible interval of $m_p$ to exceed values no larger than 0.01 M$_\odot$ of the median value with sufficient observations from CHORD, the DSA-2000 and/or ngVLA. The significance of $\dot{P}_b$, and thus the constraint on distance, will continue to improve such that its relative uncertainty scales as $T^{-5/2}$, where T is the data timespan (Damour & Taylor 1992). With ongoing observations and future prospects, we expect to gain additional, important insight from the PSR J0740+6620 system in the near future.

## ACKNOWLEDGMENTS

We acknowledge that CHIME is located on the traditional, ancestral, and unceded territory of the Syilx/Okanagan people. We thank Anne M. Archibald for useful discussion. We also thank the anonymous reference for their careful review of our work.

The NANOGrav project receives support from the National Science Foundation (NSF) Physics Frontiers Center (PFC) award number 1430284. The GBT is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio As-

<footnote n="11">¹¹ For this calculation, we assume that intrinsic variations of pulse width and brightness across frequency are negligible. Our analysis of the combined NANOGrav and CHIME/Pulsar data sets, which collectively span the projected bandwidth of CHORD, largely support this assumption.</footnote>
tronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. We thank the telescope operators and all the staff at the GBT for the essential role they played in collecting data that were used in this work.

We are grateful to the staff of the Dominion Radio Astrophysical Observatory, which is operated by the National Research Council of Canada. CHIME is funded by a grant from the Canada Foundation for Innovation (CFI) 2012 Leading Edge Fund (Project 31170) and by contributions from the provinces of British Columbia, Québec and Ontario. The CHIME Fast Radio Burst Project, which enabled development in common with the CHIME/Pulsar instrument, is funded by a grant from the CFI 2015 Innovation Fund (Project 33213) and by contributions from the provinces of British Columbia and Québec, and by the Dunlap Institute for Astronomy and Astrophysics at the University of Toronto. Additional support is provided by the Canadian Institute for Advanced Research (CIFAR), McGill University and the McGill Space Institute thanks to the Trottier Family Foundation, and the University of British Columbia. The CHIME/Pulsar instrument hardware is funded by the Natural Sciences and Engineering Research Council (NSERC) Research Tools and Instruments (RTI-1) grant EQPEQ 458893-2014.

Support for HTC is provided by NASA through the NASA Hubble Fellowship Program grant \#HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. TTP is a NANOGray PFC Postdoctoral Fellow and acknowledges support

from the MTA-ELTE Extragalactic Astrophysics Research Group, funded by the Hungarian Academy of Sciences (Magyar Tudományos Akadémia), that was used during the development of this research. Portions of this work performed at NRL is supported by ONR 6.1 basic research funding. SMR is a CIFAR Fellow. Pulsar research at UBC is supported by an NSERC Discovery Grant and by CIFAR. TD and MTL are supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the “Programme National Gravitation, Références, Astronomie, Métrologie” (PNGRAM) of L’institut national des sciences de l’Univers (INSU) du CNRS, France. ECF is supported by NASA under award number 80GSFC17M0002. WF is supported by the STEM Mountains of Excellence graduate fellowship. DCG is supported by the John I. Watters research fellowship. VMK holds the Lorne Trottier Chair in Astrophysics & Cosmology, a Distinguished James McGill Professorship and receives support from an NSERC Discovery Grant (RGPIN 228738-13) and Gerhard Herzberg Award, from an R. Howard Webster Foundation Fellowship from CIFAR, and from the Fonds de Recherche du Québec: Nature et technologies, Centre de Recherche en Astrophysique du Québec. JWM is a CITA Postdoctoral Fellow, and this work was supported by NSERC (funding reference \#CITA 490888-16).

Software: `matplotlib` (Hunter 2007), `nanopipe` (Demorest 2018), `psrarchive` (Hotan et al. 2004), `PulsePortraiture` (Pennucci et al. 2016), `tempo` (Nice et al. 2015),

REFERENCES

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101

Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJL, 896, L44, doi: 10.3847/2041-8213/ab960f

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

—. 2021b, ApJS, 252, 5

Antoniadis, J., Tauris, T. M., Özel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665

Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37, doi: 10.3847/1538-4365/aab5b0

Backer, D. C., & Hellings, R. W. 1986, Ann. Rev. Astr. Ap., 24, 537

Bergeron, P., Wesemael, F., Dufour, P., et al. 2011, ApJ, 737, 28

Beronya, D. M., Karpova, A. V., Kirichenko, A. Y., et al. 2019, MNRAS, 485, 3715

Bogdanov, S., Guillot, S., Ray, P. S., et al. 2019, ApJL, 887, L25

Chatterjee, S. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 751

CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5
Coles, W., Hobbs, G., Champion, D. J., Manchester, R. N.,
& Verbiest, J. P. W. 2011, MNRAS, 418, 561

Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints,
astro. https://arxiv.org/abs/astro-ph/0207156

Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020,
Nat. Astron., 4, 72

Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840

Demorest, P. B. 2018, nanopipe: Calibration and data
reduction pipeline for pulsar timing.
http://ascl.net/1803.004

Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts,
M. S. E., & Hessels, J. W. T. 2010, Nat., 467, 1081

Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2020,
A&A, 644, A153

DuPlain, R., Ransom, S., Demorest, P., et al. 2008, in
Society of Photo-Optical Instrumentation Engineers
(SPIE) Conference Series, Vol. 7019, Advanced Software
and Control for Astronomy II, ed. A. Bridger & N. M.
Radziwill, 70191D

Echeveste, M., Novarino, M. L., Benvenuto, O. G., & De
Vito, M. A. 2020, MNRAS, 495, 2509

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

Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021,
Refined Mass and Geometric Measurements of the
High-Mass PSR J0740+6620: Probability Density
Functions and their Credible Intervals, Zenodo,
doi: 10.5281/zenodo.4773599

Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199

Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019,
A&A, 625, L10

Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, in Bulletin
of the American Astronomical Society, Vol. 51, 255

Hoeting, J. A., D., M., & Raftery, A. E. Volinsky, C. T.
1999, Stat. Sci., 14, 382

Holberg, J. B., & Bergeron, P. 2006, AJ, 132, 1221

Hotan, A. W., van Straten, W., & Manchester, R. N. 2004,
PASA, 21, 302

Hu, H., Kramer, M., Wex, N., Champion, D. J., & Kehl,
M. S. 2020, MNRAS, 497, 3118

Hunter, J. D. 2007, Computing in Science & Engineering, 9,
90

Kowalski, P. M., & Saumon, D. 2006, ApJL, 651, L137

Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571

Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2017, ApJ,
834, 35

Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326,
274

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

Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54

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

Liu, K., Guillemot, L., Istrate, A. G., et al. 2020, MNRAS,
499, 2276

Lynch, R. S., Swiggum, J. K., Kondratiev, V. I., et al.
2018, ApJ, 859, 93

Miller, M. C., Lamb, F. K., Dittmann, A. J., & et al. 2021,
arXiv e-prints, arXiv:2105.06979

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019,
ApJL, 887, L24

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

Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo:
Pulsar timing data analysis. http://ascl.net/1509.002

Nice, D. J., & Taylor, J. H. 1995, ApJ, 441, 429

Özel, F., & Freire, P. 2016, Ann. Rev. Astr. Ap., 54, 401

Pathak, D., & Bagchi, M. 2018, ApJ, 868, 123

Pennucci, T. T. 2015, PhD thesis, University of Virginia
—. 2019, ApJ, 871, 34

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

—. 2016, Pulse Portraiture: Pulsar timing.
http://ascl.net/1606.013

Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019,
ApJL, 887, L22

Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020,
ApJL, 893, L21

Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021,
arXiv e-prints, arXiv:2105.06981

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014,
ApJ, 783, 130

Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL,
887, L21

Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv
e-prints, arXiv:2105.06980

Romani, R. W., Kandel, D., Filippenko, A. V., & et al.
2021, ApJL, 908, L46

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

Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789

Shklovskii, I. S. 1970, Sov. Astron., 13, 562

Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002,
ApJ, 581, 509

Stairs, I. H. 2003, Living Reviews in Relativity, 6, 5

Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928

Taylor, J. H. 1992, Philosophical Transactions of the Royal
Society of London Series A, 341, 117
Tremblay, P.-E., Bergeron, P., & Gianninas, A. 2011, ApJ, 730, 128

Tsokaros, A., Ruiz, M., & Shapiro, S. L. 2020, ApJ, 905, 48

Vallisneri, M., & van Haasteren, R. 2017, MNRAS, 466, 4954

van Kerkwijk, M. H., Breton, R. P., & Kulkarni, S. R. 2011, ApJ, 728, 95

van Straten, W., & Bailes, M. 2011, PASA, 28, 1

Vanderlinde, K., Liu, A., Gaensler, B., et al. 2019, in Canadian Long Range Plan for Astronony and Astrophysics White Papers, Vol. 2020, 28.

https://arxiv.org/abs/1911.01777

Verbiest, J. P. W., Lorimer, D. R., & McLaughlin, M. A. 2010, MNRAS, 405, 564

Wolff, M. T., Guillot, S., Bogdanov, S., et al. 2021, submitted to ApJL

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29