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

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), ³⁸ 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, 2021b), 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 centered at radio frequencies of 820 and 1400 MHz; these receivers possessed full bandwidths of 200 and 800 MHz, respectively. Moreover, the spectrum of each folded profile was resolved over 128 and 512 channels when using the 820 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 hr session occurred on MJD 58368 using the 820 MHz receiver during orbital phase of 0.25, and two 5 hr 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 composed of four half-cylinder reflectors and a total of 1024 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 time series. The CHIME/Pulsar back end (CHIME/Pulsar Collaboration et al. 2020) is a 10-node computing system that receives these 10 independent, complex voltage streams—at resolutions of 2.56 μs and 1024 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 1024 channels evaluated across the 400 MHz bandwidth of CHIME for each recorded profile. At the decl. of PSR J0740+6620, acquisitions typically occur with a duration of ∼33 minutes in order to record a full transit at the CHIME telescope.

Under the nominal NANOGrav program, multifrequency 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 NANOGrav session, which led to three to five multifrequency observations within a 1-week period of time during each month. See Alam et al. (2021a, 2021b) for additional discussion on the NANOGrav observing program.

Folded data were acquired with the CHIME/Pulsar coherent-dedispersion back end 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 ∼1-week periods where the CHIME correlator ceased operation for software upgrades and hardware maintenance. Such activities only occurred three to four times per year during the time span of CHIME data presented in this work.

All data acquired with the GBT and CHIME/Pulsar were processed offline using the psrchive ³⁹ 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 pipeline ⁴⁰ —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 owing 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 ratio (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 was 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.

We generated two distinct sets of times of arrival (TOAs) and TOA uncertainties (σ_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 this 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.

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), ⁴¹ 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 (FD) 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).

Zoom In Zoom Out Reset image size

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 as a result of improper packaging of timing data required for downstream, real-time processing and determination of UTC time stamps. 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 we used those values as corrections to data acquired between MJDs 58517–58543. ⁴³

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 (ϖ); 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 semimajor 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, ₁ and ₂, that describe departures from circular motion. The ELL1 parameters are related to the traditional orbital elements—eccentricity (e), argument of periastron (ω), and epoch of passage through periastron (T₀)—in the following manner:

$$\begin{eqnarray}&&e=\sqrt{{\epsilon }_{1}^{2}+{\epsilon }_{2}^{2}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\omega =\arctan \left({\epsilon }_{1}/{\epsilon }_{2}\right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{T}_{0}={T}_{\mathrm{asc}}+{P}_{{\rm{b}}}\times \omega /(2\pi ).\end{eqnarray} \tag{ 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_⊙ m_c , where T_⊙ = GM_⊙/c⁻³ = 4.925490947 μ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,

$$\begin{eqnarray}&&{f}_{m}=\displaystyle \frac{4{\pi }^{2}}{{\text{}}{T}_{\odot }}\displaystyle \frac{{x}^{3}}{{P}_{{\rm{b}}}^{2}}=\displaystyle \frac{{({m}_{{\text{}}c}\sin i)}^{3}}{{({m}_{{\text{}}p}+{m}_{{\text{}}c})}^{2}}.\end{eqnarray} \tag{ 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σ 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}}_{{\text{}}b}$) possessed sufficient statistical significance for inclusion as a degree of freedom. We interpret and discuss the implications of the ${\dot{P}}_{{\text{}}b}$ measurement in several sections below.

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 FD variations in our narrowband TOA data set using the heuristic model developed by NANOGrav Collaboration et al. (2015), where the associated time delay ${\rm{\Delta }}{t}_{\mathrm{FD}}={\sum }_{i=1}^{n}{c}_{i}\mathrm{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 FD 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).

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 front-end/back-end combinations and corresponds to the three subsets presented in Table 1 for this work. The factor F_k accounts for incorrect estimation of σ_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 it applies to simultaneously observed, multifrequency 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, 2021b) 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 multifrequency 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 ∼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 in timing models (e.g., NANOGrav Collaboration et al. 2015).

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.

