A More Precise Measurement of the Radius of PSR J0740+6620 Using Updated NICER Data

Abstract

PSR J0740+6620 is the neutron star with the highest precisely determined mass, inferred from radio observations to be $2.08\pm 0.07\,\rm M_{\odot}$. Measurements of its radius therefore hold promise to constrain the properties of the cold, catalyzed, high-density matter in neutron star cores. Previously, Miller et al. (2021) and Riley et al. (2021) reported measurements of the radius of PSR J0740+6620 based on Neutron Star Interior Composition Explorer (NICER) observations accumulated through 17 April 2020, and an exploratory analysis utilizing NICER background estimates and a data set accumulated through 28 December 2021 was presented in Salmi et al. (2022). Here we report an updated radius measurement, derived by fitting models of X-ray emission from the neutron star surface to NICER data accumulated through 21 April 2022, totaling $\sim 1.1$ Ms additional exposure compared to the data set analyzed in Miller et al. (2021) and Riley et al. (2021), and to data from X-ray Multi-Mirror (XMM-Newton) observations. We find that the equatorial circumferential radius of PSR J0740+6620 is $12.92_{-1.13}^{+2.09}$ km (68% credibility), a fractional uncertainty $\sim 83\%$ the width of that reported in Miller et al. (2021), in line with statistical expectations given the additional data. If we were to require the radius to be less than 16 km, as was done in Salmi et al. (2024), then our 68% credible region would become $R=12.76^{+1.49}_{-1.02}$ km, which is close to the headline result of Salmi et al. (2024). Our updated measurements, along with other laboratory and astrophysical constraints, imply a slightly softer equation of state than that inferred from our previous measurements.

1 Introduction

The cores of neutron stars consist of cold matter that is believed to be catalyzed to its ground state at a few times nuclear saturation density (saturation corresponds to a baryon number density $n_{s}\approx 0.16\,\rm{fm^{-3}}$, equivalent to a mass density $\rho_{s}\approx 2.7\times 10^{14}\,\rm{g\,cm^{-3}}$). Laboratories are unable to replicate the extremely high densities, and likely very high neutron-proton asymmetries, present in neutron star cores. Thus, observations of neutron stars provide unique insights about the nature of dense matter. Notably, the densities in neutron star cores partially bridge the gap between the regime $\rho\lesssim\rho_{s}$ probed by nuclear experimentation (Danielewicz et al., 2002; Adhikari et al., 2021, e.g.,) and the regime $\rho\gg\rho_{s}$ where perturbative quantum chromodynamics is currently able to make predictions about the nature of nuclear matter (e.g., Bedaque & Steiner, 2015; Hoyos et al., 2016; Ecker et al., 2017).

In recent years, observations of neutron stars have added considerably to our knowledge of the equation of state (EOS: pressure as a function of energy density) of cold matter - matter in which the thermal energies of the particles are much smaller than their Fermi energies - at high densities (e.g., Pavlov & Zavlin, 1997; Bhattacharyya et al., 2005; Miller, 2021; Özel et al., 2016b; Nättilä et al., 2017).

Three neutron stars have precisely measured masses $\sim 2\,M_{\odot}$: PSR J1614-2230, with $M=1.937\pm 0.014\,M_{\odot}$ (Demorest et al., 2010; Fonseca et al., 2016; Arzoumanian et al., 2018; Agazie et al., 2023); PSR J0348+0432, with $M=2.01\pm 0.04\,M_{\odot}$ (Antoniadis et al., 2013); and PSR J0740+6620, with $M=2.08\pm 0.07\,M_{\odot}$ (Cromartie et al., 2020; Fonseca et al., 2021), where the quoted uncertainties denote the 68% credible regions. Observations of such massive neutron stars can rule out equations of state that predict appreciably lower maximum stable masses. Additionally, observations of the gravitational wave event GW170817 have provided constraints on the tidal deformability of inspiraling neutron stars (e.g., Abbott et al., 2017, 2018; De et al., 2018; Abbott et al., 2020a), ruling out equations of state that would cause neutron stars to have relatively high radii at a given mass. Observations of the kilonovae following binary neutron star mergers may also place upper limits on the maximum stable mass of nonrotating neutron stars (e.g., Margalit & Metzger, 2017; Rezzolla et al., 2018; Pang et al., 2021).

The Neutron Star Interior Composition Explorer (NICER) can make phase- and energy-resolved measurements of the thermal X-ray pulses produced by some rotating neutron stars. These can be used to constrain the masses and radii of these stars, and thus the neutron star EOS (e.g., Watts et al., 2016; Özel et al., 2016a). The first mass and radius constraints derived from modeling NICER X-ray pulse profiles were reported in two independent analyses of the energy-resolved pulse profile of PSR J0030+0451: Miller et al. (2019) measured its mass and equatorial radius to be $1.44_{-0.14}^{+0.15}\,M_{\odot}$ and $13.02_{-1.06}^{+1.24}\,{\rm km}$, respectively, while Riley et al. (2019) measured its mass and radius to be $1.34_{-0.16}^{+0.15}\,M_{\odot}$ and $12.71_{-1.19}^{+1.14}\,{\rm km}$. Later, a pair of independent analyses of NICER and XMM-Newton data constrained the radius of PSR J0740+6620 to be $13.7_{-1.5}^{+2.6}\,{\rm km}$ (Miller et al., 2021) and $12.39_{-0.98}^{+1.30},{\rm km}$ (Riley et al., 2021), providing unique insight into the EOS, given the higher-density material present in the core of such a massive neutron star. The incorporation of XMM-Newton observations is particularly useful when analyzing NICER data on PSR J0740+6620, given the faintness of the source and the crowded field of view, because the imaging data provided by XMM-Newton constrain the phase-averaged stellar flux and the non-stellar astrophysical background (Wolff et al., 2021). Models of the spurious X-ray counts caused by particle interactions with the NICER detectors (e.g., Remillard et al., 2022) can also provide valuable information about the actual astrophysical flux, and can provide similar constraints in practice to the inclusion of XMM-Newton data (e.g., Salmi et al., 2022).

Here we present a new measurement of the equatorial radius of PSR J0740+6620, incorporating previous NICER and XMM-Newton observations as well as an additional $\sim 1.1$ Ms of NICER data. We outline our data selection and processing procedures in Section 2 and describe our methodology in Section 3. The results of our X-ray data analysis are presented in Section 4 and their implications are discussed in Section 5. We summarize our conclusions in Section 6 and report our full posteriors in Appendices B and C. The posterior samples from our analyses are available on Zenodo.

2 Data

The present analysis of the pulse profile of PSR J0740+6620 uses NICER event data collected between 21 September 2018 (ObsID 1031020101) and 21 April 2022 (ObsID 5031020445). These NICER data add roughly two years of data to the set presented in Wolff et al. (2021) and analyzed in Miller et al. (2021) and Riley et al. (2021). The net exposure time of the resulting data set (after application of the filtering procedure described below) was 2733.81 ks, in comparison to the 1602.68 ks net exposure time of the observations analyzed in Miller et al. (2021) and Riley et al. (2021).

