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

E. Fonseca 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  H. T. Cromartie Cornell Center for Astrophysics and
Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY
14853, USA  NASA Hubble Fellowship Program Einstein Postdoctoral Fellow  T. T.
Pennucci 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  P. S. Ray U.S. Naval Research
Laboratory, Washington, DC 20375, USA  A. Yu. Kirichenko Universidad Nacional
Autónoma de México, Instituto de Astronomía, AP 106, Ensenada 22800, BC,
México  Ioffe Institute, Politekhnicheskaya 26, St. Petersburg, 194021, Russia
S. M. Ransom National Radio Astronomy Observatory, 520 Edgemont Road,
Charlottesville, VA 22903, USA  P. B. Demorest  National Radio Astronomy
Observatory, P.O. Box O, Socorro, NM 87801, USA  I. H. Stairs Department of
Physics & Astronomy, University of British Columbia, 6224 Agricultural Road,
Vancouver, BC V6T 1Z1, Canada  Z. Arzoumanian  Astrophysics Science Division,
Code 662, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  L.
Guillemot 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  A. Parthasarathy Max Planck Institute for Radio Astronomy, Auf dem H
”u gel 69, D-53121 Bonn, Germany  M. Kerr U.S. Naval Research Laboratory,
Washington, DC 20375, USA  I. Cognard  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  P. T. Baker Department of Physics and
Astronomy, Widener University, Chester, PA 19013, USA  H. Blumer 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  P. R. Brook 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  M. DeCesar Department of Physics and Astronomy,
George Mason University, Fairfax, VA 22030, resident at U.S. Naval Research
Laboratory, Washington, D.C. 20375, USA  T. Dolch 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
F. A. Dong Department of Physics & Astronomy, University of British Columbia,
6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada  E. C. Ferrara
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  W. Fiore 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  N. Garver-Daniels  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  D. C. Good Department of Physics & Astronomy,
University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1,
Canada  R. Jennings Department of Astronomy, Cornell University, Ithaca, NY
14853, USA  M. L. Jones Center for Gravitation, Cosmology, and Astrophysics,
Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413,
Milwaukee, WI 53201, USA  V. M. Kaspi 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  M. T. Lam School of Physics and Astronomy, Rochester Institute of
Technology, Rochester, NY 14623, USA  Laboratory for Multiwavelength
Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA  D.
R. Lorimer 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  J. Luo Canadian
Institute for Theoretical Astrophysics, University of Toronto, 60 St. George
Street, Toronto, ON M5S 3H8, Canada  A. McEwen Center for Gravitation,
Cosmology, and Astrophysics, Department of Physics, University of Wisconsin-
Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA  J. W. McKee Canadian
Institute for Theoretical Astrophysics, University of Toronto, 60 St. George
Street, Toronto, ON M5S 3H8, Canada  M. A. McLaughlin  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  N. McMann Department of Physics and
Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235,
USA  B. W. Meyers Department of Physics & Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada  A. Naidu
University of Oxford, Sub-Department of Astrophysics, Denys Wilkinson
Building, Keble Road, Oxford, OX1 3RH, United Kingdom  C. Ng Dunlap Institute
for Astronomy & Astrophysics, University of Toronto, 50 St. George Street,
Toronto, ON M5S 3H4, Canada  D. J. Nice Department of Physics, Lafayette
College, Easton, PA 18042, USA  N. Pol Department of Physics and Astronomy,
Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA  H. A.
Radovan  Department of Physics, University of Puerto Rico, Mayagüez, PR 00681,
USA  B. Shapiro-Albert 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  C. M.
Tan 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  S. P. Tendulkar 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  J. K. Swiggum Department of Physics,
Lafayette College, Easton, PA 18042, USA  H. M. Wahl 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  W. W. Zhu National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100101, China

(Received XXX; Revised YYY; Accepted ZZZ)

###### 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
$\sim 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_{\rm
p}=2.08^{+0.07}_{-0.07}$ M⊙ (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.17}_{-0.15}$
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_{\rm p}$ for J0740+6620 by an
order of magnitude within the next few years.

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

††journal: ApJ Letters††software:  matplotlib (Hunter, 2007), nanopipe
(Demorest, 2018), psrchive (Hotan et al., 2004), PulsePortraiture (Pennucci et
al., 2016), tempo (Nice et al., 2015),