All fitted timing models exhibit a significant and consistent ${\dot{P}}_{{\text{}}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 ((${\dot{P}}_{{\text{}}b}$)_GR; e.g., Damour & Taylor 1992), differential rotation in the Galaxy ((${\dot{P}}_{{\text{}}b}$)_DR; e.g., Nice & Taylor 1995), off-plane acceleration in the Galactic gravitational potential ((${\dot{P}}_{{\text{}}b}$)_z ; Kuijken & Gilmore 1989), and apparent acceleration due to transverse motion ((${\dot{P}}_{{\text{}}b}$)_μ ; Shklovskii 1970). We assumed that the total observed variation is therefore (${\dot{P}}_{b}{)}_{\mathrm{obs}}$ = ${({\dot{P}}_{b})}_{\mathrm{GR}}$ + ${({\dot{P}}_{b})}_{\mathrm{DR}}$ + ${({\dot{P}}_{b})}_{z}$ + ${({\dot{P}}_{b})}_{\mu }$, where

$$\begin{eqnarray}&&{({\dot{P}}_{{\text{}}b})}_{\mathrm{GR}}=-\displaystyle \frac{192\pi }{5}\displaystyle \frac{{({n}_{{\text{}}b}{\text{}}{T}_{\odot })}^{5/3}}{{(1-{e}^{2})}^{7/2}}\displaystyle \frac{{m}_{{\text{}}p}{m}_{{\text{}}c}}{{m}_{\mathrm{tot}}^{1/3}}\left(1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}^{4}\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{({\dot{P}}_{{\text{}}b})}_{\mathrm{DR}}=-{P}_{{\text{}}b}\cos b\left(\displaystyle \frac{{{\rm{\Theta }}}_{0}^{2}}{{{cR}}_{0}}\right)\left(\cos l+\displaystyle \frac{{\kappa }^{2}}{{\sin }^{2}l+{\kappa }^{2}}\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{rcl}{({\dot{P}}_{{\text{}}b})}_{{\text{}}z} & = & -1.08\times {10}^{-19}\displaystyle \frac{{P}_{{\text{}}b}}{c}\\ & & \times \left(\displaystyle \frac{1.25z}{{[{z}^{2}+0.0324]}^{1/2}}+0.58z\right)\sin b,\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{({\dot{P}}_{{\text{}}b})}_{\mu }=\displaystyle \frac{{\mu }^{2}d}{c}{P}_{{\text{}}b},\end{eqnarray} \tag{ 8 }$$

and where n_b = 2π/P_b is the orbital frequency; m_tot = m_p + m_c ; l and b are the Galactic longitude and latitude, respectively; Θ₀ and R₀ 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 Θ₀ = 236.9(4.2) km s⁻¹ and R₀ = 8.178(26) ⁴⁴ 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 (${\dot{P}}_{{\text{}}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}}_{{\text{}}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}}_{{\text{}}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}}_{{\text{}}b}$, and we therefore chose to ignore those terms in subsequent calculations.

Equations (5)–(8) demonstrate that (${\dot{P}}_{{\text{}}b}$)_obs is ultimately a function of m_p , m_c , and d. ⁴⁵ To check for self-consistency in our best-fit timing models, we performed a Monte Carlo analysis of the relation (P_b )_obs=${({\dot{P}}_{b})}_{\mathrm{GR}}$ + ${({\dot{P}}_{b})}_{\mathrm{DR}}$ + ${({\dot{P}}_{b})}_{z}$ + ${({\dot{P}}_{b})}_{\mu }$ by sampling a set of values for {(${\dot{P}}_{{\text{}}b}$)_obs, m_p , m_c , R₀, Θ₀, μ}, 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 probability density function (pdf) for d obtained from measurement of ϖ = 1/d. The overlapping pdf's demonstrate that the two independent timing effects yield statistically consistent estimates of d. The width of the pdf obtained from ${\dot{P}}_{{\text{}}b}$ is dominated by its measurement uncertainty, and width contributions from sampling the {R₀, Θ₀, μ} terms are negligible. We therefore held the {R₀, Θ₀, μ} terms fixed to their best-fit values in subsequent analysis.

