A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy

The PSR J0740+6620 analysis is based on a sequence of exposures with NICER XTI (hereafter NICER) acquired in the period 2018 September 21–2020 April 17. The event data were obtained using similar filtering criteria to the previously analyzed NICER data set of PSR J0030+0451. We only used good time intervals when all 52 active detectors were collecting data but rejected all events from DetID 34, as it often shows enhanced count rates relative to the other 51 detectors. We excluded time intervals when PSR J0740+6620 was situated at an angle ≤80° from the Sun to reduce the increase in background in the lowest channels due to optical loading. We further excised events collected at low cut-off rigidity (COR_SAX values <5) to minimize particle background contamination. The resulting cleaned event list has an on-source exposure time of 1602683.761 s.

We model events registered in the PI channel subset [30, 150), corresponding to the nominal photon energy range [0.3, 1.5] keV. ³² The quoted nominal photon energy for a channel is the energy mid-point for that channel.

Below nominal energy 0.3 keV, the instrument calibration products have greater a priori uncertainty due to the sharp lower-energy threshold cutoff, as well as the presence of unrejected instrumental noise that can extend above 0.25 keV. We therefore neglect the information below channel 30 in order to reduce the risk of inferential bias. Wolff et al. (2021) detect pulsations with the highest significance when considering event data up to nominal energy ∼1.2 keV; in this Letter we include the additional information at nominal energies in the range [1.2, 1.5] keV. The number of counts generated by the PSR J0740+6620 hot regions in channels 150 and above, however, is small and diminishes relative to the counts generated by background processes (including a non-thermal component from the environment in the near-vicinity of PSR J0740+6620); we therefore neglect this higher-energy information, and focus energy resolution at nominal energies below 1.5 keV.

The NICER event data are phase-folded according to the NANOGrav radio timing solution presented in Fonseca et al. (2021a). The phase-binned count numbers are displayed in Figure 1.

Zoom In Zoom Out Reset image size

The XMM-Newton (hereafter XMM) telescope observed PSR J0740+6620 as part of a Director's Discretionary Time program in three visits: 2019 October 26 (ObsID 0851181601), 2019 October 28 (ObsID 0851181401), and 2019 November 1 (ObsID 0851181501). The EPIC instruments (pn, MOS1, and MOS2) were employed in "Full Frame" imaging mode with the "Thin" optical blocking filters in place. Due to the insufficiently fast read-out times (73.4 ms for EPIC-pn and 2.6 s for EPIC-MOS1/2), the data do not provide useful pulse timing information, so only phase-averaged spectral information is available. The XMM data were reduced using the Science Analysis Software (SAS ³³ ) using the standard set of analysis threads.

The event data were first screened for periods of strong soft proton background flaring. The resulting event lists were then further cleaned by applying the recommended PATTERN (≤12 for MOS1/2 and ≤4 for pn) and FLAG (0) filters. The final source event lists were obtained by extracting events from a circular region of radius 25'' centered on the radio timing position of PSR J0740+6620. This resulted in effective exposure times of 6.81, 17.96, and 18.7 ks for the pn, MOS1, and MOS2 instruments, respectively.

One of the ultimate aims of this Letter is to report a likelihood function of mass and (equatorial) radius that can be used for EOS posterior computation. As highlighted by Riley et al. (2018), it is desirable to define a prior PDF which is jointly flat with respect to two parameters within the prior support; these parameters can simply be mass and radius, or some deterministically related variables, such as mass and compactness. The marginal joint posterior PDF of these parameters is then proportional to the marginal likelihood function of these parameters, meaning that the marginal likelihood function can be estimated from posterior samples, for use in subsequent inferential analyses.

In this Letter we follow Riley et al. (2019) by defining a joint prior PDF of mass and radius that is flat within the prior support, which is maximally inclusive with regard to theoretical EOS predictions. The prior support is zero for R > 16 km because we are not aware of any EOS models predicting a radius higher than this limit that would be compatible with current constraints from nuclear physics, or with the constraints posed by the gravitational wave measurement of tidal deformability for the binary neutron star merger GW170817 (see, e.g., Reed et al. 2021). A difference to Riley et al. (2019) is that we define the prior support using a higher compactness limit as discussed in Section 2.2.3 (see Table 1 for the prior PDF and support). The prior support is also subject to the condition that the effective gravity lies within a bounded range at every point of the rotating oblate surface, in order to conform to bounds on the atmosphere models we condition on (see Section 2.3.2 and Table 1).

In Section 2.2.2 we introduce information about the mass in the form of a marginal PDF whose shape approximates the marginal likelihood function of mass derived in a recent radio timing analysis; we multiply this likelihood function with the jointly flat prior PDF of mass and radius described above, thereby defining an updated prior PDF for the X-ray pulse-profile modeling in this Letter. Our prior PDF for pulse-profile modeling is therefore not jointly flat with respect to mass and radius, but all shape information is likelihood-based, satisfying the requirements of Riley et al. (2018).

Radiation is transported from the oblate X-ray-emitting surface to distant static telescopes by relativistic ray-tracing, as described by Bogdanov et al. (2019b) and references therein. The gravitational mass of PSR J0740+6620, the distance of PSR J0740+6620 from Earth, and the inclination of the Earth's line of sight to the PSR J0740+6620 spin axis are all parameters of the source–receiver system that need to be specified in order to simulate an X-ray signal incident on a telescope, and can be inferred via pulsar radio timing. Note that these static model X-ray telescopes are fictitious constructs: the real telescopes are in motion relative to the pulsar (see, e.g., Riley 2019). The NICER event data pre-processing operations include the phase-folding of on-board arrival times according to a NANOGrav radio timing solution; the phase-folded events are then the events that would be registered by a telescope that is static relative to the pulsar.

The NANOGrav and CHIME/Pulsar collaborations jointly performed wideband radio timing of PSR J0740+6620 (Fonseca et al. 2021a, 2021b). A product of this radio timing work was a joint posterior PDF of the pulsar gravitational mass, Earth distance, and the inclination Earth subtends to the orbital direction ³⁴ that we can condition on as an informative prior PDF for joint NICER and XMM X-ray pulse-profile modeling. The synergy—due to degeneracy breaking—between radio timing and X-ray modeling yields a higher sensitivity to the pulsar radius and the parameters of a surface X-ray emission model. We used two sets of joint posterior PDFs provided by Fonseca et al. (2021a) as prior information over the course of our analysis: one set that was numerically estimated using a weighted least-squares solver, used in this Letter for an exploratory analysis (Section 3.2), and a second that instead used a generalized least-squares (GLS) algorithm for determining timing model parameters from wideband timing data, used in this Letter for a production analysis (Section 3.3). Fonseca et al. (2021a) report on results obtained with GLS fitting of wideband timing data, though they publicly provide both sets of PDFs as they were used in this work.