##  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, propelling 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_{\rm p}=2.14^{+0.10}_{-0.09}$ M⊙ (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
thermal 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.

Table 1: A summary of observing parameters and data sets of PSR J0740+6620 anaylzed in this work. Parameter | GBT/1400-MHz | GBT/820-MHz | CHIME
---|---|---|---
Time Range (MJD) | 56675–58945 | 56640–58944 | 58517–58975
Frequency Range (MHz) | 1151–1885 | 722–919 | 400–800
Number of Channels, raw | 512 | 128 | 1024
Number of Channels, downsampled | 64 | 64 | 32
Number of Observing Epochs | 187 | 142 | 263
Number of TOAs, narrowband | 8,208 | 5,191 | 4,862
Number of TOAs, wideband | 209 | 154 | 263
Min./Median/Max. $\sigma_{\rm TOA}$, narrowband ($\mu$s) | 0.24/2.65/32.7 | 0.16/3.16/24.9 | 0.49/5.21/62.8
Min./Median/Max. $\sigma_{\rm TOA}$, wideband ($\mu$s) | 0.05/0.39/2.04 | 0.04/0.51/3.10 | 0.14/0.97/3.28

###  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 (NANOGrav111nanograv.org), 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-$\mu$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 $\sim 33$-minute duration
in order to record a full transit at the CHIME telescope.

###  2.2 Observing Cadences

Under the nominal NANOGrav program, multi-frequency observations of PSR
J0740+6620 with the GBT typically lasted for $\sim$25 minutes per receiver and
occurred with a monthly cadence; observations with each receiver were acquired
within $\sim$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 NANOGrav 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 NANOGrav 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
$\sim$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 psrchive222http://psrchive.sourceforge.net/ suite of analysis utilities,
through the use of separate psrchive-based procedures tailored to the NANOGrav
and CHIME/Pulsar data sets. For NANOGrav data on PSR J0740+6620, we used the
nanopipe pipeline333https://github.com/demorest/nanopipe – which is employed
in ongoing NANOGrav analysis of all its sources – for the following
manipulations: mask-based excision of radio frequency interference (RFI);
calibration of flux density and polarization information 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)
$\sim\;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 over-weighting
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_{\rm 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., NANOGrav Collaboration et al., 2015). These techniques of template and
TOA generation form a standard 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 narrowband TOAs for each receiver
bandwidth and all processed, RFI-cleaned observations. Narrowband TOAs derived
from GBT data were further cleaned by removing outlier timing data using the
NANOGrav analysis procedure outlined by Alam et al. (2021a), where TOAs were
discarded from analysis based on: low profile fidelity (i.e., with a pulse
profile S/N $<$ 8); corrupted calibration data; a probabilistic classification
of outlier TOAs based on prior timing models (Vallisneri & van Haasteren,
2017); and manual TOA excision. Narrowband TOAs derived from CHIME/Pulsar data
were cleaned using a similar procedure, though only consisted 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.

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 $\mu$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.

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)444https://github.com/pennucci/PulsePortraiture, which we briefly
summarize as follows. A single, wideband TOA and instantaneous DM were
estimated for each observation using a generalized Fourier phase-gradient
algorithm 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 portraits 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 wideband data sets described in Section 2.4. For each TOA
set, we combined NANOGrav and CHIME/Pulsar TOAs and used the tempo pulsar-
timing package555https://ascl.net/1509.002 for conversion and generalized
least-squares modeling (Coles et al., 2011) in terms of parameters that are
specific to relevant physical processes as well as instrumental effects. A
summary of the timing procedure and model is provided below.

###  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.666The values of these
offsets are noted in the “-to” metadata flag in the TOA data made publicly
available in this work.

###  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 ($\varpi$); 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_{\rm 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_{\rm 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 ($\omega$), and epoch of passage
through periastron ($T_{0}$) – in the following manner:

$\displaystyle e$ $\displaystyle=\sqrt{\epsilon_{1}^{2}+\epsilon_{2}^{2}}$ (1)
$\displaystyle\omega$
$\displaystyle=\arctan{\big{(}\epsilon_{1}/\epsilon_{2}\big{)}}$ (2)
$\displaystyle T_{0}$ $\displaystyle=T_{\rm asc}+P_{\rm
b}\times\omega/(2\pi).$ (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={\rm T}_{\odot}m_{\rm c}$, where T${}_{\odot}=G{\rm
M}_{\odot}/c^{3}$ = 4.925490947 $\mu$s and $m_{\rm c}$ is the mass of the
companion star. With estimates of $m_{\rm c}$ and $\sin i$, the pulsar mass
($m_{\rm p}$) is determined by the Keplerian mass function,

$f_{m}=\frac{4\pi^{2}}{{\rm T}_{\odot}}\frac{x^{3}}{P_{\rm
b}^{2}}=\frac{(m_{\rm c}\sin i)^{3}}{(m_{\rm p}+m_{\rm c})^{2}}.$ (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}_{\rm b}$) possessed sufficient statistical significance for
inclusion as a degree of freedom. We interpret and discuss the implications of
the $\dot{P}_{\rm 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_{\rm FD}=\sum_{i=1}^{n}c_{i}\ln{(f/1\mathrm{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 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_{\rm 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., NANOGrav 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 NANOGrav 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 $\sim 1$. We therefore did not include terms that
quantify red noise in the noise model for PSR J0740+6620, as NANOGrav sets a
Bayes factor threshold of 100 for including red-noise parameters into timing
models (e.g., NANOGrav Collaboration et al., 2015).

##  4 Analysis

Table 2: A summary of timing parameters for PSR J0740+6620 when using DMX to estimate DM valuesaaWe 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. in 6.5-day bins. ParameterbbAll parameter uncertainties, listed in parentheses, denote the 68.3% (i.e., 1 $\sigma$) confidence intervals in the preceding digit(s) as obtained from tempo. | Narrowband TOAs | Wideband TOAs | Difference ($\sigma_{\rm max}$)ccDifferences between best-fit parameter estimates derived from narrowband and wideband TOAs are listed in units of the larger of the two statistical uncertainties ($\sigma_{\rm max}$).
---|---|---|---
Astrometry
Ecliptic longitude, $\lambda$ (degrees) | 103.759135333(12) | 103.759135338(13) | 0.4
Ecliptic latitude, $\beta$ (degrees) | 44.102478368(13) | 44.102478361(13) | 0.5
Proper motion in $\lambda$, $\mu_{\lambda}$ (mas yr-1) |  $-$2.735(14) |  $-$2.737(15) | 0.1
Proper motion in $\beta$, $\mu_{\beta}$ (mas yr-1) |  $-$32.48(2) |  $-$32.48(2) | 0.0
Timing parallax, $\varpi$ (mas) | 1.04(18) | 0.87(19) | 0.9
Spin
Spin frequency, $\nu_{\rm s}$ (s-1) | 346.5319964608338(3) | 346.5319964608337(3) | 0.0
Time rate of change in frequency, $\dot{\nu}_{\rm s}$ ($10^{-15}$ s-2) |  $-$1.463874(11) |  $-$1.463870(11) | 0.4
Binary Motion
Orbital period, $P_{\rm b}$ (days) | 4.76694461933(8) | 4.76694461936(8) | 0.4
Projected semi-major axis, $x$ (lt-s) | 3.97755608(10) | 3.97755607(11) | 0.0
First Laplace-Lagrange eccentricity parameter, $\epsilon_{1}$ ($10^{-6}$) |  $-$5.68(3) |  $-$5.70(3) | 0.7
Second Laplace-Lagrange eccentricity parameter, $\epsilon_{2}$ ($10^{-6}$) |  $-$1.833(18) |  $-$1.840(19) | 0.3
Epoch of ascending-node passage, $T_{\rm asc}$ | 57804.731308893(17) | 57804.731308895(18) | 0.1
Time rate of change in period, $\dot{P}_{\rm b}$ (10-12 s s-1) | 1.2(2) | 1.2(2) | 0.2
Companion mass, $m_{\rm c}$ (M⊙) | 0.251(5) | 0.253(6) | 0.3
Sine of inclination angle, $\sin i$ | 0.99909(12) | 0.99908(13) | 0.1
Profile Evolution
First coefficient of FD expansion, $c_{1}$ ($10^{-5}$) |  $-$3.2(2) | n/a | n/a
Second coefficient of FD expansion, $c_{2}$ ($10^{-5}$) |  $-$2.0(3) | n/a | n/a
Third coefficient of FD expansion, $c_{3}$ ($10^{-5}$) |  $-$1.0(2) | n/a | n/a
Fit Configuration and Statistics
Reference epoch for spin, astrometry (MJD) | 57807 | 57807 | n/a
Terrestrial clock standard | TT(BIPM2019) | TT(BIPM2019) | n/a
Solar system ephemeris | DE438 | DE438 | n/a
Barycentric time scale | TDB | TDB | n/a
Degrees of freedom | 18,041 | 1,099 | n/a
Goodness of fit, $\chi^{2}$ | 18,183 | 1,072 | n/a
Daily-averaged, weighted RMS ($\mu$s) | 0.27 | 0.28 | n/a
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 $\Delta$DM values are offset from the nominal zero
value of 14.9631 pc cm-3. 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. Figure 3: A zoomed-in view of the $\sim 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.

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 ($\sim$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 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}_{\rm
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_{\rm b}$, such as those due to: energy loss from
quadrupole-order gravitational radiation (($\dot{P}_{\rm b}$)GR; e.g., Damour
& Taylor, 1992); differential rotation in the Galaxy (($\dot{P}_{\rm
b}$)${}_{\textrm{DR}}$; e.g., Nice & Taylor, 1995); off-plane acceleration in
the Galactic gravitational potential (($\dot{P}_{\rm b}$)z; Kuijken & Gilmore,
1989); and apparent acceleration due to transverse motion (($\dot{P}_{\rm
b}$)μ; Shklovskii, 1970). We assumed that the total observed variation is
therefore ($\dot{P}_{\rm b})_{\rm obs}$ = $(\dot{P}_{\rm b})_{\rm GR}$ \+
$(\dot{P}_{\rm b})_{\rm DR}$ \+ $(\dot{P}_{\rm b})_{z}$ \+ $(\dot{P}_{\rm
b})_{\mu}$, where

$\displaystyle(\dot{P}_{\rm b})_{\rm GR}$
$\displaystyle=-\frac{192\pi}{5}\frac{(n_{\rm b}{\rm
T}_{\odot})^{5/3}}{(1-e^{2})^{7/2}}\frac{m_{\rm p}m_{\rm c}}{m_{\rm
tot}^{1/3}}\bigg{(}1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\bigg{)},$ (5)
$\displaystyle(\dot{P}_{\rm b})_{\rm DR}$ $\displaystyle=-P_{\rm b}\cos
b\bigg{(}\frac{\Theta_{0}^{2}}{cR_{0}}\bigg{)}\bigg{(}\cos
l+\frac{\kappa^{2}}{\sin^{2}l+\kappa^{2}}\bigg{)},$ (6)
$\displaystyle(\dot{P}_{\rm b})_{\rm z}$ $\displaystyle=-1.08\times
10^{-19}\frac{P_{\rm
b}}{c}\bigg{(}\frac{1.25z}{[z^{2}+0.0324]^{1/2}}+0.58z\bigg{)}\sin b,$ (7)
$\displaystyle(\dot{P}_{\rm b})_{\mu}$ $\displaystyle=\frac{\mu^{2}d}{c}P_{\rm
b},$ (8)

and where: $n_{\rm b}=2\pi/P_{\rm b}$ is the orbital frequency; $m_{\rm
tot}=m_{\rm p}+m_{\rm c}$; $l$ and $b$ are the Galactic longitude and
latitude, respectively; $\Theta_{0}$ and $R_{0}$ are the Galactrocentric
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-1 and $R_{0}$ = 8.178(26)777The 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. 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 $\sim$20% uncertainty in ($\dot{P}_{\rm b}$)obs is too
large for resolving statistically meaningful differences between model
predictions (Pathak & Bagchi, 2018).

While expected to be negligible, we included the ($\dot{P}_{\rm 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
$\dot{P}_{\rm b}$ corrections that can eventually be observed through pulsar
timing, arising from mechanisms such as mass loss 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}_{\rm b}$, and we
therefore chose to ignore those terms in subsequent calculations.

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}_{\rm b}$)obs.

Equations 5-8 demonstrate that ($\dot{P}_{\rm b}$)obs is ultimately a function
of $m_{\rm p}$, $m_{\rm c}$, and $d$.888($\dot{P}_{\rm b}$)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. To check for self-
consistency in our best-fit timing models, we performed a Monte Carlo analysis
of the relation ($\dot{P}_{\rm b})_{\rm obs}$ = $(\dot{P}_{\rm b})_{\rm GR}$
\+ $(\dot{P}_{\rm b})_{\rm DR}$ \+ $(\dot{P}_{\rm b})_{z}$ \+ $(\dot{P}_{\rm
b})_{\mu}$ by sampling a set of values for {($\dot{P}_{\rm b}$)obs, $m_{\rm
p}$, $m_{\rm 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 $\varpi=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}_{\rm 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 $\varpi$ was not corrected for the Lutz-Kelker bias in
pulsar timing-parallax measurements, which arises in low-precision
measurements of $\varpi$ and can lead to statistical underestimates of $d$
(Verbiest 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 $\varpi$, 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_{\rm p}$,
$m_{\rm c}$, $i$, and $d$}. The procedure we used for PDF and credible-
interval estimation is described as follows:

  * •

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

  * •

computed values of $\sin i$, $m_{\rm p}$, $\varpi=1/d$, and ($\dot{P}_{\rm
b}$)tot using Equations 4-8, based on the current location in the ($m_{\rm
c}$, $\cos i$, $d$) phase space;

  * •

held the corresponding values of {$\varpi$, $m_{\rm c}$, $\sin i$,
($\dot{P}_{\rm b}$)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_{\rm c}$, $\cos i$, $d$).

We ultimately obtained a three-dimensional grid of $\chi^{2}$ values computed
over uniform steps in ($m_{\rm c}$, $\cos i$, $d$), which were then mapped to
a likelihood function $p({\rm data,DMX}_{j}|m_{\rm c},\cos
i,d)\propto\exp(-\Delta\chi^{2}/2)$ where $\Delta\chi^{2}=\chi^{2}-\textrm{
min}(\chi^{2})$. The notation “DMXj” 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_{\rm c},\cos
i,d|{\rm data,DMX}_{j})$ for uniform priors on the physical parameters and
choice of DMX model.

Systematic uncertainties in the ($m_{\rm p},m_{\rm 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

$\displaystyle p(m_{\rm c},\cos i,d|{\rm data})$
$\displaystyle=\sum_{j}p(m_{\rm c},\cos i,d|{\rm data,DMX}_{j})$
$\displaystyle\phantom{{}={}}\times p({\rm DMX}_{j}|{\rm data})$ (9)

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

Figure 5: A comparison of normalized posterior PDFs of $m_{\rm 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.

Posterior PDFs and credible intervals for $m_{\rm p}$ were derived from the
original three-dimensional posterior PDF by noting that

$\displaystyle p(m_{\rm p},\cos i,d|{\rm data})$ $\displaystyle=p(m_{\rm
c},\cos i,d|{\rm data})\bigg{|}\frac{\partial m_{\rm c}}{\partial m_{\rm
p}}\bigg{|},$ (10) $\displaystyle\textrm{where }\frac{\partial m_{\rm
c}}{\partial m_{\rm p}}$ $\displaystyle=\frac{2f_{m}m_{\rm tot}}{3m_{\rm
c}^{2}\sin^{3}i-2f_{m}m_{\rm tot}},$ (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 $\sim 10^{-7}$; 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: $\sim$0.3 $\mu$s for the GBT/1400-MHz data;
$\sim$0.4 $\mu$s for GBT/820-MHz; and $\sim$1 $\mu$s for CHIME/Pulsar data.
The factor of $\sim$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 NANOGrav data
set. Therefore, while slightly less precise than its NANOGrav 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_{\rm
c}=0.28(2)$ M⊙ and $\sin i=0.9989(8)$ that correspond to $m_{\rm 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
$\sim$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 $\sim$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 $\sim
8\times 10^{-4}$ pc cm-3, while the same measure in the CHIME/Pulsar data set
is $\sim 2\times 10^{-4}$ pc cm-3. The median DMX measurement uncertainty in
the CHIME/Pulsar era has comparable precision to that obtained with GLOW by
Donner et al. (2020), $\sim 4\times 10^{-5}$ pc cm-3. 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

Table 3: Credible intervals for the Shapiro-delay and distance parameters in the J0740+6620 system.  Parameter |  Cromartie et al. (2020)aaValues 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. | DMX=3.0 | DMX=6.5 | Hybrid-DMX | Model-Averaged
---|---|---|---|---|---
Pulsar mass, $m_{\rm p}$ (M⊙) | $2.14^{+0.10}_{-0.09}$ | $2.08^{+0.07}_{-0.07}$ | $2.07^{+0.07}_{-0.06}$ | $2.10^{+0.07}_{-0.07}$ | $2.08^{+0.07}_{-0.07}$
Companion mass, $m_{\rm c}$ (M⊙) | $0.260^{+0.008}_{-0.007}$ | $0.253^{+0.006}_{-0.005}$ | $0.252^{+0.005}_{-0.005}$ | $0.254^{+0.006}_{-0.006}$ | $0.253^{+0.006}_{-0.005}$
System inclination, $i$ (degrees) | $87.38^{+0.2}_{-0.2}$ | $87.53^{+0.17}_{-0.18}$ | $87.61^{+0.16}_{-0.16}$ | $87.53^{+0.17}_{-0.18}$ | $87.56^{+0.17}_{-0.18}$
Distance, $d$ (kpc) |  $<4.4$bbNeither $\dot{P}_{\rm 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. | $1.16^{+0.17}_{-0.15}$ | $1.11^{+0.17}_{-0.15}$ | $1.15^{+0.17}_{-0.15}$ | $1.14^{+0.17}_{-0.15}$

We performed the $\chi^{2}$-grid and PDF analyses outlined in Section 4.3 for
all three DMX timing solutions derived from the combined NANOGrav+CHIME/Pulsar
data set for PSR J0740+6620. A summary of the credible intervals derived from
all posterior PDFs is presented 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.

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.

Our 68.3% credible interval of the model-averaged $m_{\rm
p}=2.08^{+0.07}_{-0.07}\textrm{ M}_{\odot}$ is reduced by $\sim$30% in
comparison to the estimate made by Cromartie et al. (2020). The model-averaged
$m_{\rm 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_{\rm 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_{\rm
c}=0.253^{+0.006}_{-0.005}\;{\rm M}_{\odot}$, 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_{\rm c}$ is in improved
agreement with correlation parameters that define the $m_{\rm c}-P_{\rm 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_{\rm c}\sim 0.25\;{\rm M}_{\odot}$. 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_{\rm 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 $\sim 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 $\varpi=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 $\varpi=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,
instead of a finite-term Fourier expansion that is more appropriate for low-
inclination systems. The combination 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).

Figure 7: Color–magnitude diagram with white dwarf evolutionary sequences,
evaluated at the $r^{\prime}$ and $i^{\prime}$ 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⊙, 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.

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)
concluded, based on their derived magnitudes and colors of the optical
counterpart and the range of contemporaneous 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 measurements, are currently not
possible as the $r^{\prime}$ and $i^{\prime}$ apparent magnitudes of the
white-dwarf companion to PSR J0740+6620 are $\sim 26$ (Beronya et al., 2019).
These observed magnitudes are well above the limiting $g^{\prime}$ magnitude
of $\sim 19$ for Gaia Data Release
2.999https://www.cosmos.esa.int/web/gaia/dr2

Figure 7 shows a comparison of the absolute magnitude and color of the optical
source with cooling predictions for white dwarfs with hydrogen and helium
atmospheres, and over a range of
masses101010http://www.astro.umontreal.ca/$\sim$bergeron/CoolingModels
(Holberg & Bergeron, 2006; Kowalski & Saumon, 2006; Tremblay et al., 2011;
Bergeron et al., 2011). The range of prior constraints on $d$ led to a large
range in absolute magnitude, shown as the blue-shaded region in Figure 7, that
encompasses cooling tracks of white dwarfs with $0.2-1\;{\rm 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 $\sim 0.2-0.3\;{\rm M}_{\odot}$ remain consistent
with observations. This result is in agreement with the initial estimate of
$m_{\rm c}$ made by Cromartie et al. (2020), and with our updated $m_{\rm
c}=0.253^{+0.006}_{-0.005}\textrm{ M}_{\odot}$. The improved synergy in
optical- and radio-based estimates of intrinsic 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
measurement of $\dot{P}_{\rm 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,
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., $\sim$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_{\rm
p}=2.08^{+0.07}_{-0.07}$ M⊙, $m_{\rm c}=0.253^{+0.006}_{-0.005}$ M⊙ and
$i=87.56^{+0.17\circ}_{-0.18}$, consistent with the estimates made by
Cromartie et al. (2020). Our constraint on $m_{\rm p}$ is improved by
$\sim$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}_{\rm b}$ which we argued is due to apparent acceleration from
transverse motion. We found that the model-averaged $d=1.14^{+0.17}_{-0.15}$
kpc when using the observed $\dot{P}_{\rm 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_{\rm TOA}\propto S_{\rm sys}/\sqrt{t_{\rm obs}\Delta f}$, where:
$S_{\rm sys}$ is the system-equivalent flux density of the observatory;
$t_{\rm obs}$ is the observation timespan, and $\Delta f$ is the receiver
bandwidth.111111For 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. The expected
improvement of TOA precision between CHIME/Pulsar (with $S_{\rm 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_{\rm sys}=9$ Jy and $\Delta f=$ 1200 MHz (Vanderlinde et al., 2019). For
observations of the same $t_{\rm obs}$, the band/time-averaged TOA uncertainty
obtained with CHORD is $\sigma_{\rm TOA}\sim 0.1$ $\mu$s when using
CHIME/Pulsar values listed in Table 1. The DSA-2000 and ngVLA have similar
planned bandwidths, but a factor of $\sim$ 5 smaller $S_{\rm sys}$ than that
of CHORD. With CHORD projected to be operational by 2025, it is likely that
daily TOA uncertainties of 0.1 $\mu$s and lower will be regularly achieved
with next-generation observatories as soon as several years from now.

A reduction in $\sigma_{\rm 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_{\rm p}$ to exceed values no
larger than 0.01 M⊙ of the median value with sufficient observations from
CHORD, the DSA-2000 and/or ngVLA. The significance of $\dot{P}_{\rm 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.

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 referee 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 Astronomy 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 NANOGrav 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).

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJL, 896, L44, doi: 10.3847/2041-8213/ab960f
  * Alam et al. (2021a) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021a, ApJS, 252, 4
  * Alam et al. (2021b) —. 2021b, ApJS, 252, 5
  * Antoniadis et al. (2016) Antoniadis, J., Tauris, T. M., Özel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448
  * Arzoumanian et al. (2018) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37, doi: 10.3847/1538-4365/aab5b0
  * Backer & Hellings (1986) Backer, D. C., & Hellings, R. W. 1986, Ann. Rev. Astr. Ap., 24, 537
  * Bergeron et al. (2011) Bergeron, P., Wesemael, F., Dufour, P., et al. 2011, ApJ, 737, 28
  * Beronya et al. (2019) Beronya, D. M., Karpova, A. V., Kirichenko, A. Y., et al. 2019, MNRAS, 485, 3715
  * Bogdanov et al. (2019) Bogdanov, S., Guillot, S., Ray, P. S., et al. 2019, ApJL, 887, L25
  * Chatterjee (2018) 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 et al. (2021) CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5
  * Coles et al. (2011) Coles, W., Hobbs, G., Champion, D. J., Manchester, R. N., & Verbiest, J. P. W. 2011, MNRAS, 418, 561
  * Cordes & Lazio (2002) Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0207156
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nat. Astron., 4, 72
  * Damour & Taylor (1992) Damour, T., & Taylor, J. H. 1992, Phys. Rev. D, 45, 1840
  * Demorest (2018) Demorest, P. B. 2018, nanopipe: Calibration and data reduction pipeline for pulsar timing.  http://ascl.net/1803.004
  * Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nat., 467, 1081
  * Donner et al. (2020) Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2020, A&A, 644, A153
  * DuPlain et al. (2008) 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 et al. (2020) Echeveste, M., Novarino, M. L., Benvenuto, O. G., & De Vito, M. A. 2020, MNRAS, 495, 2509
  * Ellis et al. (2020) 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 et al. (2021) 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 & Wex (2010) Freire, P. C. C., & Wex, N. 2010, MNRAS, 409, 199
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10
  * Hallinan et al. (2019) Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 255
  * Hoeting et al. (1999) Hoeting, J. A., D., M., & Raftery, A. E. Volinsky, C. T. 1999, Stat. Sci., 14, 382
  * Holberg & Bergeron (2006) Holberg, J. B., & Bergeron, P. 2006, AJ, 132, 1221
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Hu et al. (2020) Hu, H., Kramer, M., Wex, N., Champion, D. J., & Kehl, M. S. 2020, MNRAS, 497, 3118
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Kowalski & Saumon (2006) Kowalski, P. M., & Saumon, D. 2006, ApJL, 651, L137
  * Kuijken & Gilmore (1989) Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571
  * Lam et al. (2017) Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2017, ApJ, 834, 35
  * Lange et al. (2001) Lange, C., Camilo, F., Wex, N., et al. 2001, MNRAS, 326, 274
  * Lentati et al. (2014) Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004
  * Linares et al. (2018) Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54
  * Liu et al. (2014) Liu, K., Desvignes, G., Cognard, I., et al. 2014, MNRAS, 443, 3752
  * Liu et al. (2020) Liu, K., Guillemot, L., Istrate, A. G., et al. 2020, MNRAS, 499, 2276
  * Lynch et al. (2018) Lynch, R. S., Swiggum, J. K., Kondratiev, V. I., et al. 2018, ApJ, 859, 93
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., & et al. 2021, arXiv e-prints, arXiv:2105.06979
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24
  * NANOGrav Collaboration et al. (2015) NANOGrav Collaboration, Arzoumanian, Z., Brazier, A., et al. 2015, ApJ, 813, 65
  * Nice et al. (2015) Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo: Pulsar timing data analysis.  http://ascl.net/1509.002
  * Nice & Taylor (1995) Nice, D. J., & Taylor, J. H. 1995, ApJ, 441, 429
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Ann. Rev. Astr. Ap., 54, 401
  * Pathak & Bagchi (2018) Pathak, D., & Bagchi, M. 2018, ApJ, 868, 123
  * Pennucci (2015) Pennucci, T. T. 2015, PhD thesis, University of Virginia
  * Pennucci (2019) —. 2019, ApJ, 871, 34
  * Pennucci et al. (2014) Pennucci, T. T., Demorest, P. B., & Ransom, S. M. 2014, ApJ, 790, 93
  * Pennucci et al. (2016) —. 2016, Pulse Portraiture: Pulsar timing.  http://ascl.net/1606.013
  * Raaijmakers et al. (2019) Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, ApJL, 887, L22
  * Raaijmakers et al. (2020) Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020, ApJL, 893, L21
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, arXiv e-prints, arXiv:2105.06981
  * Reid et al. (2014) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv e-prints, arXiv:2105.06980
  * Romani et al. (2021) Romani, R. W., Kandel, D., Filippenko, A. V., & et al. 2021, ApJL, 908, L46
  * Shannon & Cordes (2010) Shannon, R. M., & Cordes, J. M. 2010, ApJ, 725, 1607
  * Shapiro (1964) Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789
  * Shklovskii (1970) Shklovskii, I. S. 1970, Sov. Astron., 13, 562
  * Splaver et al. (2002) Splaver, E. M., Nice, D. J., Arzoumanian, Z., et al. 2002, ApJ, 581, 509
  * Stairs (2003) Stairs, I. H. 2003, Living Reviews in Relativity, 6, 5
  * Tauris & Savonije (1999) Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928
  * Taylor (1992) Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 117
  * Tremblay et al. (2011) Tremblay, P.-E., Bergeron, P., & Gianninas, A. 2011, ApJ, 730, 128
  * Tsokaros et al. (2020) Tsokaros, A., Ruiz, M., & Shapiro, S. L. 2020, ApJ, 905, 48
  * Vallisneri & van Haasteren (2017) Vallisneri, M., & van Haasteren, R. 2017, MNRAS, 466, 4954
  * van Kerkwijk et al. (2011) van Kerkwijk, M. H., Breton, R. P., & Kulkarni, S. R. 2011, ApJ, 728, 95
  * van Straten & Bailes (2011) van Straten, W., & Bailes, M. 2011, PASA, 28, 1
  * Vanderlinde et al. (2019) 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 et al. (2010) Verbiest, J. P. W., Lorimer, D. R., & McLaughlin, M. A. 2010, MNRAS, 405, 564
  * Wolff et al. (2021) Wolff, M. T., Guillot, S., Bogdanov, S., et al. 2021, submitted to ApJL
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29