We utilized NICER event data processed using the xti20210707 calibration release and HEASoft version 6.30.1 to extract events (Nasa High Energy Astrophysics Science Archive Research Center (2014), Heasarc). NICER event data were filtered to minimize the background signal resulting from high-energy particles interacting with the NICER X-ray timing instrument (XTI). We rejected events obtained at low cosmic-ray cut-off rigidities (COR_SAX $<$ 2). We excluded all events from NICER XTI detector 34 because it exhibited anomalously high levels of noise compared to the others, and in some time intervals we excluded data recorded by detector 14 for the same reason. We did not explicitly exclude data based on the angle between NICER’s pointing direction and the Sun. Rather, we limited our event extractions to time intervals where the maximum “undershoot rate” (the detector reset rate) was below a set threshold. This limited contamination by the accumulation of optical photons (“optical loading”), which leads to spectral broadening of low-energy electronic noise that would otherwise be confined below $\approx 0.2$ keV, and negatively affects both the detector gain and spectral resolution. Holding all other data reduction choices constant, we compared the significances with which pulsations were detected (as quantified by the $Z^{2}_{2}$ test; see Buccheri et al. 1983) for data selected by requiring maximum undershoot rate limits of 50, 100, 150, and 200 cps (counts per second). The 100 cps undershoot rate limit produced the most significant detection and was therefore used in our subsequent analysis. Unlike in Wolff et al. (2021), we did not employ the good time interval (GTI) sorting method of Guillot et al. (2019), which might have biased the resulting radius measurements by up to $\sim 10\%$ (under simplified assumptions; Essick, 2022), but in practice by $\mathrel{\raise 1.29167pt\hbox{$<$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}1\%$ (see also Salmi et al., 2022). Because our data selection procedure is based on the optical-photon-induced detector reset rate, we expect any biases potentially introduced by our consideration of pulsation detection significance to be even smaller than the $\lesssim 1\%$ bias in previous analyses.

We analyzed data from NICER pulse invariant (PI) channels 30 through 123 inclusive, which correspond to a calibrated photon energy range of 0.3 to 1.24 keV. We folded the waveform observed in each NICER energy channel on the spin period of the pulsar to produce a pulse profile at each energy for subsequent analysis. In our analyses incorporating XMM-Newton data, we used the same XMM-Newton data sets and channel ranges as Miller et al. (2021): channels 57 through 299 (inclusive) for the European Photon Imaging Camera (EPIC) pn instrument and channels 20 through 99 (inclusive) of the EPIC MOS1 and MOS2 instruments.

3 Methods

Our approach to modeling the X-ray data is similar in most aspects to the procedures outlined in Miller et al. (2019) and Miller et al. (2021). We review our models of emission from the stellar surface in Section 3.1 and our pulse profile models in Section 3.2. We then describe our parameter estimation procedure in Section 3.3, which is generally similar to that described in Miller et al. (2021).

3.1 Modeling X-Ray Emission from the Stellar Surface

PSR J0740+6620 emits soft X-rays that are modulated on the rotational period of the neutron star (Wolff et al., 2021). The accepted model for soft X-ray emission in old, non-accreting pulsars such as PSR J0740+6620 is that particles with large ($\gg 100$) Lorentz factors, produced by magnetospheric pair cascades (e.g., Timokhin & Harding, 2015; Tolman et al., 2022), penetrate deep into the neutron star’s atmosphere, where their interaction with the atmosphere heats it (Tsai, 1974; Harding & Muslimov, 2001, 2002; Bogdanov et al., 2021). The dependence of the specific intensity of the radiation emerging from the atmosphere on the angle between the direction of propagation and the local normal to the surface depends on the atmospheric composition, which we expect to be pure hydrogen or, less probably, pure helium (see Section 3.1.1); the ionization state (see 3.1.2); and the strength and orientation of the stellar magnetic field, which we assume has a negligible effect on the specific intensity emerging from the atmosphere of PSR J0740+6620 (see 3.1.3).

3.1.1 Upper Atmosphere Composition

Calculations in Alcock & Illarionov (1980) suggest that the strong surface gravity of neutron stars causes the lightest elements present in their atmospheres to rise to the surface within seconds or minutes. Because PSR J0740+6620 is expected to have accreted a substantial quantity of matter from its binary companion to reach its current spin frequency, the prevalence of hydrogen in the outer layers of the envelopes of stars like its companion suggests that PSR J0740+6620’s upper atmosphere is likely to consist of pure hydrogen (see, e.g., Blaes et al., 1992; Wijngaarden et al., 2019, 2020). Although some neutron stars in X-ray binaries are thought to have accreted little or no hydrogen (see, e.g., 4U 1820-30, Strohmayer & Brown, 2002), such binaries have very short periods (685 s for 4U 1820-30, Stella et al., 1987), implying that the neutron star and its companion are very close to one another and that the companion’s envelope has been completely stripped of hydrogen. The orbital period of the system containing PSR J0740+6620 is 4.77 days (Cromartie et al., 2020), suggesting that the pulsar’s companion had a significant hydrogen envelope when its envelope was accreting onto PSR J0740+6620. There is some evidence that the companion of PSR J0740+6620 may now be a cool helium-atmosphere white dwarf (e.g., Beronya et al., 2019; Echeveste et al., 2020). Even if little hydrogen was accreted onto PSR J0740+6620, spallation may have created enough hydrogen to dominate the atmospheric composition (Bildsten et al., 1992; Chang & Bildsten, 2004; Randhawa et al., 2019). For the sake of completeness, we performed analyses using helium (Ho & Heinke, 2009) as well as hydrogen model atmospheres, but consider them less likely to accurately represent the surface composition of PSR J0740+6620 a priori.

3.1.2 Ionization Fraction

The analyses of PSR J0740+6620 presented in Miller et al. (2021), Riley et al. (2021), and Salmi et al. (2022) inferred that the emitting regions of the stellar surface had temperatures $k_{b}T_{\rm eff}\gtrsim 0.07\,{\rm keV}$. In spite of these high temperatures, atmospheric densities may be sufficient to create an appreciable neutral hydrogen fraction via recombination. We therefore performed two sets of analyses using tabulated atmospheric data generated by the NSX code, one assuming complete ionization (Ho & Lai, 2001) and a second allowing for partial ionization (Ho & Heinke, 2009). We report results from both sets of analyses, although we consider the fully ionized models to be more reliable because of systematic uncertainties in the partially ionized hydrogen atmosphere models due to the incomplete coverage of the temperature, density, and energy ranges relevant to neutron star atmospheres in the opacity tables that are currently available (Bogdanov et al., 2021).

3.1.3 Surface Magnetic Fields

Sufficiently strong magnetic fields can significantly affect radiative transfer of X-rays through neutron star atmospheres (e.g., Ho & Lai, 2003; Ho et al., 2003). Using the measured spin period ($P=0.00289\,{\rm s}$) and period derivative ($\dot{P}=1.219\times 10^{-20}$) of PSR J0740+6620 (Cromartie et al., 2020) along with Equation (12) of Contopoulos & Spitkovsky (2006) assuming $\alpha=0$ (i.e., that the closed field line region extends to the light cylinder) and a magnetic field inclination of $\pi/2$, we estimate a surface magnetic field strength of $B\approx 3\times 10^{8}\,{\rm G}$ if the field configuration is a centered dipole. Miller et al. (2021) inferred that the centers of emitting regions were separated on the surface by $\approx 2$ radians in the best-fit model and were thus not antipodal. The chord length between the two spots was therefore $\approx\sin{(1)}\approx 0.84$ times the diameter, suggesting a surface magnetic field strength $\approx 0.84^{-3}\approx 1.7$ times larger than that of a centered dipole, i.e., a surface magnetic field strength $B\approx 5\times 10^{8}\,{\rm G}$. This field strength would have only a small effect on atomic structure (Lai, 2001; Potekhin, 2014) and the corresponding electron cyclotron energy is only $\approx 6\,{\rm eV}$, much smaller than the typical thermal energies we infer for particles in the atmosphere of PSR J0740+6620. Therefore, we use nonmagnetic atmosphere models throughout the present work. However, we caution that we have not rigorously tested this assumption, which would require the construction of grids of model atmospheres with magnetic fields covering the range $B\sim 10^{9}-10^{10}\,{\rm G}$, an especially challenging task due to the comparable strengths of Coulomb and magnetic effects.

3.1.4 X-ray Propagation from Surface to Detector

We model the propagation of light in the vicinity of the pulsar using the oblate Schwarzschild approximation (e.g., Miller & Lamb, 1998; AlGendy & Morsink, 2014; Bogdanov et al., 2019; Silva et al., 2021), which takes into account the oblate shape and Doppler boosts resulting from the rapidity of the spins of millisecond pulsars, but treats the exterior spacetime as Schwarzschild; the oblate Schwarzschild approximation typically results in fractional systematic errors $\mathrel{\raise 1.29167pt\hbox{$<$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}0.1\%$ for pulsars rotating at frequencies of a few hundred Hz (Bogdanov et al., 2019; Silva et al., 2021). We account for the absorption of X-rays in the interstellar medium using the TBabs model (Wilms et al., 2000). We convert our model photon spectra reaching the telescope to model count spectra, as would be registered by the NICER detectors, using a version of the NICER response matrix constructed to match that used over the course of a set of GTIs comprising each data set (see Section 2). The XMM-Newton data are treated similarly, using response matrices generated by the XMM-Newton Scientific Analysis System tools “arfgen” and “rmfgen” relevant to the individual CCDs where the pulsar source image fell during the observations for the EPIC-pn and EPIC-MOS instruments (Gabriel et al., 2004; SAS development team, 2014).

3.2 Pulse Profile Modeling

Our approach to modeling the NICER pulse profile of PSR J0740+6620 is the same as that described in Miller et al. (2021), which built upon the approach discussed in Miller et al. (2019). Detailed descriptions of our pulse profile generation procedure and tests of its accuracy are given in Bogdanov et al. (2019) and Bogdanov et al. (2021). The X-ray pulse profile of PSR J0740+6620 has two clear peaks separated by $\sim 0.4$ in phase (Wolff et al., 2021), and is therefore incompatible with emission from a single circular emitting region of uniform temperature (hereafter we use the term “spot” to refer to X-ray-emitting regions of uniform temperature). While modeling the pulse profile of PSR J0030+0451 necessitated oval (Miller et al., 2019) or crescent (Riley et al., 2019) spots, Miller et al. (2021) and Riley et al. (2021) found that emission from two circular spots, although not antipodal ones, was sufficient to describe the X-ray pulse profile of PSR J0740+6620, with no evidence for more complicated spot geometries; however, we have not rigorously explored this possibility with the present data set.

3.2.1 Model Parameters

Parameter	Definition	Assumed Prior
$c^{2}R_{e}/(GM)$	Inverse stellar compactness at equator	3.2 to 8.0
$M$	Gravitational mass	$\exp{[-(M-2.08\,M_{\odot})^{2}(0.09\,M_{\odot})^{-2}/2]}$
$\theta_{\rm c,1}$	Colatitude of spot 1 center	0 to $\pi$ radians
$\Delta\theta_{1}$	Angular radius of spot 1	0 to 3 radians
$kT_{\rm eff,1}$	Effective temperature of spot 1	0.011 to 0.5 keV
$\theta_{\rm c,2}$	Colatitude of spot 2 center	0 to $\pi$ radians
$\Delta\theta_{2}$	Angular radius of spot 2	0 to 3 radians
$kT_{\rm eff,2}$	Effective temperature of spot 2	0.011 to 0.5 keV
$\Delta\phi$	Longitudinal offset between spots 1 and 2	0 to 1 cycles
$\theta_{\rm obs}$	Observer inclination to stellar rotation axis	1.44 to 1.62 radians
$N_{H}$	Neutral H column density	0 to $2\times 10^{21}\,{\rm cm^{-2}}$
$D$	Distance	$\exp{[-(D-1.136\,{\rm kpc})^{2}(0.20\,{\rm kpc})^{-2}/2]}$, $D\geq 1.136\,{\rm kpc}$
		$\exp{[-(D-1.136\,{\rm kpc})^{2}(0.18\,{\rm kpc})^{-2}/2]}$, $D\leq 1.136\,{\rm kpc}$

Table 1: The primary model parameters, and their corresponding priors, for the pulse profile models considered in this work. Except where noted, the prior is uniform across the given range and zero elsewhere. We place no additional restrictions on the potential values for each parameter.

Having specialized our models to a pair of circular spots, we construct them as follows: each spot ($i=1,2$) is characterized by an effective temperature ($kT_{{\rm eff},i}$), where $k$ is Boltzmann’s constant; an angular radius ($\Delta\theta_{i}$); and a colatitude ($\theta_{{\rm c},i}$). We allow the two spots to have an arbitrary phase offset, measured in rotational cycles from one spot center to the other ($\Delta\phi$). Furthermore, we characterize the star using its gravitational mass ($M$) and inverse stellar compactness ($c^{2}R_{e}/(GM)$), where $R_{e}$ is its equatorial circumferential radius. We also model the observer inclination to the pulsar spin axis ($\theta_{\rm obs}$), the neutral hydrogen column density between NICER and PSR J0740+6620 ($N_{H}$), and the distance from NICER to PSR J0740+6620 ($D$). We allow arbitrary overlaps or lack thereof between spots. This is accomplished by labeling the spots “1” and “2,” and assigning any overlapping regions the effective temperature of the lower-numbered spot. This allows, for example, the creation of a single crescent-shaped spot by covering a hotter spot with another so cool that the latter emits virtually no radiation.

The full set of parameters used in our analyses is listed in Table 1, along with the priors we assume for each parameter. For most parameters we use a flat prior, defined as having a uniform nonzero probability density between an upper and lower bound (inclusive) and zero probability density outside of those bounds. Because PSR J0740+6620 is in a binary system, we are able to incorporate additional informative priors. Specifically, Fonseca et al. (2021) reported based on Shapiro delay measurements that the mass of PSR J0740+6620 was $2.08_{-0.069}^{+0.072}\,M_{\odot}$, that the distance to the pulsar was $1.136_{-0.152}^{+0.174}\,\rm{kpc}$, and that the angle between our line of sight and the orbital axis of the system was $87.5^{\circ}\pm 0.17^{\circ}$, all at 68% credibility.

We use these radio-timing-derived measurements to place priors on the mass, distance, and observer inclination in our analysis. We implement a Gaussian prior on the gravitational mass of the pulsar, with a median of $2.08\,M_{\odot}$ and standard deviation of $0.09\,M_{\odot}$, where we have linearly added an estimate of the systematic error in the Shapiro delay measurement (E. Fonseca, personal communication). We implement an asymmetric quasi-Gaussian distance prior, which has a probability distribution that falls off at different rates above and below the mode, according to the functional form listed in Table 1. This distribution is based on the results presented in Fonseca et al. (2021) and includes an additional, linearly added estimate of the systematic error of 0.03 kpc (E. Fonseca, personal communication).

Radio observations are able to constrain the angle between the observer’s line of sight and the orbital axis of the binary system, but our analysis depends on the angle between the observer’s line of sight and the rotational axis of PSR J0740+6620. Accretion from their companion stars is thought to be the primary mechanism for the spin-up of millisecond pulsars such as PSR J0740+6620 (e.g., Bhattacharya & van den Heuvel, 1991). Accordingly, the spin axis of the pulsar should gradually align with the orbital axis of the system. However, it is unlikely that the alignment is perfect, and thus the pulsar spin axis may be tilted a few degrees with respect to the orbital axis of the binary system. Thus, we have adopted a flat prior on the angle between the observer’s line of sight and pulsar spin axis centered at $87.5^{\circ}$, which is the best estimate of the orbital inclination of the binary, with a width of $5^{\circ}$.

3.2.2 Pulse Profile Models

Given values for the set of model parameters described in Section 3.2.1, a stellar rotation period (which is known to high precision from radio observations, e.g., Cromartie et al. 2020; Fonseca et al. 2021), a neutron star atmosphere model, and a specified observation time, we can generate a spots-only pulse profile. However to compare these model pulse profiles with data, it is necessary to account for non-spot (or ‘background’) contributions to the observed X-ray counts, which are not modulated on the spin period of the pulsar. To compare with NICER data it is also necessary to select a starting rotational phase for the pulse profile. We are then able to calculate the log likelihood of the available data sets given a particular model. When we incorporate information from various XMM-Newton instruments, the overall log likelihood is simply the sum of the log likelihoods for the data given the model for each data set, i.e.,

$$\begin{split}\log{\mathcal{L}_{\rm tot}}=\log{\mathcal{L}_{\rm NICER}}+\log{\mathcal{L}_{\rm XMM-pn}}\\ +\log{\mathcal{L}_{\rm XMM-MOS1}}+\log{\mathcal{L}_{\rm XMM-MOS2}}\end{split}$$ (1)

for analyses of NICER data along with XMM-Newton imaging data from the pn, MOS1, and MOS2 instruments, or $\log{\mathcal{L}_{\rm tot}}=\log{\mathcal{L}_{\rm NICER}}$ for analyses of NICER data alone. The straightforward summation of the terms in Equation (1) assumes that the different terms are uncorrelated, which is valid for data from NICER and the different instruments mounted on XMM-Newton (Turner et al., 2001; Strüder et al., 2001).

We compute the log likelihood corresponding to each data set by summing the Poisson log likelihoods of the data given the model for each energy channel ($E_{j}$) and phase ($\phi_{k}$). Given an observed number of counts $d_{jk}$, which is a non-negative integer, and a predicted number of counts $m_{jk}$, which is a positive real number, the Poisson likelihood of the data given the model is $m_{jk}^{d_{jk}}e^{-m_{jk}}/d_{jk}!$. However, the factor $1/d_{jk}!$ is shared by all models, and can thus be neglected. The log likelihood we use is therefore $\log{\mathcal{L}}=\sum_{\phi_{k}}\sum_{E_{j}}(d_{jk}\log{m_{jk}}-m_{jk})$, where NICER data are analyzed over 32 phases and XMM-Newton data that we use are resolved only in energy, with a single time-averaged phase.

As previously mentioned, observed pulsar X-ray pulse profiles may include contributions not modulated on the spin period of the pulsar, from sources such as particle background, optical loading, and both resolved and unresolved sources in the telescope field of view, among others. We analytically marginalize over a phase-independent background parameter for each NICER energy channel, as described in section 3.4 of Miller et al. (2019). Concerning the XMM-Newton background, we take previously measured XMM-Newton blank-sky background spectra to be a Poisson realization of the actual XMM-Newton background, and calculate the probability of the data given the spot counts and the distribution of possible background counts, as described in Section 3.4.2 of Miller et al. (2021).

3.3 Bayesian Parameter Estimation

We adopt a hybrid sampling approach that utilizes an initial suite of nested sampling (Skilling, 2004) analyses followed by a suite of Markov chain Monte Carlo (MCMC) analyses for the purpose of refining our posterior inferences. To wit, each of the initial nested sampling analyses, which we carry out using the MultiNest package and its Python bindings (Feroz et al., 2009; Buchner, 2016a), is able to sample the parameter space globally, identify multiple modes if they exist, and approximate the Bayesian evidence for each mode. The evidence ($\mathcal{Z}$), which is the average of the likelihood $\mathcal{L}$ over the normalized prior $\uppi$, is given by the following integral over the model parameters $\mathbf{\theta}$

$$\mathcal{Z}=\int\mathcal{L}(\mathbf{\theta})\uppi(\mathbf{\theta})d\mathbf{\theta},$$ (2)

and can be used to compare models, where the model with the highest evidence is preferred.

Basic nested sampling algorithms such as MultiNest, while often relatively fast, produce a very limited number of posterior samples during each analysis and provide neither theoretical nor practical convergence guarantees. Therefore, following each nested sampling analysis we used the posterior probability distribution inferred from that analysis to initialize an MCMC analysis, which we continued until the posterior distribution appeared to be stationary and we were able to obtain approximately $10^{6}$ independent posterior samples from the stationary distribution.

3.3.1 Nested Sampling Analyses

We began our analysis using the publicly available MultiNest algorithm (Feroz et al., 2009; Buchner, 2016a), which begins by randomly sampling a number of points (a user-specified number of “live points,” $N_{\rm live}$) from the prior. Subsequently, MultiNest sequentially replaces the lowest-likelihood point with a higher-likelihood point. Proposals for these higher-likelihood points are drawn from within the region bounded by an approximation of the isolikelihood surface corresponding to the likelihood of the point to be replaced. A larger number of live points will lead to more thorough sampling of the parameter space, all other parameters being equal. MultiNest constructs approximate isolikelihood surfaces using multiple hyper-ellipsoids constructed to envelop the set of live points. If these minimal bounding hyper-ellipsoids contain a volume smaller than the expectation value of the remaining prior volume, the bounding ellipsoids are expanded until the overlap-corrected enclosed volume matches the expected remaining prior volume (Feroz et al., 2009, Section 5.2). However, if the true isolikelihood surface is not described accurately by the (possibly expanded) hyper-ellipsoids encompassing a given set of live points, MultiNest can draw samples from a biased region of parameter space, resulting in biased estimates of the posterior and Bayesian evidence (e.g., Buchner, 2016b; Nelson et al., 2020; Buchner, 2023; Lemos et al., 2023; Dittmann, 2024). The MultiNest algorithm attempts to ameliorate this shortcoming through an ‘efficiency’ parameter, which is the inverse of a factor used to further increase the hypervolume of the MultiNest bounding ellipsoids (for example, an efficiency value of $\epsilon=0.1$ multiplies by 10 the target volume for the expanded hyper-ellipsoids). However, this procedure is still not guaranteed to encompass the intended isolikelihood surface.

Previous analyses of NICER data have found that MultiNest rarely produces converged posterior estimates. For example, Miller et al. (2021) reported that a MultiNest analysis of an earlier PSR J0740+6620 data set using $N_{\rm live}=1000$ and $\epsilon=0.01$ systematically underestimated the stellar radius and its uncertainty. Analyzing the same data set, Riley et al. (2021) found that the settings $\epsilon=0.1$ and $N_{\rm live}=4\times 10^{4}$ underestimated the width of the inferred radius posterior compared to a run with the same efficiency and twice as many live points. Moreover, Miller et al. (2021) showed, by performing MultiNest inferences and comparing the number of parameters inferred to fall within the $\pm(1,2,3)\sigma$ credible intervals with statistical expectations, that MultiNest systematically underestimated the width of credible intervals, although less significantly when more live points and lower sampling efficiencies were used. Analyzing NICER data for PSR J0030+0451, Vinciguerra et al. (2024) found that analyses using $\epsilon=0.3$ and $N_{\rm live}=10^{3}$ often overestimated the median stellar radius and underestimated the width of the posterior distribution compared to higher-resolution analyses using $\epsilon=0.1,~{}N_{\rm live}=10^{3}$ and $\epsilon=0.3,~{}N_{\rm live}=10^{4}$; in some cases the posteriors estimated by lower-resolution analyses strongly excluded the median from higher-resolution analyses (Vinciguerra et al., 2024). Despite these shortcomings, MultiNest is in our experience able to provide a suitable starting point for subsequent MCMC analyses using orders of magnitude fewer computational resources than would initializing an MCMC analysis using uniform samples drawn from the prior.

Based on these considerations and a series of convergence tests, our default MultiNest settings were $\epsilon=0.01$ and $N_{\rm live}=4096$. In each analysis we enabled MultiNest to evaluate multiple modes. Because our parameter inferences do not rely on the capacity of MultiNest to estimate posteriors with high accuracy, these settings were sufficient for our preliminary nested sampling analysis.

3.3.2 Markov Chain Monte Carlo Analyses

Following the completion of each nested sampling analysis, we used the posterior probability distribution estimates it produced to draw initial walker positions for an MCMC analysis using the emcee package (Foreman-Mackey et al., 2013) and the affine-invariant ‘stretch’ proposal of Goodman & Weare (2010). Because the proposal distribution utilized used in this analysis satisfies detailed balance, the distribution of walker positions will converge to and provide samples from the stationary posterior distribution. In principle the convergence of each chain of walkers to the stationary distribution may take an arbitrarily long time if the walkers are initialized very far from equilibrium, which is why we use a prelimiary nested sampling analysis to generate initial walker positions that more accurately approximate the posterior distribution.

Each analysis used 4096 walkers, for which we drew initial positions by re-sampling a Gaussian kernel density estimate of the corresponding MultiNest-derived posterior distributions. We performed MCMC sampling until we had collected $\sim 10^{6}$ effective samples (in practice $\sim 10^{7}$ samples, because our analyses had proposal acceptance fractions of $\sim 0.1$ and autocorrelation times of $\sim 10$ iterations). We judged the sampling to be converged when over the aforementioned $\sim 10^{7}$ MCMC iterations we observed no secular variation in the 1st, 16th, 50th, 84th, or 99th percentiles.

4 Pulse Profile Modeling Results

(a) Lunate spot configuration

(b) Dual-spot configuration

Figure 1: Examples of the two spot configurations identified during our analysis of the updated data set for PSR J0740+6620 able to produce good fits to the NICER data. The black circle denotes the colatitude of the observer line of sight. Dashed lines indicate features on the side of each sphere opposite the observer. To reiterate, in our models regions that overlap emit with the temperature of the ‘lowest-number’ spot; accordingly, emission from the configuration shown in panel (a), for which the temperature of ‘spot 1’ tends towards the lower end of our prior, effectively comes from the protruding crescent of ‘spot 2.’

We first present in Section 4.1 the different modes identified during the initial stage of our analysis, assess the thoroughness of that preliminary nested sampling stage, and report approximate evidences for each atmosphere model and data set analyzed. Subsequently, we report in Section 4.2 the mass and radius inferences from each data set and atmosphere model combination. We assess the goodness of fit of our models in Section 4.3, and compare our results with those of Salmi et al. (2024) in Section 4.4.

4.1 Inferred Emitting Region Geometries

Dataset	Atmosphere	Lunate spot	Dual spot	$\Delta\log{Z}$
NICER	$\rm H_{full}$	$\checkmark$	-	-
NICER	$\rm H_{partial}$	$\checkmark$	$\checkmark$	$-0.11$
NICER	$\rm He_{full}$	$\checkmark$	-	$1.33$
NICER+XMM	$\rm H_{full}$	-	$\checkmark$	-
NICER+XMM	$\rm H_{partial}$	-	$\checkmark$	$1.0$
NICER+XMM	$\rm He_{full}$	-	$\checkmark$	$2.08$

Table 2: We present the inferred spot configuration(s) and log Bayes factors for each data set and each atmosphere model considered in this work. In the second column, subscripts indicate whether each atmosphere model assumed full ionization or allowed for partial ionization. Check marks ($\checkmark$) indicate which spot configurations provided good fits to each data set under the assumptions of a particular atmosphere model. For each data set (NICER and NICER+XMM), we report the evidence for each model (including all identified modes) relative to the evidence obtained using that data set assuming the fully ionized hydrogen atmosphere, which we consider to be the most astrophysically relevant atmosphere (see Section 3).

Miller et al. (2021) identified two emitting region geometries that were capable of reproducing the observed pulse profile and spectrum of PSR J0740+6620: one was a crescent (or lunate) configuration formed by an extremely cool circular spot eclipsing a hotter spot; the other consisted of two fairly small spots with similar temperatures, well-separated in azimuth. Both spot configurations are also allowed by the priors used in this work (Table 1).

Figure 1 shows ilustrative schematic examples of both configurations. Miller et al. (2021) found that the configuration with two fairly small spots was favored by a Bayes factor of $\sim 3000$ over the lunate configuration. When only the updated NICER data are considered, we find that either spot configuration can provide good fits, depending on the assumed atmospheric composition and whether or not we also consider the XMM-Newton data.

Table 2 shows the spot configurations that provided good fits to each data set under each assumed atmosphere model. We also include in Table 2, for each data set, the log evidence for each model relative to the fit assuming a fully ionized hydrogen atmosphere to the same data set. For example, the logarithm of the Bayesian evidence for our fully ionized helium fit to only the NICER data was $\approx 1.33$ larger than the logarithm of the Bayesian evidence for the fully ionized hydrogen fit to the same data set. Given the sensitivity of the evidence to priors, and the potential for appreciable systematic errors in evidences estimated by MultiNest, no atmosphere model is significantly favored over another.

When only the available NICER data are considered, the preferred mode depends on the assumed atmospheric composition and conditions: while the fully ionized hydrogen and helium atmosphere models overwhelmingly favor the lunate configuration, the partially ionized hydrogen atmosphere models favor both configurations. However, when the XMM data—which constrain the stellar count rate and non-stellar background—are included, the data prefer models with two smaller, near-antipodal spots.

Figure 2: A comparison of the $M$, $R_{e}$, and joint $M-R_{e}$ posterior probability distributions for PSR J0740+6620 inferred by analyzing the NICER data only (left panels) and the NICER and XMM-Newton data analyzed jointly (right panels). The results shown here assume a fully ionized, purely hydrogen atmosphere, which we consider the most likely atmosphere for the reasons explained in Section 3.1. The solid curve in each bottom-left panel shows the one-dimensional probability distribution for the stellar mass inferred from these data sets; the mass prior derived from Fonseca et al. (2021) is shown by the dashed gray line. The color shading in the two-dimensional $M$–$R_{e}$ plots illustrates the shapes of the minimum-area credible regions; the overplotted white lines show the $1\sigma$ and $2\sigma$ contours outside which $\sim 31.7\%$ and $\sim 4.6\%$ of the posterior mass of each distribution lies. The solid curve in the bottom-right panel on the left shows the probability density of the equatorial radius inferred by analyzing only the NICER data. The solid curve in the bottom-right panel on the right shows the probability density inferred by jointly analyzing the NICER and XMM-Newton data; the black dotted curve in this panel shows the probability density inferred by analyzing the NICER data alone.

Before attempting to sample any of the posteriors more thoroughly, we first verified that our initial nested sampling analyses were thorough enough to identify the most significant mode, and that each did not find one mode or the other only by chance. The first test to which we subjected our nested sampling results was to evaluate the best fit from an analysis that identified only one mode using the data and atmospheric model of analyses which found other modes. Taking the best fits from our fully ionized hydrogen atmosphere analyses as an example, when we apply the (lunate) best fit parameters from our analysis of the NICER data alone to the joint NICER and XMM-Newton data set, we find a log likelihood $\sim 400$ lower than for the (dual-spot) best joint fit to the NICER and XMM-Newton data. Similarly, the log likelihood of the best fit found analyzing the NICER data alone assuming a helium atmosphere was $\sim 600$ lower than the log likelihood of the best fit found in the joint analysis of NICER and XMM-Newton data. This analysis thus suggests that the XMM-Newton data strongly disfavor the lunate spot configuration. For comparison, the model that best fits the NICER data only, using the partially ionized hydrogen atmosphere tables, which is the smaller near-antipodal spot configuration, gives a NICER+XMM log likelihood only $\sim 35$ lower than the best fit for this model found in the joint NICER and XMM-Newton analysis.

Based on the comparisons above, and additional tests discussed in Appendix A, we believe that our MultiNest analyses have sampled parameter space thoroughly enough to identify the relevant modes in each analysis. Because we are able to sample each posterior further using MCMC techniques, we do not need to rely on the often-questionable convergence of the posterior estimates reported by MultiNest. Instead, we require only that the different modes identified by MultiNest be assigned roughly the correct weights. In practice, the two modes present in $\Delta\phi$ for the smaller-spot configurations (see Appendix C), which are a result of the near reflection symmetry due to the observer’s inclination being so close to the rotational equator of the star, have very little effect on the radius or mass posteriors. Thus, the only capability of MultiNest on which we rely is its ability to differentiate between the lunate and near-antipodal configuration for each data set and atmosphere model, a capability our analyses have demonstrated.

4.2 Mass and Radius Constraints

After we completed our preliminary nested sampling analyses, we initialized a series of MCMC analyses to thoroughly sample each posterior. We discuss the results in this section.

Unsurprisingly, given the faintness of PSR J0740+6620, we are able to infer little about the stellar mass, although we systematically infer masses very slightly lower than the median value of our prior (e.g., $2.06\,M_{\odot}$ rather than $2.08\,M_{\odot}$). We present the credible regions for the equatorial radius of PSR J0740+6620 derived from each of our analyses in Table 3, the posterior probability distributions of the mass and equatorial radius derived from our analyses assuming a fully ionized hydrogen atmosphere in Figure 2, and a comparison of the present constraints on the equatorial radius to those derived in Miller et al. (2021) in Figure 3.

We find that, when the constraints provided by XMM-Newton data are included, the radii we infer from analyses utilizing differing neutron star atmosphere models are consistent with one another: the $-1\sigma$ values range only from $11.76\,\rm km$ to $11.82\,\rm km$ and the $+1\sigma$ values range from $15.01\,\rm km$ to $15.52\,\rm km$. In contrast, when we analyze the NICER data alone, we find that the constraints on the stellar radius vary much more between our analyses. This is largely due to the faintness of PSR J0740+6620 and the presence of other sources within the field of view of NICER when it is observing PSR J0740+6620, as shown in Wolff et al. (2021). Because the number of phase-dependent counts in the pulse profile of PSR J0740+6620 is much smaller than the number of phase-independent counts, the NICER data alone allow a wide range of emitting region geometries and values of the stellar compactness, which limits the precision of the resulting constraints on the radius of the pulsar. When the XMM-Newton imaging data are included, the stellar and background fluxes are better constrained.

Table 3: Summary of Inferred Equatorial Radii

Dataset	Atmosphere Model	$-2\sigma$	$-1\sigma$	median	$+1\sigma$	$+2\sigma$
NICER	$\rm{H_{full}}$	9.60	10.54	12.48	16.54	21.84
NICER	$\rm{H_{partial}}$	9.72	10.59	11.95	14.46	18.32
NICER	$\rm{He_{full}}$	9.62	10.56	12.72	17.39	22.55
NICER+XMM	$\rm{H_{full}}$	10.99	11.79	12.92	15.01	18.57
NICER+XMM	$\rm{H_{partial}}$	10.94	11.76	12.99	15.36	19.75
NICER+XMM	$\rm He_{full}$	10.98	11.82	13.07	15.52	20.50

Note. — A comparison of the $-2\sigma$, $-1\sigma$, median, $+1\sigma$, and $+2\sigma$ constraints on the equatorial radius of PSR J0740+6620 (measured in kilometers), using different data sets and atmospheric models. When we include only the NICER data, the inferred stellar radius and especially its upper limit depend on the model atmosphere assumed. However, when we also include the XMM-Newton data, which constrain the spectra of the stellar flux and the non-stellar background, we find more consistent results across different assumed model atmospheres.

A general trend illustrated in Table 3 and Figure 3 is that in each analysis, when the XMM-Newton data are included the lower limits on the stellar radius increase. Because the XMM-Newton observations constrain the number of non-stellar counts in the star’s pulse profile the inferred fractional modulation in the profile increases when the XMM-Newton data are included, placing a strong upper bound on the stellar compactness. The XMM-Newton data also constrain the total flux from the stellar surface and its spectrum, and when we assume the atmosphere is fully ionized hydrogen or helium we find that, compared to fitting the NICER data alone, including the XMM-Newton data places more stringent upper limits on the stellar radius. As shown in Appendices B and C, the maximum-likelihood radius is also slightly reduced. When we instead assume the atmosphere is partially ionized hydrogen, we find that the upper limits on the radius are more lenient when the XMM-Newton data are included, as previously found in Miller et al. (2021), which used these same model atmospheres.

Figure 3: Radius posteriors obtained when only the NICER data are used (top panel) and when the NICER and XMM-Newton data are both used (bottom panel). The solid blue lines show the posterior for a fully ionized hydrogen atmosphere whereas the dashed orange lines are for a partially ionized hydrogen atmosphere. The dotted green lines show the posterior reported by Miller et al. (2021), who assumed a partially ionized hydrogen atmosphere. The results obtained here using only the NICER data depend on the atmospheric model, because different spot patterns are favored by the different models (see Figure 1). However, when both the NICER and XMM-Newton data are used, the posterior is almost independent of the atmospheric model assumed.

If we compare our results to those presented in Miller et al. (2021), temporarily restricting ourselves to the model atmospheres used in that work, we find a $68\%$ credible interval for the equatorial radius that is narrower (only $83\%$ as wide) when the NICER and XMM-Newton data are jointly fit. On the other hand, when we fit the NICER data alone we find no improvement over Miller et al. (2021). We have not rigorously determined why this is the case. One clear difference is that our current analyses have identified support for lunate modes in our NICER-only analyses, which tend to support larger-radius configurations (see Tables 2 and 3, and Appendices B and C). Miller et al. (2021) did not find evidence for these modes. We find that including the XMM-Newton data excludes models with those modes, possibly due to the information about the total stellar X-ray flux provided by these data. This may explain why the analyses that include the XMM-Newton data find lower upper limits on the equatorial radius. A possible but less likely explanation is the different data selection choices made in this work.

We find that fits that assume a fully ionized helium atmosphere infer slightly larger values of the stellar radius than the corresponding fits that assume either a partial or fully ionized hydrogen atmosphere. Similar though much stronger trends have been noted in analyses of the PSR J0030+0451 NICER data (Miller et al., 2019; Salmi et al., 2023). These trends are likely due to differences between the emergent spectra and the beaming patterns produced by the two model atmospheres, which are larger for neutron stars that emit X-rays from larger fractions of their surfaces. The previous analyses of the NICER data on PSR J0740+6620 assuming an ionized helium atmosphere that were reported in Riley et al. (2021) and Salmi et al. (2023) did not show any appreciable differences between the results for the helium and hydrogen atmospheres used in these works. However, both works assumed an a priori upper limit on the stellar radius of $16\,\rm km$, which reduced the variation of the upper limits they derived from their posteriors.

4.3 Adequacy of the Fits

Figure 4: The value of $\chi_{i}=(m_{i}-d_{i})/\sqrt{m_{i}}$ (where $m_{i}$ is the predicted number of model counts and $d_{i}$ is the number of observed counts in phase-channel bin $i$) for each of the 3008 NICER phase-energy bins considered in this study (32 bins in rotational phase and 94 energy channels), for the energy-resolved pulse profile model that best fits the NICER and XMM-Newton data, assuming a fully ionized hydrogen atmosphere. No patterns in the values of $\chi$ are evident, as expected of a good fit. For this fit, $\chi^{2}/{\rm dof}=2818.8/2901$, where $\chi^{2}=\sum\chi^{2}_{i}$. The probability of finding $\chi^{2}\geq 2818.8$ for 2901 degrees of freedom is $\sim 86$% if the model is correct.

As in Miller et al. (2019) and Miller et al. (2021), we performed several tests to assess the adequacy of the fits we obtained. Here we consider only the maximum-likelihood parameters from our joint analysis of the NICER and XMM-Newton data, assuming a fully ionized hydrogen atmosphere.

We first performed a standard $\chi^{2}$ test comparing our best-fit phase- and energy-resolved pulse profile model the NICER data and found the ratio of the resulting $\chi^{2}$ to the number of degrees of freedom of 2818.8/2901. The probability of finding a value of $\chi^{2}$ this high or higher for this many degrees of freedom using the correct model is 86%. Figure 4 shows the resulting signed residuals in each phase and energy bin. There is no evidence for clustering or systematic trends. Although this test cannot show that the model being used is correct, this value of $\chi^{2}$ shows that this model is not obviously deficient.

We have also checked the ability of our analysis to reproduce the bolometric pulse profile. This is a nontrivial test, because our Bayesian analysis fits an energy-resolved pulse profile model to the energy- and phase-resolved NICER data. Because the $\chi^{2}$ statistic is nonlinear, it is possible to significantly over- or under-estimate the total flux at a particular phase, even though the pulse-phase–energy-channel residuals are negligible, causing the residuals obtained by comparing the model bolometric pulse profile to the observed bolometric profile to be significant. We find that the $\chi^{2}$ per degree of freedom obtained by comparing this bolometric pulse profile model to the observed bolometric pulse profile is $28.8/27$. The probability of finding a value of $\chi^{2}$ this high or higher given the correct model is $\sim 37\%$. We show the bolometric fit in Figure 5, which also shows the size of the background relative to the X-ray emission from the pulsar.

Figure 5: Top: A comparison of the bolometric NICER pulse profile (in 32 bins, plotted in blue) to the bolometric pulse profile from our best fit to the combined NICER and XMM-Newton data set, shown in orange. The dashed green line plots the unmodulated background in our best fit, which was added to the pulse profile generated by the two hot spots, as described in Section 3.2.2. Bottom: bolometric residuals, $\chi$, as a function of phase. For our best fit we find a bolometric $\chi^{2}$/dof of 28.8/27, which has a probability of 37% if the model is correct.

Figure 6 compares the observed energy spectra predicted by our model of a fully ionized hydrogen atmosphere that best fits the combined NICER and XMM-Newton data with the spectra observed by the MOS1, MOS2, and pn detectors on XMM-Newton in the energy channels used in our analysis. Because very few counts were recorded in each energy channel, we evaluated the quality of the fit by generating a set of $10^{5}$ synthetic data sets by Poisson sampling the spectrum predicted by our best-fit model in each channel of each instrument. We then computed the log likelihood for each synthetic data set and instrument, and compared it to the log likelihood of the actual data given our model. We found the total log likelihood of the real data was at the 78th percentile, that the log likelihood of the real pn data alone was at the 99th percentile, and that the log likelihoods of the real data for MOS1 and MOS2 were at the 15th and 4th percentiles, respectively. The overall fit to the XMM data is therefore good, the fit to the pn data is anomalously good, and the fits to the MOS1 and MOS2 data are relatively bad. Miller et al. (2021) observed the same trend, although to a lesser extent.

Figure 6: Comparison of the observed energy spectra predicted by our model of a fully ionized hydrogen atmosphere that best fits the combined NICER and XMM-Newton data with the spectra observed by the MOS1, MOS2, and pn detectors on XMM-Newton in the energy channels used in our analysis. Each panel shows the spectrum observed by the indicated XMM-Newton detector in blue as a histogram and the spectrum predicted by the best-fit model in orange. The model provides acceptable fits to all three observed spectra.

4.4 Differences Between Our Results and Those of Salmi et al.(2024)

An independent analysis of the data analyzed in this paper is presented in Salmi et al. (2024), which reports a radius measurement of $12.49^{+1.28}_{-0.88}$ km (68% credibility). This radius measurement overlaps with our measurement of $12.92_{-1.13}^{+2.09}$ km, particularly at the $-1\sigma$ point in the posterior distribution, which is especially important for constraining the dense matter EOS, because larger radii are disfavored by the current gravitational wave observations of binary neutron star mergers (e.g., Abbott et al., 2017; De et al., 2018). However, substantial discrepancies remain between the median values and uncertainties in these two measurements. As we demonstrate below, these discrepancies are due largely to the different methods used to sample the posterior, the different priors assumed for the equatorial radius (Salmi et al. 2024 imposed a hard upper limit of 16 km on the equatorial radius), and the different methods used to model the XMM-Newton data. For reference, the results of both our preliminary nested sampling analyses and our converged MCMC analyses are provided in Appendix A.

4.4.1 Analyses of the NICER Data Alone

For our initial, preliminary analysis of the data we used MultiNest, which Salmi et al. (2024) use exclusively for their parameter estimation. We can therefore make a nearly apples-to-apples comparison of some of our preliminary results with their final results in some cases. As we now discuss, our results are consistent with those of Salmi et al. (2024) when we use the same sampling methodology, consider only the NICER data, and employ the same radius prior.

Consider first analyses that assume a fully ionized hydrogen atmosphere. Our preliminary, MultiNest analysis of only the NICER data gave an estimate for the equatorial radius of $12.15^{+2.3}_{-1.51}~{}{\rm km}$; the corresponding absolute and fractional widths of the $\pm 1\sigma$ interval are $3.81$ km and 31%. The final radius estimate quoted in Salmi et al. (2024) for this same case is slightly smaller, namely, $11.29^{+1.13}_{-0.81}$ km, while the absolute and fractional widths of the $\pm 1\sigma$ interval are almost a factor of two smaller, namely, 1.94 km and 17%. If we now discard all posterior samples with $R_{e}>16\,\rm km$ to mimic the radius prior assumed in Salmi et al. (2024), our radius estimate becomes $11.98^{+1.9}_{-1.4}$ km, and the absolute and fractional widths of the $\pm 1\sigma$ interval become 3.28 km and 27%. We note, however, that the lunate spot mode was enormously favored over the dual-spot configuration for this data set and atmosphere model in our analysis, and that the analysis of Salmi et al. (2024) was precluded from finding this mode by their constraint that each spot have an angular radius less than $\pi/2$ and that the spots were not allowed to overlap.

Consider now the preliminary, MultiNest results we obtained assuming a partially ionized hydrogen atmosphere and the final results Salmi et al. (2024) obtained assuming a fully ionized hydrogen atmosphere. If we again discard from our analysis all samples with $R_{e}>16$ km, our radius estimate becomes $11.47^{+1.2}_{-0.91}$ km, and the absolute and fractional widths of the $\pm 1\sigma$ interval become 3.28 km and 27%. If we then also discard all samples with spot angular radii larger than $0.4$ radians, which effectively eliminates the lunate spot mode, we obtain a radius estimate of $11.36^{+0.97}_{-0.81}$ km and the absolute and fractional widths of the $\pm 1\sigma$ interval become 1.78 km and 16%, in agreement with the result reported by Salmi et al. (2024).

4.4.2 Analyses of the NICER and XMM-Newton Data

Our preliminary joint analysis of the NICER and XMM-Newton data using MultiNest and assuming a fully ionized hydrogen atmosphere gave a radius estimate of $12.70^{+1.48}_{-0.97}$ km. Using a prior as above to exclude radii greater than 16 km reduced the estimated radius to $12.66^{+1.37}_{-0.94}$ km. Because the lower bounds on the radius provide the most information about the EOS, it is useful to compare the value of the radius at the $-1\sigma$ point in the posterior distribution obtained by various analyses. In our preliminary analysis we obtained, whereas Salmi et al. (2024) found 11.61 km when analyzing the same data. Based on the above comparison of analyses that use only the NICER data, this difference, as well as some of the difference in the reported widths of the credible regions, appears to be due largely to differences in the treatment of the XMM-Newton data. As one example, we treated the distribution of blank-sky counts in each energy channel as a realization of a Poisson process, whereas Salmi et al. (2024) followed Riley et al. (2021) in assigning a flat prior for the distribution of blank-sky counts in each channel with bounds chosen to be a few times larger than the statistical uncertainty in the number of blank-sky counts, clipped to maintain positivity. Although notable, the differences in the assumed radius prior and in the treatments of the XMM-Newton data are insufficient to explain the difference between the width of our radius posteriors and those of Salmi et al. (2024).

Another significant difference between our analysis and that of Salmi et al. (2024) is the statistical sampling method employed. Although MultiNest is often suitable for exploring the parameter space, numerous studies have shown that it can be unreliable when used for parameter estimation, as discussed in Section 3.3.1. Appendix A shows that using MultiNest for parameter estimation tends to underestimate the widths of radius credible regions, and that these are typically biased towards low radii.

To derive converged results, we performed MCMC sampling, using our initial MultiNest results only to select the initial positions of the walkers. Over the course of the MCMC sampling, the radius credible interval expanded from the initial, MultiNest-derived interval of $12.70^{+1.48}_{-0.97}$ km to $12.92^{+2.09}_{-1.13}$ km. In section 4.3 of Salmi et al. (2024), the authors report an additional MultiNest analysis that used a lower sampling efficiency ($\epsilon=10^{-3}$) than the efficiency they used to obtain the headline result featured in the abstract of their paper. A lower efficiency allows MultiNest to better sample the parameter space. The additional analysis gave a radius of $12.55^{+1.37}_{-0.92}$, in better agreement with our result than the headline result, despite the different treatments of the XMM data and radius prior in these two analyses.

The headline result of Salmi et al. (2024) ($R_{e}\approx 12.49^{+1.28}_{-0.88}$ km) is only very slightly narrower than the value reported in Riley et al. (2021) ($R_{e}\approx 12.39^{+1.30}_{-0.98}$ km), despite the inclusion of about 50% more data. Adopting instead the results presented in Section 4.3 of Salmi et al. (2024), which used a MultiNest sampling efficiency that is expected to better sample the parameter space, their radius estimate using the present data is $12.55^{+1.37}_{-0.92}$ km, very slightly broader than the result reported in Riley et al. (2021). In contrast, our radius estimate improved from $13.71^{+2.61}_{-1.50}$ km in Miller et al. (2021) to $12.92^{+2.09}_{-1.13}$ km in the present work. Overall, the differences between the results obtained in different analyses have diminished as NICER has accumulated more data.

5 Discussion

As discussed in Section 2, we have analyzed an additional $1.1$ Ms of NICER data on PSR J0740+6620, in addition to the NICER and XMM-Newton data analyzed in Miller et al. (2021). Analysis of these NICER and XMM-Newton data improves the precision of the constraint on the radius of this neutron star from $13.71^{+2.61}_{-1.50}$ km to $12.92^{+2.09}_{-1.13}$ km, reducing the fractional uncertainty in the radius from $\sim 30\%$ to $\sim 25\%$. This more precise measurement of the radius incrementally improves the constraints on the EOS of dense matter. We expect the EOS to be further constrained by analysis of the additional data on PSR J0740+6620 that is currently being collected. Although measurements of the radius of PSR J0740+6620 made using NICER still depend somewhat on the procedure used to analyze the data (see Section 4.4), the consequences of using different procedures have become less significant as the amount of NICER data has increased.

5.1 Equation of State Constraints

Our current radius measurement is more precise than that of Miller et al. (2021), which was based on the data collected through 2020 (see Miller et al. (2021)). For example, the value of the equatorial radius at the $+1\sigma$ point in the posterior distribution has decreased from 16.32 km to 14.34 km, and the radius width at $\pm 1\sigma$ is now 3.22 km, which is 78% of the 4.12 km width reported in Miller et al. (2021). This improvement is almost exactly what is expected under the assumption that the radius uncertainty scales as the inverse square root of the exposure time, because $\sqrt{1602.68~{}{\rm ks}/2733.81~{}{\rm ks}}\approx 0.77$.

However, this decrease has a relatively minor effect on our EOS inferences, because in making these inferences we also include other results that bear on the properties of cold, catalyzed matter at densities above the saturation density of nuclear matter. In particular, the upper bound placed on the tidal deformability of neutron stars with moderate masses by the observations of GW170817 (Abbott et al., 2017, 2018; De et al., 2018) disfavors large stellar radii, so equatorial radii $\gtrsim 14$ km already had low posterior weight. Our new constraints on the radius of PSR J0740+6620 can be viewed as being in better agreement with current tidal deformability measurements.

In this section we update our EOS constraints based on our new measurement of the radius of PSR J0740+6620. Our method is described in detail elsewhere (Miller et al., 2019, 2020, 2021). In brief:

1. We assume that we know the EOS below a threshold density, which we take to be half the saturation density of nuclear matter, because this is roughly the density of the transition between the solid crust and the fluid core (Hebeler et al., 2013). We use the QHC19 EOS of Baym et al. (2019) below this density, but the mass, radius, tidal deformability, and other stellar properties we derive are insensitive to the low-density EOS.

2. We then extend the EOS to higher densities. There are many frameworks for such extensions, but for the most direct comparison with previous work we use the Gaussian process approach that was introduced in this context by Landry & Essick (2019). We use the same set of 100,000 realizations of a high-density EOS that were used in Miller et al. (2021), so any differences between the results we obtain here and our previous results are due to the new measurement rather than to different samples in the EOS space. Our particular Gaussian process sampling of possible equations of state does not explicitly take into account the possibility of phase transitions, but initial studies suggest that phase transitions (or more generally, complex variations of the sound speed with density) are neither favored nor disfavored by current data (Essick et al., 2023; Mroczek et al., 2023a, b). As discussed in Miller et al. (2021), this prior on the sound speed tends towards the speed of light at very high densities.

3. The EOS constraints presented here also impose a Gaussian prior on the nuclear symmetry energy at nuclear saturation density $S=32\pm 2$ MeV, (Tsang et al., 2012), take into account the existence of the three high-mass pulsars (Antoniadis et al., 2013; Arzoumanian et al., 2018; Fonseca et al., 2021), and make use of the tidal deformability posteriors from the gravitational wave observations of GW170817 (Abbott et al., 2017, 2018; De et al., 2018) and GW190425 (Abbott et al., 2020b).

Figure 7: The constraints we infer on the square of the sound speed in dense matter using the results obtained in this paper. At a given baryonic number density (shown here in units of the nuclear saturation density $n_{s}=0.16~{}\rm fm^{-3}$), the lower line shows the 5th percentile of the squared sound speed, the upper line shows the 95th percentile, and the middle line shows the median value of the squared sound speed. The dotted green lines show the constraints inferred by Miller et al. (2021) whereas the solid blue lines show the constraints obtained using the results presented in this paper. The dotted black lines show the prior 5th and 95th percentiles of our EOS samples. We find a slightly lower sound speed at densities below $n\lesssim 3n_{s}\,(\log_{10}{n/n_{s}}\lesssim 0.5)$ and slightly tighter overall constraints on the sound speed in the density range from $n\sim 3-5n_{s}~{}(\log_{10}{n/n_{s}}\sim 0.5-0.7)$ than those found by Miller et al. (2021).

Figure 8: The posterior distributions (shaded regions) and the 25th, 50th, and 75th percentiles of these distributions (horizontal lines) for the maximum stable gravitational mass of a nonrotating neutron star (top panel) and the circumferential radius of a nonrotating $1.4\,M_{\odot}$ neutron star (bottom panel) inferred from the current constraints on the EOS of neutron star matter. The results reported by Miller et al. (2021) are shown in green, the constraints based on the updated measurements presented in this work are shown in blue, and the projected constraints that might follow from analysis of an additional $\sim 2.7\,{\rm Ms}$ of NICER data on PSR J0740+6620 are shown in red.

Table 4: Updated Maximum Mass and Fiducial Radius

Quantity	Data set	$-1\sigma$	median	$+1\sigma$
$M_{\rm TOV}(M_{\odot})$	Miller et al. (2021)	2.08	2.23	2.47
$M_{\rm TOV}(M_{\odot})$	This work	2.08	2.22	2.44
$M_{\rm TOV}(M_{\odot})$	Projected	2.08	2.19	2.37
$R_{e}(1.4~{}M_{\odot})({\rm km})$	Miller et al. (2021)	12.17	12.63	13.11
$R_{e}(1.4~{}M_{\odot})({\rm km})$	This work	12.09	12.57	13.06
$R_{e}(1.4~{}M_{\odot})({\rm km})$	Projected	12.07	12.46	12.85

Note. — Comparison of the $-1\sigma$, median, and $+1\sigma$ points in the posterior distributions of the maximum gravitational mass of a nonrotating neutron star ($M_{\rm TOV}$) and the equatorial circumferential radius of a fiducial nonrotating $1.4~{}M_{\odot}$ neutron star ($R_{e}$) inferred by Miller et al. (2021) using the PSR J0740+6620 data then available and the values we infer here after including an additional $1.1$ Ms of NICER data. The improvements achieved by using the larger present data set are relatively small because other information, notably that provided by the gravitational waves produced by the double neutron star coalescence event GW170817, plays an important role in both analyses and by itself disfavors the high radii that are more disfavored by the present NICER data set than by the earlier NICER data set. The table also shows a projection of the constraints on $M_{\rm TOV}$ and $R_{e}(1.4\,M_{\odot})$ that could potentially be achieved using twice the presently analyzed NICER data on PSR J0740+6620.

Figure 7 shows how the constraints on the square of the sound speed in dense matter are revised when the additional NICER results obtained in this paper are included in the analysis. The details of this analysis and the priors we used can be found in Section 5.3 and Figure 13 of Miller et al. (2021). Table 4 and Figure 8 show how the inferred maximum gravitational mass of a nonrotating neutron star ($M_{\rm TOV}$) and the circumferential radius of a nonrotating neutron star of fiducial gravitational mass $1.4~{}M_{\odot}$ are revised when our current results are included. All the inferences in Table 4 include the current measurements of the masses of the three high-mass pulsars, the gravitational wave-inferred limits on the tidal deformability of neutron stars, and the current experimental constraints on the nuclear symmetry energy (see section 5.3 of Miller et al. 2021 for details).

We find evidence for a slightly softer EOS and consequent reduction of the maximum mass and fiducial radius, but the changes are relatively small. The high values for the radius of PSR J0740+6620 that are more disfavored by the present NICER data than by the previous data are also disfavored by the gravitational wave observations. However, as the slight improvements in Figure 7 suggest, the NICER data has begun to exclude high stellar radii independent of the constraints provided by gravitational wave data. We have included, in Table 4 and Figure 8, a forecast of how the constraints on the equatorial radius of a 1.4 solar mass neutron star and the maximum mass of a nonrotating neutron star might tighten using twice the present data on PSR J0740+6620, or equivalent measurements of another pulsar. However, we caution that the precise values of the constraints on the equatorial radius of a 1.4 solar mass neutron star and the maximum mass of a nonrotating neutron star depend more strongly on the details of a given EOS parameterization than do the EOS constraints themselves, as shown in Table 4 of Miller et al. (2021). Thus, additional NICER observations will continue to improve our knowledge of the properties of dense matter; for example, our projections suggest that doubling the current exposure time on PSR J0740 would by itself reduce the uncertainty of the radius of a fiducial $1.4\,M_{\odot}$ pulsar to $\sim 80\%$ of its present value, at least based on the set of EOS that we currently consider.

5.2 Future Prospects

The X-ray faintness of PSR J0740+6620 makes measurement of its radius using NICER challenging. Based on the data analyzed here, which yield an uncertainty in the radius of $\pm\sim 12\%$, and assuming the uncertainty is proportional to the inverse of the square root of the observing duration, more than 12 Ms of NICER data would be required to achieve an uncertainty of $\pm\sim 5\%$. Although it is unlikely that NICER will be able to collect such a large quantity of data on PSR J0740+6620, even including the data from $\sim 2-3$ Ms of additional observations could appreciably improve the constraints on the equation of state of dense matter, as Table 4 shows.

Future X-ray timing missions such as STROBE-X (Ray et al., 2019), which would have an effective area $\sim 2~{}\rm m^{2}$ at 1.5 keV, or eXTP (Zhang et al., 2019) should be able to provide estimates of the radii of even faint pulsars like PSR J0740+6620 with a precision of a few percent or better. Given the paucity of counts observed by XMM-Newton and the apparent discrepancy between the XMM-Newton pn and MOS1/MOS2 data (see Section 4.3), additional imaging observations of PSR J0740+6620 may further improve our knowledge of the flux from the pulsar and the background flux, and therefore decrease the limit on the radius of PSR J0740+6620.

6 Conclusions

The present NICER observations, which include approximately 1.1 Ms of data in addition to the $\sim\,1.6$ Ms of data that was previously available, have improved the constraints on the radius of PSR J0740+6620 from $13.71^{+2.61}_{-1.50}$ km (Miller et al., 2021) to $12.92_{-1.13}^{+2.09}$ km. These additional data have also reduced the effects of using different analysis methods. Being able to obtain results that are insensitive to the analysis methods used is particularly important when analyzing data on faint sources such as PSR J0740+6620 (see, e.g., Section 4.4 and Salmi et al., 2024). Although the improvements in estimates of the mass and radius of PSR J0740+6620 presented here are incremental, these new estimates provide further information about the properties of matter at supranuclear densities. They imply that the EOS of cold, dense matter is slightly softer than was implied by the previously available data. Future NICER observations will provide even more information about the EOS.

Software

emcee (Foreman-Mackey et al., 2013), MultiNest (Feroz et al., 2009), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Virtanen et al., 2020), PyMultiNest (Buchner, 2016a), PGF/Ti$k$Z (Tantau, 2013), PINT (Luo et al., 2019, 2021), XMM-Newton SAS (Gabriel et al., 2004; SAS development team, 2014), HEASoft (Nasa High Energy Astrophysics Science Archive Research Center (2014), Heasarc)

Appendix A Assessing the Thoroughness of our Nested Sampling Analyses

Table 5: The equatorial radii inferred by our various analyses.

Dataset	Atmosphere Model	Nested Sampling	Nested Sampling $R_{e}<16\,\rm km$	Nested Sampling $R_{e}<16\,\rm km$ $\Delta\theta_{i}<0.4$	MCMC	MCMC $R_{e}<16\,\rm km$	MCMC $R_{e}<16\,\rm km$ $\Delta\theta_{i}<0.4$
NICER+XMM	$\rm{H_{full}}$	$12.704_{-0.970}^{+1.484}$	$12.661_{-0.942}^{+1.368}$	$12.661_{-0.942}^{+1.368}$	$12.922_{-1.129}^{+2.089}$	$12.756_{-1.020}^{+1.495}$	$12.756_{-1.021}^{+1.495}$
NICER+XMM	$\rm{H_{partial}}$	$12.754_{-1.039}^{+1.656}$	$12.692_{-1.004}^{+1.454}$	$12.692_{-1.004}^{+1.454}$	$12.992_{-1.228}^{+2.364}$	$12.760_{-1.074}^{+1.534}$	$12.760_{-1.074}^{+1.534}$
NICER+XMM	$\rm He_{full}$	$12.773_{-1.048}^{+1.657}$	$12.705_{-1.008}^{+1.442}$	$12.705_{-1.007}^{+1.442}$	$13.067_{-1.248}^{+2.455}$	$12.812_{-1.076}^{+1.560}$	$12.814_{-1.073}^{+1.562}$
NICER	$\rm{H_{full}}$	$12.150_{-1.507}^{+2.305}$	$11.984_{-1.394}^{+1.881}$	N/A	$12.475_{-1.938}^{+4.062}$	$11.856_{-1.481}^{+2.302}$	N/A
NICER	$\rm{H_{partial}}$	$11.479_{-0.922}^{+1.232}$	$11.465_{-0.914}^{+1.216}$	$11.365_{-0.811}^{+0.970}$	$11.946_{-1.361}^{+2.512}$	$11.781_{-1.247}^{+1.946}$	$11.804_{-1.191}^{+1.826}$
NICER	$\rm He_{full}$	$12.552_{-1.799}^{+2.816}$	$12.211_{-1.567}^{+2.039}$	N/A	$12.723_{-2.159}^{+4.666}$	$11.877_{-1.514}^{+2.351}$	N/A

Note. — The median and $68\%$ credible regions for the equatorial radius (in km) derived from both our preliminary nested sampling analyses (using MultiNest), and from our MCMC analysis after convergence (using emcee). For the sake of comparison with other works, such as Riley et al. (2021) and Salmi et al. (2024), we also report the credible regions resulting from each analysis after imposing an upper limit of $16\,\rm km$ on the radius. Because some of our analyses identified lunate spot configurations which were neither significant after incorporating XMM-Newton data nor allowed by the priors used in Salmi et al. (2024), we also include the radius constraints after eliminating these solutions by requiring the angular diameters of both spots to be less than 0.4 radians. As this table shows, our preliminary, MultiNest-derived results are typically biased towards smaller radii, and often have erroneously small uncertainties.

In Section 4.1 we found that the best-fitting parameter values from our analysis of NICER data alone, using fully ionized hydrogen and helium atmosphere models, which produce a lunate spot configuration, were drastically disfavored by differences in log likelihood of $>400$ when the XMM-Newton data were included. On the other hand, the best-fitting parameters from our NICER-only analysis using neutron star atmosphere models consisting of hydrogen that could be partially ionized, which describe a mode with separate and much smaller emitting regions, gave a log likelihood $\sim 35$ smaller than the lunate spot configuration. While this first analysis strongly suggested that the lunate configurations that are capable of fitting the NICER data alone are disfavored when the XMM-Newton data are included, the possibility remained that comparably good fits exist for that or other spot configurations that differ from the maximum likelihood configurations identified during the analyses of NICER data alone. To check whether these modes were explored by our analyses, we re-analyzed the nested samples produced by MultiNest using a tempered likelihood function $L^{\prime}=L^{\beta}$ for $\beta\in[1,0]$, such that $\beta=1$ generates samples from the posterior $\beta=0$ generates samples from the prior, using the procedure suggested in Wolfe et al. (2023) to initialize parallel-tempered MCMC runs using nested sampling analyses.

To make things explicit, nested sampling determines a set of samples $i$ and the prior volumes $w_{i}$ associated with each likelihood isocontour. Given a likelihood function $L$, these prior volumes can be used to calculate the evidence as $Z\approx\sum w_{i}L_{i}$, after which the points can be used as weighted posterior samples with weights $p_{i}=w_{i}L_{i}/Z.$ Because the prior volume of each sample is independent of the likelihood function itself as long as the order between points is preserved, we are free to re-weight samples according to $p_{i,\beta}=L_{i}^{\beta}w_{i}/Z_{\beta}$. For example, using $\beta\sim 1/100$, the differences observed earlier between points in different modes can be smoothed over. Thus, plotting effective pseudo-posteriors with $0\leq\beta\leq 1$ can inform us whether certain regions of parameter space were favored over the prior, and explored during sampling, before the highest-likelihood regions dominated the sampling.

Figure 9: A demonstration, from our analyses using a fully ionized helium atmosphere model, that MultiNest explored subdominant modes in each analysis thoroughly enough to assign them appreciable mass at intermediate stages of the analysis. We accomplish this by tempering the likelihood function to place less weight on samples acquired at later times in the analysis. The top panels present effective posteriors from the analysis of NICER data alone, and the bottom panels show results from the joint analysis of NICER and XMM-Newton data. The lightest yellow lines plot the inferred posterior, while the others plot effective posteriors, tempered to varying degrees. The left panels present probability distributions for the angular radius of the first spot, and the right panels present distributions for the separation between spots; in both analyses, the subdominant mode was thoroughly explored before being abandoned for higher-likelihood ground.

Figure 9 illustrates the aforementioned analysis for two parameters, namely, the angular radius of the first spot and the phase offset between the two spots, for fits using fully ionized hydrogen atmospheres to the NICER data alone and to both the NICER and XMM-Newton data. Considering first the fits to only the NICER data, the results show that the mode with a large first spot, and spots separated by $\sim 0$ rotational cycles, clearly dominates the posterior. However, after tempering the likelihood by a factor of $\beta\lesssim 0.06$, we observe substantial pseudo-posterior mass at small spot angular radii and at intermediate phased offsets, suggesting that these families of modes were explored early in the analysis before the far-higher likelihood points belonging the lunate mode were identified. Similarly, while the posterior from the joint analysis of NICER and XMM data is comprised almost entirely by configurations with small, nearly antipodal spots, before honing in on these much higher-likelihood solutions the analysis thoroughly explored the crescent-shaped configurations which better suit the NICER data without the flux constraints from the data provided by XMM-Newton.

Although this analysis indicates that MultiNest was able to identity the various spot configurations capable of describing the present NICER and XMM-Newton data analyzed here, with the settings used in this work, it was unable to produce accurate posterior inferences. Table 5 lists the $68\%$ uncertainty and median values of the equatorial radius for each of our preliminary nested sampling analyses and our converged MCMC analyses. In every case, MultiNest estimated a posterior that is both erroneously narrow and biased towards low radii. We also list in Table 5 the radius credible regions obtained after placing various cuts on our posterior samples to simulate the effects of imposing the priors used in Salmi et al. (2024). The posterior distributions resulting from these constraints agree better with the MultiNest-derived results reported in Salmi et al. (2024).

Appendix B Posterior Distributions from Analysis of Only NICER Data

Table 6 lists the median, $\pm 1\sigma$, and $\pm 2\sigma$ points in the posterior distributions obtained by fitting our models to only the NICER data, assuming a fully ionized hydrogen atmosphere, and the resulting maximum likelihood values for each of the parameters in these models. We display the complete corner plot of the posteriors from these same analyses in Figure 10.

Table 6: Fits to Only Nicer Data

Parameter	median	$-1\sigma$	$+1\sigma$	$-2\sigma$	$+2\sigma$	Maximum Likelihood
$R_{e}\,\rm(km)$	12.475	10.537	16.538	9.602	21.840	11.731
$GM/c^{2}R_{e}$	0.245	0.185	0.289	0.140	0.308	0.259
$M\,(M_{\odot})$	2.066	1.972	2.161	1.878	2.252	2.058
$\theta_{\rm c1}\rm(rad)$	1.331	0.384	2.632	0.149	2.972	1.680
$\Delta\theta_{1}\,\rm(rad)$	2.096	1.266	2.485	0.552	2.671	1.159
$kT_{\rm eff,1}\,\rm(keV)$	0.020	0.014	0.028	0.011	0.036	0.021
$\theta_{\rm c2}\rm(rad)$	1.718	0.629	2.677	0.193	2.980	2.907
$\Delta\theta_{2}\,\rm(rad)$	1.167	0.719	1.993	0.490	2.570	0.617
$kT_{\rm eff,2}\,\rm(keV)$	0.098	0.088	0.109	0.078	0.122	0.107
$\Delta\phi_{2}\,\rm(cycles)$	0.003	0.000	1.000	0.000	1.000	1.000
$\theta_{\rm obs}\,\rm(rad)$	1.528	1.467	1.589	1.444	1.610	1.471
$N_{H}\,(10^{20}\,\rm cm^{-2})$	1.056	0.273	2.524	0.037	4.182	0.063
$d\,\rm(kpc)$	1.228	1.039	1.423	0.861	1.619	1.116

Note. — A comparison of the $-2\sigma$, $-1\sigma$, median, $+1\sigma$, and $+2\sigma$, and maximum likelihood values inferred from our analysis of NICER data for PSR J0740+6620. These results were derived under the assumption of a fully ionized pure hydrogen atmosphere.

Figure 10: Posterior probability density distributions from our analysis of NICER data alone, using a fully ionized hydrogen atmosphere model. In the one-dimensional plots of the distance to and gravitational mass of the pulsar, dashed lines indicate the priors that we applied.

Appendix C Posterior Distributions from Analysis of both NICER and XMM-Newton Data

Table 7 lists the median, $\pm 1\sigma$, and $\pm 2\sigma$ points in the posterior distributions obtained by fitting our models to both the NICER data and XMM-Newton data, assuming a fully ionized hydrogen atmosphere, and the resulting maximum likelihood values for each of the parameters in these models. We display the complete corner plot of the posteriors from these same analyses in Figure 11.

Table 7: Fits to NICER and XMM-Newton Data

Parameter	median	$-1\sigma$	$+1\sigma$	$-2\sigma$	$+2\sigma$	Maximum Likelihood
$R_{e}\,\rm(km)$	12.922	11.793	15.010	10.992	18.574	11.541
$GM/c^{2}R_{e}$	0.236	0.204	0.257	0.165	0.272	0.250
$M\,(M_{\odot})$	2.066	1.976	2.155	1.886	2.246	1.958
$\theta_{\rm c1}\rm(rad)$	1.573	0.946	2.222	0.641	2.546	1.387
$\Delta\theta_{1}\,\rm(rad)$	0.101	0.072	0.140	0.051	0.199	0.092
$kT_{\rm eff,1}\,\rm(keV)$	0.095	0.085	0.104	0.076	0.114	0.106
$\theta_{\rm c2}\rm(rad)$	1.578	0.957	2.220	0.640	2.531	1.980
$\Delta\theta_{2}\,\rm(rad)$	0.100	0.073	0.138	0.052	0.193	0.112
$kT_{\rm eff,2}\,\rm(keV)$	0.095	0.085	0.104	0.077	0.114	0.099
$\Delta\phi_{2}\,\rm(cycles)$	0.553	0.424	0.576	0.413	0.587	0.428
$\theta_{\rm obs}\,\rm(rad)$	1.529	1.469	1.587	1.444	1.610	1.539
$N_{H}\,(10^{20}\,\rm cm^{-2})$	0.953	0.255	2.209	0.034	3.871	0.023
$d\,\rm(kpc)$	1.205	1.020	1.397	0.844	1.588	1.297

Note. — A comparison of the $-2\sigma$, $-1\sigma$, median, $+1\sigma$, and $+2\sigma$, and maximum likelihood values inferred from our joint analysis of NICER and XMM-Newton data for PSR J0740+6620. These results were derived under the assumption of a fully ionized pure hydrogen atmosphere.

Figure 11: Posterior probability density distributions from our joint analysis of NICER and XMM-Newton data, using a fully ionized hydrogen atmosphere model. In the one-dimensional plots of the distance to and gravitational mass of the pulsar, dashed lines indicate the priors that we applied.