Systematic error in the parameter estimates by NANOGrav and CHIME/Pulsar is dominated by sensitivity to the choice of the dispersion-measure variability (DMX) model. Posterior PDFs were computed for three DMX variants, and the systematic error implied by the posterior variation was substantially smaller than the formal posterior spread conditional on any one DMX model. We average (marginalize) the posterior PDFs over the three DMX models with a uniform weighting. ³⁵

The prior PDFs we implement for every model that conditions on joint NANOGrav and CHIME/Pulsar wideband radio timing are as follows. The distance is separable from the pulsar mass and the Earth inclination to the orbital axis, but the mass and (cosine of the) inclination are correlated. For the distance, we implement the NANOGrav × CHIME/Pulsar measurement by first marginalizing the PDFs over the DMX models and then reweighting from the flat distance prior conditioned on by NANOGrav and CHIME/Pulsar to a physical distance prior in the direction of PSR J0740+6620 following the exemplar treatment of distance information by Igoshev et al. (2016). ³⁶ The physical distance prior, however, remains relatively flat in the context of the likelihood function—meaning the likelihood function is dominant in the measurement—and the modification to the original DMX-averaged PDF is entirely minor. Additional distance likelihood information derives from the Shklovskii effect which effectively truncates the PDF, putting an upper limit on the distance (Shklovskii 1970); such truncation, however, occurs well into the tail of the NANOGrav × CHIME/Pulsar distance distribution, so it is also an unimportant detail. Finally, we approximate the NANOGrav × CHIME/Pulsar marginal PDF of distance as a skewed and truncated Gaussian distribution (see Table 1 for details). We display the marginal distance PDF variants referenced above in Figure 2.

Zoom In Zoom Out Reset image size

For the mass and the (cosine of the) orbital inclination we implement a multi-variate Gaussian distribution with the covariance matrix of the DMX-averaged PDF. Implicit in this PDF is a prior PDF defined by the radio timing collaboration. The marginal prior PDF of the pulsar mass is not flat: the prior PDF of the cosine of the orbital inclination and of the mass of the white dwarf companion of PSR J0740+6620 is jointly flat, and these variables map deterministically to the pulsar mass through the binary mass function and (precise) radio timing of the system. We do not modify the joint prior PDF of the pulsar mass and the inclination despite the likelihood function of the pulsar mass (marginalized over all other radio timing parameters) being formally desirable for EOS parameter estimation (see Section 2.2.1 and Riley et al. 2018). In this case the posteriors are sufficiently dominated by the likelihood function to ignore the small structural modifications that would result from tweaking relatively diffuse priors. Moreover, changing the pulsar mass prior PDF would change the prior PDF of the companion, which may be undesirable.

Ultimately, handling the detailed structure ³⁷ of these prior PDFs—i.e., beyond the location and marginal spread of the parameters—is not important because the NICER likelihood function is not sufficiently sensitive to some combination of these radio-timing parameters and its native parameters (such as equatorial radius) for posterior inferences to change to any discernable degree. That is, the changes will be difficult to resolve from sampling noise and implementation error for instance (Higson et al. 2018, 2019), and relative to model-to-model systematic variations that, when marginalized over, further broaden the posteriors of shared parameters of interest. Also note that the distance only enters in the likelihood function in combination with the effective-area scaling parameter defined in Section 2.4 that operates on all X-ray instrument response models because of global calibration error; the distance to PSR J0740+6620 could therefore be combined with this absolute scaling parameter, further diminishing the importance of treating fine details in the prior PDF of the distance.

The original version of the X-PSI package used to model PSR J0030+0451 (Riley et al. 2019) via relativistic ray-tracing assumed that no more than one ray connected the telescope to each point on the stellar surface. However, photons emitted from a star with compactness greater than GM/Rc² = 0.284 can have a deflection angle larger than π, which makes multiple images of small regions at the "back" of the star possible. This was not an issue for the analysis of PSR J0030+0451, because the 95% credible region only allowed compactness values in the range of GM/Rc² ≤ 0.171. Additionally, the independent analysis by Miller et al. (2019) allowed for the possibility of multiply imaged regions on the star but still led to a credible range on compactness that is too small to allow for multiple images.

The radio observations of Shapiro delay in the PSR J0740+6620 binary system by NANOGrav and CHIME/Pulsar (Fonseca et al. 2021a) lead to the marginal 68% credible interval of M = 2.08 ± 0.07 M_⊙. For reference, a small selection of EOSs that cover a realistic range of stiffness are shown in Figure 3. The EOSs shown include a set of three EOSs (HLPS soft, intermediate, and stiff) constructed by Hebeler et al. (2013) that span a range of stiffness allowed by nuclear experiments, as well as the A18+δ(v)+UIX* EOS (Akmal et al. 1998) (abbreviated to APR). Each of the four curves shows the values of the equatorial compactness and mass for stars rotating at the rate of 346 Hz, computed using the code rns (Stergioulas & Friedman 1995). Figure 3 shows that for this set of EOSs, the large mass values (>2.0 M_⊙) lead to a significant range of models that have large enough compactness (i.e., are above the horizontal line at 0.284) that multiple images of some part of the star are possible. Because the stars are oblate, the compactness at the spin poles is larger than the equatorial compactness, so the range of multiply imaged stars extends to slightly lower values of compactness; however, the change is not perceptible on this figure.

Zoom In Zoom Out Reset image size

Due to the expectation that PSR J0740+6620 could be very compact, the X-PSI code was extended to allow multiple rays from any point to reach the telescope. This improvement to X-PSI as well as details of code validation are discussed in Bogdanov et al. (2021). X-PSI sums over the primary and the visible higher-order images of any regions on a star that lie in a multiply imaged region. The X-PSI package can typically detect up to the quaternary, quinary, or senary order depending on resolution settings. In practice only secondary images, and potentially the tertiary images for some configurations, might be important; but omission of secondary images can lead to large errors of ${ \mathcal O }(10 \% )$ in the light-curve calculation (see, e.g., the X-PSI documentation ³⁸ ).

The X-ray signal is attenuated to some degree by the intervening interstellar medium along the line of sight to the pulsar. In all models, the attenuation physics is parameterized solely by the neutral hydrogen column density N_H, relative to which the abundances of all other attenuating gaseous elements, dust, and grains are fixed by the state-of-the-art TBabs model (Wilms et al. 2000, updated in 2016). We implement this attenuation model as a one-dimensional lookup table with respect to energy at a fiducial column density, and then raise the attenuation factor to the power of the ratio of the column density to the fiducial value.