Zoom In Zoom Out Reset image size

The pdf derived from ϖ was not corrected for the Lutz–Kelker bias in pulsar timing parallax measurements, which arises in low-precision measurements of ϖ 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 ϖ, and we therefore ignore this bias in subsequent analysis.

For all timing models presented here, we used the χ²-grid method developed by Splaver et al. (2002) for determining posterior pdf's and robust credible intervals for {m_p , m_c , i, and d}. The procedure we used for pdf and credible-interval estimation is described as follows:

• 1.  
  selected a set of values for ${m}_{c},\cos i=\sqrt{1-{\sin }^{2}i}$, and d, each from uniform distributions;
• 2.  
  computed values of sin i, m_p , ϖ = 1/d, and (${\dot{P}}_{{\text{}}b}$)_tot using Equations (4)–(8), based on the current location in the (m_c , cos i, d) phase space;
• 3.  
  held the corresponding values of {ϖ, m_c, $\sin i$, (${\dot{P}}_{{\text{}}b}$)_tot} fixed in the timing model;
• 4.  
  refitted the timing model, allowing all timing model parameters not defined on the grid to remain unconstrained during the fit;
• 5.  
  recorded the best-fit χ² value, and repeated the above steps for different values of (m_c, cos i, d).

We ultimately obtained a three-dimensional grid of χ² values computed over uniform steps in (m_c , cos i, d), which were then mapped to a likelihood function $p(\mathrm{data},\ {\mathrm{DMX}}_{j}| {m}_{c},\cos i,d)\propto \exp (-{\rm{\Delta }}{\chi }^{2}/2)$, where ${\rm{\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's theorem to compute the posterior distribution $p({m}_{{\rm{c}}},\cos i,d| \mathrm{data},\ {\mathrm{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 as a result of choices in modeling DM variations. We used Bayesian model averaging (e.g., Hoeting et al. 1999) to obtain pdf's 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 pdf's; in the case of the three-dimensional posterior distribution described above, the model-averaged pdf is

$$\begin{eqnarray}\begin{array}{rcl}p({m}_{c},\cos i,d| \mathrm{data}) & = & \displaystyle \sum _{j}p({m}_{c},\cos i,d| \mathrm{data},\ {\mathrm{DMX}}_{j})\\ & & \times p({\mathrm{DMX}}_{j}| \mathrm{data}),\end{array}\end{eqnarray} \tag{ 9 }$$

where p(DMX_j ∣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 pdf's. We used Equation (9) to compute model-averaged posterior pdf's and their corresponding credible intervals for all four gridded parameters. The DMX_j and model-averaged posterior pdf's for m_p are shown in Figure 5.

Zoom In Zoom Out Reset image size

Posterior pdf's and credible intervals for m_p were derived from the original three-dimensional posterior pdf by noting that

$$\begin{eqnarray}&&p({m}_{p},\cos i,d| \mathrm{data})=p({m}_{c},\cos i,d| \mathrm{data})\left|\displaystyle \frac{\partial {m}_{c}}{\partial {m}_{p}}\right|,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\mathrm{where}\displaystyle \frac{\partial {m}_{c}}{\partial {m}_{p}}=\displaystyle \frac{2{f}_{m}{m}_{\mathrm{tot}}}{3{m}_{c}^{2}{\sin }^{3}i-2{f}_{m}{m}_{\mathrm{tot}}},\end{eqnarray} \tag{ 11 }$$

and where Equation (11) is determined by the Keplerian mass function (f_m ). Using these original and translated pdf's, we computed credible intervals by marginalizing over the relevant subset of parameters, and we also calculated two-dimensional posterior pdf's for all possible parameter pairs. We ignored the uncertainty in f_m for all pdf calculations owing 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.

Table 3. Credible Intervals for the Shapiro Delay and Distance Parameters in the J0740+6620 System

Parameter	Cromartie et al. (2020) a	DMX = 3.0	DMX = 6.5	Hybrid-DMX	Model-averaged
Pulsar mass, mp (M⊙)	${2.14}_{-0.09}^{+0.10}$	${2.08}_{-0.07}^{+0.07}$	${2.07}_{-0.06}^{+0.07}$	${2.10}_{-0.07}^{+0.07}$	${2.08}_{-0.07}^{+0.07}$
Companion mass, mc (M⊙)	${0.260}_{-0.007}^{+0.008}$	${0.253}_{-0.005}^{+0.006}$	${0.252}_{-0.005}^{+0.005}$	${0.254}_{-0.006}^{+0.006}$	${0.253}_{-0.005}^{+0.006}$
System inclination, i (deg)	${87.38}_{-0.2}^{+0.2}$	${87.53}_{-0.18}^{+0.17}$	${87.61}_{-0.16}^{+0.16}$	${87.53}_{-0.18}^{+0.17}$	${87.56}_{-0.18}^{+0.17}$
Distance, d (kpc)	<4.4 b	${1.16}_{-0.15}^{+0.17}$	${1.11}_{-0.15}^{+0.17}$	${1.15}_{-0.15}^{+0.17}$	${1.14}_{-0.15}^{+0.17}$

Notes.

^a Values initially determined by Cromartie et al. (2020) are listed for comparison. The χ² gridding procedure used by Cromartie et al. (2020) only considered one DMX model, with maximum bin size of 6.5 days. ^b Neither ${\dot{P}}_{{\text{}}b}$ nor timing parallax—and thus the distance—was significantly constrained with the data set presented by Cromartie et al. (2020). We instead quote here the 2σ upper limit obtained from their reported constraint.

Download table as:  ASCIITypeset image

All best-fit timing models developed for this work yielded consistent weighted rms residuals: ∼0.3 μs for GBT/1400 MHz data, ∼0.4 μs for GBT/820 MHz data, 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. 2020).

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 owing 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 yr 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.

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 (β ≈ 44°; 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 time series for PSR J0740+6620 will be the subject of future works.

We performed the χ²-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 pdf's is presented in Table 3. The model-averaged pdf's and credible intervals are shown in Figure 6.

Zoom In Zoom Out Reset image size

For all four parameters defined on each χ² grid, the three posterior pdf's 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 owing 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}_{\odot }$ 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.005}^{+0.006}\ {M}_{\odot }$, remains largely consistent with the prediction from expected correlations between companion masses and orbital sizes that arise as a result of 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 ∼ 10⁻³ or lower), which was noted by Cromartie et al. (2020) to yield an expected m_c ∼ 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.

We obtained a constrained, model-averaged distance of $d={1.14}_{-0.15}^{+0.17}$ kpc from the combined NANOGrav and CHIME/Pulsar TOA data set when using the χ²-grid method described in Section 4.3. This updated distance estimate is consistent with the distance of d ≈ 0.9 kpc derived from the observed mean DM, placement within the Milky Way, 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 ≈ 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σ level with the marginal estimate of d derived from ϖ = 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.1}^{+0.2}$ kpc made from initial NANOGrav timing of PSR J0740+6620 (derived from ϖ = 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 suboptimal 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) 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 ultracool, helium-atmosphere white dwarf with an effective temperature <3500 K and cooling age >5 Gyr. However, Beronya et al. (2019) were not able to constrain the mass of the companion owing to the significant distance ambiguity. Observations with the Gaia satellite, designed for astrometric measurements, 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' magnitude of ∼19 for Gaia Data Release 2. ⁴⁶

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 masses ⁴⁷ (Holberg & Bergeron 2006; Kowalski & Saumon 2006; Bergeron et al. 2011; Tremblay 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 M_⊙. 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_⊙ remain consistent with observations. 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 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).

Zoom In Zoom Out Reset image size

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 modeling (Echeveste et al. 2020). The increasing time span of the PSR J0740+6620 data set will further improve the measurement of ${\dot{P}}_{{\text{}}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.