We report Bayesian inference of the mass radius and hot Xray emitting region propertiesusing data from the Neutron Star Interior Composition ExploreR NICERfor the brightest rotationpowered millisecond Xray pulsar PSR J04374715. Our modeling is conditional on informative tight priors on mass distance and binary inclination obtained from radio pulsar timing using the Parkes Pulsar Timing Array PPTA Reardon et al. and we use NICER background models to constrain the nonsource background crosschecking with data from XMMNewton. We assume two distinct hot emitting regions and various parameterized hot region geometries that are defined in terms of overlapping circles while simplified these capture many of the possibilities suggested by detailed modeling of return current heating. For the preferred model identified by our analysis we infer a mass of M  1.418  0.037 M SUBSUB largely informed by the PPTA mass prior and an equatorial radius of inlineformula mmlmath overflowscrollmmlmiRmmlmimmlmommlmommlmsubsupmmlmrowmmlmn11.36mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.63mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.95mmlmnmmlmrowmmlmsubsupmmlmath inlineformula km each reported as the posterior credible interval bounded by the 16 and 84 quantiles. This radius favors softer dense matter equations of state and is highly consistent with constraints derived from gravitational wave measurements of neutron star binary mergers. The hot regions are inferred to be nonantipodal and hence inconsistent with a pure centered dipole magnetic field.
['Neutron stars', 'Pulsars', 'Rotation powered pulsars', 'Millisecond pulsars', 'Neutron star cores', 'Nuclear astrophysics', 'Magnetic fields', '1108', '1306', '1408', '1062', '1107', '1129', '994', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Solar and Stellar Astrophysics', 'Nuclear Theory']
A NICER View of the Nearest and Brightest Millisecond Pulsar PSR J04374715
{'A NICER VIEW OF THE NEAREST AND BRIGHTEST MILLISECOND PULSAR: PSR J0437 -4715': "Devarshi Choudhury, 1 Tuomo Salmi, 1 Serena Vinciguerra, 1 Thomas E. Riley, 1 Yves Kini, 1 Anna L. Watts, 1 Bas Dorsman, 1 Slavko Bogdanov, 2 Sebastien Guillot, 3, 4 Paul S. Ray, 5 Daniel J. Reardon, 6, 7 Ronald A. Remillard, 8 Anna V. Bilous, 9 Daniela Huppenkothen, 10 James M. Lattimer, 11 Nathan Rutherford, 12 Zaven Arzoumanian, 13 Keith C. Gendreau, 13 Sharon M. Morsink, 14 and Wynn C. G. Ho 15  \n1 Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090GE Amsterdam, the Netherlands 2 Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA 3 IRAP, CNRS, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France 4 Universit'e de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France. 5 Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375, USA 6 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia 7 OzGrav: The Australian Research Council Centre of Excellence for Gravitational Wave Discovery, Hawthorn VIC 3122, Australia 8 MIT Kavli Institute for Astrophysics & Space Research, MIT, 70 Vassar St., Cambridge, MA 02139, USA 9 ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA Dwingeloo, The Netherlands 10 SRON Netherlands Institute for Space Research, Niels Bohrlaan 4, 2333 CA Leiden, the Netherlands 11 Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA 12 Department of Physics and Astronomy, University of New Hampshire, Durham, New Hampshire 03824, USA 13 X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA 14 Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB, T6G 2E1, Canada 15 Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA  \n(Received; Revised; Accepted; Published)", 'ABSTRACT': 'Wereport Bayesian inference of the mass, radius and hot X-ray emitting region properties - using data from the Neutron Star Interior Composition ExploreR (NICER) - for the brightest rotation-powered millisecond X-ray pulsar PSR J0437 -4715. Our modeling is conditional on informative tight priors on mass, distance and binary inclination obtained from radio pulsar timing using the Parkes Pulsar Timing Array (PPTA) (Reardon et al. 2024), and we use NICER background models to constrain the non-source background, cross-checking with data from XMM-Newton. We assume two distinct hot emitting regions, and various parameterized hot region geometries that are defined in terms of overlapping circles; while simplified, these capture many of the possibilities suggested by detailed modeling of return current heating. For the preferred model identified by our analysis we infer a mass of M = 1 . 418 ± 0 . 037 M ⊙ (largely informed by the PPTA mass prior) and an equatorial radius of R = 11 . 36 +0 . 95 -0 . 63 km, each reported as the posterior credible interval bounded by the 16% and 84% quantiles. This radius favors softer dense matter equations of state and is highly consistent with constraints derived from gravitational wave measurements of neutron star binary mergers. The hot regions are inferred to be non-antipodal, and hence inconsistent with a pure centered dipole magnetic field.  \nKeywords: dense matter - equation of state - pulsars: general - pulsars: individual (PSR J0437 -4715) - stars: neutron - X-rays: stars', '1. INTRODUCTION': "NICER, the Neutron Star Interior Composition Explorer (Gendreau et al. 2016), collects high time resolu-  \nCorresponding author: D. Choudhury  \nD.Choudhury@uva.nl  \ntion X-ray spectral-timing data of rotation-powered millisecond pulsars (MSPs) to build up high quality phase and spectrally resolved pulse profiles (Bogdanov et al. 2019a). Relativistic effects due to the neutron star's rapid spin and strong gravitational field affect the radiation emitted from the hot magnetic polar caps of the pulsar in a well-understood fashion (see Bogdanov et al. 2019b, and references therein). Via Pulse Profile Mod-  \neling (PPM) we can thus infer the pulsar's mass and radius as well as geometric properties such as the size and location of the emitting polar caps (Watts 2019; Bogdanov et al. 2021). Mass and radius measurements constrain the Equation of State (EoS) of the supranuclear dense matter in neutron star cores (see e.g. Lattimer & Prakash 2007), while the hot polar cap properties map the magnetic field structure (e.g. Bilous et al. 2019).  \nTo date the NICER collaboration have reported PPM results for two MSPs: PSR J0030+0451 (Miller et al. 2019; Riley et al. 2019; Salmi et al. 2023; Vinciguerra et al. 2024) and the heavy pulsar PSR J0740+6620 (Miller et al. 2021; Riley et al. 2021; Salmi et al. 2022, 2023, 2024; Dittmann et al. 2024). Radius uncertainties for both are at present ∼ ± 10 % (68% credible interval), with this expected to improve as more photons are accumulated (Lo et al. 2013; Psaltis et al. 2014). PSR J0030+0451 is the second brightest source in the NICER PPM target list, but is an isolated pulsar for which there is no prior information on the mass. PSR J0740+6620 is in a binary and has a wellconstrained mass derived from radio pulsar timing (Fonseca et al. 2021) that is used as a prior in the pulse profile modeling of NICER data. With a high neutron star mass ∼ 2 . 1 M ⊙ , it is a very interesting source in terms of potential constraints on dense matter. However, being more distant, it is fainter, which makes it more challenging to obtain sufficient photons for tight constraints on radius. NICER's results for these two sources have already started to place interesting constraints on both dense matter (see, e.g., Miller et al. 2021; Raaijmakers et al. 2021; Biswas 2022; Annala et al. 2023; Tak'atsy et al. 2023) and pulsar magnetic field geometry and emission mechanisms (see, e.g., Bilous et al. 2019; Chen et al. 2020; Kalapotharakos et al. 2021; Carrasco et al. 2023). NICER PPM analysis has been supported by a programme of testing and simulation (Bogdanov et al. 2019b, 2021; Vinciguerra et al. 2023; Choudhury et al. 2024).  \nThe rotation powered pulsar PSR J0437 -4715, which is the focus of this Letter, was discovered by Johnston et al. (1993) and is the closest observed MSP and the brightest source in the NICER PPM target list. Its brightness and proximity allows for precise measurements of astrometric and orbital parameters from radio pulsar timing, and it is possible to obtain a tightlyconstrained distance value of 156 . 98 ± 0 . 15 pc (Reardon et al. 2024). It has a spin frequency of 174 Hz and forms a binary system with a 0.2 M ⊙ helium-core white dwarf (Bailyn 1993). The binary system has an inclination of 137 . 5 · , which is favorable for measuring the Shapiro delay (Shapiro 1964) in arriving radio pulses, which in turn provides us with an independent mass constraint. Reardon et al. (2024) report a tight constraint on the mass of PSR J0437 -4715 with a value of M = 1 . 418 ± 0 . 044 M ⊙ , which is used to inform the models in this work.  \nThe Letter is organized as follows: In Section 2, we introduce the different data sets and associated background estimates for PSR J0437 -4715. In Section 3 we summarize the various aspects of building a pulse profile model and introduce the unique aspects required for modeling PSR J0437 -4715, building upon the formalism described in previous NICER analyses using the pulse profile modeling code X-PSI (Riley et al. 2023). Section 4 provides details of the sampling process. We document our exploratory analyses in Section 5 where we test the effects of various modeling choices. In Section 6, we provide the results for our final model and showcase our headline results. We discuss the physical implications of the inferred model parameters by our headline model in Section 7 and conclude in Section 8.", '2.1. NICER XTI': "For this work, observations of PSR J0437 -4715 conducted by the NICER X-ray Timing Instrument (XTI) have been used to generate two different datasets. The datasets vary in the pre-processing techniques applied which are detailed in the subsections below. First we describe the common features present in both datasets.  \nA bright Seyfert II Active Galactic Nuclei (AGN), RX J0437.4 -4711 (Halpern & Marshall 1996) is located within 4 ' .18 of our target (See e.g., Figure 8 of Bogdanov et al. 2019a). At this distance from the telescope aim point, ∼ 8% of the AGN flux is detected, therefore providing a non-negligible contribution to the total astrophysical background. Bogdanov et al. (2019a) developed an optimization technique to minimize background contamination from nearby point sources by finding the position that maximize the target's Signal-to-Noise Ratio (SNR, as defined in Equation 3 of Bogdanov et al. 2019a). For PSR J0437 -4715, the optimal pointing is 1 ' .15 south-west of the pulsar, which improved the SNR by 16% compared to directly pointing at the pulsar. With this pointing, ∼ 2 . 25% of the AGN flux (i.e., ∼ 0 . 2 counts s -1 ) falls in the Field of View (FoV) of NICER and contributes to the background 1 .  \nFor our analysis, we only consider those events registered in the Pulse Invariant 2 (PI) channel subset [30, 300), nominally corresponding to the photon energy range 0.3-3.0 keV. The sharp low-energy threshold cutoff in the instrument calibration products coupled with the presence of increased electronic noise induced by optical light at energies below 0.3 keV can lead to biased inferences. To prevent this, we ignore the events reg-  \nFigure 1. Phase-folded PSR J0437 -4715 (3C50) event data shown over two rotational cycles. We use 32 phase intervals (bins) per cycle and the count numbers in bins separated by one cycle are identical for a given channel(s). The sum over all phase-channel pairs yields the total counts. The top panel displays the waveform obtained by summing the channels in the contiguous subset [30, 300). The bottom panel displays the phase-channel-resolved count numbers for the same channel subset, where the color bar represents counts per channel per two cycles.  \n<!-- image -->  \nistered below channel 30 altogether. The contribution of PSR J0437 -4715 hot regions to the counts registered above channel 300 are negligible compared to the counts generated by background processes, and are therefore also neglected.  \nAll the registered events accounted for are phasefolded based on the ephemeris obtained from the PPTA radio timing solutions (Reardon et al. 2024), using the PINT 3 photonphase tool. A corresponding pulse profile thus generated (using methods described in Section 2.1.2) is shown in Figure 1. Beside a prominent pulse, the profile also features a diffuse hump around phases 0.5-0.7 which Bogdanov (2013) suggests is a consequence of non-antipodal hot regions on the surface. This possibility is explored in our models.  \n2.1.1. Delta dataset, space weather background estimate, and AGN spectrum  \nThe 'Delta' dataset (consisting of all NICER pointings in the time range 2017 July 6 to 2021 July 31). was constructed using HEASoft version 6.29c and CALDB  \nversion xti20210707 . We first applied the standard NICER L2 processing, which applies basic filters ensuring that NICER is pointed at the target, is outside the South Atlantic Anomaly (SAA) region, and not pointed within 15 · of the Earth limb or within 30 · of the bright Earth, as well as standard cuts on the overshoot rate to remove time intervals with high particle background. As this analysis is sensitive to the background rate and optical loading, we applied further cuts including: (1) cutoff rigidity ( COR SAX ) > 1 . 5 GeV/ c , (2) mask DET ID 34, (3) planetary K index ( K p ) < 5, (4) undershoot rate < 200 counts s -1 . Lastly, we removed time intervals based on count rates, first by filtering out times where DET ID 14 had a count rate greater than 1.0 counts s -1 in 8.0 s bins, then by filtering on total X-ray count rate (0.25-8.0 keV) above 6.0 counts s -1 in 2.0 s bins. This catches times of high background that are not caught by proxies such as K p or high overshoot rate. The final good time in this dataset is 2.32 Ms.  \nFor this dataset, we generated a background estimate using an internal version of the 'Space Weather' background estimator tool 4 . This model uses a library of blank sky observations taken over a range of conditions of space weather ( K p ) and the magnetic cutoff rigidity ( COR SAX ) which are combined to match the conditions of the observation in question. To this, a soft component representing the optical loading contamination is added based on the sun angle distribution in the observation. This background estimate represents the blank sky and does not include any contribution from contaminating sources in the field of view or non-spot emission from the pulsar.  \nFor this dataset, we also made an estimate of the contribution to the spectrum from the AGN. To do this, we created a new response matrix using the same method above but specifying the AGN position instead of the pulsar. The effective area for the AGN is a factor of 3050 less than for the pulsar, depending on the individual detector misalignments. We constructed a model spectrum of the AGN from an absorbed double powerlaw fit to XMM-MOS1 data 5 ( N H = 1 . 9(5) × 10 20 cm -26 , Γ 1 = 2 . 9(1), N 1 = 1 . 7(2) × 10 -3 photons/keV/cm 2 /s, Γ 2 = 1 . 4(1), N 2 = 6 . 6(17) × 10 -4 ) photons/keV/cm 2 /s. We also extracted spectra of the AGN from 45 NICER ObsIDs where the AGN was the target. We fit one long observation (ObsID 4060180109, 8.9 ks of good exposure) to an absorbed powerlaw ( N H = 1 . 8(4) × 10 20 cm -2 , Γ = 2 . 42(5), N = 1 . 84(6) × 10 -3 photons/keV/cm 2 /s) and measured the flux. Using the measured XMM and  \n4 https://heasarc.gsfc.nasa.gov/docs/nicer/tools/ nicer bkg est tools.html  \n5 We used the XMM MOS1 imaging data from ObsID 0603460101, acquired on 2009-12-15 at 19:41:41 for a total filtered exposure of 119.1 ks.  \nFigure 2. Count rate spectrum of the Delta dataset, along with the Space Weather background estimate and a representative range of the attenuated AGN spectra. Green is the estimated AGN spectrum for NICER based on XMMNewton observations from 2009 and the red band shows the range of fluxes observed in NICER pointings at the AGN.  \n<!-- image -->  \nNICER spectra, we constructed an estimated counts spectrum for the AGN. According to this, the AGN contributes about 0.2 counts s -1 in the 0.3-3 keV band, which is noticeably less than the pulsar at about 1 counts s -1 . However the AGN is known to vary, so this estimate is only representative. We do not specifically know the flux state of the AGN at the time of each of the NICER observations of the pulsar. Figure 2 shows the observed count rate spectrum from the Delta dataset, along with the Space Weather background estimate and and the attenuated AGN spectra. The NICER spectrum is shown as a range covering the lowest to highest fluxes observed in the individual NICER ObsIDs.", '2.1.2. 3C50 dataset, background model, and AGN spectrum': "The Delta dataset is processed to minimize background contamination and provide a space weather estimate. However, this estimate does not have any quantifiable statistical and systematic uncertainty associated with it. For our analysis using the space weather estimate, we therefore only use it to place a conservative lower limit on the background (see Section 4).  \nA substantial portion of the non-source background in a typical NICER observation originates from the local environment of the telescope. The NICER team has developed several approaches for modeling this timedependent background emission. One is the 3C50 model (based on a library of blank fields, Remillard et al. 2022), which makes use of the NICER event data for the source of interest to produce an estimate of the background level as a function of energy.  \nThe 3C50 dataset is based on a particular list of Good Time Intervals (GTIs) that pass filtering tests designed to diagnose when the 3C50 background model is per-  \nforming optimally. The 3C50 model is defined in Remillard et al. (2022), and the analysis and filtering procedures are described in detail in Salmi et al. (2022). The filtering steps reduce the usable exposure time, but this method provides access to quantified uncertainties in the net background spectrum, which is dominated by systematic effects that are otherwise difficult to determine. The investigation of seven rotation-powered X-ray pulsars which have NICER archival exposures of 1-3 Ms revealed that the root-mean-square (rms) deviations in the net source intensity (0.3-2.0 keV), when the data is binned at 200-400 ks, is very close to 0.018 counts s -1 (1 σ ) for each pulsar, regardless of the average pulsar intensity (Salmi et al. 2022). Given the defensible assumption that these pulsars are stable in X-ray intensity, the common rms values are interpreted as a measure of the systematic uncertainty for the 3C50 background model, for filtered exposure times of 200 ks or longer. Furthermore, it was found that the rms deviations in net intensity, in different energy bands, are proportional to the average background rate in the same bands. This implies that the deviations are tied to the precise level of the average background spectrum, and it enables estimation of the background uncertainty in any energy band. The 3C50 filtering effort therefore constrains pulse fractions, to some degree, in the effort to physically model the pulsar emission. One caveat, here, is the possible presence of additional 'background' components, in the form of extra X-ray sources in the NICER FoV. Such considerations have to be made for each target, and they are particularly important in the case of NICER investigations of PSR J0437 -4715.  \nThe procedures to define GTIs and to filter results for the 3C50 dataset (Salmi et al. 2022) are briefly summarized, as follows. Exposure intervals for a given target are first determined with the HEASoft tool nimaketime , and the times are masked by ± 30 s at any sunshine on/off boundary. The remaining continuous intervals, each of duration T i , are selected if T i ≥ 200 s, and they are divided into N equal pieces if T i > 450 s, where N = int ( T i / 300 + 0 . 5). These steps provide a list of GTIs with an average exposure ∼ 300 s and exposure range 200-450 s.  \nFor each GTI, the raw NICER spectrum is extracted, a background prediction is obtained with the 3C50 model, and then six filtering tests are performed on the results. Two filters apply to parameters needed for the background model: nz < 220 counts s -1 , where nz is the noise rate normalized to 50 Focal Plane Modules (FPMs) at 0.0-0.25 keV, and ibg < 0 . 2 counts s -1 , where ibg is the normalized rate of good events at 15-18 keV (beyond the effective area limits of the optics). The next two filters (as described in Remillard et al. 2022) operate on the net background-subtracted spectrum, S 0 net < 0 . 15 counts s -1 and hbg net < 0 . 05 counts s -1 , where the corresponding energy bands are 0.2 to 0.3 keV for S 0 and 13-15 keV for hbg . Finally, when the net pulsar  \nspectrum shows extremely faint or undetectable intensity above 2.0 keV, we add filters in the C band (2-4 keV) and D band (4-12 keV): C net < 0 . 1 counts s -1 and D net < 0 . 3 counts s -1 .  \nThe dataset for PSR J0437 -4715 considered in this letter consists of all NICER pointings in the time range 2017 July 6 to 2021 October 11. The raw exposure ( nimaketime , with geometric filters, but no rate filters) is 2.736 Ms, and the yield is reduced to 1.328 Ms after all of these filtering steps are applied (therefore, lower than the final good time of the Delta dataset). Therefore the total number of counts in the final 3C50 dataset (2095698 counts) is lower than that of the Delta dataset (3863743 counts). In addition, 1.310 Ms of these GTIs uniformly utilize 50 FPMs (excluding only 14 and 34), and we sidestep the need to normalize the input data and vary the response files, for each GTI, by adopting the 50-FPM subset (4442 GTIs; 1.310 Ms) as the 3C50 dataset (see also Section 3.3). In the 3C50 dataset, the background-subtracted average intensity at 0.3-2.0 keV is 1.079 counts s -1 . Pertinent to filtering in the C and D bands, the net intensities there are 0.020 and -0.004 counts s -1 , respectively.  \nAs noted previously, the NICER observations of PSR J0437-4715 contain a known AGN, RX J0437.4 -4711, near the edge of the FoV. Contributions from this source can be regarded as an additional background component, with its own uncertainty. NICER conducted separate, on-axis observations of RX J0437.4 -4711 during 2017 July 10-20, 2017 October 8-10, 2021 December 1026, and 2022 February 28 - March 19. The mean spectrum for the AGN was determined for each epoch, using the 3C50 background model and level 2 filtering applied to each GTI (Remillard et al. 2022). The average intensity (0.3-3.0 keV) is 6.26 counts s -1 , and the rms deviation between the four epochs is 25%. On the other hand the response calculator for NICER ( nicerarf ) determines that the AGN offset position in the primary observations of pulsar PSR J0437 -4715, integrated over the GTIs in the PSR J0437 -4715 dataset, has an effective area reduced by a factor of 42. The expected contribution from the AGN in the pulsar dataset is then represented by the mean AGN spectrum, rescaled by a factor of 42. The uncertainty range for this component is conservatively estimated as -50% to +100% of the values in this spectrum, based on NICER and XMM-Newton observations.", '2.2. XMM-Newton EPIC': "The XMM-Newton European Photon Imaging Camera (EPIC) event dataset used in this analysis was extracted from two archival observations of PSR J0437 -4715 acquired on 2002 October 09 (ObsID 0112320201) and 2009 December 15 (ObsID 0603460101), with total unfiltered on-source exposures of 69.4 and 129.4 ks, respectively. In this work, we only make use of the 'Full Frame' imaging mode data from the EPIC MOS1 and MOS2  \ninstruments (Turner et al. 2001). Although lacking high time resolution, they provide reliable low background phase-averaged source spectra desirable for the analyses described below. While the EPIC pn data acquired in 'Timing' mode allows the study of the pulsed emission, it is not used here due to the considerably larger uncertainties in the calibration of the instrument in this observing mode compared to MOS1 and MOS2 7 .  \nThe source event list extraction of the EPIC MOS data were conducted with the Science Analysis Software (SAS 8 ) version xmmsas 20211130 0941 . Following the same filtering procedure as in Bogdanov (2013), the event data were screened for instances of high particle background count rates and the recommended PATTERN ( ≤ 12 for MOS1/2) and FLAG (0) filters were applied. This resulted in 182.5 and 182.8 ks of clean effective exposure for the MOS1 and MOS2, respectively. The X-ray count rates of the pulsar in the cleaned event data from the two epochs are consistent within uncertainties, as expected from the steady thermal emission from rotation-powered MSPs. We extracted source and background counts from circular regions of radius 32 '' and 105 '' , respectively, to produce spectra in the form of counts per instrument channel. The SAS command backscale was then used to calculate the area of source and background regions used to extract the spectral files 9 . The value of this parameter is written into the header of the SPECTRUM table of the file as keyword BACKSCAL.", '3. PULSE PROFILE MODELING USING X-PSI': 'We make use of the X-ray Pulse Simulation and Inference (X-PSI) software package (Riley et al. 2023) 10 to develop models and perform Bayesian parameter estimation by analyzing the X-ray emission observed from PSR J0437 -4715. The core modeling procedure is extensively described in Riley et al. (2019, 2021); Bogdanov et al. (2019b, 2021). We summarize the process in this section and detail the novel components required in this work. Our analyses span X-PSI versions v0.7.0 and v2.2.2 , and the exact versions used in each of the inference runs can be found in the output files in the Zenodo repository 11 . The analysis files can be found on the Zenodo repository including the data products, the numeric model files, scripts, and Jupyter notebooks in the Python language to produce the results and figures  \n7 For details, see the XMM-Newton calibration documentation at https://www.cosmos.esa.int/web/xmm-newton/calibration  \n8 The XMM-Newton SAS is developed and maintained by the Science Operations Centre at the European Space Astronomy Centre and the Survey Science Centre at the University of Leicester.  \n9  \nSee https://www.cosmos.esa.int/web/xmm-newton/  \nsas-threads  \n- 10 https://github.com/xpsi-group/xpsi  \nof this article using the X-PSI framework, and posterior sample files.', '3.1. Radiation from source to telescope': "X-PSI performs relativistic ray-tracing of the neutron star surface emission using the Oblate Schwarzschild plus Doppler approximation (Poutanen & Gierli'nski 2003; Poutanen & Beloborodov 2006; Cadeau et al. 2007; Morsink et al. 2007; AlGendy & Morsink 2014; Bogdanov et al. 2019b). We define a flat prior density distribution on the equatorial radius ( R eq ), which is then shaped by limits imposed on the compactness and surface effective gravity as described in section 2.2.1 of Riley et al. (2021).  \nRadio-frequency timing of PSR J0437 -4715 provides precise measurements for some of our model parameters. The pulsar mass ( M ), distance ( D ), and orbital inclination angle ( i ) can be measured with much higher precision with this technique than with the Xray observations. We therefore inform our models using the latest pulsar timing analyses from the PPTA, described in Reardon et al. (2024). Specifically, we adopt the PPTA measurements of M , D , and i as Gaussiandistributed prior Probability Density Functions (PDFs) for our Bayesian inference.  \nThe pulsar distance is measured through the Shklovskii effect (Shklovskii 1970) on the orbital period ( P b ) of the pulsar, which induces an apparent time derivative, ˙ P b . The measured ˙ P b is corrected for contributions from differential acceleration in the Galaxy (circular and Galactic potential) and for the small contribution due to gravitational wave emission from the system. The resulting distance constraint is D = 156 . 98 ± 0 . 15 pc (Reardon et al. 2024), which is consistent with independent (but more uncertain) distances derived from other methods.  \nThe orbital inclination is derived from a measurement of annual-orbital parallax (Kopeikin 1995). This effect is now particularly useful because precise independent measurement of i breaks the strong covariance between the range (i.e., companion mass) and shape (sin i ) parameters that describe the Shapiro delay. Therefore, the measured Shapiro delay primarily informs the companion mass ( M c ), and a precise constraint for M then follows from the binary mass function. From these effects, Reardon et al. (2024) derive i = 137 . 506 ± 0 . 016 degrees and M p = 1 . 418 ± 0 . 044 M ⊙ (where the uncertainties represent the standard 68% confidence interval). The parameters can be considered uncorrelated.  \nFor our final fiducial results, we use priors corresponding to the values quoted above, which were derived from the recent PPTA third data release (PPTA-DR3) and its noise analysis (Zic et al. 2023; Reardon et al. 2023). However, model selection for the hot spot geometry (Section 5) was conducted using a slightly different prior for these parameters, which was derived from an earlier generation of PPTA data, described as  \nPPTA-DR2.5 in Reardon et al. (2024). The values used for these earlier analyses were D = 156 . 69 ± 0 . 30 pc, i = 137 . 485 ± 0 . 038 degrees, and M p = 1 . 411 ± 0 . 056 M ⊙ , which are highly consistent with the updated PPTADR3 measurements.  \nThe interstellar medium along the line of sight of the pulsar can get ionized by photoabsorption of the emergent soft X-rays thereby attenuating the signal generated by the source. The attenuation factor at any given energy is dependent on the neutral hydrogen column density N H , which forms a parameter in our models. The relative abundances from Wilms et al. (2000, updated in 2016) are adopted to compute tables, using the tbnew 12 model, relating the attenuation factor to energy for a given N H . For further details, we refer the reader to section 2.4 of Bogdanov et al. (2021).  \nVarious estimates are available for the N H along the line of sight of PSR J0437 -4715. Employing the relation between Dispersion Measure (DM) and N H from He et al. (2013), along with the DM measurement reported by Reardon et al. (2024), yields N H ≈ 0 . 8 × 10 20 cm -2 . The HEASARC mapping tool 13 also arrives at the same value. Estimates based on 3D E (B - V) extinction maps (E (B - V) ∼ 0 . 002 ± 0 . 014 using Stilism 14 ) (Lallement et al. 2014; Capitanio et al. 2017), coupled with the relation between extinction and N H provided by Foight et al. (2016), yield N H ≈ 0 . 2 × 10 20 cm -2 . Finally, spectral fits of UV and X-ray data to PSR J0437 -4715 by Gonz'alez-Caniulef et al. (2019) yield a range of 0 . 7 -2 . 4 × 10 20 cm -2 for the N H .  \nBased on these results, we initially tested different models with a wide and uniform prior U (0 . 01 , 5 . 0) × 10 20 cm -2 , and always inferred rather low values. Therefore, for our later runs, we opted for a narrower prior of U (0 . 004 , 2 . 0) × 10 20 cm -2 to better sample the lower tail.", '3.2. Surface emission': 'Since it is essentially impossible to spatially resolve a neutron star, we have to rely on phenomenological models of the surface radiation field, described in this section, to reproduce the observed X-ray event data. These models include sufficient complexity such that they replicate our basic understanding of possible magnetic field configurations, wherein the magnetospheric return currents deposit energy on the surface and heat it up (see e.g. Harding & Muslimov 2001). Note that internal heating by magnetic field decay is not significant for old and low field pulsars such as PSR J0437 -4715 and other MSPs (Gonzalez & Reisenegger 2010).', '3.2.1. Hot spot configurations': '12 https://pulsar.sternwarte.uni-erlangen.de/wilms/research/ tbabs/ 13 https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl 14 https://stilism.obspm.fr  \nFollowing the prescription laid down in Riley et al. (2019), we model the hot regions using (potentially compound) spherical caps on the neutron star surface. We start with simplistic models and increase complexity based on model performance, also checking for complexities to which the likelihood function is insensitive and is hence unhelpful.  \nAll our models involve two hot regions that are barred from overlapping with each other and are referred to as primary and secondary spots. Schematics of our various models can be found in figure 1 of Vinciguerra et al. (2023) . We report our analysis for the following increasingly complex models 15,16 :  \n- · ST-U : The Single Temperature -Unshared ( ST-U ) model consists of two disjoint simplyconnected spherical caps of uniform effective temperature for the atmosphere. In this particular model, the spot with the lower colatitude is designated as the primary. This condition aside, unshared here implies that the parameters of one spot are independent of the other.\n- · ST+PDT : Alongside the primary ST spot, we have a secondary Protruding Dual Temperature ( PDT ) spot consisting of two spherical overlapping emitting components. In the overlapping portion, emission of only one of the components is accounted for, and that is referred to as being the superseding component.\n- · CST+PDT : In this model, we have a Concentric Single Temperature ( CST ) spot as our primary, consisting of two spherical components, one emitting and one masking , thus forming a ring. The masking component (referred to as the omitting component in Vinciguerra et al. (2023, 2024)) is ensured to be concentric within the emitting component. The secondary spot here is a PDT .  \nAlthough further complexities can be introduced, we end our model sequence with CST+PDT given its ability to sufficiently explain the data (as explained in Section 6) and the computational expense that would be incurred upon thoroughly exploring more complex models.', '3.2.2. Atmosphere': '15 In our preliminary analyses on an initial dataset not reported in this letter, we also tested ST-S , CST-U , ST+PST (also tested on the Delta dataset), and PST-U configurations not included in this list. ST-S was utterly incapable of explaining the data, resulting in large residual structures. CST-U and ST+PST only led to very nominal model performance improvements over ST-U . PST-U was plagued by multiple degenerate modes and the sampler could not converge in a reasonable timescale.  \n16 We refer the reader to section 2.5 of Riley et al. (2019) and to sections 2.2 and 2.3.4 of Vinciguerra et al. (2023) for the full prior definitions of these models.  \nAll emitting hot region components incorporate a geometrically thin, fully-ionized hydrogen atmosphere model to describe the specific intensity of the local comoving effective temperature field. The intensity is interpolated over lookup tables of precomputed values, generated using NSX (Ho & Lai 2001), for a given effective temperature, surface gravity, photon energy and cosine of the emission angle (see Section 2.4.1. of Riley et al. (2019) for more details). Compared to Riley et al. (2019), this work utilizes an extended version of the lookup tables, as was done in all the later publications that used X-PSI. The limitations and uncertainties that this choice of atmosphere model implies are discussed in detail in Bogdanov et al. (2021); Salmi et al. (2023).', '3.3. Instrument response models': "For developing the NICER instrument response models, we utilize tailored Ancillary Response Files (ARFs) and Redistribution Matrix Files (RMFs) to generate an off-axis response matrix corresponding to a particular dataset. The off-axis pointing induces projection effects on the instrument's effective area, which is nearly energy independent, reducing it by an average factor of ∼ 0 . 92 over all channels.  \nFor the Delta dataset, we created custom response matrices using the nicerarf and nicerrmf tools, specifying the pulsar position rather than the XTI pointing direction. This produces an ARF that is corrected for the vignetting and time-variable number of active detectors. The product of the ARF and RMF is generated as a single response matrix (RSP) file. The 3C50 dataset uniformly selects the same 50 FPMs for each GTI (Section 2.1.2). A single pair of response files, computed with nicerarf and nicerrmf , are used for this dataset, combining 50 of the 52 FPMs, while excluding FPMs 14 and 34.  \nThe XMM-Newton EPIC MOS1 and MOS2 RMF and ARF products were generated using the rmfgen and arfgen tools in SAS. These resulting response files, along with the source and background spectral files, for the two archival observations were merged into a single set using a weighted average based on the relative exposures of the two observations.  \nTo capture the ∼ ± 10 % calibration uncertainty associated with the effective area (and assuming that this is energy-independent), we define an absolute energyindependent instrument-scaling factor α XTI that operates on the response matrix. In some of the earlier Delta runs, we combined this with the pulsar distance D [kpc] into a single flux-scaling variable β [kpc -2 ] = α XTI D -2 that we are sensitive to. This choice reduces the number of free parameters to be sampled and helps save on computational resources. This parameterization is further detailed in Vinciguerra et al. (2023). For the later runs, we fix the distance altogether considering its highly informative prior and only sample α XTI .  \nTable 1. Default and headline-result X-PSI settings 17 used in this work. The default settings are used unless explicitly stated otherwise. Some of the headline values are smaller than the default values to compensate for computational expense incurred by improved sampling of the most complex model in our analysis.  \nWhen we jointly fit for NICER and XMM-Newton data, we introduce another energy-independent effective area term α MOS that is applied to the responses of the XMM-Newton cameras, MOS1 and MOS2. α XTI and α MOS share a correlated prior PDF as explained in sections 2.4 and 4.2 of Riley et al. (2021). We apply a shared scaling uncertainty factor of 10% and telescopespecific uncertainty factors of 3% leading to a global 10.4% uncertainty in the overall scaling factors.", '3.4. X-PSI Settings': 'Modeling in X-PSI requires specification of certain resolution settings related to the hot region cell mesh discretisation, ray emanation, energy and temporal hyperslice discretisation (for stellar rotation), and multiple imaging. We specify the default values of the X-PSI settings 17 that we use for all runs in Table 1, unless explicitly stated otherwise 18 .', '4.1. Likelihood function and background': 'In this section, we briefly describe the likelihood function formulation for the different background constraints as applied to various models in this letter. The origin and consequences of various assumptions and choices briefly mentioned here are elaborated in Riley et al. (2019), Riley et al. (2021) and Salmi et al. (2022).  \nOur framework involves a background model aimed at capturing any background emission external to the neutron star. These are expected to be registered as homogeneous Poissonian events, which we attempt to quantify by defining a set of statistically independent  \nbackground count-rate (nuisance) variables { E [ b ] } , operating one per channel. This entails the notion that background emission processes, unlike the neutron star signal, will not coherently sum up with the rotation period of the pulsar upon phase-folding and will hence be subsumed by the phase-invariant background model.  \nOur likelihood function is therefore conditional on our regular model parameters and the background countrate parameters. To save on computational expense, we numerically marginalize this likelihood function over these nuisance parameters (see equations 4 and 9 in Riley et al. 2019). We specify channel-separable flat prior PDFs between a given range, i.e., E [ b ] ∼ U ( L , U ), where L and U are lower and upper bounds for the prior support at a given channel. When our analyses involve multiple instruments, we have separate background models (and associated prior supports). Here, we denote the NICER background model as E [ b N ] ∼ U ( L N , U N ) and the XMM-Newton background model as E [ b X ] ∼ U ( L X , U X ).  \nIn this work, for the NICER-only analyses, when we do not have any specific information on the upper background limit, the upper bound of the prior support is left unspecified (see Appendix B.2.3.2 of Riley 2019, for more details). When we use the Delta dataset and the space weather estimate, we specifically define only a lower limit as {L N } := 0 . 9 ×{B N , SW } , where {B N , SW } is the set of NICER background count rates from the space weather estimate.  \nWhen using the 3C50 dataset and background estimate in combination with either of the AGN spectra (Delta or 3C50), we apply them in three different ways. We can choose to specify only a lower instrument background limit 19 as {L N } := max(0 , {B N , 3C50 } -n l { σ B N , 3C50 } ), where {B N , 3C50 } is the set of NICER background count rates from the 3C50 instrument background estimate, { σ B N , 3C50 } is the total standard deviation 20 , and n l is the degree of conservatism acting on the lower standard deviation. We can also specify only a lower limit consisting of the instrument background estimate and an AGN spectrum as {L N } := max(0 , {B N , 3C50 } -n l { σ B N , 3C50 } +0 . 25 ×{B N , AGN } ), where {B N , AGN } is the set of AGN count-rate observed by NICER. We sometimes specify an upper limit on the background consisting of the instrument background estimate and an AGN spectrum as {U N } := {B N , 3C50 } + n u { σ B N , 3C50 } + 2 × {B N , AGN } , where n u is the degree of conservatism acting on the upper standard deviation. Considering the variability of the AGN, the above range  \n20 The total standard deviation here is square root of the quadratic sum of the systematic and statistical errors (See Section 3.4 and Appendix A.2 of Salmi et al. 2022, for more details.)  \nTable 2. Default and headline-result MultiNest settings used in this work. The default settings are used unless explicitly stated otherwise.  \nNote -Reducing the values of sampling efficiency and evidence tolerance, and increasing the number of live points leads to enhanced accuracy in the calculation, albeit at an increased computational expense. The multi-mode variant locks in a select number of live points onto identified modes, the number of which depends on its prior mass, in order to better explore them. This, however, slightly compromises the precision of the final results, assuming the same settings otherwise.  \nof 0.25-2 × the AGN estimate is chosen to encapsulate the observed historic maximum and minimum of the AGN X-ray flux (see Figure 2 and the discussion in Section 2.1.2).  \nWhen involving XMM, the corresponding lower and upper instrument background limits are set as {L X } := max(0 , {B X }-n √ {B X } ) and {U X } := {B X } + n √ {B X } , where {B X } is the set of XMM-Newton background count numbers based on blank-sky observation, and n = 4 is the chosen degree of conservatism (as was done in Riley et al. (2021); Salmi et al. (2022)).  \nThe exact applications of the different lower and upper limits in the different models are detailed in Sections 5 and 6.', '4.2. Posterior computation': 'Once we construct the likelihood function and define our prior PDFs using X-PSI, we employ MultiNest (Feroz et al. 2009), interfaced via PyMultiNest (Buchner et al. 2014a). MultiNest uses a Bayesian inference technique known as nested sampling (Skilling 2004) that targets estimation of the evidence. This requires exploration of the parameter space, which in turn yields posterior distributions as a byproduct. Table 2 summarizes the most relevant settings 21 pertaining to our analyses and their default values that we use for most of our exploratory runs, unless explicitly stated otherwise. Our two CST+PDT headline runs use a much higher resolution of 20,000 live points for thorough exploration of the parameter space.', '5. EXPLORATORY ANALYSES': 'In our analyses, we performed a number of exploratory runs using the ST-U and ST+PDT models and the PPTADR2.5 radio priors described in Section 3.1. We tested how the different dataset choices, spot patterns, background models, and runtime settings affect our inferences. In this section, we document the results of these exploratory runs. Thus, having identified the important factors, and via model comparison and quality checks, we arrive upon CST+PDT as our preferred model for PSR J0437 -4715, the results of which are described in Section 6. We performed production quality runs for the CST+PDT model which constitute our headline result. We summarize all the runs reported in this letter in Table 3.', '5.1.1. Effect of datasets': "The ST-U model was initially run on the Delta dataset. Given that this is the simplest and most computationally inexpensive model in our list, we also tested the model on the 3C50 dataset, once that was made available. Besides bench-marking the effects of switching to a new dataset, we also use this model to test the possible implementations of the associated background constraints.  \nTo assess model consistency between the two datasets, we opt to compare the runs involving the fewest differences between them. The inference runs on the Delta dataset were performed both with and without the space weather estimate. However in the first case, our model also included the X-PSI elsewhere component, which allowed for the rest of the star (outside of the hot regions) to emit blackbody radiation with a singular uniform effective temperature. The runs on the 3C50 dataset were performed for a variety of cases, namely: without any background constraints; with only lower instrument-background constraints; with only lower background constraints involving both instrument and AGN estimates; and finally with both lower and upper background constraints. The details of these limits are specified in Section 4.1. The closest comparable runs between the Delta and 3C50 datasets are therefore the ones involving no background limits whatsoever.  \nBefore comparing the two runs, we note the following differences between them. The Delta run uses all of our default MultiNest settings while all the 3C50 runs for this model use lower resolution settings consisting of half the number of live points, i.e., LP 2000 and a SE of 0.8, to save on computational resources for testing purposes. During some early tests on an initial dataset 22 , we found that the inferred results and model performance were  \nTable 3. Summary of datasets used and background implementations for all the runs listed in Sections 5 and 6.  \nvirtually identical upon increasing the live points from 2000 to 4000. As for the X-PSI resolution settings, the Delta run uses a lower number for max sqrt num cells of 64 but a higher number for num leaves of 100, compared to the default values defined in Section 3.4 used by the 3C50 run. The effect of these changes on the likelihood value for a given parameter vector is expected to be minimal. However, the variation in value is not consistent for different vectors, and is not necessarily always higher or lower for one of the resolution settings. Being an earlier dataset, the Delta run incorporated the β parameter combining distance and effective area priors, while the 3C50 run keeps the distance fixed.  \nThe derived marginalized 1-D 68% credible interval (CI 68 ) 23 for radius using the Delta dataset is 10 . 39 +0 . 23 -0 . 27 km (Figure 17), and using the 3C50 dataset is 10 . 01 +0 . 13 -0 . 16 km (Figure 3), thus overlapping with each other. The posterior shifts towards lower values and gets narrower for the 3C50 run. The shift is likely a consequence of the 3C50 data being cleaner and involving different observing periods. Also considering the differences in resolution settings, and accompanying differences in some of the other inferred parameters, particularly N H , the change in inferred radius is not necessarily surprising. As for the tightening of the intervals for the 3C50 run,  \nwhile it is the opposite of what is typically expected upon a reduction in counts (see, e.g., Lo et al. 2013; Psaltis et al. 2014; Salmi et al. 2022), it could be explained as a consequence of a combination of one or many of the following factors - the 3C50 dataset being cleaner, the distance not being fixed for the Delta run, and/or the decrement in MultiNest resolution for 3C50 simply not having sampled the parameter space sufficiently.  \nFor most model parameters (see Figures 17 and 19 24 ), we similarly find an overlap of the CI 68 . Only certain model parameters and derived quantities such as compactness, neutral hydrogen column density, primary spot colatitude and radius, and secondary spot phase do not. Nevertheless, considering the stringent statistical constraints on these parameters, and the general proximity of their median values (except N H ), we do not find it to be particularly of concern. The inferred N H is significantly higher for the 3C50 run which could be explained by the difference in the shape of the spectra between the two datasets. The lower channels of the 3C50 dataset have a particularly low number of counts compared to the Delta dataset, which can be attributed to stronger interstellar attenuation. Both datasets indicate a surface emission pattern consisting of both spots near the north pole with a large primary spot encompassing the pole and a small secondary spot next to it, of slightly lower colatitude and nearly anti-phased.  \nWe perform a graphical posterior predictive check to inspect the model's ability to fit the data. This is done by taking an ensemble of posterior samples and generating (Poisson) standardized residuals over the full phase-channel intervals. Upon visual inspection of the residuals, we find no distinctly identifiable systematic structures in either case 25 , suggesting that this model is generally capable of explaining either dataset when the background is unconstrained, and yields results that are consistent with each other.", '5.1.2. Effect of NICER background constraints': "In the absence of any background constraints, for both datasets, the inferred background is very high and constitutes the bulk of the total spectral count rate across all channels compared to the contribution of the spots' emission (see, e.g., left panel of Figure 4).  \nLet us first consider only the cases where we impose a lower limit on the background. For the Delta run involving the space weather estimate and elsewhere component, we find no significant differences for the  \n24 We provide the full set of corner plots in Appendix B, which include the median and CI 68 , for all parameters of all models discussed in this letter. Any parameter overlap referred to in text without quoted numbers can be seen in these plots.  \n25 See residual generated by CST+PDT model in Section 6.1 for an example of a clean residual plot.  \nposterior distributions of any of the model parameters when compared to the run without any background constraint, with the CI 68 overlapping well with each other for all parameters. The elsewhere component is inferred as having a broad posterior in very low temperatures (log 10 ( T else [K]) = 5 . 32 ± 0 . 20) which therefore does not have any major contribution in the observed energy range. The inferred background is still high, although this time, the associated standard deviations are slightly broader below channel 40. The residuals are virtually indistinguishable compared to those of the background unconstrained run.  \nFor the 3C50 dataset, we tested for two scenarios when only imposing a lower background limit: the first involving only the 3C50 instrument background estimate; the second involving an estimate of the AGN contribution in combination with the first. For the instrument-only case, we also test for n l ∈ [1 , 4] (see Section 4 for definition of n l ). In all cases, there is good overlap for all parameters, typically at the CI 68 level, except the size and location of the primary spot (see Figure 18). For these parameters, the n l = 1 case slightly differs from the rest, but is still very much in the vicinity, especially accounting for their extremely narrow posteriors (for all n l ).  \nThe deviations of the posteriors do not exhibit any unidirectional trend for any of the parameters when going from n l = 1 to 4. This is likely due to insufficient sampling resolution 26 , and the strong influence of one or a very few channels when using the raw 3C50 background spectrum, as was found in Salmi et al. (2022). They all infer roughly similar background contributions, barring slight variations in the breadth of the associated standard deviations. The inferred backgrounds also follow the same trend as the run without any constraints, being very high.  \nGoing forward we choose n l = 3 (and for n u too, wherever applicable), for which we also tested the effect of smoothing the background. Similar to Salmi et al. (2022), we used a Savitzky-Golay filter with a window length of 17 (number of coefficients) and a fifth order polynomial. We find that the inferred radius of 9 . 89 +0 . 14 -0 . 15 km now matches more closely with that of the background unconstrained case (Figure 3), compared to a radius value of 9 . 47 +0 . 07 -0 . 10 km for the unsmoothed version. The same trend is observed for almost all model parameters. The residuals are nearly identical and the backgrounds inferred are very similar with slightly broader intervals for the smoothed run.  \nIn the second scenario, most of the 3C50 testing was done using an AGN spectrum before it was filtered using the 3C50 technique, and therefore involved a higher count rate. The limits on the AGN were implemented  \nFigure 3. 1-D and 2-D marginalized posteriors of mass, radius and compactness inferred by ST-U for the different implementations of background constraints on the 3C50 dataset. The shaded region in the 1-D marginalized posterior denotes the 68.3% credible interval for the very last model listed in the legend. The contours in the 2-D marginalized posterior for a given model denote the 68.3%, 95.4%, and 99.7% credible intervals, and the posterior for the last model in the legend is shaded. The dash-dotted lines represent priors for these parameters (also for all other corner plots in this Letter). In the legend, 3C50 refers to the dataset used, BKG refers to the background and the two numbers that follow correspond to n l and n u ( n u =X indicates an absence of an upper limit), ∆-AGN refers to the inclusion of the Delta AGN spectra (scaled to 3C50 exposure time) when defining lower and upper background limits. The inferences differ vastly upon inclusion of an upper background limit.  \n<!-- image -->  \nas described in Section 4.1. We raise the lower background limit by adding the lower AGN limit onto the 3C50 instrument background (smoothed, n l = 3; used hereon). This run again infers a somewhat similar radius of 9 . 51 +0 . 11 -0 . 11 km. The slight drop, given the tight constraints for all the runs however, implies that the 68% credible interval of this run does not overlap with the corresponding interval of the lower instrumentbackground-only and the free background cases (Figure 3).  \nSpot region parameters aside, the other parameters are in good agreement between the runs involving no background limits or only a lower limit (see Figure 19).  \nThe posterior distributions of most hot region parameters are nearly completely detached between this run and the background-unconstrained run, while some of them show some degree of overlap with the instrumentbackground-only run. All the parameters can still be said to be very close to each other, again considering the stringent statistical constraints. Again, the residuals and inferred backgrounds are highly similar between them. Overall, the presence of only lower background limits do not seem to affect the results much.  \nThe results change drastically upon the introduction of an upper background limit. The radius drops steeply to 6 . 53 +0 . 05 -0 . 09 km, making the star far more compact as can be seen in Figure 3. We also see a significant increase in the value of N H and a slight decrease in the value of α XTI , compared to the other runs discussed so far, all of which contribute to a lower number of counts. The inferred geometries are also very different.  \nThe inferred background is now bound to be within the specified limits and we find that it is particularly low in the channels ≲ 50, where it is close to the lower limits (Figure 4). Looking at the residuals for this run (Figure 5), we find prominent phase-dependent structures where the model alternatively either over-predicts or under-predicts the data. This indicates a clear model deficiency for ST-U in explaining the data in the presence of an upper background constraint.", '5.1.3. Effect of XMM-Newton constraints': 'We attempted joint fits for NICER and XMM-Newton only on the Delta dataset for the ST-U model. The addition of XMM-Newton imposes lower and upper background constraints on our models, as explained in Section 4.1. We do not add any separate NICER background constraints for this run. We find that this has similar effects to the imposition of an upper NICER background limit. We infer a low radius value of 6 . 83 ± 0 . 07 km resulting in a highly compact star. The effective areas for both telescope instruments drop to unrealistic values, requiring a scaling factor of 0 . 55 +0 . 008 -0 . 007 for NICER and 0 . 38 +0 . 007 -0 . 006 for XMM-Newton. The residuals also show similar phase-dependent residual structures but even more pronounced this time around, once again showing that ST-U is unable to explain the data in the presence of upper background constraints.', '5.2.1. Effect of datasets': "Since ST-U is our simplest model and is clearly unable to explain data in the presence of upper background constraints, we also check for model consistencies between the datasets using ST+PDT , one of our more capable models in this regard (as demonstrated in the coming paragraphs). ST+PDT appeared to be promising based on tests on earlier datasets. We directly tested the Delta dataset inclusive of space weather es-  \nFigure 4. Inferred background and spot contributions by ST-U in the absence of any background constraints (left), and in the presence of both lower and upper background constraints (right) involving the instrument background and the Delta AGN spectrum (scaled to 3C50 exposure time), as applied to the 3C50 dataset. The inferred background constitutes the bulk of the observed data in the absence of an upper background constraint. For more details of the Figure elements and notations, see caption of Figure 5 in Vinciguerra et al. (2024)  \n<!-- image -->  \ntimates on this model. The closest matching run with the 3C50 dataset is one consisting of a lower instrumentbackground-only constraint.  \nDatasets and background limits aside, the two runs also differ in their resolution settings. The Delta run uses higher than default MultiNest resolution settings of LP 8000, SE 0.1, motivated by tests on an initial dataset. The 3C50 runs on the other hand use all of the default settings except for a lower LP of 2000 to save on computational resources for testing purposes. The Delta run also uses slightly higher X-PSI resolution for min sqrt num cells of 16 and num leaves of 100 whereas the 3C50 runs stick to the default values.  \nThe Delta run infers a radius of 11 . 14 +0 . 46 -0 . 40 km (Figure 20) while the corresponding 3C50 run infers a value of 10 . 33 +0 . 41 -0 . 39 km (Figure 6), thus overlapping marginally with each other. This trend of reduction in radius for the 3C50 dataset is therefore somewhat consistent with what was observed for ST-U , but more pronounced. Unlike ST-U , the N H values are very similar for the two runs. The differences in radii can therefore be attributed to a combination of the variations in the other parameters, the resolution differences and the inherent differences between the two datasets.  \nThe spot geometry parameters are generally similar to each other between the two runs (see Figures 20 and 21), but are quite different from the patterns inferred by  \nST-U 27 . This model infers a large ST spot encompassing the north pole and a smaller PDT spot near the equator. Both Delta and 3C50 runs infer overlapping CI 68 for the ST angular size and temperature. The ST colatitude credible intervals overlap at a 95% level. For the Delta run, we find two adjacent but disjoint posterior modes for the phase of the ST spot. The primary mode contains a much larger posterior mass such that the secondary mode only shows up in the 95% credible interval. The secondary mode overlaps well with the inferred values of the 3C50 run.  \nThe PDT spot consists of a very small and high temperature superseding component compared to the larger and low temperature ceding component. For the Delta run, all hot spot parameters (except the ceding component's colatitude) also exhibit two posterior modes, which are close to each other. Similar to the ST phase, one of the modes has a much larger posterior mass than the other. Between the two modes, some of the parameters share the 99% credible interval while it is disjoint for others. The secondary modes for all these parameters overlap significantly with the corresponding inferred values of the 3C50 run. The ceding component's colatitude 68% credible intervals also overlap for the two runs.  \n27 All models attempted consisting of only single temperature regions infer roughly similar geometries to ST-U .  \nFigure 6. Mass, radius and compactness posteriors inferred by the ST+PDT model configuration for the different implementations of background constraints on the 3C50 dataset. The inclusion of an upper background limit tightens the radius posterior and favors a more compact solution.  \n<!-- image -->  \nFigure 5. Top panel: 3C50 data. Middle panel: Posterior expected count numbers generated by the ST-U model when imposing both lower and upper background limits. Bottom panel: Poisson standardized residuals ( c ik and d ik represent the model-generated and data counts resolved over phases i and channels k ). Prominently visible systematic structures indicate poor model performance.  \n<!-- image -->  \nThe secondary modes of all the parameters, wherever present, extend nominally towards lower radii. This could indicate that the dominant primary mode of the Delta run is either not identified by the 3C50 run given its lower resolution settings, or that the cleaner 3C50 dataset does not favor such solutions. We find the latter to be more likely given the similar posterior widths of the two modes and their close proximity to each other, which MultiNest should therefore have identified, especially considering the larger posterior mass of the primary mode recognized for the Delta run.  \nBoth runs infer similarly high background count rates as was seen for ST-U in the absence of any upper constraint. Neither run shows any distinct residual structure suggesting that model output resembles the data well. The results obtained using either dataset can also be considered to be generally in agreement with each other.", '5.2.2. Effect of NICER background constraints': "The ST+PDT model was tested on the 3C50 dataset for three cases: without any background constraints, with  \nonly lower instrument background constraints, and with both lower and upper constraints involving the 3C50 AGN estimate. Note that here we are utilizing the actual 3C50 AGN spectrum instead of scaling the Delta AGN spectrum to 3C50 exposure times as was done during the inference runs using ST-U .  \nThe radius posteriors for all three runs overlap well with each other, as shown in Figure 6. We obtain a radius of 10 . 13 +0 . 33 -0 . 32 km for the free background run, and 9 . 76 ± 0 . 23 km when we impose both limits. The posterior is therefore much tighter in the presence of an upper background limit, with a preference for more compact solutions.  \nThe other non-spot parameters are also consistent between the three runs, particularly for the runs without upper background limits (Figure 21). The run with upper limits infers slightly higher values for α XTI and lower values for N H which could partially account for the lower radii. The main differences for the upper background constrained runs show up in the spot geometries. In the presence of upper background limits, the ST is inferred to be smaller than the PDT and is quite offset from the north pole. The PDT , instead of being near the equator, now encompasses the south pole, thus contributing significantly to the unpulsed emission. The only similarity between this run and ones without upper background  \nlimits is that the PDT ceding component is much larger with respect to the superseding component. The ceding component inferred is larger than that inferred by the other two runs. The superseding component's angular radius contours, although not consistent with the other two runs, can be said to infer similar values given the narrow posteriors. As for the spot temperatures, the upper background constrained run infers a much hotter ST region and a slightly cooler PDT superseding component. The PDT ceding component temperature is similar between all three runs.  \nThe inferred background for the runs without upper constraints are largely similar. With the imposition of the upper limit, the inferred background seems to sit about halfway between the two limits as seen in Figure 7. The residuals have no apparent features present that would indicate poor model performance.", '6. PRODUCTION ANALYSES - CST+PDT': 'Having performed our exploratory analyses, we arrive at CST+PDT as our final model in this work, based on model evidence comparisons (as explained later in this section). For the headline result version of this model, we provide all the inferred posterior and maximum likelihood values of the model parameters in Table 4. Having established the consistency between datasets, we only analyze the 3C50 dataset for this model because it is cleaner, and allows us to place both lower and upper background constraints.  \nFor this model, in addition to default MultiNest resolution runs, we also conducted high resolution runs (simply referred to as high-res runs hereon) using LP 20,000 albeit having reduced certain X-PSI resolution settings ( sqrt num cells = 18, max sqrt num cells = 32, num leaves = 32, num energies = 64). This choice enables us to more thoroughly explore the parameter space while balancing out computational expense 28 . For these high-res runs exclusively, we utilized the latest PPTA-DR3 priors on the mass, inclination, and distance.', '6.1. Effect of NICER background constraints': "The high-res runs were conducted for two cases only, namely, with no background constraint and with both lower and upper background constraints using the 3C50 AGN estimate. We have also tested the application of only a lower instrument background constraint for the default resolution runs.  \nIn Figure 8, we show that for the high-res runs, in the absence of any background constraints, we infer a radius  \n28 While this choice affects the likelihood calculations, the deviations are at a level that is not expected to have a major effect on the overall shape of the posterior modes for such a model based on the tests conducted by Choudhury et al. (2024).  \nof 11 . 01 +0 . 51 -0 . 48 km; and with both lower and upper limits in place, we infer a radius of 11 . 36 + 0 . 95 -0 . 63 km , our headline result for this letter. The two radius intervals are therefore well in agreement with each other, with the latter consisting of a slightly higher median and a longer tail towards higher radii, but otherwise encompassing the former. The narrower posterior for the run without any background information could be attributed to the presence of solutions with much higher likelihood values whose inferred backgrounds are far higher than those allowed by our upper background constraint.  \nThe default resolution run with only a lower limit inferred a radius of 10 . 88 ± 0 . 56 km, thus overlapping with the aforementioned values. The lower median and the overall posterior being shifted towards lower radii is likely due to the use of PPTA-DR2.5 priors for this run, along with some minor effects of using different resolution settings. We deduce this based on the high degree of overlap between this run and the backgroundunconstrained default resolution PPTA-DR2.5 run, for which we infer a radius of 10 . 91 ± 0 . 58 km. If using lower MultiNest resolution settings were the main reason for the lower radius inferred, we would expect stochastic processes to vary the results more between these runs. The differing X-PSI resolution settings are expected to have only a minor influence on the inferred radius.  \nSince CST+PDT encompasses the ST+PDT model, we also compare the inferences between these two models to see how addition of a masking component to the primary spot affects our results. The radius inferred by the highres run using no background constraints does not overlap with the corresponding run of ST+PDT at a CI 68 level. However, given their general proximity, the differences in resolution between these runs, the different radio priors used, and the observed trend of moving towards higher radii for the PPTA-DR3 radio priors, these solutions can be considered to be in agreement with each other. The same cannot be argued, however, when we compare the corresponding cases using both lower and upper background constraints due to the much larger differences between the inferred radii, and the very narrow radius posterior inferred by ST+PDT .  \nThe other non-spot parameters (Figure 9) overlap well between the two high-res runs at the CI 68 level 29 . The run with both background limits in place infers slightly higher values for α XTI and lower values for N H . For ST+PDT , similar trends seen for these two parameters partially explained the lower radii values inferred compared to the corresponding background-unconstrained case. For this model however, we see an increment in  \n29 We abstain from further describing the results of the run with only a lower background limit, due to its use of the older radio priors, default resolution settings, and its overall overlap with the corresponding run with no background constraints. For more details on the results of the default resolution runs, Figure 22  \nFigure 7. Inferred background and spot contributions by ST+PDT in the absence of background constraints (left), and in the presence of both lower and upper 3C50 background limits (right) involving the instrument background and the 3C50 AGN spectrum, as applied to the 3C50 dataset. Similarly to ST-U, this model also infers a very high background contribution in the absence of upper limits.  \n<!-- image -->  \nthe inferred radii values instead, due to the difference in the PDT spot geometry (described in the coming paragraphs). For CST+PDT , the differences in α XTI , N H and radii, compared to the background-unconstrained case, all compensate for the missing unpulsed emission upon imposition of an upper background limit.  \nThe main differences show up in some of the hot region parameters (Figure 10). The CST hot regions are still somewhat similar between the two runs. The angular radius of the emitting component is about the same in both cases. The size of the masking component inferred, however, is different. When using no background constraints, the posterior peaks at very small angular radii for the masking component followed by a gentle drop-off, resulting in a CI 68 of 0 . 06 +0 . 05 -0 . 03 radians. The maximum likelihood solution has a value of 0.037 radians, making the CST spot somewhat degenerate with a ST spot. In the presence of both background limits, a larger hole is preferred, with a CI 68 of 0 . 14 +0 . 10 -0 . 08 radians for the masking component. Both CST s have similar locations where they encompass the north pole (see Figure 11 showing the projection of the spot geometry for the maximum likelihood vector of the high-res run consisting of both background limits). They also have comparable temperatures, even though they do not overlap at the CI 68 level given their stringent statistical constraints.  \nThe PDT hot regions are comparatively more dissimilar between the two high-res runs, particularly in terms of their location. For the background-unconstrained  \nrun, the colatitude posteriors of both the PDT components consist of a main peak (comprising the CI 68 ), that is just south of the equator, and a plateau-like extension towards higher colatitudes. However for the run with both background limits, both PDT components are located half-way between the equator and the south pole, almost at the observer inclination angle, with a tighter posterior mode for their colatitudes. This mode is encompassed within the plateau-like extension of the background-unconstrained run. The PDT regions are also well-separated in phase between the two runs, with a difference of 0.38 cycles between the medians of the superseding component's phases, and consisting of very narrow posteriors in each case.  \nThe similarities lie in the pattern of sizes and temperatures of the PDT components between the two runs. In both cases, we infer a very small and high temperature superseding component, and a comparatively larger and cooler ceding component, similar to the findings of the ST+PDT run without an upper background constraint. For both the PDT components, the run with lower and upper background limits infers slightly larger angular radii such that the CI 68 do not overlap with the background-unconstrained inferred values. The former also infers a slightly hotter superseding component, again with no CI 68 overlap. The two runs only strictly agree on the temperature of the ceding component.  \nSimilar to ST-U and ST+PDT , the background-unconstrained run infers a very high background contribution, far be-  \nFigure 8. Mass, radius and compactness posteriors inferred by the high MultiNest resolution CST+PDT runs, both in the absence of any background constraints, and in the presence of lower and upper background constraints involving the instrument background and 3C50 AGN spectrum. The 2-D marginalized posteriors of the latter model, that constitutes our headline result, are shaded in blue. Both runs infer highly consistent radius posteriors, with the backgroundconstrained run consisting of a distribution tail that extends slightly more towards higher radii.  \n<!-- image -->  \nyond the limits we consider, particularly in the lower energy channels (Figure 12). Upon imposition of the lower and upper constraints, the inferred background is rather high (within the bounds) for some of the lowest channels ( ∼ 30 -40) where it pushes against the upper limit, and then slowly drops and gets rather low for the channel range ∼ 80 -150, pushing against the lower limit. The background-constrained model compensates for the lower overall inferred background (compared to the unconstrained case) by inferring the position of the PDT spot such that it is never eclipsed, which contributes significantly to the unpulsed emission, and lowers the pulse amplitude.  \nThe residuals show no distinctly identifiable systematic structures for either high-res run (Figure 13), indicating the model is generally capable of explaining the data, both in the presence and absence of background constraints.  \nIn order to gauge improvement in model performance with respect to the corresponding runs of ST+PDT , we first compare the default resolution CST+PDT runs  \nagainst ST+PDT , both of which use PPTA-DR2.5 radio priors, to strictly isolate the performance gain by the model alone. Then we compare the default resolution CST+PDT run against the high-res runs to understand the difference in performance in light of the latest PPTADR3 radio priors 30 .  \nWhen imposing no background constraints, the default resolution CST+PDT run has a higher log-evidence 31 than the ST+PDT run by 3.889 units, which is not a very large difference based on the scale provided by Kass & Raftery (1995). However, given the larger model space of CST+PDT which encompasses ST+PDT , we select the former as our headline model. The high-res run only nominally improves the log-evidence value by another 0.257 units. Similar trends are also observed for the runs with both lower and upper background limits in place. The default resolution CST+PDT run is higher by 1.395 units of log-evidence compared to the ST+PDT runs, and the high-res run improves the log-evidence by a further 1.230 units.", '6.2. Effect of XMM-Newton constraints': "We also tested the CST+PDT joint fits using NICER and XMM-Newton. We place both upper and lower background constraints on the NICER data, with XMMNewton providing effectively an additional constraint on the background. This run was performed using the default resolution settings and the PPTA-DR2.5 radio priors. The computational expense of a high resolution joint NICER/XMM-Newton run was estimated to be prohibitive given available resources; taking this into account and the additional cross-calibration uncertainty introduced by the joint fit, we treat this lower resolution run only as a cross-check (particularly since our high resolution NICER-only headline result run already includes background constraints).  \nIn Figure 14, we show the inferred mass, radius, and compactness for this run, compared against the default resolution NICER-only run with both lower and upper background constraints. The NICER+XMM run infers a radius CI 68 of 10 . 04 +0 . 67 -0 . 57 km, which overlaps with the results of the default resolution NICER-only run including both background limits, but does not with the results of the corresponding high-res run. The median also shifts to a lower value with respect to the NICER-only run for this model. Given the observed trend of moving towards higher radii upon using the latest PPTA-DR3 radio priors, we expect a similar shift upwards for this run, especially considering the drop in the inferred mass,  \nFigure 9. Radius, inclination, instrument scaling and N H posteriors inferred by the high MultiNest resolution CST+PDT runs for the two different implementations of background constraints on the 3C50 dataset.  \n<!-- image -->  \nwith the median of the mass lying below the lower 1σ limit of the new mass prior. The application of the latest prior thus ought to reduce the offset with the high-res NICER-only run with both background limits in place.  \nThe compactness posterior exhibits a slightly bimodal structure although such modality does not appear separately in the marginalized 1D mass and radius posteriors. The less compact mode overlaps better with the posterior of the NICER-only run. The correlation between radius and compactness results in this more compact mode dragging the radius posterior towards lower values.  \nThe other non-spot parameters (Figure 23) have very good agreement between the two runs. We obtain a much tighter constraint on the NICER instrument scal-  \ning factor with a CI 68 of 1 . 04 +0 . 02 -0 . 03 , and also a very tight constraint on the XMM-Newton instrument scaling factor with a CI 68 of 0 . 99 +0 . 01 -0 . 01 . The spot geometries are somewhat similar between this run and the NICER-only run with both background limits in place. They more closely resemble each other compared to the runs involving no upper background constraint.  \nOnce again, we infer a CST region encompassing the north pole, although this time the center is slightly more offset from the pole compared the NICER-only runs. The colatitude credible intervals only overlap at the 95% level with the default resolution NICER-only run with both background limits. The angular radius of the emitting components agree at a CI 68 level between these two runs, with the XMM run slightly smaller angular radii.  \nFigure 10. Radius and spot geometry posteriors inferred by the high MultiNest resolution CST+PDT runs for the two different implementations of background constraints on the 3C50 dataset.  \n<!-- image -->  \nFigure 11. Projections of the inferred maximum likelihood (ML) geometry by the high MultiNest resolution CST+PDT run involving both lower and upper background constraints as viewed from the observer inclination (left) and from the north pole (right). Hot spots in the opposite hemispheres are indicated using dimmed colors. Lensing effects are not shown in these projections. The ML geometry consists of a ring encompassing the north pole, and a two-temperature spot in the southern hemisphere.  \n<!-- image -->  \nFigure 12. Inferred background and spot contributions by the high MultiNest CST+PDT runs in the absence of background constraints (left), and in the presence of both lower and upper 3C50 background limits (right) involving the instrument background and the 3C50 AGN spectrum, as applied to the 3C50 dataset. Similar to the exploratory models, we once again infer a very high background contribution in the absence of upper limits.  \n<!-- image -->  \nFigure 13. Top panel: 3C50 data. Middle panel: Posterior expected count numbers generated by the high resolution MultiNest CST+PDT run when imposing both lower and upper background limits. Bottom panel: Poisson standardized residuals showing no apparent systematic structures, indicating that the model is capable of explaining the data well.  \n<!-- image -->  \nThe CI 68 of the angular radius of the masking component strongly overlap with each other, with the XMM run posterior peaking at lower values. The temperatures inferred for the CST region are almost identical.  \nThe colatitude and angular radius of the PDT spot's superseding component are nearly identical between the two runs. The corresponding temperatures agree at a CI 68 level, with the XMM run preferring slightly higher values. The credible intervals of their phases only overlap at the 95% level, but can be considered to be roughly similar accounting for their stringent statistical constraints. As for the PDT ceding component, the colatitude and azimuthal offset credible intervals only overlap at the 95% level, with the XMM run inferring the spot  \nto be a bit closer to the south pole. They both have almost the same sizes, with the XMM run preferring slightly larger angular radii. They also share roughly similar temperatures, overlapping at a CI 68 level, and the posterior of the XMM run moves to higher temperatures.  \nQuite a few of these geometric parameters also show correlations with the compactness, particularly the colatitude of the CST region, the size and temperature of its emitting component, the size of its masking component, and the colatitude and temperatures of both the PDT components. The size of the CST emitting and masking components exhibit especially strong correlations. The various combinations of these parameters and the star's  \nFigure 15. Mass, radius and compactness posteriors inferred by CST+PDT when the models are not informed by radio priors and do not involve any background constraints. In the absence of radio priors, the model prefers much smaller masses.  \n<!-- image -->  \nFigure 14. Mass, radius and compactness posteriors inferred by the CST+PDT for the NICER-only run (default resolution) involving both lower and upper background limits, and for the joint NICER-XMM run. The joint NICER-XMM run exhibits a bimodal structure, not seen in any of the previous NICER-only analyses. The less compact mode agrees more with the background-constrained NICER-only analysis for this model.  \n<!-- image -->  \nradius could result in the formation of such a bimodal structure in compactness.  \nFor the joint NICER and XMM-Newton run, the inferred NICER background in the lower channels lies almost halfway in between the NICER background limits, whereas the NICER-only run was almost pushing up against the upper limit. For the remaining channels, they appear to be consistent with each other. The inferred background for XMM-Newton MOS1 and MOS2 are highly consistent with the input background data. The residuals for this run also do not exhibit any identifiable structures that would hint at a model deficiency.", '6.3. Effect of radio constraints': 'The change from PPTA-DR2.5 to the latest PPTADR3 radio priors only has relatively minor effects on our overall inferences. We have previously shown that the radius slightly increases upon using the latest priors. We also tested the importance of informing our models with these priors in the context of PSR J0437 -4715, by setting up inference runs where we do not use the radio measurements that can only be obtained for binary systems. We attempted these runs both without any background constraints, and with lower and upper limits in  \nplace. In both cases, we set the mass to have a uniform prior between 1 and 3 M ⊙ , and a uniform prior over the cosine of the inclination angle between -1 and 1. We still utilize the distance information although this time we do not fix it, but use the actual prior distribution to set up the β parameter. We slightly lowered the X-PSI resolution ( max sqrt num cells = 64, num energies = 64) to save on computational expense, especially since accurate parameter estimation is not the primary goal of this run.  \nWith the background constraints in place, the run did not converge in a reasonable time scale. The acceptance fraction reported by MultiNest kept dropping very rapidly as the run progressed. This is not too surprising since the background constraints do not actually reduce the prior hypervolume. Instead, they disfavor many solutions that would previously have had reasonable likelihood values, thus complicating the likelihood surface and making it more challenging for the sampler to correctly identify the relevant modes. We suspect that broadening the prior ranges in light of such highly informative data and in the presence of background constraints likely resulted in an extremely jagged likelihood surface consisting of multiple highly peaked modes, that then turned out to be far too computationally expensive to resolve.  \nWhen not using any background constraints, the results are vastly different and do not conform with the radio observations. We obtain a very low mass of 1 . 074 +0 . 098 -0 . 053 M ⊙ with the posterior peaking right at the lower edge of our prior bounds (Figure 15). Given the extent of the inclination prior space, we infer two separate degenerate modes for inclination caused by reflection symmetry in combination with spot colatitudes. One of these modes is about halfway between the south pole and the equator, much like the radio observations. However, given the tight constraints on both the radio measurement and the inferred posterior mode for inclination, we cannot claim to have accurately recovered the radio observation. The radius that we infer is also much lower than what we have previously obtained using CST+PDT with a CI 68 of 9 . 42 +0 . 95 -0 . 78 km.  \nWith no constraints in place, the inferred background is once again relatively high with respect to the limits we consider, although with a much larger spread in possible background contributions. The residuals show no discernible structures, indicating that this model can reproduce the data. This run also outperforms the highres run with no background constraints by 9.228 units in log-evidence. This indicates that while such solutions are far more preferable for explaining the data, our confidence in the radio measurements implies that such physically-viable solutions do not actually represent the true situation. Therefore, for pulse profile modeling of PSR J0437 -4715, the radio measurements play a vital role in our inference runs.', '7. DISCUSSION': "The preferred model that emerges from our analysis, taking into account both upper and lower limits on the background and using NICER data and background estimates, has a mass of M = 1 . 418 ± 0 . 037 M ⊙ (largely informed by the PPTA-DR3 mass prior) and a radius of R = 11 . 36 +0 . 95 -0 . 63 km (68% credible intervals). The associated hot spot geometry has a large single temperature ring encompassing the north pole, and a dual temperature spot (consisting of a small high temperature component overlapping a larger low temperature component) between the equator and the south pole. The observer views the star from the southern hemisphere at an inclination angle of 137 . 506 ± 0 . 019 degrees. This inclination angle does not match the 16 · inferred from modeling the position angle of the linearly polarized part of PSR J0437 -4715's radio emission (Gil & Krawczyk 1997), however see Section 7.4 for a discussion of the limitations of the radio data analysis.", '7.1. Modeling assumptions': "Our modeling is conditional upon the tight mass prior obtained from radio timing. In the analysis of PSR J0740+6620 by Riley et al. (2021), the impact of the mass prior provided by radio timing was assessed; without it, no constraining mass posterior was recov-  \nered (the posterior changed very little from the unconstrained mass prior) - it would not have been concluded independently that it was a high mass pulsar. For PSR J0437 -4715, when we drop the informative mass prior, we find that the recovered posterior does differ substantially from the broad prior but it disfavors the known mass very strongly (at least when not imposing any background constraints). Unfortunately due to computational cost we have not yet been able to determine whether the imposition of background constraints would change this. Interestingly, for PSR J0030+0451, which does not have a mass prior since it is not in a binary, simulations have demonstrated that it is nevertheless possible to recover the injected mass (Vinciguerra et al. 2023). The factors determining reliance on an informative mass prior are clearly something that merits further detailed investigation, since several of the NICER PPM targets lack mass priors.  \nIn our modeling we have made use of available background constraints: when analyzing NICER data alone, we used both 3C50 instrument background estimates and an estimate of the AGN contribution to put lower and upper limits on the background. This does make a difference to our results: without background constraints we would have inferred a slightly lower radius of 11 . 01 +0 . 51 -0 . 48 km. Including XMM-Newton data in joint modeling provides a useful additional check: there the solution is somewhat bimodal, but results are consistent with the NICER-only analysis. The fact that background constraints have an effect shows the importance of developing good instrumental background models and also monitoring variable sources in the field of view.  \nThe AGN is variable on a range of timescales. Halpern & Marshall (1996) observed it for 20 days with EUVE in the 0.1-0.2 keV band and found that it varied unpredictably with as short as a 5 hour doubling time and a total range of a factor of 4 during the full span. We also looked at a 110 ks XMM-Newton observation from 2009 December 15 (ObsID 0603460101) and again saw rapid variability of about 50% on 10 ks timescales. This means that only fully simultaneous observations would be sufficient to accurately measure the AGN contribution to the observed NICER spectrum, but this is unfeasible for the large number of short observations that are a consequence of NICER's visibility constraints on the ISS.  \nOur primary result uses a fully ionized hydrogen atmosphere. As shown in Salmi et al. (2023), the atmosphere choices, such as hydrogen versus helium composition, have the potential to significantly affect the inferred neutron star parameters, including radius. However, not all choices are equally probable or able to explain the data. Hydrogen composition is considered most likely due to the rapid sinking of heavier elements (Alcock & Illarionov 1980; Hameury et al. 1983; Brown et al. 2002; Zavlin & Pavlov 2002). The tight mass and inclination priors for PSR J0437 -4715 may also help to reduce the  \nsensitivity to the atmosphere choices, similar to the case of PSR J0740+6620 (Salmi et al. 2023).  \nFinally, the present analysis neglected the possible contribution of the non-thermal emission to the phaseresolved signal from PSR J0437 -4715. It has been known for some time that this pulsar exhibits nonthermal emission, which had been detected in Chandra data as a power-law tail (Zavlin et al. 2002), then confirmed in the XMM-Newton phase-averaged spectrum (Bogdanov et al. 2006; Bogdanov 2013). One possible explanation for the origin of this non-thermal emission is inverse Compton scattering of the thermal X-ray photons by relativistic electrons, located above the polar caps (Bogdanov et al. 2006; Bogdanov 2013). When fitted with a power-law, this non-thermal emission dominates over the thermal emission above ∼ 3 keV. However, the XMM-Newton data only could not precisely constrain the slope of the power-law because it was largely correlated with the choices of model for the emission at lower energies (below 3 keV). Observations with NuSTAR in the 3-20 keV energy range provided a more robust estimate of power-law (Γ ∼ 1 . 60 ± 0 . 25, 90% confidence interval) when using the NuSTAR alone, i.e., free of correlation with the low-energy emission (Guillot et al. 2016).  \nThe NuSTAR observation also revealed that the hard X-ray, non-thermal, emission was also pulsed at the frequency of the pulsar (a detection at 3 . 7 σ 32 Guillot et al. 2016). While this pulsed hard X-ray emission is present above 3 keV, i.e., outside of the energy range of the NICER data analysed in the present letter, its contribution could extend to lower energies and account for a small fraction (of the order of a few percent) of the phase-resolved signal from the pulsar. However, since the mechanism responsible for the non-thermal emission remains unknown, it is difficult to assert the energy range where one might expect the pulsed nonthermal emission to be present. For example, inverse Compton scattering of thermal X-ray photons by a hot plasma of electrons above the polar caps could result in pulsed hard X-ray emission, but it would mostly be above ∼ 3 keV (see Figure 2 in Bogdanov et al. 2006, for Comptonized blackbody emission). In the absence of a firm identification of the origin and energy extent of the pulsed non-thermal emission, we have for now neglected any potential contribution to the phase-resolved signal.", '7.2. Radius measurement in context': "In this section we compare our results for the radius of a ∼ 1 . 4 M ⊙ neutron star with those obtained from  \nother studies or techniques that do not involve the specification of an EoS model.  \nThere are several previous studies that have attempted to place constraints on the radius of PSR J0437-4715. Some early work attempting to fit radiative profiles and constrain polar cap properties and the mass to radius ratio (before a precision mass measurement was available) was carried out by Pavlov & Zavlin (1997); Zavlin & Pavlov (1998). Once a mass measurement was available, Bogdanov (2013) modeled the pulse profile of the XMM-EPIC dataset, assuming two single temperature circular hot spots and concluding that R > 11 . 1 km (at the 3 σ confidence level) for the then current radio timing mass measurement of 1.76 M ⊙ . Using the current mass measurement would have reduced the lower limit on the radius. Our lower limit ( R > 10 . 5 km at the 68% level) is smaller than the values from Bogdanov (2013), however given the improvement in data quality and modeling, we do not find this difference surprising.  \nMore recently Gonz'alez-Caniulef et al. (2019) carried out joint modeling of phase-averaged far ultraviolet (FUV) and soft X-ray data for PSR J0437 -4715, using data from the Hubble Space Telescope (HST) and Rontgensatellit (ROSAT, in the 0.1-0.3 keV range), respectively. They inferred a radius of 13 . 6 +0 . 9 -0 . 8 km, higher than the value inferred in this current study. That analysis modeled the ∼ 10 5 K emission from the rest of the surface (outside the hot spots) using atmosphere models compatible with that range of temperatures, for hydrogen, helium and iron compositions. Using fixed values for the pulsar mass and distance and priors on the amount of extinction E ( B -V ) for the FUV data, the authors favored a hydrogen composition of the atmosphere, although a helium composition was not fully excluded. In the final results of Gonz'alez-Caniulef et al. (2019), the radius was obtained with a model that included the phase-average contributions of the hot spots to the spectra, assuming a blackbody emission for these and using parameter values from the work of Guillot et al. (2016). Although the spots have minor contributions to the total phase-averaged spectrum below 0.3 keV, this approach ensured that this contribution would not bias the spectral fits of the rest of the surface at ∼ 10 5 K.  \nHowever, an update to the HST calibration revealed that the absolute flux of the HST Advanced Camera for Surveys (ACS) Solar Blind Channel (SBC) was incorrect. Specifically, the zero-point flux is actually lower than previously thought by ∼ 30%, thus overestimating the measured FUV fluxes 33 . These data, comprising the bulk of the FUV data points in Gonz'alezCaniulef et al. (2019) likely biased the radius measure-  \n33 See https://archive.stsci.edu/contents/newsletters/ october-2019/sbc-absolute-flux-calibration-revised-by-30  \nment, possibly explaining the difference in radius with the one in the present work. In an effort to mitigate this bias, a re-analysis of the FUV + X-ray spectrum of PSR J0437 -4715 was performed (Stammler et al., in prep) excluding the ACS-SBC data, and only using data from the HST Space Telescope Imaging Spectrograph FUV Multi-Anode MicroChannel Array for which the calibration has not been revised. In addition, this new analysis used an updated, much tighter prior than previously on the extinction E ( B -V ), obtained from the GAIA datasets (Vergely et al. 2022). The full priors 34 on the mass and on the source distance were also used in that recent analysis, instead of fixed values. Finally, the phase-averaged contribution of the hot spots was modeled realistically with atmosphere models, instead of assuming blackbody emission. Overall, the new, updated, radius of PSR J0437 -4715 measured from the FUV and soft X-ray data, is 12 . 3 ± 0 . 9 km (Stammler et al., in prep.), which is compatible with the radius presented in this work.  \nWe can also compare the mass and radius inferred for PSR J0437 -4715 to the masses and radii inferred for other NICER sources. For the heavy pulsar PSR J0740+6620, the current best X-PSI result 35 is M = 2 . 073 +0 . 069 -0 . 069 M ⊙ and R = 12 . 49 +1 . 28 -0 . 88 km (Salmi et al. 2024). This means that ∆ R = R 2 . 0 -R 1 . 4 = 1 . 13 +1 . 59 -1 . 08 km. The implications of this change in radius over this mass range are discussed in more detail in Section 7.3.  \nFor PSR J0030+0451 the picture is more complex (Vinciguerra et al. 2024): the current best NICERonly (with no background constraints) results are M = 1 . 37 ± 0 . 17 M ⊙ , R = 13 . 11 ± 1 . 30 km for the ST+PST model, consistent with the earlier results of Riley et al. (2019) 36 . However when taking into account XMMNewton data to provide constraints on background, two different solutions are preferred, with masses and radii of [ M = 1 . 40 +0 . 13 -0 . 12 M ⊙ , R = 11 . 71 +0 . 88 -0 . 83 km] and [ M = 1 . 70 +0 . 18 -0 . 19 M ⊙ , R = 14 . 44 +0 . 88 -1 . 05 km], depending on the assumed surface pattern model (Vinciguerra et al. 2024). The mass and radius for the first mode identified in the joint NICER-XMM analysis are extremely consistent with our findings for PSR J0437 -4715; the NICERonly solution (which also has a similar mass but a higher radius) less so. The higher mass, higher radius NICERXMM mode of PSR J0030+0451 (although favored according to the evidence calculation) is also harder to reconcile with our new findings.  \n34 In that analysis, the priors used were: E ( B -V ) = 0 . 005 ± 0 . 003, M = 1 . 44 ± 0 . 07 M ⊙ and d = 157 ± 0 . 2 pc  \n35 See Dittmann et al. (2024) for analysis using an independent PPM pipeline for this source.  \n36 See also Miller et al. (2019) for analysis using an independent PPM pipeline for this source.  \nFigure 16. Inferred mass and radius posteriors for the different capable models explored in this paper, compared to mass-radius curves from EoS models. The shaded contours are the 68% and 95% credible regions for the highres CST+PDT runs, and the remaining contours represent the 68% credible regions for the default resolution CST+PDT and ST+PDT runs. The colored dash-dot lines are example massradius curves generated using one of the families of EoS models explored in Rutherford et al. (2024). It employs the Baym-Pethick-Sutherland crust EoS (Baym et al. 1971) at densities less than 0.5 times the nuclear saturation density n s ; the N 3 LO χ EFT band from Keller et al. (2023) at densities 0.5-1.5 n s ; and a parametrized piecewise polytropic (PP) model (Hebeler et al. 2013) at higher densities. The parameters for the mass-radius relations shown here have been chosen to span the range of stiffness allowed by this model that is also compatible with the existence of ∼ 2 . 0 M ⊙ neutron stars.  \n<!-- image -->  \nNote that while rotation does increase the equatorial radius of a neutron star (and the resulting oblateness is incorporated into our ray-tracing models), the change for a 1.4 M ⊙ neutron star as one goes from zero rotation to 174 Hz is at most ∼ 0 . 2 km (for the stiffest EoS with the largest radii), see Figure 5 of Raaijmakers et al. (2019). For lower radii, such as we infer for PSR J0437 -4715, the difference is much smaller than our inferred credible intervals.", '7.3. Equation of State implications': "A full EoS analysis of the implications of our mass and radius for PSR J0437 -4715 is presented in a companion letter by Rutherford et al. (2024). In this Section we therefore highlight only a few key points.  \nThe radius that we infer for PSR J0437 -4715 favors softer EoS models in the vicinity of the nuclear saturation density n s . Figure 16 shows the comparison to several example mass-radius curves from the EoS models explored in Rutherford et al. (2024).  \nThe radius inferred for PSR J0437 -4715 is also more consistent with the results that have emerged to date from gravitational wave measurements of tidal deformabilities. Using EoS insensitive relations, Abbott et al. (2018) reported masses and radii for the two neutron stars involved in the binary merger GW170817 of [ M 1 = 1 . 15 -1 . 36 M ⊙ , R 1 = 10 . 8 +2 . 0 -1 . 7 km] and [ M 2 = 1 . 36 -1 . 62 M ⊙ , R 2 = 10 . 7 +2 . 1 -1 . 5 km] (90% credible intervals). When using an EoS model their inferred masses and radii were [ M 1 = 1 . 18 -1 . 36 M ⊙ , R 1 = 11 . 9 +1 . 4 -1 . 4 km] and [ M 2 = 1 . 36 -1 . 58 M ⊙ , R 2 = 11 . 9 +1 . 4 -1 . 4 km]. Meanwhile Raaijmakers et al. (2021), using different EoS models, inferred R 1 . 4 = 11 . 51 +1 . 51 -1 . 47 km (PP model, 95% credible intervals) and R 1 . 4 = 11 . 18 +1 . 33 -1 . 51 km (CS model) when using only the tidal deformabilities from GW170817 and GW190425.  \nOur results appear less consistent with the neutron star radii that have been inferred from the results of the parity-violating electron scattering experiment PREX2 (Adhikari et al. 2021). PREX measures the neutron skin thickness of 208 Pb, which constrains the density dependence of the symmetry energy, with the results of Adhikari et al. (2021) implying a very stiff symmetry energy. Given the expected correlation between symmetry energy and neutron star radius, Reed et al. (2021) have argued that this would imply R 1 . 4 ≳ 13 . 25 km. However the PREX-2 uncertainties are broad (Reed et al. 2021) and there are also apparent tensions with more recent measurements of the neutron skin thickness of 48 Ca taken by the CREX experiment (Adhikari et al. 2022) and the Large Hadron Collider (LHC) inference of the neutron skin thickness of 208 Pb (Giacalone et al. 2023). In fact, if one treats the results from the PREX and CREX experiments with equal weight, one finds that only nuclear interactions with symmetry energy parameters J = 32 . 2 ± 1 . 7 MeV and L = 52 . 7 ± 13 . 2 MeV (Lattimer 2023) can accommodate both experiments simultaneously ( J and L being the symmetry energy and its slope at nuclear saturation density). In addition, the LHC measurement implies L ≃ 70 ± 40 MeV (Sun et al. 2024). In comparison, nuclear mass fitting with a wide variety of both non-relativistic and relativistic nuclear interactions yields J = 32 . 5 ± 2 . 1 MeV and L = 57 . 7 ± 20 . 7 MeV. Moreover, combining Chiral Effective Field Theory ( χ EFT) neutron matter calculations with estimates of the nuclear saturation window ( n s = 0 . 155 ± 0 . 010 fm -3 and symmetric matter binding energy at saturation B = 16 . 1 ± 0 . 3 MeV) gives J = 32 . 0 ± 1 . 1 MeV and L = 51 . 9 ± 7 . 9 MeV. These estimates are in remarkable agreement with one another. Neutron star models constructed with those  \ninteractions that best fit both PREX and CREX data imply R 1 . 4 = 11 . 6 ± 1 . 0 km (Lattimer 2023). Thus, nuclear experiment and theory are very consistent with our results for the radius of PSR J0437 -4715.  \nThe radius measurement of PSR J0437 -4715, combined with that for PSR J0740+6620 (Salmi et al. 2024), implies that the radius difference between 2.0M ⊙ and 1.4M ⊙ stars (a quantity whose significance was first discussed by Drischler et al. 2021), ∆ R = R 2 . 0 -R 1 . 4 ≃ 1 . 13 +1 . 59 -1 . 08 km. The central value for ∆ R is extremely large (although the uncertainties are also large). For comparison, Sun et al. (2024)'s compilation of 313 nonrelativistic and relativistic equations of state that satisfy M max ≥ 2 . 0M ⊙ , where M max is the maximum possible neutron star mass, have one EoS model with ∆ R ≥ 1 . 0 km, and only 4% have ∆ R ≥ 0 . 5 km. Drischler et al. (2021) states that if ∆ R > 0 km, standard extrapolations would predict some unusual stiffening (including, perhaps, a quarkyonic or crossover transition to quark matter) should occur below 2 . 6 -2 . 8 n s . Furthermore, if ∆ R turns out to be greater than 0.5 km, then one should expect that this stiffening would not only occur for densities ≤ 2 n s , but also would be accompanied by a very high value for M max , one that could be in conflict with expectations from the neutron star merger GW170817 and in particular modeling of its electromagnetic counterpart ( M max ≲ 2 . 3M ⊙ , Margalit & Metzger 2017; Rezzolla et al. 2018; Shibata et al. 2019; Lattimer 2021).", '7.4. Surface heating configuration': "The favored hot spot configuration identified by our analysis, shown in Figure 11, is not consistent with a centered dipole magnetic field but could potentially be explained by an offset dipole magnetic field or a field with quadrupole components. Via magnetospheric modeling, Lockhart et al. (2019) (building on the work by Gralla et al. 2017) infer similar surface emission patterns for PSR J0437 -4715, consisting of a quadrudipolar magnetic field. While there are some dissimilarities to our inferred geometry, these likely stem from the more restrictive fitting procedure, with more fixed model neutron star parameters, used in their analysis.  \nThese findings are also consistent with the picture from radio observations: PSR J0437 -4715 has complex radio pulse morphology that is challenging to model assuming a dipolar magnetic field and plausible radio emission microphysics (Bhat et al. 2014). The swing of position angle (PA) of linearly polarized radio emission does not conform to the prescriptions of a widely used phenomenological rotating vector model, which again relies on the dipole field assumption (Manchester & Johnston 1995; Navarro et al. 1997). Some complexities of the PA behavior may arise from mixing between two orthogonally polarized modes. However a detailed study of PSR J0437 -4715 's single-pulse emission, which resolved the modes, showed that mode mixing can not  \nexplain all peculiarities of PA swing (Osglyph[suppress]lowski et al. 2014). It has been supposed therefore that radio emission from PSR J0437 -4715 originates either close to surface where higher-order multipoles are present (Osglyph[suppress]lowski et al. 2014), or closer to the light cylinder where largescale magnetic field is modified by magnetospheric currents (Gil & Krawczyk 1997).  \nThe fact that the field configuration appears to be more complex than a simple centered dipole is consistent with findings for the other two NICER MSPs analyzed to date. Initial analysis of PSR J0030+0451 NICER data indicated a highly non-dipolar field structure and complex emitting region shapes (Riley et al. 2019; Bilous et al. 2019; Miller et al. 2019). More recent joint analysis of NICER and XMM data suggests that a more antipodal geometry could be possible, but with the hot regions having a more complex dual temperature distribution (Vinciguerra et al. 2024). The heavy pulsar PSR J0740+6620 also has a definitively non-antipodal hot region geometry: Salmi et al. (2022, 2024) found that the offset angle from the antipode is above 25 · with 84% probability.  \nNeutron stars are initially formed with complex magnetic field structures as a result of supernova explosions (Ardeljan et al. 2005; Obergaulinger & Aloy 2017). These fields are further modified over time by various diffusive processes (Reisenegger 2009; Vigan'o et al. 2013; Mitchell et al. 2015; Gourgouliatos & Hollerbach 2018; Igoshev et al. 2023). It remains uncertain whether a complex multipolar field can persist for long periods in isolated neutron stars, but the field can also be altered during the accretion-induced spin-up phase that produces rotation-powered millisecond pulsars (Romani 1990; Melatos & Phinney 2001; Payne & Melatos 2004). Although quantitative modeling of hot spots is still pending, non-antipodal hot spots are expected as a natural outcome of accretion-induced magnetic burial (Suvorov & Melatos 2020). Accurate modeling of hot spot shapes and locations must also incorporate the coupled magneto-thermal evolution of the crust: in isolated neutron stars, hot spots can be displaced from their antipodal positions by past heat deposition events, even in a dipole-dominated magnetic field, provided that the crustal transport properties are properly considered (De Grandis et al. 2020).", '8. CONCLUSIONS': "In this letter, we have presented NICER's first mass, radius and geometry inferences for the bright ∼ 1 . 4 M ⊙ pulsar PSR J0437 -4715, using pulse profile modeling of NICER XTI X-ray spectral-timing data, with priors on mass, distance and inclination from PPTA radio pulsar timing data. For our preferred model, we find an inferred mass M = 1 . 418 ± 0 . 037 M ⊙ (dominated by the radio prior) and radius of R = 11 . 36 +0 . 95 -0 . 63 km. This radius favors softer dense matter equations of state and is highly consistent with constraints derived from gravita-  \nwave measurements of neutron star binary mergers.  \nThe information from the radio prior is critical - our (limited) tests indicate that without it we would not have recovered the correct mass. Our results are also dependent on having good estimates for the astrophysical background, including a contribution from a bright AGN in the NICER field of view. Better monitoring of the AGN may enable tighter background constraints in the future, and NICER continues to build up a larger dataset on PSR J0437 -4715 that will also enable tighter constraints on radius.  \nPPM also allows us to constrain the geometry of the X-ray emitting hotspots on the stellar surface, which are assumed to be the heated magnetic poles. The configuration that we find is not consistent with a centered dipole magnetic field, but could be consistent with either an offset dipole or the presence of quadrupole components. This is consistent with the picture from radio and γ -ray modeling of this pulsar, which has long pointed to a complex magnetic field structure.  \nThis work was supported in part by NASA through the NICER mission and the Astrophysics Explorers Program. D.C., T.S., S.V., T.E.R., Y.K., A.L.W. and B.D. acknowledge support from ERC Consolidator Grant No. 865768 AEONS (PI: Watts). S.B. acknowledges funding from NASA grants 80NSSC20K0275 and 80NSSC22K0728. N.R. was supported by NASA Grant No. 80NSSC22K0092. S.G. acknowledges the support of the CNES. S.M. acknowledges support from NSERC Discovery Grant RGPIN-2019-06077. W.C.G.H. acknowledges support through grant 80NSSC23K0078 from NASA. Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project numbers CE170100004 and CE230100016. The use of the national computer facilities in this research was subsidized by NWO Domain Science. Part of the work was carried out on the HELIOS cluster including dedicated nodes funded via the aforementioned ERC CoG. Astrophysics research at the Naval Research Laboratory is supported by the NASA Astrophysics Explorer Program. We acknowledge extensive use of NASA's Astrophysics Data System (ADS) Bibliographic Services and the ArXiv.  \nSoftware: X-PSI (Riley et al. 2023), GNU Scientific Library ( GSL ; Gough 2009), HEASoft (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014), MPI for Python (Dalc'ın et al. 2008), MultiNest (Ferozetal.2009), PyMultiNest (Buchneretal.2014b), nestcheck (Higson2018), GetDist (Lewis2019), Jupyter (Perez&Granger2007;Kluyveretal.2016), astropy (AstropyCollaborationetal.2013,2018,2022), scipy (Virtanenetal.2020;Gommersetal.2024), matplotlib (Hunter 2007), numpy (Harris et al. 2020), python (Van Rossum &  \nDrake 2009) and Cython (Behnel et al. 2011). Software citation information aggregated using The Software  \nCitation Station (Wagg & Broekgaarden 2024; Wagg & Broekgaarden 2024).", 'REFERENCES': "Shibata, M., Zhou, E., Kiuchi, K., & Fujibayashi, S. 2019, PhRvD, 100, 023015  \nShklovskii, I. S. 1970, Soviet Ast., 13, 562  \nSkilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395-405  \n- Sun, B., Bhattiprolu, S., & Lattimer, J. M. 2024, PhRvC, 109, 055801  \nSuvorov, A. G., & Melatos, A. 2020, MNRAS, 499, 3243 Tak'atsy, J., Kov'acs, P., Wolf, G., & Schaffner-Bielich, J. 2023, PhRvD, 108, 043002  \nTurner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27  \nVan Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)  \nVergely, J. L., Lallement, R., & Cox, N. L. J. 2022, A&A, 664, A174  \n- Vigan'o, D., Rea, N., Pons, J. A., et al. 2013, MNRAS, 434, 123", 'APPENDIX A. CST+PDT RESULTS TABLE': 'Table 4 . Summary table for the headline high-res CST+PDT run described in Section 6', 'B. SUPPLEMENTARY CORNER PLOTS OF POSTERIORS INFERRED FOR DIFFERENT MODEL PARAMETERS BY EXPLORATORY RUNS': 'Figure 17. 1-D and 2-D marginalized posteriors inferred by ST-U when applied to the Delta (∆) dataset, using either no background, Space Weather (SW) + elsewhere component, or informed by XMM-Newton (data and background).  \n<!-- image -->  \nFigure 18. Posteriors obtained by ST-U as applied to the 3C50 dataset, testing the implementation of different n l values when applying a lower 3C50 instrument-only background constraint.  \n<!-- image -->  \nFigure 19. Posteriors obtained by ST-U as applied to the 3C50 dataset, testing the different implementations of background constraints.  \n<!-- image -->  \nFigure 20. Posteriors obtained by ST+PDT as applied to the Delta dataset with lower background constraints imposed by the space weather estimate.  \n<!-- image -->  \nFigure 21. Posteriors obtained by ST+PDT as applied to the the 3C50 dataset, testing the different implementations of the 3C50 background constraints.  \n<!-- image -->  \nFigure 22. Posteriors obtained by CST+PDT (default resolution run) as applied to the the 3C50 dataset, testing the different implementations of the 3C50 background constraints.  \n<!-- image -->  \nFigure 23. Posteriors obtained by CST+PDT for the NICER-only (default resolution run) involving both lower and upper background limits, and for the joint NICER-XMM run.  \n<!-- image -->'}