The column density along the line of sight to PSR J0740+6620 can be estimated using several techniques. The HEASARC neutral hydrogen map tool ³⁹ yields N_H ≈ 3.5 × 10²⁰ cm⁻² given the most modern map (HI4PI Collaboration et al. 2016). ⁴⁰

Using the relation between dispersion measure (DM) and N_H from He et al. (2013) together with the DM measurement reported by Cromartie et al. (2020), N_H ≈ 4.5 × 10²⁰ cm⁻². Finally, estimates based on 3D $E\left(B-V\right)$ extinction maps (Lallement et al. 2018), together with the relation between extinction and N_H from Foight et al. (2016), yield values of N_H ≈ 4.5 × 10²⁰ cm⁻².

Given that N_H could be as large as 4.5 × 10²⁰ cm⁻² based on these estimates, and given the large uncertainties in the relations between N_H and other quantities such as the DM, a conservative prior of N_H ∼ U(0, 10²¹) cm⁻² is warranted.

The surface effective temperature field of a rotation-powered MSP is one of the primary sources (if not the primary source) of uncertainty a priori regarding the physical processes that generate the X-ray event data. The image of a neutron star cannot be resolved with any current X-ray telescope, so all our knowledge about the surface temperature field comes from models. These models heavily rely on the assumptions about magnetic field configuration, which is a crucial part in calculating heating by space-like or return current in some sub-region of the open field line footprints (Kalapotharakos et al. 2021) or anisotropic thermal conductivity of the outer star layers (e.g., De Grandis et al. 2020; Kondratyev et al. 2020). There is, therefore, an essentially unknown degree of complexity in the structure of the temperature field.

A given likelihood function is insensitive to complexity beyond some degree. We focus on a simple model of the surface temperature field: two disjoint hot regions that are simply connected spherical caps (when projected onto the unit sphere) wherein the effective temperature of the atmosphere is uniform. This model may be identified as ST-U in Riley et al. (2019). The phase-folded NICER pulse profile is suggestive of two phase-separated hot regions which may or may not be disjoint.

We could begin with an even simpler model that restricts the hot regions to be related via antipodal reflection symmetry (ST-S in Riley et al. 2019). Computing a posterior conditional on such a simple surface hot region model is justifiable, and it also delivers a lower-dimensional target distribution to check for egregious model implementation error (as reasoned by Riley et al. 2019). However, for the analysis reported in this Letter, we immediately explored the ST-U model because antipodal reflection symmetry is physically unrealistic and the ST-U model is sufficiently simple and inexpensive to condition on; nevertheless, if one is interested in the question of whether we are sensitive to deviations from this symmetry, we open-source the entire analysis package for this Letter, and the online X-PSI documentation offers guidance on how to condition on the ST-S model.

Nodes can be added to the model space by incrementing the complexity of the surface hot regions, following Riley et al. (2019), wherein symmetries between hot regions are broken ⁴¹ and temperature components are added. ⁴² In this Letter we introduce a new binary flag for the atmosphere composition (hydrogen or helium as discussed in Section 2.3.2) within the (single-temperature) hot regions. However, we consider only one higher-complexity model than ST-U. We opt to do this because subject to limited (computational) resources, ST-U performs well and is ultimately determined to be sufficiently complex.

We now define prior PDFs that are shared by all models. Every hot region has one effective temperature component. The prior PDF of that effective temperature is flat in its logarithm, diffuse, ⁴³ and is separable from the joint prior PDF of all other model parameters. Every hot region also has a colatitude coordinate and an azimuth coordinate (i.e., phase shift) at the center of a constituent spherical cap with a finite effective temperature component (Riley et al. 2019). The prior PDF of these coordinates is non-trivial because it is not separable from the prior PDF of other parameters when there are two surface hot regions in the model, as we proceed to explain.

In order to eliminate a trivial form of hot region exchange-degeneracy, ⁴⁴ we define the prior support by imposing that one hot region object in the X-PSI model has a center colatitude that is always less than or equal to that of the companion hot region. We also impose that the hot regions are disjoint, meaning that they cannot overlap in the prior support, because otherwise the model cannot always be characterized by the number of disjoint hot regions; the non-overlapping condition is a function of the center colatitudes, the center azimuths, and the hot region angular radii, meaning that the joint prior is not separable, and the marginal prior PDFs are modulated relative to the PDFs that were used to construct the form of the joint PDF in six dimensions. For instance, the joint prior PDF of the azimuthal coordinates exhibits rarefaction where the azimuths are close in value.

We now construct the prior PDF specifically for the ST-U model (Section 3.3) that is the focus of this Letter. First, as a construction tool, define a flat PDF of the cosine of each hot region center colatitude, with support being the interval $\cos ({\rm{\Theta }})\in [-1,1];$ such a PDF is founded on the argument of isotropy, a priori, of the direction to Earth (see below). Second, define the PDF of each hot region center azimuth ⁴⁵ to be flat on the interval 2πϕ ∈ [0,2π] radians. Third, define a flat PDF of each angular radius on the interval [0, π/2]. To calculate the marginal prior PDF of one of these parameters, marginalize over the other five parameters subject to the prior support condition of non-overlapping hot regions. For instance, the prior PDF of a hot region angular radius is given by marginalizing over the angular radius of the companion hot region, and the hot region center colatitudes and azimuths. Marginalization yields a marginal prior PDF that deviates in form from the initially flat PDF of the parameter; the support of the PDF remains unchanged, however.

Our hot region models are phenomenological and, to a degree, agnostic, despite being influenced by temperature fields implied by pulsar magnetospheric simulations. Our choice here is contrary to Riley et al. (2019), wherein the prior PDF of the colatitude at the center of a hot region in isolation—meaning one hot region on the surface—was flat. A flat PDF of colatitude, combined with a flat prior in azimuth, yields an anisotropic probability density field on the unit sphere, weighted toward polar regions. One reason this might be a sub-optimal assumption is because of X-ray selection effects: pulsed emission from two hot regions will exhibit a larger amplitude if the hot regions are closer to the equatorial zone than a polar zone, for equatorial observers as is the case assumed for PSR J0740+6620 based on the NANOGrav × CHIME/Pulsar wideband radio timing measurements (Section 2.2.2). However, such arguments are tenuous, taking no heed of radio selection effects and magnetospheric physics, for instance.

The specific intensity is determined from lookup tables generated using the nsx atmosphere code assuming either fully ionized hydrogen or helium (Ho & Lai 2001) or partially ionized hydrogen (Ho & Heinke 2009). ⁴⁶ In this work, we consider only fully ionized models since the limited parameter ranges of existing opacity tables (Iglesias & Rogers 1996; Badnell et al. 2005; Colgan et al. 2016) reduce the accuracy of atmosphere models employing these tables (see Bogdanov et al. 2021 and Miller et al. 2021 for further discussion and comparisons).

The surface exterior to the hot regions does not explicitly radiate in any model we condition on. The signal that would be generated by a cooler exterior surface can be robustly subsumed in the phase-invariant count rate terms described in Section 2.5.1, provided that the angular scale of the hot regions (whose images are explicitly integrated over) is small and that exterior emission is soft and thus dominated by other NICER backgrounds in the channels with low nominal photon energies. Explicitly including radiation from the exterior surface would add a very weak informative mode of dependence on the mass, radius, and other parameters of interest, and thus the likelihood function is assumed insensitive to its exact treatment. Explicit treatment would also require handling of the ionization state of the cooler atmosphere.

Let d_N be the NICER count matrix over phase interval–channel pairs. We now define variables that the NICER likelihood is a function of: let s be a vector of parameters, collected from Section 2.2, of which the incident signal generated by the pulsar is a function; let nicer denote a (parameterized) model for the response of the instrument in response to incident radiation; and let $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$ be a set of statistically independent nuisance variables, one per detector channel that contains event data to be modeled, defined as the phase-invariant expected count rate. The NICER event data are phase-resolved, so there are multiple random variates—distributed in phase—per background random variable. The background model with variables $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$ is free-form as discussed by Riley et al. (2019). The set of these variables is designed to handle complex channel-to-channel variations in the expected phase-invariant count rate. A physical background model should need (far) fewer random variables for the underlying background-generating process to capture these complexities. The likelihood function is denoted by $p({d}_{{\rm{N}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$, nicer).

We numerically marginalize over the variables $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$ to yield a marginal likelihood function that is combined with a joint prior PDF to define a target distribution to draw samples from. The count rate variables have separable flat prior PDFs that are strictly improper because we do not explicitly define upper bounds on the prior support of each variable (Riley et al. 2019); the posterior is considered integrable, however, so these ill-defined prior PDFs do not result in posterior pathologies. The separable prior PDF of the count rate variables is overly diffuse, with extremely high prior-predictive complexity, such that inferences will be insensitive to minor changes to the function: ${\mathbb{E}}[{b}_{{\rm{N}}}]\sim U(0,{ \mathcal U })$, where the upper bound ${ \mathcal U }$ of the prior support is left unspecified. ⁵⁰

The marginalization operation is separable over channels, yielding a product of one-dimensional integrals. We need to perform fast numerical marginalization over the $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$ in order to compress the dimensionality of the sampling space. The simpler the form of the integrands—the functions of the $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$—the more straightforward fast numerical marginalization is. If the prior PDF $p({\mathbb{E}}[{b}_{{\rm{N}}}])$ is flat, the integrand has a single global maximum, and would be highly Gaussian if the peak of the conditional likelihood function p(${d}_{{\rm{N}}}| s,{\mathbb{E}}[{b}_{{\rm{N}}}]$, nicer) lies within the support of $p({\mathbb{E}}[{b}_{{\rm{N}}}])$.

As discussed by Riley et al. (2019), a major open question regards the total expected spectral signal that is attributable to the surface hot regions in reality. ⁵¹ In lieu of a physical background model—which as discussed above is difficult to formulate for NICER—independent information conditional on data acquired with an imaging X-ray telescope such as XMM can be fundamentally valuable for deriving robust inferences if the exposure time is sufficiently long; it will, however, be substantially shorter than the order megasecond exposures of the near-dedicated NICER telescope.

Note that in order to infer the contribution from contaminating sources to the NICER event data, ⁵² we would need to separate out the $\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}$ into an environmental background model—with some informative prior PDF including space weather contributions—and some model for the contribution from contaminating (point) sources in the field that cannot be resolved from the point-spread function (PSF) of PSR J0740+6620. However, for the principal purpose of constraining the physical properties of the pulsar, we are uninterested in the distinction between NICER backgrounds. Moreover, we do not have the statistical power to distinguish the contributions if the priors for the components are all rather diffuse. In other words, if some of the priors are informative, we do not have the statistical power to gain much information a posteriori.

Information about the total background contribution to the NICER event data—including contaminating (point) sources in the field (see Bogdanov et al. 2019a for a breakdown of components)—is encoded in the XMM spectroscopic imaging event data. Due to the relatively brief exposure times of the observations, very few background counts are available for a reliable background estimate. Thus, we obtained representative background estimates with higher photon statistics from the blank-sky event files provided by the XMM-Newton Science Operations Centre. ⁵³ The blank-sky images were filtered in the same manner as the PSR J0740+6620 field images and the background was extracted from the same location on the detector as the pulsar. The resulting background spectrum was then rescaled so that the exposure times and BACKSCAL factors matched those of the PSR J0740+6620 exposures.

Leveraging spatial resolving power, we can derive posterior inferences about the signal from the pulsar in isolation with little confusion about the expected contributions from the pulsar and background sources. ⁵⁴ When this information about the pulsar signal is injected into modeling of NICER observations—either explicitly as a prior or as a likelihood factor—the contribution from the pulsar to the NICER event data is informed; here we work toward a likelihood function factor for the XMM telescope.

Let d_X be the XMM time-integrated count vector (over detector channels) from a region ${ \mathcal S }$ of a CCD of a camera, that is some sufficiently large subset of the support of the PSF of the PSR J0740+6620 while optimizing signal-to-noise. ⁵⁵ We now define variables that the XMM likelihood is a function of. The vector of parameters of the incident signal generated by the pulsar remains as s, shared with the NICER likelihood function. Let xmm denote a (parameterized) model for the response of an XMM camera in response to incident radiation. Once more, in lieu of a physical model for the underlying background-generating process, let us define one statistically independent random nuisance variable per detector channel: the expected count rate. ⁵⁶ These variables are collectively denoted $\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}$. The likelihood function for each XMM camera is denoted by p(${d}_{{\rm{X}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}$, xmm). We numerically marginalize over these model variables in the same vein as for the NICER background-marginalized likelihood function; however, to constrain the signal from PSR J0740+6620 we are in need of a more informative prior PDF of $\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}$.

To form the prior PDF of $\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}$ for an XMM camera, we consider a set of independent Poisson-random variates $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$ defined as time-integrated astrophysical (sky) background count numbers in detector channels. The events constituting $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$ are extracted from a blank-sky ⁵⁷ region ${ \mathcal B }$ of the camera CCD array that is disjoint from region ${ \mathcal S }$ associated with PSR J0740+6620. The XMM cameras are composed of multiple side-by-side CCD detectors. For each camera, it is preferable to choose the region ${ \mathcal B }$ on the same detector as ${ \mathcal S }$ because the different CCD detectors have slightly different responses to incident radiation and therefore exhibit slightly different astrophysical (sky) backgrounds. For blank-sky exposures, the conditional sampling distribution in the space of the data is $p(\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}| \left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\})$, where the variables $\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}$ are per-channel expected count numbers.

To constrain the background in the XMM images of PSR J0740+6620, we use blank-sky estimates generated using the XMM-Newton SAS tools. The blank-sky exposures are much longer in duration than the exposures to PSR J0740+6620 by almost two orders of magnitude, allowing us to constrain the astrophysical (sky) background count rate in the on-source exposures more tightly (in the absence of systematic error) than possible from blank-sky extraction regions from the same CCD during the shorter on-source exposure. For each XMM camera the respective region ${ \mathcal B }$ is not only localized to the same CCD as the PSR J0740+6620 PSF, but is the same region of the CCD. However, the blank-sky estimates are based on an ensemble of pointings over the sky. Our cross-check of the sky background in these reference exposures against the sky background in the vicinity of PSR J0740+6620 did not yield evidence of systematic difference.

We now formally derive the constraints on the expected number of background counts per detector channel registered within the source region ${ \mathcal S }$ of an XMM camera, $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$, conditional on the counts registered in region ${ \mathcal B }$ over the ensemble of blank-sky exposures. A prior PDF of the variables $\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}$ is needed. In the same vein as for the NICER background variables, we define a separable prior PDF

$$\begin{eqnarray}&&{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\sim U(0,{\mathscr{V}}),\end{eqnarray} \tag{ 1 }$$

where the upper-bound ${\mathscr{V}}$ of the prior support is left unspecified. This trivial model over-fits the background data, with posterior PDF

$$\begin{eqnarray}&&p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})\propto p(\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}| \left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\});\end{eqnarray} \tag{ 2 }$$

the vector of random variates $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$ is coincident in value with the maximum a posteriori vector in the space of $\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}$.

As is the case for the NICER background marginalization operation described in Section 2.5.1, if the PDF $p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ is flat, the marginalization is more straightforward to execute. For example, we could form a flat prior PDF spanning the highest-density x% posterior credible interval in ${\mathbb{E}}[{b}_{{\rm{N}}}]$, given $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$. We form a flat prior PDF in a simpler way. The PDF $p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ has the approximate structure of a truncated Gaussian with deviation dependent on the number of counts ${{\mathscr{B}}}_{{\rm{X}}}$. We let the lower and upper bounds of the support respectively be ${\mathscr{L}}:= \max \left(0,{{\mathscr{B}}}_{{\rm{X}}}-n\sqrt{{{\mathscr{B}}}_{{\rm{X}}}}\right)$ and ${\mathscr{U}}:= {{\mathscr{B}}}_{{\rm{X}}}+n\sqrt{{{\mathscr{B}}}_{{\rm{X}}}}$, where n is a setting that controls the degree of conservatism. For some channels the number of counts is low and the PDF $p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ deviates substantially from being a truncated Gaussian, but we nevertheless adopt the same procedure—the lower-bound is pushed to zero or far into the lower tail of the distribution, but the upper bound remains sensible. For some channels, however, ${{\mathscr{B}}}_{{\rm{X}}}=0$. In these cases, we simply set ${\mathscr{L}}:= 0$ and then iterate upwards in channel number to locate the next finite value of ${{\mathscr{B}}}_{{\rm{X}}}+n\sqrt{{{\mathscr{B}}}_{{\rm{X}}}}$. ⁵⁸ We choose n = 4. With the flat PDF of ${\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]$ defined, we transform variables to derive the prior PDF $p({\mathbb{E}}[{b}_{{\rm{X}}}]| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ according to

$$\begin{eqnarray}&&{\mathbb{E}}[{b}_{{\rm{X}}}]=\displaystyle \frac{{A}_{{ \mathcal S }}}{{A}_{{ \mathcal B }}}\displaystyle \frac{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]}{{T}_{{\mathscr{B}}}},\end{eqnarray} \tag{ 3 }$$

where ${T}_{{\mathscr{B}}}$ is the blank-sky exposure time needed to transform to a count rate, and ${A}_{{ \mathcal S }}/{A}_{{ \mathcal B }}$ is the ratio of the areas of the extraction regions ${ \mathcal S }$ (encompassing the PSR J0740+6620 PSF) and region ${ \mathcal B }$.

In many respects this flat PDF of ${\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]$ is a remarkably conservative choice for the prior constraint on $p({\mathbb{E}}[{b}_{{\rm{X}}}]| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$. A reason it is not conservative is the assumption of zero prior mass above an upper count-rate limit ${\mathscr{U}}$ somewhere in the upper tail of the true probability PDF $p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$, especially because the count-number data d_X for each XMM camera are moderately consistent with being generated by background processes. The upper limit ${\mathscr{U}}$ can be decreased so that $p({\mathbb{E}}[{b}_{{\rm{X}}}]| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ is more informative at the risk of bias. On the other hand, ${\mathscr{U}}$ can be increased, which naturally weakens the constraining power but can be justified as a safety precaution to capture a contribution such as a power-law component originating from the magnetosphere of PSR J0740+6620 that is not explicitly modeled. ⁵⁹ Alternatively, we could justify a higher upper limit in terms of systematic error due to blank-sky estimates derived from an ensemble of exposures over the sky and variation in time of the XMM camera response matrices between the PSR J0740+6620 exposure and blank-sky exposure ensemble. The setting of n = 4 seems like a reasonable—albeit arbitrary—choice to balance information loss versus bias; we probe posterior sensitivity to this hyperparameter and conclude our inferences are insensitive to its value (see Section 3). We display the XMM pn camera count number spectrum in Figure 4 together with the blank-sky-derived prior PDF of the expected count-rate variables $\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\};$ we also display the NICER count-number spectrum for data quality comparison.

Zoom In Zoom Out Reset image size

The likelihood function given the NICER and XMM data sets may be written as

$$\begin{eqnarray}&&\begin{array}{l}p\left({d}_{{\rm{N}}},{d}_{{\rm{X}}},\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}| s,\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\},\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}\right)\\ \ =\ p\left({d}_{{\rm{N}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\},{\rm\small{NICER}}\right)p\left({d}_{{\rm{X}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\},{\rm\small{XMM}}\right)\\ \ \,\times \ p\left(\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}| \left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where we combine the expected background count rate variables over the XMM cameras into $\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\}$, we combine the count numbers in the PSR J0740+6620 exposures into d_X, and we combine the blank-sky count numbers into $\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}$. Introducing the flat prior densities from Equation (1), approximating the posterior PDF $p(\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}| \left\{{{\mathscr{B}}}_{{\rm{X}}}\right\})$ as flat and bounded as described in Section 2.5.2, and marginalizing over all expected background count rate variables yields the background-marginalized likelihood function

$$\begin{eqnarray}&&\begin{array}{l}p\left({d}_{{\rm{N}}},{d}_{{\rm{X}}},\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}| s\right)\propto {\displaystyle \int }_{\{0\}}^{\{{ \mathcal U }\}}p\left({d}_{{\rm{N}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\},{\rm\small{NICER}}\right)\\ \ \ \times \ d\left\{{\mathbb{E}}[{b}_{{\rm{N}}}]\right\}{\displaystyle \int }_{\{{\mathscr{L}}\}}^{\{{\mathscr{U}}\}}p\left({d}_{{\rm{X}}}| s,\left\{{\mathbb{E}}[{b}_{{\rm{X}}}]\right\},{\rm\small{XMM}}\right)d\left\{{\mathbb{E}}[{{\mathscr{B}}}_{{\rm{X}}}]\right\}.\end{array}\end{eqnarray} \tag{ 5 }$$

This function is fed as a callback to a sampling process, together with a joint prior PDF callback for the pulsar signal parameters s and parameters associated with the nicer and xmm instrument response models.

The informative joint NANOGrav and CHIME/Pulsar prior is critical for deriving a useful constraint on the radius of PSR J0740+6620. For comparison, we compute a posterior conditional on a fully ionized hydrogen atmosphere and a diffuse, separable prior PDF of the mass, the distance, and the cosine of inclination angle. The mass prior is such that the joint prior PDF of mass and radius is flat within the prior support (see Table 1). The distance prior PDF is adopted from Igoshev et al. (2016) and displayed in Figure 2, with support D ∈ [0.1, 10.0] kpc. The prior PDF of the cosine of the inclination angle is isotropic, meaning flat with support $\cos i\in [0,1]$.

Table 1. Summary Table for ST-U nsx Fully Ionized Hydrogen Hot Regions Plugged into the NICER × XMM Likelihood Function, Conditional on the NANOGrav × CHIME/Pulsar Prior PDF

Parameter	Description	Prior PDF (Density and Support)	${\widehat{\mathrm{CI}}}_{68 \% }$	${\widehat{D}}_{\mathrm{KL}}$	$\widehat{\mathrm{ML}}$
P (ms)	coordinate spin period	P = 2.8857, a fixed	⋯	⋯	⋯
M (M⊙)	gravitational mass	$M,\cos (i)\sim N({{\boldsymbol{\mu }}}^{\star },{{\boldsymbol{\Sigma }}}^{\star })$	${2.072}_{-0.066}^{+0.067}$	0.01	2.070
	joint prior PDF N( μ ⋆, Σ⋆)	μ ⋆ = [2.082, 0.0427]⊤
		${{\boldsymbol{\Sigma }}}^{\star }=\left[\begin{array}{c}{0.0703}^{2}\quad {0.0131}^{2}\\ {0.0131}^{2}\quad {0.00304}^{2}\end{array}\right]$
Req (km)	coordinate equatorial radius	Req ∼ U(3rg(1),16) b	${12.39}_{-0.98}^{+1.30}$	0.58	11.02
	compactness condition c	Rpolar/rg(M) > 3
	effective gravity condition d	$13.7\leqslant {\mathrm{log}}_{10}{g}_{\mathrm{eff}}(\theta )\leqslant 15.0$, ∀ θ
Θp (rad)	p region center colatitude	$\cos ({{\rm{\Theta }}}_{p})\sim U(-1,1)$	${1.35}_{-0.39}^{+0.46}$	0.25	1.622
Θs (rad)	s region center colatitude	$\cos ({{\rm{\Theta }}}_{s})\sim U(-1,1)$	${1.89}_{-0.46}^{+0.40}$	0.24	2.303
ϕp (cycles)	p region initial phase e	ϕp ∼ U(−0.5, 0.5), wrapped f	bimodal	3.52	0.185
ϕs (cycles)	s region initial phase g	ϕs ∼ U(−0.5, 0.5), wrapped	bimodal	3.51	0.243
ζp (rad)	p region angular radius	ζp ∼ U(0, π/2)	${0.147}_{-0.041}^{+0.070}$	2.18	0.093
ζs (rad)	s region angular radius	ζs ∼ U(0, π/2)	${0.146}_{-0.042}^{+0.071}$	2.15	0.127
	no region-exchange degeneracy	Θs ≥ Θp
	non-overlapping hot regions	function of (Θp , Θs , ϕp , ϕs , ζp , ζs )
${\mathrm{log}}_{10}\left({{ \mathcal T }}_{p}\,[{\rm{K}}]\right)$	p region nsx effective temperature	${\mathrm{log}}_{10}\left({{ \mathcal T }}_{p}\right)\sim U(5.1,6.8)$, nsx limits	${5.988}_{-0.059}^{+0.048}$	2.95	6.080
${\mathrm{log}}_{10}\left({{ \mathcal T }}_{s}\,[{\rm{K}}]\right)$	s region nsx effective temperature	${\mathrm{log}}_{10}\left({{ \mathcal T }}_{s}\right)\sim U(5.1,6.8)$, nsx limits	${5.992}_{-0.058}^{+0.047}$	2.98	6.058
$\cos (i)$	cosine Earth inclination to spin axis	$M,\cos (i)\sim N({{\boldsymbol{\mu }}}^{\star },{{\boldsymbol{\Sigma }}}^{\star })$	${0.0424}_{-0.0029}^{+0.0029}$	0.01	0.044
D (kpc)	Earth distance	D ∼ skewnorm(1.7, 1.0, 0.23) h	${1.21}_{-0.15}^{+0.15}$	0.10	0.995
NH [1020 cm−2]	interstellar neutral H column density	NH ∼ U(0,10)	${1.587}_{-1.092}^{+1.953}$	0.91	0.216
αNICER	NICER effective-area scaling	αNICER, αXMM ∼ N( μ , Σ)	${1.026}_{-0.137}^{+0.136}$	0.03	1.111
αXMM	XMM effective-area scaling	αNICER, αXMM ∼ N( μ , Σ)	${0.93}_{-0.13}^{+0.14}$	0.17	0.638
	joint prior PDF N( μ , Σ)	μ = [1.0, 1.0]⊤
		${\boldsymbol{\Sigma }}=\left[\begin{array}{c}{0.150}^{2}\quad {0.106}^{2}\\ {0.106}^{2}\quad {0.150}^{2}\end{array}\right]$
	Sampling process information
	number of free parameters: i 15
	number of processes: j 1
	number of live points: 4 × 104
	hypervolume expansion factor: 0.1−1
	termination condition: 10−1
	evidence: k $\widehat{\mathrm{ln}{ \mathcal Z }}=-20714.61\pm 0.02$
	number of core l hours: 28320
	likelihood evaluations: 23710136
	nested replacements: 1250386
	effective sample size: m 420281

Notes. The description in this caption is largely adopted from Riley et al. (2019) for consistency. We provide: (i) the parameters that constitute the sampling space, with symbols, units, and short descriptions; (ii) any notable derived or fixed parameters; (iii) the joint prior distribution, including hard truncation bounds and constraint equations that define the hyperboundary of the support; (iv) one-dimensional (marginal) 68.3% credible interval estimates symmetric in posterior mass about the median (${\widehat{\mathrm{CI}}}_{68 \% }$); (v) Kullback–Leibler (KL)-divergence estimates in bits (${\widehat{D}}_{\mathrm{KL}}$) representing prior-to-posterior information gain (see the appendix of Riley et al. 2019 for a high-level description of the divergence); (vi) the parameter vector ($\widehat{\mathrm{ML}}$) estimated to maximize the background-marginalized likelihood function, corresponding to a nested sample. Note that, strictly, the target of nested sampling is not to generate a maximum likelihood estimator—it is a by-product of the sampling process for evidence estimation. Moreover, there is degeneracy in the likelihood function and high-likelihood solutions with remarkably different parameter values—such as a radius near or above the posterior median—may be retrieved from the public sample information. Constraint equations in terms of two or more parameters result in marginal distributions that are not equivalent to those inverse-sampled.

^a Cromartie et al. (2020); Wolff et al. (2021). ^b The function r_g(M) $:=$ GM/c² denotes the gravitational radius with dimensions of length. ^c The coordinate polar radius of the source 2-surface, R_polar(M, R_eq, Ω), is a quasi-universal function adopted from AlGendy & Morsink (2014), where Ω $:=$ 2π/P is the coordinate angular rotation frequency. ^d The range of effective gravity from the equator (minimum gravity) to the pole (maximum gravity) must lie within nsx limits. A quasi-universal function is adopted from AlGendy & Morsink (2014) for effective gravity g_eff(θ; M, R_eq, Ω) in units of cm s⁻² in the table, where Ω $:=$ 2π/P is the coordinate angular rotation frequency. ^e With respect to the meridian on which Earth lies. ^f The periodic boundary is admitted and handled by MultiNest. However, this is an unnecessary measure because we straightforwardly define the mapping from the native sampling space to the space of a phase parameter ϕ such that the likelihood function maxima are not in the vicinity of this boundary. ^g With respect to the meridian on which the Earth antipode lies. ^h Specifically, the PDF defined as scipy.stats.skewnorm.pdf(D, 1.7, loc = 1.002, scale = 0.227), truncated to the interval D ∈ [0,1.7] kpc. ⁱ In the sampling space; the number of background count rate variables is equal to the number of channels defined by the NICER and XMM data sets. ^j The mode-separation MultiNest variant was deactivated, meaning that isolated modes are not evolved independently and nested sampling threads contact multiple modes. In principle this also allows us to combine the processes in a post-processing phase using nestcheck (Higson 2018), if more than one process is available for a given posterior; the posteriors derived in the production analysis are high-resolution, so we neglect combining repeat processes. ^k Defined as the prior predictive probability $p({d}_{{\rm{N}}},{d}_{{\rm{X}}},\left\{{{\mathscr{B}}}_{{\rm{X}}}\right\}\,| \,{\mathtt{ST}}-{\mathtt{U}})$. Note, however, that in order to complete the reported evidence for comparison to models other than those defined in this work, upper-bounds for the NICER background parameters need to be specified. ^l Intel^®Xeon^® E5-2697A v4. ^m The effective sample size estimator invoked, following DNest4 (Brewer & Foreman-Mackey 2016, https://github.com/eggplantbren/DNest4), is the perplexity measure$\begin{eqnarray*}&&\widehat{\mathrm{ESS}}:= \exp \left(-\sum _{i}^{I}{w}_{i}\mathrm{log}{w}_{i}\right),\end{eqnarray*}$ where the {w_i }_i=1,...,I are the sample weights (e.g., Martino et al. 2017).

Download table as:  ASCIITypeset image

The radius posterior conditional on the diffuse prior is much broader than when we condition on the joint NANOGrav and CHIME/Pulsar prior, and there is no independent indication from the pulse-profile modeling for a high mass (see Figure 5). Once the models are conditioned on the joint NANOGrav and CHIME/Pulsar prior PDF, we gain very little additional information a posteriori from the X-ray likelihood function about the pulsar mass, distance, and inclination angle with respect to the spin axis. We can verify this by examining the marginal posterior distributions in comparison to the respective prior distributions, which is summarized for each parameter by the KL-divergence estimate (Kullback & Leibler 1951).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (7 panels)

Tarball of figure set images

The NICER likelihood function is sensitive to the basic configuration of the surface hot regions and their temperatures, despite the relatively small number of pulsed counts (those above the phase-invariant background). The XMM likelihood function is less informative both due to the lack of phase information, and because the events are sparse for all three EPIC cameras and moderately consistent with the expected background signal derived from blank-sky exposures. However, the XMM likelihood function acts to reduce the NICER posterior mode volume substantially, affecting the inferred radius and geometry. The reason for this is because the XMM likelihood is sensitive to the combined phase-averaged signal from the hot regions being too bright. Therefore, given the NICER likelihood function, we constrain the contribution to the unpulsed portion of the pulse profile that must be generated by the hot regions rather than the backgrounds. Posterior figures demonstrating this are reserved for Section 3.3.

The atmospheric composition of PSR J0740+6620 is not known a priori. We therefore compared ST-U posteriors for hydrogen and helium atmospheres assuming full ionization—see the discussion in Section 2.3.2—in the second online figure of the set associated with Figure 5. The marginal radius posteriors were indistinguishable, although there were some small changes in the properties of the hot regions. However, given the apparent lack of sensitivity to atmospheric composition, inferences reported hereafter are conditioned on a fully ionized hydrogen atmosphere—we do not need to marginalize over the binary atmosphere parameter. Both hot regions are inferred to have effective temperatures T ≈ 10⁶ K, at which partial ionization effects should be small.

The Riley et al. (2019) analysis of PSR J0030+0451 reported a number of hot region models that provided an adequate description of the data according to their performance measures (largely adopted here). These included ST-U and variants in which one of the hot regions was permitted increasingly complex forms, including rings and crescents. The inferred radius changed as model complexity increased, but evidence calculations showed a substantial improvement in model performance. As a result of this, we deemed the ST+PST model—in which one hot region is a single temperature spherical cap and the other is, a posteriori, a crescent—to be superior for PSR J0030+0451.

For PSR J0740+6620 the ST-U model also provides an adequate description of the data. In order to assess whether additional complexity is useful we also condition on the ST+PST model. The posterior configuration and properties of the hot regions conditional on this more complex model (which includes the possibility of hot regions that are simply spherical caps) did not differ in an important way from the configuration inferred from ST-U: the hot region for which more complexity was permissible exhibited degeneracy a posteriori—we were not sensitive to the existence of additional emission structure beyond that of a simple spherical cap, and the evidence estimates are consistent. There was therefore no extended crescent as inferred for PSR J0030+0451; the likelihood function degeneracy included some subset of possible crescent structures—those on smaller angular scales (see Riley et al. 2019 for a discussion about hot region structure degeneracy)—which may be of interest to pulsar modelers. The inferred radius changed very little (see the third online figure of the set associated with Figure 5), and there was no increase in model performance. For this reason we hereafter report inferences exclusively for the ST-U model.

As described in Section 2.5.2, the XMM background is free-form, but each variable (one per channel) has a prior with compact support. For each variable, a flat prior PDF is defined whose width is controlled by a hyperparameter n. For the headline inferences reported in this Letter we used n = 4, having tested sensitivity to varying n in the range n ∈ [0.01, 8]. In the limit that n tends to zero, the background information would be treated as a point estimate of the XMM background. The posterior distribution of the radius is insensitive to n being varied through the range n ∈ [0.01, 4]; see the fourth online figure of the set associated with Figure 5. It broadens slightly for n = 8 because fainter combined signals from the hot regions have greater background-marginalized likelihoods, yielding additional posterior weight for higher-radius configurations that reduce the unpulsed component while conserving the pulsed component to satisfy the NICER event data. However, the value n = 8 is arguably too conservative even when considering potential systematic error.

The X-PSI likelihood function has a number of resolution settings, most notably settings that control the discretization of the computational domain for computation of signals (pulse profiles) incident on telescopes. The photon-specific flux signal we require is an integral over a distant observer's sky of the photon specific intensity from the hot regions, yielding a two-dimensional function of time (rotational phase) and photon energy. The level of discretization with respect to four variables in the domain of the incident photon specific intensity field generally controls the computational expense of likelihood evaluation. These variables are the number of rotational phases and energies the photon-specific flux signal is computed at, the number of hot region surface elements, and the number of rays calculated. See the X-PSI documentation ⁶³ for additional information.

We tested posterior sensitivity to increasing the discretization degrees for these variables by recomputing a posterior PDF with a new nested sampling process. We found that doubling these discretization degrees does not yield a change in the posterior PDF that is clearly distinguishable from Monte Carlo sampling noise; ⁶⁴ see the fifth online figure of the set associated with Figure 5.

For a fixed bounding hypervolume expansion factor of 0.1⁻¹, the posterior PDFs were insensitive to doubling live-point number from 10³ up to 4 × 10³. Following sampler comparisons within the NICER collaboration, we then increased resolution to 4 × 10⁴ live points, leading to broadening of the radius posterior; see the sixth and seventh online figures in the set associated with Figure 5. Increasing the sampling resolution by using 8 × 10⁴ live points led to some further broadening, but doubled an already large computational cost. Given the computational resources required for posterior computation with such a large number of live points, we were not able to rigorously prove convergence with live-point number. We decided to adopt 4 × 10⁴ live points for the production analysis; all information necessary to reproduce and improve upon our posterior computation is available in open-source repositories.

Such posterior mode broadening with increased nested sampling resolution is naturally expected because nested sampling algorithms approximate hypersurfaces in parameter space of constant likelihood; these approximations improve with sampling resolution but their sufficiency is difficult to prove for non-trivial likelihood functions encountered in real problems and when subject to resource limitations. It is desirable to transform away nonlinear modal degeneracies so that an approximation conforms more efficiently to structure in the sampling space; however, this can in practice be an intractable task for a given problem. Moreover, when sampling from a target distribution with two or more modes of commensurate posterior mass, the live-point resolution is roughly split between the modes, reducing the resolution of a given mode due to the approximations alluded to above. For PSR J0740+6620, posterior bimodality arises due to the near-equatorial inclination of the source, leading to two competitive geometric configurations of the hot regions (see Section 3.3, where we discuss this further).

The NICER background is difficult to estimate directly, but there are two tools available. Figure 15 shows the NICER background estimated using the "space weather" model (Gendreau 2020) and the "3C50" model (Remillard et al. 2021); see also Bogdanov et al. (2019a). The former models background due to the space weather environment, which varies as NICER moves through different geomagnetic latitudes, and depends on solar activity. The "3C50" model is empirical, taking into account two different types of particle-induced events, the cosmic X-ray background, and a soft X-ray noise component related to operation in sunlight. We also render a set of NICER background estimates for comparison: using joint NICER and XMM posterior samples, we display in Figure 15 the NICER background spectrum that maximizes the conditional likelihood function, yielding a band that is a proxy for background variable posterior mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

The background spectra displayed in Figure 15 exceed the space weather model at low energies. This excess is not unreasonable given the presence of multiple other point sources in the NICER field of view. ⁶⁷ In higher channels, the background spectra appear to agree well with the space weather model. There is a possible small systematic under-prediction compared to the space weather model in channels 80−100; the level of agreement is satisfactory given the current uncertainties on the background modeling. The level of agreement with the 3C50 model, which exceeds the space weather model at low energies and is consistent with it at higher energies, also appears good. Neither background model includes off-axis X-ray sources that might contaminate the target spectrum, and therefore inferred backgrounds that exceed the two estimates are not in principle a problem.

Once the uncertainties on the two NICER background models are understood more fully, we anticipate being able to use them to reduce systematic error in the pulse-profile modeling. The current radius measurement conditional on the XMM data set—which yields an indirect NICER background constraint—will therefore be superseded. Efforts are also ongoing to quantify the level of background due to any off-axis sources in the field of view.

In Figure 16 we display the XMM pn background spectra that maximize the conditional likelihood function given nested samples from the joint NICER and XMM posterior, together with supplementary information. Graphical checking of the spectra against the blank-sky estimate and the event data does not reveal any problems a posteriori.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images