# A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation

Thomas E. Riley, Anna L. Watts, Slavko Bogdanov, Paul S. Ray, Renee M. Ludlam, Sebastien Guillot, Zaven Arzoumanian, Charles L. Baker, Anna V. Bilous, Deepto Chakrabarty, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Thomas E. Riley, Anna L. Watts, Slavko Bogdanov, Paul S. Ray, Renee M. Ludlam, et al. A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation. The Astrophysical Journal Letters, 2019, 887 (1), pp.L21. 10.3847/2041-8213/ab481c hal-02440907

HAL Id: hal-02440907

https://hal.science/hal-02440907v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

CrossMark

# A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation

T. E. Riley<sup>1</sup>, A. L. Watts<sup>1</sup>, S. Bogdanov<sup>2</sup>, P. S. Ray<sup>3</sup>, R. M. Ludlam<sup>4,16</sup>, S. Guillot<sup>5,6</sup>, Z. Arzoumanian<sup>7</sup>, C. L. Baker<sup>8</sup>, A. V. Bilous<sup>1</sup>, D. Chakrabarty<sup>9</sup>, K. C. Gendreau<sup>7</sup>, A. K. Harding<sup>10</sup>, W. C. G. Ho<sup>11,12</sup>, J. M. Lattimer<sup>13</sup>, S. M. Morsink<sup>14</sup>, and T. E. Strohmayer<sup>15</sup>

*Received 2019 July 12; revised 2019 September 24; accepted 2019 September 25; published 2019 December 12*

## Abstract

We report on Bayesian parameter estimation of the mass and equatorial radius of the millisecond pulsar PSR J0030+0451, conditional on pulse-profile modeling of Neutron Star Interior Composition Explorer X-ray spectral-timing event data. We perform relativistic ray-tracing of thermal emission from hot regions of the pulsar's surface. We assume two distinct hot regions based on two clear pulsed components in the phase-folded pulse-profile data; we explore a number of forms (morphologies and topologies) for each hot region, inferring their parameters in addition to the stellar mass and radius. For the family of models considered, the evidence (prior predictive probability of the data) strongly favors a model that permits both hot regions to be located in the same rotational hemisphere. Models wherein both hot regions are assumed to be simply connected circular single-temperature spots, in particular those where the spots are assumed to be reflection-symmetric with respect to the stellar origin, are strongly disfavored. For the inferred configuration, one hot region subtends an angular extent of only a few degrees (in spherical coordinates with origin at the stellar center) and we are insensitive to other structural details; the second hot region is far more azimuthally extended in the form of a narrow arc, thus requiring a larger number of parameters to describe. The inferred mass *M* and equatorial radius *R*<sub>eq</sub> are, respectively, 1.34<sup>+0.15</sup><sub>-0.16</sub> *M*<sub>☉</sub> and 12.71<sup>+1.14</sup><sub>-1.19</sub> km, while the compactness *GM*/*R*<sub>eq</sub>c<sup>2</sup> = 0.156<sup>+0.008</sup><sub>-0.010</sub> is more tightly constrained; the credible interval bounds reported here are approximately the 16% and 84% quantiles in marginal posterior mass.

*Unified Astronomy Thesaurus concepts: Matter density (1014); Rotation powered pulsars (1408); Millisecond pulsars (1062); Pulsars (1306); X-ray stars (1823); Neutron stars (1108); Neutron star cores (1107); Nuclear astrophysics (1129)*

*Supporting material: figure sets*

## 1. Introduction

Neutron star (NS) cores are thought to harbor nucleonic matter under extreme conditions: high in density, neutron rich, and potentially strange. Stable states of strange matter may either be bound in the form of hyperons, or deconfined as a mixture of up, down, and strange quarks. The density in the stellar core may reach up to several times the nuclear saturation density $\rho_{\text{sat}} = 2.8 \times 10^{14} \text{ g cm}^{-3}$, a range for which we cannot yet calculate the state of nuclear matter from first principles.<sup>17</sup> Instead, theorists develop phenomenological models of particle

interactions and phase transitions, which must be tested by experiment and observation. Heavy ion collision experiments explore the high-temperature and lower-density parts of the nuclear matter phase diagram; but NSs are unique laboratories for the study of strong and weak force physics in cold, dense matter (for recent reviews see Lattimer & Prakash 2016; Oertel et al. 2017; Baym et al. 2018).

The particle interactions on a microphysical scale emerge macroscopically as an equation of state (EOS)—in the context of cold dense matter, a relationship between pressure and (energy) density. The EOS forms part of the relativistic stellar structure equations that enable us, given a central density and a spin rate, to compute model NSs (e.g., Hartle 1967; Hartle & Thorne 1968). An EOS function thus maps to a sequence of stable global spacetime solutions, each controlled in the exterior domain by parameters such as the mass and equatorial radius (at low orders in a small dimensionless spin parameter;

<footnote n="1"><sup>1</sup>Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090GE Amsterdam, The Netherlands; T.E.Riley@uva.nl, A.L.Watts@uva.nl</footnote>

<footnote n="2"><sup>2</sup>Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="3"><sup>3</sup>Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="4"><sup>4</sup>Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="5"><sup>5</sup>IRAP, CNRS, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="6"><sup>6</sup>Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France</footnote>

<footnote n="7"><sup>7</sup>X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA</footnote>

<footnote n="8"><sup>8</sup>Applied Engineering and Technology Directorate, NASA Goddard Space Flight Center, Code 592, Greenbelt, MD 20771, USA</footnote>

<footnote n="9"><sup>9</sup>MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="10"><sup>10</sup>Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="11"><sup>11</sup>Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA</footnote>

<footnote n="12"><sup>12</sup>Mathematical Sciences, Physics and Astronomy, and STAG Research Centre, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="13"><sup>13</sup>Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA</footnote>

<footnote n="14"><sup>14</sup>Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB T6G 2E1, Canada</footnote>

<footnote n="15"><sup>15</sup>Astrophysics Science Division and Joint Space-Science Institute, NASA's Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="16"><sup>16</sup>Einstein Fellow.</footnote>

<footnote n="17"><sup>17</sup>Although note that calculations at sub-saturation densities (see, e.g., Hebeler et al. 2010), and perturbative quantum chromodynamics (QCD) modeling (valid for densities that are several times higher than the maximum expected in NS cores; see Kurkela et al. 2014) can impose some constraints on the low- and high-density limits of NS core parameter space.</footnote>
Hartle 1967). In this work we constrain the total mass (sometimes referred to as the gravitational mass) and the equatorial radius of the star, respectively defined as the mass and coordinate radius in the Schwarzschild metric. If we can statistically estimate the masses and radii of a set of stars whose central densities span some sufficiently broad range, we can in principle map out the EOS and hence make inferential statements about the microphysics (Özel & Psaltis 2009; Steiner et al. 2010, 2013; Nättilä et al. 2016; Özel et al. 2016; Raithel et al. 2017; Raaijmakers et al. 2018; Riley et al. 2018; Greif et al. 2019).

The strongest statistical constraints<sup>18</sup> on NS masses are derived by timing radio pulsars in (compact) binaries, and rely on our well-established understanding of relativistic orbital dynamics. Every EOS function (corresponding to a parameter vector for a parameterized model) permits stable spacetime solutions with a maximum gravitational mass—associated with a specific central density—beyond which no stable solutions can exist. High-mass NSs with tight constraints can therefore effectively exclude<sup>19</sup> a subset of EOS function space, barring strong contention with future analyses of independently acquired data. The most informative pulsars in this regard are PSR J0348+0432, with mass 2.01 ± 0.04 M<sub>☉</sub> (Antoniadis et al. 2013, where the mass is derived by combining pulsar timing and models of the white dwarf companion), and PSR J1614−2230 with mass 1.908 ± 0.016 M<sub>☉</sub> (Arzoumanian et al. 2018, where the mass comes from Shapiro delay estimation).<sup>20</sup> More recently, Cromartie et al. (2019) have reported a higher—but at present more uncertain—mass of 2.14<sup>+0.10</sup><sub>-0.09</sub> M<sub>☉</sub> for PSR J0740+6620. We note, however, that the Cromartie et al. (2019) measurement is not subject to the systematic uncertainty that should be added to the formal uncertainty on the mass of PSR J0348+0432 due to the latter’s dependence on theoretical models of white dwarf evolution.

The radio pulsar timing of compact binaries has yet to deliver a radius constraint, although this is feasible and indeed anticipated, via moment of inertia estimation (Kramer & Wex 2009). There are, however, constraints on radius via X-ray spectral modeling of transiently accreting and bursting NSs (see, e.g., Steiner et al. 2013; Özel et al. 2016; Nättilä et al. 2017; Shaw et al. 2018; Baillot d'Etivaux et al. 2019); we refer the reader to Miller (2013) and Özel & Freire (2016) for detailed reviews that include an explanation of these X-ray modeling techniques and associated uncertainties. NS mass and tidal deformability estimates are now also being reported based on the first binary NS merger gravitational wave event, GW170817 (Abbott et al. 2017, 2018; De et al. 2018). These can be translated, usually by means of universal relations or EOS model assumptions, into constraints on mass and radius (see, e.g., Abbott et al. 2018; Annala et al. 2018; De et al. 2018; Most et al. 2018; Tews et al. 2018b). Generally, assuming that both NSs have the same EOS, it is found that their radii are nearly equal and (for the 68% credible interval) have the common value 11.9 ± 1.1 km.

NASA's *Neutron Star Interior Composition Explorer* (NICER; Gendreau et al. 2016), a soft X-ray telescope installed on the *International Space Station* (ISS) in 2017, was developed in part to estimate masses and radii of NSs using pulse-profile modeling of nearby rotation-powered millisecond pulsars (MSPs). Pulse-profile modeling is a technique that probes (approximations to) general relativistic effects on thermal emission from hot regions on the stellar surface (Bogdanov et al. 2019b); these effects are, predominantly, local radiation beaming<sup>21</sup> due to the bulk motion of material on the rotationally deformed surface,<sup>22</sup> and subsequent ray propagation on the exterior spacetime. Ray propagation includes the canonical *bending* of light, gravitational redshift, and the increasingly small imprints of rotational metric deformation: frame-dragging, a finite mass quadrupole moment, and higher-order (mass and current) multipole moments. As the star spins, the flux and spectrum of X-ray emission registered by a distant observer is modulated in a periodic manner: we can determine the rotational phase evolution of pulsars precisely and build up a pulse profile (X-ray counts per rotational phase bin per detector channel) by phase-folding X-ray events according to an ephemeris.<sup>23</sup> The mapping of surface emission into the pulse profile detected by a distant observer, via relativistic ray-tracing through the spacetime of a rapidly rotating (and hence oblate) star, is well understood (Pechenick et al. 1983; Miller & Lamb 1998; Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006; Cadeau et al. 2007; Morsink et al. 2007; Bauböck et al. 2013; AlGendy & Morsink 2014; Psaltis & Özel 2014; Nättilä & Pihajoki 2018; Vincent et al. 2018). Thus, given a model for the surface emission (e.g., a geometrically thin atmosphere of some chemical and ionic composition together with a local comoving effective temperature field as a function of surface coordinates), one calculates the expected pulse profile for a given exterior spacetime solution and a given instrument. By coupling such light-curve models to statistical sampling software via efficient software implementations, we can use Bayesian inference to derive posterior probability distributions for spacetime parameters such as mass and equatorial radius directly from pulse-profile data.

For the pulse-profile modeling technique to deliver tight constraints on mass and radius, rapid spin (≥100 Hz) is desirable (Psaltis et al. 2014; Miller & Lamb 2015; Stevens et al. 2016), and one needs high-quality phase- and energy-resolved pulse profiles with time resolution ≤10 µs and a large number of photons. The precise number of photons needed to deliver constraints on mass and radius at levels of a few percent—and by extension tight constraints on EOS models—depends on the geometry of a given source, but is ~10⁶ pulsed photons (Lo et al. 2013; Psaltis et al. 2014; Miller & Lamb 2015). For the brightest of the rotation-powered MSPs targeted by NICER, it is feasible to collect sufficient data with observation times ~1 Ms. The hot regions on rotation-powered MSPs in theory arise as

<footnote n="18"><sup>18</sup> Or measurements, or estimates; in any case, this means some probabilistic measure that is a function of, or otherwise pertains to, model parameters. See Section 2.3 for the probabilistic measures that we consider in this present work.</footnote>

<footnote n="19"><sup>19</sup> In a Bayesian context, by truncating the mass likelihood function only far in the tails, leading to a finite but small marginal posterior density for EOSs that do have substantially smaller maximum supported masses.</footnote>

<footnote n="20"><sup>20</sup> Note that PSR J1614−2230 was initially reported as having mass 1.97 ± 0.04 (Demorest et al. 2010); the inferences have since been updated via analysis of newly acquired data.</footnote>

<footnote n="21"><sup>21</sup> Note that local effective gravity in local comoving frames (instantaneously inertial during rotation) also enters calculation of atmospheric beaming of radiation emergent from the local comoving atmosphere.</footnote>

<footnote n="22"><sup>22</sup> Where for statistical applications the surface is either self-consistently computed via matching to a numerical interior solution to the field equations, or is embedded via a quasi-universal relation in an ambient spacetime solution (for an overview see, e.g., Riley et al. 2018, and references therein).</footnote>

<footnote n="23"><sup>23</sup> Note a key difference to the X-ray spectral modeling mentioned two paragraphs earlier: pulse-profile modeling involves phase-resolved spectroscopy: spectral modeling is phase-averaged, and does not fully leverage the temporal dimension of information provided by the star's rotation.</footnote>
magnetospheric currents—including return currents—deposit energy in the surface layers of the star; the resulting surface radiation field is a priori highly uncertain (Harding & Muslimov 2001; Gralla et al. 2017; Bauböck et al. 2019). NICER pulse-profile modeling can therefore also help to constrain the characteristics of the hot regions.

In this Letter we undertake pulse-profile modeling of NICER X-ray Timing Instrument (XTI) observations of the rotation-powered MSP PSR J0030+0451. Discovered as a radio pulsar by Lommen et al. (2000) and then identified as an X-ray pulsar (Becker et al. 2000), PSR J0030+0451 has a spin frequency of 205 Hz and lies at a distance of $325 \pm 9$ pc (Arzoumanian et al. 2018). There are no independent prior constraints on either mass or radius. Our analysis uses the X-ray Pulsation Simulation and Inference package (X-PSI$^{24}$ v0.1; Riley & Watts 2019). X-PSI is a software package for Bayesian modeling of astrophysical X-ray pulsations generated by the rotating, radiating surfaces of relativistic compact stars. X-PSI couples X-ray pulsation likelihood functionality to open-source statistical sampling software for use on high-performance computing systems; we apply nested sampling (Skilling 2006) in our analysis (refer to Appendix A). The work presented here is based on usage of a 500,000 core-hour grant on the Dutch national supercomputer Cartesius.$^{25}$

Section 2 outlines the modeling choices and introduces details specific to this analysis, including issues associated with the PSR J0030+0451 surface radiation field parameterization, instrument response, prior definition and implementation, and the consequences for computational efficiency of posterior sampling. In particular we restrict this analysis to models with two distinct hot regions with various structures. While our choices are physically motivated, it is important to emphasize that our inferences are conditional upon these choices. Posterior$^{26}$ inferences for the models are presented in Section 3, including the inferred posterior probability distributions for the spacetime parameters (mass, radius, and their combination into compactness), surface radiation field parameters (e.g., heating distribution and resulting temperature field), and instrument parameters. The Bayesian evidence for each model is reported, and we also summarize the computational resources required for each parameter estimation run. We also compare our inference to predictions derived via earlier study of PSR J0030+0451 using XMM-Newton (Bogdanov & Grindlay 2009). Sections 2 and 3 are long and detailed; in Section 2.1 we provide a brief overview of some of the key aspects from those sections, to help orient the reader. In Section 4 we discuss the implications of our results for our understanding of dense matter, pulsar emission mechanisms, and stellar evolution. We conclude with a discussion of future work for PSR J0030+0451: variations in the model that should be considered, tests and cross-checks, and the potential for improving the constraints for this source via longer observations or more in-depth analysis.

## 2. Modeling Procedure

### 2.1. Executive Summary of Modeling Procedure and Inferences

This Letter (and its companions) are the first pulse-profile modeling analyses to emerge from the NICER mission. We

have therefore provided (in Sections 2 and 3) a very detailed description of the methodology, the flow of the analysis, and the results. Since this is lengthy, we summarize the key aspects in this subsection.

The pulse-profile modeling technique requires us to define a model for the data-generating process, incorporating the physics that we initially assume to be most important. For a given choice of parameters, this model can be used to generate synthetic pulse-profile data sets. The model in part defines the likelihood function (the probability of the data as a function of parameters); the model also defines the prior probability distribution of parameters entering in the likelihood function. The posterior probability distribution of the model parameters (conditional on the observed data) is then sampled during the inference process. We must also pay attention to model complexity, in order to keep the computational load tractable.

In this Letter we assume that there are two separate hot regions on the stellar surface; this choice was motivated by the presence of two distinct pulses in the observed (phase-folded) pulse-profile. However, we considered a number of different possible configurations for the shapes and temperature functions of the hot regions: circular spots, annuli (rings, both centered and off-centered), and crescents; with one or two temperature components. These choices were motivated by contemporary theories for pulsar surface heating distributions as a result of magnetospheric return currents. We tested configurations where we insisted that the two hot regions were antipodal and identical; and where the hot regions were completely independent and potentially non-antipodal. We then assumed a geometrically thin fully ionized hydrogen atmosphere model (NSX, Ho & Heinke 2009), which characterizes the beaming and spectrum of the emergent thermal radiation (see, e.g., Zavlin et al. 1996, for details).

To propagate the emergent radiation toward the observer via relativistic ray-tracing, we use the Oblate Schwarzschild plus Doppler approximation of Morsink et al. (2007) for the NS spacetime. This is sufficiently accurate for our analysis, given the rotation rate of PSR J0030+0451. We define a joint prior distribution of mass and radius (the key parameters specifying the spacetime) that facilitates the subsequent inference of EOS model parameters (Riley et al. 2018; Raaijmakers et al. 2019). For the distance to the source, we use the (Gaussian-distributed) value inferred from radio observations (Arzoumanian et al. 2018) as a prior in our modeling. We then need to model the instrument response matrix (which includes both the effective area and the way in which incident photons of a given energy are assigned to specific detector energy channels). We develop a parameterized model that includes both energy-independent and energy-dependent components. The former attempts to capture absolute calibration uncertainty; for the latter we base our parameterization on residuals derived from NICER observations of the Crab Nebula and pulsar (Ludlam et al. 2018). We also assume a non-source background component, which we treat as a rotational phase-independent channel-by-channel contribution, rather than invoking a specific physics-driven spectral model. There are no prior constraints on either observer inclination or interstellar absorption for PSR J0030+0451, so we adopt a wide and diffuse prior for both parameters.

Note that Miller et al. (2019) have made an independent analysis of the same data set using different modeling choices and methodology. The choices we have made in this Letter

<footnote>$^{24}$ https://github.com/ThomasEdwardRiley/xpsi</footnote>

<footnote>$^{25}$ https://userinfo.surfsara.nl/systems/cartesius</footnote>

<footnote>$^{26}$ Appendix A provides an overview of the methodology used for posterior computation, and of the format used to present the posterior information.</footnote>
differ in several regards from those made by Miller et al. (2019); some of the most notable differences are in the models for the hot region configurations and the instrumental response, and the specification of the prior on distance.

During the inference analysis reported in Section 3 we considered a sequence of increasingly complex models for the shape and temperature function of the hot regions. All of the other aspects of the modeling described above are shared between models. Model assessment and comparison then enabled the identification of a superior configuration. We use a combination of performance measures: the evidence (the prior predictive probability of the data); graphical posterior predictive checking (to verify whether or not a model generates synthetic data without obvious residual systematic structure in comparison to the real data); visualization of the combined signals from the hot regions; Kullback-Leibler (KL) divergences (a measure of the parameter-by-parameter information gain of the posterior over the prior); background-marginalized likelihood functions (useful in combination with evidence to assess whether additional model complexity is helpful); model tractability (posterior computational accuracy being higher for less-complex models); and cross-checking of the inferred background against earlier analysis of PSR J0030+0451 with *XMM-Newton*.

Before beginning our analysis, we had mapped out an initial route through the model space of different heating configurations. This was modified as we progressed, informed by the results of each stage. We began with the simplest model, with single-temperature circular spots. Having the spots be antipodal and identical was quickly ruled out due to large residuals between model and data. Relaxing the requirement that the spots be identical and antipodal largely resolved this issue. We then moved to a more complex model where each hot region consisted of a circular spot—a core—and a surrounding annulus with an independently determined temperature. This model was superior to the simpler one, based on the evidence, but appeared to be overly complex: one hot region was dominated by a small hot circular spot, with negligible emission in the NICER waveband from the cooler annulus; for the other, emission was dominated by a hot annulus, with a much cooler core making almost no contribution. Simplifying the model such that one hot region was a single-temperature circular spot and the other a single-temperature annulus (with a centered, non-emitting core) produced congruent inferences at lower computational cost. At this point (after assessing the contribution from this component, and our remaining computational resources), we elected to restrict the model for one of the hot regions to be a single-temperature circular spot. The other was restricted to a single temperature, but we increased the complexity of the shape, testing two additional models: an annulus with an off-centered non-emitting core, and a crescent.

The superior configuration to emerge from this sequence of models, in terms of the performance measures listed above, was the final one: one hot region a small circular spot (sufficiently small that we would be insensitive to shape changes), and the other was an extended thin crescent. The results that we report in the abstract for the mass and radius of PSR J0030+0451 are those associated with this configuration.

## 2.2. Data Pre-processing

The NICER XTI data set *d* associated with the rotation-powered MSP PSR J0030+0451 is necessarily a product of

pre-processing, and its curation is largely described in Bogdanov et al. (2019a). However, there are a number of details that are specific to each data set and so we record a summary of these details below. The raw XTI event data are publicly available via the NICER archive (accessible via HEASARC$^{27}$). The processed data set *d* may be found in the persistent repository of Bogdanov et al. (2019c).

For this analysis we consider only the detector channel$^{28}$ subset [25, 300]—meaning channels 25 through 299 inclusive —nominally corresponding to 0.25–3 keV. Below channel 25, there is increased “optical loading” contamination (electronic noise due to ambient light), and there is greater uncertainty in the detector readout triggering efficiency for valid X-ray events. Above channel 300 the soft thermal emission from PSR J0030+0451 becomes negligible relative to the non-source background.

For PSR J0030+0451 we use 49 out of the 52 active detectors (excluding the three detectors that are frequently “hot,” i.e., prone to excess electronic noise) and only NICER pointings subtending an angle >80° to the direction of the Sun. In addition to the standard filtering criteria applied to XTI events, we excluded times where the planetary K-index $K_P \ge 5$, and excluded regions in the NICER orbit where the cutoff rigidity (the COR_SAX parameter) was less than 1.5. We further screened for occurrences of elevated background and eliminated all 16 s time intervals in which the 0.25–8.00 keV count rate exceeded 3 counts s$^{-1}$, yielding an integrated exposure time of $T_{\text{exp.}} = \sum_{\ell} \Delta t_{\ell} = 1,936,864$ s, where each $\Delta t_{\ell}$ is the time interval for the $\ell$th exposure.

Harnessing the flux and timing stability of PSR J0030+0451, the entire event data set in each detector channel is phase-folded coherently according to the best available radio pulsar timing solution from Arzoumanian et al. (2018) and using two approaches to check for consistency: using the TEMPO2 (Hobbs et al. 2006) photons plug-in and the PINT$^{29}$ photonphase tool. Differences of $\lesssim 1$ µs as a phase-offset are observed, but are deemed negligible for the present analysis as such an offset corresponds to $\lesssim 0.02\%$ of the total phase, or less than 0.7% of a bin width. The resulting folded event list, obtained with the PINT photonphase tool, is summarized as count data in a set of rotational phase intervals and detector channels. We display in Figure 1 the count data for channels in the interval [25, 300), nominally corresponding to the energy range 0.25–3.00 keV. We reserve the remaining details of data-space definition for Section 2.4.3, wherein we formalize the likelihood function applied in this specific work.

## 2.3. Generative Modeling

We begin the description of our modeling process by outlining useful mathematical objects that are assigned statistical meaning, and which can interface with open-source computational machinery. In Section 2.4, we build upon this conceptual groundwork and assign astrophysical meaning to such objects.

<footnote>$^{27}$ https://heasarc.gsfc.nasa.gov/</footnote>

<footnote>$^{28}$ Conventionally termed a pulse-invariant (PI) channel, or alternatively a charge pulse amplitude (or height) channel. Note the distinction between channel, which is only indicative of photon energy in units of 10 eV, and incident photon energy intervals that map probabilistically to the channel set via the redistribution matrix, a standard calibration product provided by the instrument team.</footnote>

<footnote>$^{29}$ https://github.com/nanograv/PINT</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Phase-folded PSR J0030+0451 event data split over two rotational cycles for clarity: we use 32 phase intervals (bins) per cycle and the count numbers in bins separated by one cycle—in a given channel—are identical. The total number of counts is given by the sum over all phase-channel pairs. The top panel displays the pulse-profile summed over the contiguous subset of channels [25, 300]. The bottom panel displays the phase-channel resolved count numbers for channel subset [25, 300]. For likelihood function evaluation (see Section 2.4.3) we group all event data registered in a given channel into phase intervals spanning a single rotational cycle. Moreover, we do not indicate the count-number noise in the top panel to avoid confusion: the observed events are viewed as fixed random variates that do not have errors, and whose parameterized joint sampling distribution is to be modeled (see Section 2.3).</figcaption></figure>

There exists an inherent degree of freedom as to the precise definition of a model. We consider a generative model for the X-ray data set d curated for PSR J0030+0451. Moreover, we opt to define a generative model as the union of the following components: (i) a data space D in which a data set d ∈ D exists as a fixed vector of numbers; (ii) an abstract model space M of elements, such that each element completely specifies a joint sampling distribution on space D (i.e., the sampling distribution can be evaluated conditional on the element); (iii) a joint probability distribution defined on space M that is not conditional on data vector d. Such a model is considered as generative because one can define a Bayesian joint distribution on the joint space of D and M: the data vector and a model vector are both interpreted as (finite-dimensional) multivariate random variables, and in order to simulate (or generate) data sets, one can jointly draw random variates to populate the vector elements.

The model space M is in general a discrete-continuous mixed space, meaning that its elements form a discrete-continuous mixed set; it follows that one could in principle identify each element as a model. Ultimately, however, one can distinguish continuous subsets of M as models, which is perhaps more common. To formalize the framework that we are working in further, let us define a discrete-continuous mixed space M as a union ⋃m∈F Mm, where F ⊂ N is the space of a discrete parameter m—a flag or label—and where each Mm may be considered as a model within M. Let each model Mm

define a continuous space Θm ⊆ ℝn_m such that M defines a space given by the union of Cartesian products

$$ \Theta := \bigcup_{m \in \mathcal{F}} (\{m\} \times \Theta_m). \qquad (1) $$

The number of dimensions of the space Θm is given by n_m; a parameter vector θ ∈ Θm ⊆ ℝn_m then has n_m elements. We thus consider each space in the union M as a model: unless explicitly stated otherwise, we hereafter consider a model as an element in a discrete set, where each element has an associated continuous parameter space Θm (among other constructs, e.g., probability measures such as a joint prior density distribution).

Working within the scope of a given model Mm with continuous parameter space Θm, the prior support Sm ⊆ Θm is a (compact) subset of ℝn_m on which the joint prior density is finite. The target distribution is the joint posterior density distribution, denoted by π(θ | d, Mm), which is related to the joint prior density distribution, p(θ | Mm), via the probability identity (Bayes' theorem)

$$ \pi(\theta | d, M_m)p(d | M_m) = p(\theta, d | M_m) \\ = p(d | \theta, M_m)p(\theta | M_m), \quad (2) $$

where p(d | θ, Mm) is the likelihood function for model Mm. The likelihood function is the sampling distribution on the space D evaluated at the fixed data vector d, as a function of θ, conditional on model Mm; the normalization of the posterior is
the prior predictive probability<fnref n="30" /> $p(d | \mathcal{M}_m)$ of the data conditional on model $\mathcal{M}_m$. We document the techniques implemented for posterior computation in Appendix A. It is necessarily the case that for construction of a generative model, a joint prior distribution injects a finite quantity of information about the elements of the model space. From a Bayesian perspective, models and parameters are random variables, whose joint prior distribution strictly encodes the information available before acquisition of—and computation given—data. Hereafter we do not need to typeset the symbol $\Theta$ (with or without a subscript) to represent a space, and thus we explicitly free it for a different use.

The model space that we consider is some mixture of $\mathcal{M}$-complete and $\mathcal{M}$-open (Vehtari & Ojanen 2012): we do not believe the true data-generating process exists in $\mathcal{M}$, nor do we believe that an element in $\mathcal{M}$ is the closest approximation (to the true data-generating process) that is achievable and tractable (in the future) given more resources (i.e., both cognition and computing time). However, we proceed as though we can imagine a true data-generating process to physically exist—i.e., is plausibly the product of simple physical laws and initial conditions—and that parameterized approximations can in principle approach the real process arbitrarily closely, at least in terms of predictive performance. Moreover, we view the model space for this specific work as effectively the best available to us given current resources, and consider it plausible that it is sufficiently rich to predict features in the procured X-ray event data. A subset of continuous parameters are shared by all models in the discrete set: these are considered to be of an (effectively) fundamental physical nature, such as the distance to a source system and parameters controlling the exterior spacetime solution of a compact star in general relativistic gravity. On the other hand, parameters are also defined that are of an overtly phenomenological nature and are intended to approximate reality only insofar as: (i) predicting—via interplay with the (almost) fundamental parameters—features observed in data, and (ii) supporting future development of physics models that are considered more realistic but crucially remain tractable for statistical falsification.

When we increment model complexity, it is intended to be in an intuitive and natural manner: we generally do so by breaking a form of symmetry, forming nested relationships between models (see Section 2.5 and then Section 2.5.7). Nevertheless, the form of the simplest model and the increments in complexity are ultimately of a subjective nature and exhibit a degree of arbitrariness—different practitioners would have defined their model spaces differently.

The discrete set of models that we condition on was determined partly out of consideration for available resources. Given that posterior computation cost generally increases with prior predictive complexity,<fnref n="31" /> we consider it justified to organize a set of increasingly costly problems wherein each is an extension of the former. The salient advantage in this respect is facilitating robustness: by monitoring and analyzing sampling processes operating on gradually more difficult problems, tractability can be gauged for subsequent problems given available computing resources. Moreover, if a pair of models form an exact or approximate nested relationship then

we expect the posterior parameter inferences to be consistent with that of the simpler model if the additional complexity is unhelpful.

Aside from resource management, and conditioned on the assumption that posterior computation is sufficiently accurate, one can pose the question of how much complexity is useful to capture structure in observational data. In a Bayesian framework one can in principle estimate prior predictive probabilities of data conditional on a model. It is generally argued that Occam's razor is inherent to prior predictive probability integrals: predictive complexity is penalized if predictions are not *expected* to be at least as commensurate with the data as a simpler (or nested) model. The interpretation of prior predictive probabilities is often fraught with problems, principally sensitivity to prior definition. In this work the joint prior distribution defined for the continuous parameters of each model is not rigorously chosen according to an information-theoretic criterion, nor to accurately quantify belief for overtly phenomenological parameters. However, the prior choices are viewed as being weakly informative<fnref n="32" /> for most parameters of interest in the absence of existing constraints, and are viewed as being consistent between models.

As acknowledged above, a widely held view is that evidence estimation<fnref n="33" /> does not solve the problem of model comparison. In order to evaluate model performance we thus employ both a form of graphical posterior-predictive checking, and prior predictive probabilities that hereafter we will refer to simply as *evidences*—a less accurate but canonical descriptor of $p(d | \mathcal{M})$. If the evidence increases, it is generally accurate to conclude that additional complexity is warranted; if evidence does not increase, however, graphical posterior predictive checking on local modes is useful for determining whether or not facets of the higher-complexity model are a promising avenue to pursue in model development—i.e., if the likelihood function maxima are larger.

## 2.4. Overarching Definitions

In this section we describe model aspects that are generally shared between all models in the discrete set $\mathcal{M}$. These model facets are in some cases described in detail elsewhere in the NICER literature (Bogdanov et al. 2019a, 2019b, 2019c) and thus we are brief where possible. Due to the large number of symbols required to describe the models in this Letter, the symbols used to describe geometric variables shared with Bogdanov et al. (2019b) are different; Table 1 provides symbol translation from Bogdanov et al. (2019b), the theory in which underpins the present work.

### 2.4.1. Source

The pulsed sources are assumed to be thermally emitting, rotating hot surface regions of PSR J0030+0451.

**Parameterization.** The exterior spacetime solution is approximated as follows: we embed in each temporal hyperslice of an ambient Schwarzschild spacetime, a (quasi-universal) oblate two-surface, such that the geometric center coincides with the origin of the Schwarzschild coordinate chart (Morsink et al. 2007). The coordinate equatorial radius is

<footnote n="30">³⁰ The scalar expectation of the likelihood function with respect to the prior $p(\theta | \mathcal{M}_m)$, commonly referred to as the evidence or fully marginal likelihood.</footnote>

<footnote n="31">³¹ Which is easily proven for (non-dynamic) nested sampling processes (see Appendix A for an outline of the sampling procedure).</footnote>

<footnote n="32">³² Also known as *vague* or *diffuse*.</footnote>

<footnote n="33">³³ For calculation of Bayes' factors.</footnote>
<table><caption>**Table 1**
Translation of Symbols for Angle Variables Typeset in *both* this Work and Bogdanov et al. (2019b)</caption><thead><tr><th>Symbol</th><th>Description</th><th>Bogdanov et al. (2019b)</th></tr></thead><tbody><tr><td><i>i</i></td><td>Earth inclination to pulsar rotation axis</td><td>ζ</td></tr><tr><td>Θ</td><td>colatitude of center of a circular hot spot<sup>a</sup></td><td>θ<sub>c</sub></td></tr><tr><td>ζ, ψ</td><td>angular radius of a circular hot spot<sup>b</sup></td><td>θ<sub>spot</sub></td></tr><tr><td>φ</td><td>pulsar rotational phase<sup>c</sup></td><td>φ</td></tr></tbody></table>

**Notes.**

<sup>a</sup> We also use the symbol Θ, with subscripts, to denote the colatitudes of hot regions whose shapes are more complex than circular spots.

<sup>b</sup> We also use these symbols to parameterize hot regions with more complex shapes, such as rings with outer angular radius ζ and inner angular radius ψ.

<sup>c</sup> We also use the symbol φ, with subscripts, to denote the azimuthal coordinates of hot regions.

denoted by $R_{eq}$, and the circumference of the equator is $2\pi R_{eq}$; the total mass in the ambient spacetime is denoted by $M$. The polar axis of the Schwarzschild chart is defined as the pulsar rotational axis. Rotational deformation of the metric away from spherical symmetry is neglected; the current dipole and mass quadrupole moments of the exterior metric enter at first- and second-order in dimensionless angular velocity $\bar{\Omega}$ (Hartle 1967) but are sufficiently small—regardless of the EOS—that they can be neglected in pursuit of a tractable likelihood function (Morsink et al. 2007), especially during initial modeling and in a context where we do not expect the likelihood function to be sensitive to the finiteness of these moments. The perturbations, at constant baryon number, to both the total mass and circumferential radius of a nonrotating star are second-order in $\bar{\Omega}$ (Hartle 1967), and are thus small for a spin of 205 Hz (see also Section 4.2) but are implicitly accounted for. The perturbation to the polar coordinate radius is second-order in $\bar{\Omega}$, and the surface oblateness is controlled by an EOS-insensitive constraint equation (as is the effective gravity along the surface; Morsink et al. 2007). These small changes in the shape of the stellar cross-section induce tilt to the surface, which affects the rays that connect spacetime events at the surface to an observer; together with the change in the projected surface area of a tilted surface, and the change in effective gravity, the effect on light-curves manifests at first-order in $\bar{\Omega}$. Crucially, there is a performance floor for light-curve integration demarcated by a spherically symmetric exterior spacetime solution: embedding an oblate surface in such an ambient spacetime results in negligible increase in computation time per call to a light-curve integrator. We therefore are not concerned about quantifying the difference in our statistical inferences due to inclusion of oblateness over a spherical surface—more resources are required to quantify this rigorously than to simply account for oblateness.

A distant, static, and fictitious<sup>34</sup> instrument (see Section 2.4.2) is located at radial coordinate $D$, and subtends colatitude—hereafter termed inclination—denoted by $i$. Interstellar light-matter interaction is described by absorption within a column of material. The attenuation factor is parameterized solely by the column density $N_H$ of neutral hydrogen and we assume relative abundances for the interstellar medium from Wilms et al. (2000); we implement the tbnew<sup>35</sup> model to precompute a set of lookup tables for attenuation as a function of photon energy.

In each model the spatial dependence of the surface radiation field is of a phenomenological nature: the aim is to introduce sufficient complexity so as to represent the basic notion of pulsar surface heating due to energy deposition by magnetospheric currents (in the vicinity of the magnetic poles). We are largely ignorant of spatial structure in the surface radiation field because the star is not spatially resolved; moreover, it is intractable for us to consider more self-consistent numerical models of the surface radiation field, in part due to the expense of statistical computation. In the simplest case the radiation field is constructed by filling two closed simply connected regions on the surface—which do not mutually overlap—with radiating material;<sup>36</sup> these regions may be interpreted to each result from magnetospheric polar cap heating. We only compute a radiative signal from these hot regions, and therefore in the context of each of our models, *hot region* can be viewed as synonymous with *radiating region*.

For all models a geometrically thin (and thus plane-parallel) fully ionized hydrogen NSX atmosphere is invoked for the radiating material (Ho & Lai 2001; Ho & Heinke 2009). The radiation field is precomputed and represented as a lookup table for cubic polynomial interpolation of specific intensity, $I_E/k_B T_{eff}^3$, with respect to four variables defined in a surface local comoving frame: effective temperature, $T_{eff}$; effective gravity; photon energy, $E/k_B T_{eff}$; and the cosine of the ray zenith angle (to surface normal). A quasi-universal relation for surface effective gravity<sup>37</sup> is adopted from AlGendy & Morsink (2014) in order to evaluate local radiation field intensities. We do not explicitly compute emission from the stellar surface exterior to the closed regions (i.e., as a function of source parameters controlling the exterior radiation field),<sup>38</sup> but the phase-invariant background model that we invoke in all cases

<footnote n="36"><sup>36</sup> Note that although a cooler radiating hydrogen atmosphere should exist globally over the stellar surface (as observed for PSR J0437-4715, Durant et al. 2012; Guillot et al. 2016; Gonzalez-Caniulef et al. 2019, and PSR J2124-3358, Rangelov et al. 2017), we make no explicit reference to it when defining our likelihood function—i.e., we do not compute any radiative signal from the atmosphere exterior to the closed regions. The atmosphere cannot be globally uniform because local heating by magnetospheric currents will affect the local temperature and ionization degree; effective gravity also varies due to rotation. Reference to the global atmosphere is implicit due to the fluid properties required for containment of hot material in the closed regions.</footnote>

<footnote n="37"><sup>37</sup> Equatorially reflection-symmetric.</footnote>

<footnote n="38"><sup>38</sup> Such computation is supported by X-PSI with specialization to ensure (almost) exact areas as described above (see also Bogdanov et al. 2019b), but requires a choice of surface radiation field (e.g., atmosphere ionization degree and chemical composition). We therefore opt to capture the *non-pulsed* fraction of emission from the stellar surface exterior to the regions via our default background treatment. If evidence for unmodeled soft pulsed emission arises a posteriori one could then consider explicit computation of such emission.</footnote>

<footnote n="34"><sup>34</sup> A notion borne from the nature of the event-data pre-processing (see, e.g., Riley & Watts 2019, for details pertaining to X-PSI).</footnote>

<footnote n="35"><sup>35</sup> https://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/</footnote>
can capture non-pulsed components of surface and off-surface emission (see Section 2.4.3).

Let us hereafter refer to the geometric configuration of the infinitesimal radiating surface elements simply as the *shape* of a hot region—including both exterior and interior boundaries. In more complex models each hot region is constructed using additional shape parameters and in some cases a second temperature component. We consider the shape and temperature of a hot region to result from the interaction of *two* closed regions—hereafter *members*. The members partially or wholly<fnref n="39" /> overlap, and one member takes precedence when evaluating local radiation intensities along rays (null geodesics) that connect spacetime events on the rotating surface to a distant observer. We term one member as *ceding*, and the other as *superseding*. As discussed in the appendix of Bogdanov et al. (2019b), the X-PSI implementation (referred to as the AMS code in Bogdanov et al. (2019b)) of a radiating region is specialized for fast likelihood function evaluation when said function is a callback for sampling processes: while numerical approximations are necessary in general, it is relatively inexpensive to ensure that the proper area of each finite-element (discretized) radiating region is computed to a precision that (almost) exactly matches their mathematical definition.<fnref n="40" /> When two members overlap to form a hot region (i) the area of a discretely represented superseding member is (almost) exact; and (ii) the area the discretely represented, *non-superseded subset* of a ceding member is also (almost) exact. The subset that is not superseded can itself be simply connected or non-simply connected depending on the model (where the set of configurations assigned finite prior density is model dependent). For the precise details of the hot regions, refer forward to Sections 2.5.1 through 2.5.7.

We opt for two disjoint hot regions: the two distinct pulses visible in the phase-folded event data (Figure 1) are suggestive of two such regions being widely separated. Initially, we impose parity in the complexity of each hot region—i.e., an equal number of shape parameters and temperature components. We then consider models in which the hot regions have unequal complexities because it becomes clear that increasing the complexity of a particular hot region is unwarranted. In all cases we define the support of the joint prior so as to exclude limiting configurations in which the hot regions overlap; the reason being that extension of scope to such configurations requires specification and implementation<fnref n="41" /> of additional logical conditions for a complete order of precedence in local radiation intensity evaluation. For a subset of models we impose antipodal symmetry of the hot regions, in order to crudely represent a heating distribution that is consistent with symmetry in the physical mechanisms driving surface X-ray emission, such as a dominantly or perfectly centered-dipolar field configuration. A magnetic field with finite higher-order structure is viewed as a closer approximation of physical reality; we represent this case crudely by breaking antipodal symmetry and defining additional parameters for a secondary hot region that are not derived from parameters of a primary hot region. However, self-consistent coupling of the magnetosphere to the surface radiation field is beyond the scope of this

work, partly due to the associated increase in complexity of efficient model implementation for posterior computation.

We consider models with three or more mutually disjoint and separated radiating regions to be a logical extension of the model space if it is deemed that incrementally increasing the complexity of only two such regions is yielding insufficient advancement in posterior predictive performance—i.e., is not satisfactorily capturing observed structure in data for the resource expenditure—when approaching or extrapolating to the limit of what is considered computationally tractable by a group executing posterior computation. A salient advantage of such an approach is that it is more exhaustive with ideas for two hot regions and incrementally breaks symmetries; an apparent disadvantage is that some small set of conceivable closely related models with equal (continuous parameter) dimensionality are not applied.<fnref n="42" /> It is, however, necessary to be selective—inherent to which is subjectivity and arbitrariness.

The above choices for surface radiation field configuration are somewhat consistent with the notion of the source being a rotation-powered X-ray pulsar with two relatively small hot regions that are disjoint. Therefore, we consider the proposed model space to be logically structured and a reasonable representation of widely held conceptions of such stars that are yet to be falsified statistically.

Whether or not such a model space is tractable given algorithm properties and computing resources is highly sensitive to the choice of parameterization for posterior computation, especially in phenomenological contexts. If degeneracies plague the problem at hand, it may be considered an indicator that either (i) the model is simply ill-defined, leading to forms of invariance of the parameterized sampling distribution on the space of the data; (ii) the model is needlessly or at least *unhelpfully* complex for describing observations, because despite (physical) parameters having a finite effect on forward data-generation, one is ultimately insensitive to such model structure. Usually this equates—at least in part—to transforming away nonlinear likelihood function degeneracies where possible. A number of sophisticated open-source sampling software packages efficiently handle linear degeneracies, even in multi-modal contexts, but nonlinear degeneracy remains fiendish: certain sampling algorithms can perform accurately,<fnref n="43" /> but coupled with an expensive (numerical) likelihood function and moderately high-dimensional sampling spaces, still require massive computing resources.

As we highlight in Section 2.5 (where we provide more precise definitions of surface radiation field structure), the choice of parameterization of a (largely phenomenological) hot region plays a crucial role in sampling-space definition. In Appendix A we summarize the techniques adopted for posterior computation: we opt to perform nested sampling (Skilling 2006). The natural space for nested sampling is usually that of a unit hypercube, which maps to an equal-dimensional physical parameter space according to an inverse transformation derived from a joint (prior) probability distribution on the physical space (e.g., Feroz et al. 2009). Given that nested sampling algorithms tend to operate in such a native space, a parameterization that approaches optimality involves

<footnote n="39">³⁹ Such that one member is a superset of the other.</footnote>

<footnote n="40">⁴⁰ The overall numerical accuracy remains implementation dependent.</footnote>

<footnote n="41">⁴¹ Efficient finite-element representations.</footnote>

<footnote n="42">⁴² For example, configurations in which three single-temperature regions are disjoint.</footnote>

<footnote n="43">⁴³ At reduced efficiency relative to simpler contexts—see the MULTINEST sampler cited in Appendix A.</footnote>
both the physical parameterization and the inverse transformation from the native to physical space, and is such that continuous likelihood function degeneracies, if existent, manifest effectively linearly in the native space.

Usually a joint prior distribution is chosen to be weakly informative—or “flat”—in the context of the likelihood function, and in some cases is defined as an absolutely flat density function with respect to a joint space. It follows that the mapping from the parameter space to the native sampling space will then approximately preserve the linear degeneracy of a posterior mode. An exception to this occurs if the boundary of the support of the joint prior satisfies some set of non-trivial constraint equations.<fnref n="44" />

*Priors.* We define the joint prior density distribution $p(M, R_{eq})$ to be jointly flat with compact support: a prophylactic choice that eases future use of samples on the $(M, R_{eq})$-subspace for computing an approximative marginal likelihood function, which in turn can be used for estimation of interior source-matter properties—principally EOS parameters (see Riley et al. 2018; Riley & Watts 2019). We also choose the boundary of prior support to be close to maximally inclusive in regards to theoretical EOS predictions: we impose hard bounds $M \in [1, 3] M_{\odot}$, and impose that $R_{eq} \in [3 r_g, 16]$ km, where $r_g = r_g(M) = GM/c^2$ is the gravitational radius.

Although we allow the boundary of the prior support to extend down to the photon sphere of the ambient spacetime solution, when computing pulse profiles we only integrate over the primary images (along rays with angular deflection $\le\pi$) of radiating elements subtended on the sky of the instrument. For a spherical star of radius $R < 3.52 r_g(M)$, multiple images of parts of the star will be visible (Pechenick et al. 1983; Bogdanov et al. 2019b) requiring that light from the primary, secondary, and higher-order images be included—at additional computational expense—in an exact calculation of the flux; this issue, and how it pertains to oblate stars, is discussed in more detail elsewhere (Bogdanov et al. 2019b). From a computational statistics perspective, one could view the inclusion of one or more higher-order images as a modeling refinement to be made if, a posteriori, a rotating star is favored to be sufficiently compact: e.g., a substantial fraction of posterior mass lies at $R_{eq}/r_g(M) \lesssim 3.6$. The parameter inference reported in this work favors much less compact stars a posteriori (refer forward to Figure 19), so multiple imaging is not deemed important. However, when images are neglected, the issue of choosing the most appropriate support for a joint prior distribution of $M$ and $R_{eq}$ remains an open problem for statistical modeling.

A typical likelihood function for pulse-profile modeling will express many modes of dependence on the compactness $r_g(M)/R_{eq}$, and will generally be more sensitive to this combination than to $M$ (or $R_{eq}$) individually. We ensure the mapping from the parameter space to the native space preserves such linear degeneracy between $M$ and $R_{eq}$. In Appendix B we provide implementation details for the joint density $p(M, R_{eq})$.

Finally, there exists a constraint on the distance of the PSR J0030+0451 system (Arzoumanian et al. 2018) that we adopt—in approximation—as an informative prior.

### 2.4.2. Instrument

In defining a generative model, the data space is constructed by phase-folding X-ray events in each detector channel and grouping those events into a uniform set of phase intervals (bins) to curate a set count numbers, typically with cardinality $O(10^3)$. The conditional joint sampling distribution of these count numbers is always constructed in terms of a phase-energy-resolved signal that is generated during a single rotation of the source. Together with an appropriate nuisance background model (see Section 2.4.3), it follows that the instrument in such an analysis represents the temporal-mean operation of all detectors collectively in response to the incident radiation field from the source during the observation time intervals. We reserve a more elaborate discussion on these modeling facets for Section 2.4.3 (and also refer the reader to Riley & Watts 2019).

For every model we invoke the instrument response model: the on-axis $v1.02$ ancillary response function (ARF) and an updated version of the $v1.02$ redistribution matrix file (RMF; private communication from James Steiner, see Hamaguchi et al. 2019, for details of the updates) to generate a reference (or nominal) response matrix $R^*$ derived from microphysical knowledge. Let detector channels increment with row number $i$, and energy intervals increment with column number $j$, such that an element of the reference response matrix is denoted by $R_{ij}^*$. We use this reference matrix as a basis for a parameterized family of response matrices, and aim to compute (for each model) a joint posterior density distribution of continuous source parameters and continuous instrument (response matrix) parameters. It is well founded to parameterize the instrument because despite its synthetic nature, we do not consider its microphysical operation to be sufficiently *known*;<fnref n="45" /> nevertheless, we define far fewer parameters for the instrument than for the source. While the following model ensures that operational uncertainty is included, the continuum of response models and the associated prior density distribution does not attempt to rigorously represent uncertainty in microphysical knowledge; given a close approximation to the radiation field incident on the telescope, we would expect the model to be conservative in terms of *prior* predictive performance.

We parameterize the response matrix using a calibration product derived from observations of a calibration source. For this work we use instrumental residuals derived from NICER observations of the Crab. These residuals are derived using the observations and following the procedure outlined in Ludlam et al. (2018), modified to use the appropriate number of detectors, ARF, RMF, and Sun-angle cut consistent with the PSR J0030+0451 data set (refer to Section 2.2). We acknowledge that the Crab is a remarkably different source to PSR J0030+0451: the expected operation of NICER (and X-ray instruments in general) in response to incident radiation fields is a function of its properties. The Crab exhibits a very different spectrum to PSR J0030+0451, being harder, more absorbed, and subject to astrophysical features; the Crab is also an extended source, not a point source, and is much brighter than the rotation-powered MSPs targeted by NICER.

The calibration product is a channel-by-channel vector $\mathcal{R}$ of ratios of observed Crab count numbers to count numbers derived

<footnote n="44">⁴⁴ An example of trivial constraint equations here are those that generate a (hyper-)rectangular support boundary in parameter space.</footnote>

<footnote n="45">⁴⁵ In-flight astrophysical calibration sources, for instance, are in practice far brighter than science targets, and operation is conditional on the radiation field incident on the detectors.</footnote>
using a theoretical incident spectrum and the reference response matrix $\mathcal{R}^*$. Let the vector elements $\mathcal{R}_i := C_i/(\mathcal{R}^* \cdot \mathcal{F})_i$ in words, the element-wise division of an observed count vector $C$ from some calibration source, by a vector $\mathcal{R}^* \cdot \mathcal{F}$ where $\mathcal{F}$ is a vector of photon fluences (in the set of energy intervals inherent to the definition of $\mathcal{R}^*$) computed given some theoretical model of said calibration source during the calibration observations. Note that the calibration product is derived from observations of a single chromatic source and thus is not resolved over elements of the matrix, only over the set of channels; therefore we apply the ratio $\mathcal{R}_i$ for the ith channel to all elements $\mathcal{R}_{ij}^*$. The theoretical model is uncertain in the lowest 10 channels that we consider: the Crab is highly absorbed so that there is less data at low energies, and the telescopes used to generate the reference spectra for the residuals also perform poorly in this regime. In this work we therefore assume that $\mathcal{R}_i$ for $i \in [25, 35)$ is equal to $\mathcal{R}_{35}$.

We choose to construct the response matrix as a continuous three-parameter family, where the parameters are denoted by NICER $\alpha$, NICER $\beta$, and NICER $\gamma$ where possible, but reduced to the aliases $\alpha, \beta$, and $\gamma$, respectively, for clarity of mathematical expressions. First we give the definition, and then we offer an interpretation in words. The parameterized matrix is defined as

$$ \begin{aligned} \mathcal{R}_{ij}(\alpha, \beta, \gamma; \mathcal{R}, \mathcal{R}^*) &:= (1 - \beta)\gamma\mathcal{R}_{ij}^* + \beta\alpha\mathcal{R}_i\mathcal{R}_{ij}^* \\ &\equiv \gamma\mathcal{R}_{ij}^* + \beta\frac{[\alpha\mathcal{R}_i\mathcal{R}_{ij}^* - \gamma\mathcal{R}_{ij}^*]}{\text{calibration shift}} \\ &\equiv \mathcal{R}_{ij}^*[\gamma + \beta(\alpha\mathcal{R}_i - \gamma)]. \end{aligned} \quad (3) $$

In general an additional condition must be invoked: $\mathcal{R}_{ij}(\alpha, \beta, \gamma) := 0$ if $[(1 - \beta)\gamma\mathcal{R}_{ij}^* + \beta\alpha\mathcal{R}_i\mathcal{R}_{ij}^*] \le 0$ for $(\alpha, \beta, \gamma) \in S_{(\alpha, \beta, \gamma)}$, where $S_{(\alpha, \beta, \gamma)}$ is the joint prior support on the $(\alpha, \beta, \gamma)$-subspace.<fnref n="46" /> We note, however, that this parameterization exhibits a finite degree of degeneracy.

The parameter $\alpha$ scales the calibration vector $\mathcal{R}$, and manifests to target the assumption that the vectors $C$ and $\mathcal{F}$ are known: in regards to $C$, the assumption that the temporal-mean operation of the instrument between calibration and science observations is invariant (including, e.g., pointing vector relative to source line of sight and the flux-dependent effects); and in regards to $\mathcal{F}$, the assumption that the expectation of the incident radiation field from the Crab during the calibration observations is known. Thus, the parameter $\alpha$ represents the product of dimensionless element-invariant scaling factors applied, respectively, to vectors $C$ and $\mathcal{F}$. The support of $\alpha$ is such that $S_\alpha = \{\alpha: \alpha \in \mathbb{R}^{>0} \wedge 1 - \epsilon_\alpha \le \alpha \le 1 + \epsilon_\alpha\}$ where $0 < \epsilon_\alpha < 1$; the prior density function allocates mass mostly to the near vicinity of unity.

The parameter $\gamma$ is an element-invariant scaling factor applied to the reference matrix $\mathcal{R}^*$ alone. The reason that we choose $\alpha \ne \gamma$ is that $\gamma$ operates on $\mathcal{R}^*$, and thus because $\mathcal{R}$ is a function of $\mathcal{R}^*$ (via the element-wise vector division written above), $\gamma$ cancels. The support of $\gamma$ is such that $S_\gamma = \{\gamma: \gamma \in \mathbb{R}^{>0} \wedge 1 - \epsilon_\gamma \le \gamma \le 1 + \epsilon_\gamma\}$ where $0 < \epsilon_\gamma < 1$; the

prior density function allocates mass mostly to the near vicinity of unity.

The parameter $\beta$ is a weighting factor between matrices with elements $\gamma\mathcal{R}_{ij}^*$ and $\alpha\mathcal{R}_i\mathcal{R}_{ij}^*$—i.e., the element-invariant coefficient of the element-dependent calibration shift away from the $\gamma\mathcal{R}^*$ matrix. Thus $\beta$ may be interpreted as the degree to which the calibration vector $\mathcal{R}$ is used to modify the reference matrix $\mathcal{R}^*$. Note that the prior support for a weighting parameter such as $\beta$ is the unit interval $\beta \in [0, 1]$. We consider the limit $\beta \to 0$ as a useful safeguard against erroneous calibration—e.g., artifacts may be introduced by invoking a calibration source to which the instrument responds appreciably differently than it does to PSR J0030+0451. As $\beta \to 0$, $\mathcal{R}_{ij}(\alpha, \beta, \gamma) \to \gamma\mathcal{R}_{ij}^*$, capturing a simple element-invariant scaling of $\mathcal{R}^*$. On the other hand, in the limit $\beta \to 1$, we have $\mathcal{R}_{ij}(\alpha, \beta, \gamma) \to \alpha\mathcal{R}_i\mathcal{R}_{ij}^*$.

*Priors.* To demonstrate the properties of the family of response matrices defined above, we (i) define a joint probability density distribution $p(\alpha, \beta, \gamma)$ that could be plausibly viewed as a weak prior on its compact support given that the instrument is an artificial system that has been closely studied, and (ii) then visualize the distribution in terms of derived properties (such as total effective area as a function on energy interval). Let $\alpha \sim N(1, 0.1)$ truncated such that $\alpha \in [0.5, 1.5]$; $\beta \sim U(0, 1)$; $\gamma \sim N(1, 0.1)$, truncated such that $\gamma \in [0.5, 1.5]$. In Figure 2 we display the prior distribution of the total on-axis effective area as a function of energy.

### 2.4.3. Likelihood Functions and Background

In the X-PSI documentation (Riley & Watts 2019) we offer a more complete overview of the supported class of generative models than is appropriate for this work—instead we adapt the X-PSI documentation to provide a summary. A generative model for raw on-board event data is eschewed by subsuming a non-parameterized portion of the modeling within a data pre-processing phase. In this work, a pulsar radio timing solution is invoked to transform events into a simpler time domain: a fictitious instrument, which is static (or Eulerian) and distant in the (Schwarzschild) spacetime of the source, is implicitly constructed to register events against the elapsed natural number rotations of the star, which is a clock related to the Schwarzschild coordinate time simply by an affine transformation. In all models we condition on the phase-resolved specific flux signal (incident on the instrument) generated by precisely one rotation of the star: accurate computation of such a signal, even when invoking spacetime spherical symmetry, is approaching the limit of what we consider tractable at present in terms of likelihood function callback cost for sampling processes in $O(10)$-dimensional spaces.

In the context of a model for the joint probability distribution of observed events (see Riley & Watts 2019),<fnref n="47" /> phase-folding said events is equivalent to computing one (average) rotational pulse and replicating it over the many rotational cycles in order to evaluate the likelihood function. It follows that, in this limit, no information is lost by transforming events to the unit interval because the underlying information content in the model is not summarized for comparison to data. In each instrument channel we choose to group events into a set of uniform-width phase intervals (bins) that are subsets of the unit interval, and define the data space as $\mathcal{D} := \mathbb{N}^{I \times K}$, where $I \in \mathbb{N}$

<footnote n="46">⁴⁶ Alternatively, one might define the joint prior support such that if for $(\alpha, \beta) = (\alpha', \beta', \gamma')$ the inequality is true for any $(i, j)$, then the joint density is locally zero at $(\alpha', \beta', \gamma')$.</footnote>

<footnote n="47">⁴⁷ Leading to an unbinned likelihood function.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** In the left panel we display marginal conditional prior probability density distributions of total on-axis effective area A as a function of energy (integrated over uniform intervals of width $5 \times 10^{-3}$ keV), summed over the contiguous channel subset [25, 300). Note that the total NICER effective area at energies $\gtrsim 3$ keV in particular is far greater than shown here. The prior $p(A; E)$ is represented by the blue contours. At energy E a set of (one-dimensional) highest-density credible intervals are estimated for $A = A(\alpha, \beta, \gamma)$; the credible intervals are connected as a function of energy, such that, e.g., the second-darkest band encodes the energy-dependence of the estimated 68.3% highest-density credible interval. We give the marginal prior distributions for the parameters $\alpha, \beta$, and $\gamma$ in Table 2. Note that the posterior information shown is not that of the joint distribution of effective areas over energy intervals: the effective areas are coupled by a functional form with three parameters. The range of the energy intervals is determined based on the RMF of the reference matrix $\mathcal{R}^*$ and on the curated data set in the contiguous set of channels [25, 300), where these channels are summed over. The solid curve is that of the reference matrix $\mathcal{R}^*$ (with $\gamma = 1$); the dashed-dotted curve is that of the calibrated matrix $\mathcal{R}_{ij} = \mathcal{R}_i \mathcal{R}_{ij}^*$ (with $\alpha = 1$). In the right panel we display as black points the elements of the vector $\mathcal{R}$ of multipliers supplied for instrument parameterization based on the Crab as a calibration source. The prior $p(\eta, i)$, where $i$ enumerates channels and $\eta$ is defined in Equation (3), is represented by the blue contours; the probabilistic information is otherwise congruent in nature to that described for the left panel.</figcaption></figure>

is the number of channels over which the folded events are distributed, and $K \in \mathbb{N}$ is the number of phase intervals; we choose $K = 32$ and $I = 275$ corresponding to channel subset [25, 300). Information loss is an inherent consequence of compression of events into a smaller set of summary quantities, but the phase resolution is sufficiently high here to mitigate our concern about the use of a binned likelihood function instead of an unbinned likelihood function. The conditional joint sampling distribution on the space $\mathcal{D}$ is assumed to be purely Poissonian and separable over channels: the NICER instrument exhibits sufficiently high-resolution event-timing capabilities for the Poissonian nature of the incident radiation field to be effectively conserved as an event arrival process in the on-board time domain (see Riley & Watts 2019, for a more explicit set of arguments pertaining to this matter).

Let the likelihood function for phase-folded and binned events be defined by (see also Miller & Lamb 2015)

$$L(\boldsymbol{\theta}, \mathbf{B}) := p(\{d_i\}_{i=1,...,I} | \boldsymbol{\theta}, \{B_i\}) = \prod_{i,k} p(d_{ik} | \boldsymbol{\theta}, B_i), \quad (4)$$

where $i \in [1, I]$ enumerates channels of the instrument; each $d_i$ is a data vector associated with the $i$th channel, constituted by count numbers $\{d_{ik}\}_{k=1,...,K}$ where $k$ enumerates phase intervals $\phi_k \subset [0, 1]$; $B_i$ is the (nuisance) background count-rate parameter in the $i$th channel; and $\boldsymbol{\theta}$ are the continuous source parameters that constitute the sampling space, and on which source expected count numbers $s_{ik}(\boldsymbol{\theta})$ are dependent. Each $B_i$ is defined as the expectation of a homogeneous—i.e., time

invariant—Poisson arrival process: up to a known constant,

$$\ln L(\boldsymbol{\theta}, \mathbf{B}) = -2 \sum_{i,k} d_{ik} \ln \left[ s_{ik}(\boldsymbol{\theta}) + B_i \left( \sum_{\ell} \Delta t_{\ell} \right) \int_{\phi_k} d\phi \right] - s_{ik}(\boldsymbol{\theta}) - B_i \left( \sum_{\ell} \Delta t_{\ell} \right) \int_{\phi_k} d\phi. \quad (5)$$

Such background treatment is the default for NICER parameter estimation work and X-PSI (for implementation details see Appendix B of Riley & Watts 2019), and is based on Miller & Lamb (2015). In reality, the statistical properties of backgrounds do not exhibit time invariance, but exhibit long-term variation over the $\mathcal{O}(1)$ year observation epoch of MSPs target by the NICER mission, especially due to factors such as dynamical space weather (Bogdanov et al. 2019a). However, any background emission processes—and dynamical emission processes in the local vicinity of PSR J0030+0451—that are not harmonically coupled to the surface X-ray emission will *decohere* over the unit phase interval whose boundary is periodic. It follows that event phase-folding enables invocation of phase-invariant channel-by-channel (background) count-rate terms.

The background below $\sim 3$ keV for NICER observations of PSR J0030+0451 observations consists of (Bogdanov et al. 2019a) (i) cosmic energetic particle events and diffuse X-ray emission over the $\sim 30 \text{ arcmin}^2$ field of view (Arzoumanian et al. 2014); (ii) many nearby X-ray point sources in the field
that make a small total contribution relative to the targeted MSP; and (iii) solar system contamination, including optical loading (pointing Sun-angle dependent), and high-energy non-cosmic particles and radiation. Considering a proper subset of detector channels and filtering background events during the pre-processing phase acts to reduce background contribution, but some subset of background events survive and must be modeled (Bogdanov et al. 2019a).

The source terms $s_{ik}(\theta)$ in Equation (5) are then derived as follows. Let $F(\phi, E; \theta)$ denote the incident specific photon flux from the source as a function of rotational phase $\phi$. The function $F(\phi, E; \theta)$ is evaluated numerically at a regular discrete set of points in the joint space of energy and phase, as an approximating two-dimensional integral over the solid angle of the image of the source subtended on the sky of a point in the vicinity of the distant static instrument; given the discrete representation, a continuous representation is constructed via spline interpolation in X-PSI.

Let the symbol $R_{ij}(\hat{r})$ denote a temporal-mean point-source response matrix invoked for the $l$th observing interval, which is dependent on the radial coordinate unit vector $\hat{r}$ in the Schwarzschild chart at the location of the instrument—i.e., the pointing of the telescope relative to the source. The response matrix with elements

$$R_{ij}(\hat{r}) = \sum_{l} \frac{\Delta t_{l}}{(\sum_{l} \Delta t_{l})} R_{ij}(\hat{r}) := R_{ij}(\alpha, \beta, \gamma) \quad (6)$$

is the exposure-time-weighted mean response matrix that is modeled as the matrix defined in Equation (3) and Section 2.4.2.

The source contribution to the expected number of counts is given by

$$\begin{aligned} s_{ik}(\theta) &= \sum_{\ell,j} R_{\ell ij}(\hat{r}) \Delta t_{\ell} \int_{\phi_k \times E_j} F(\phi, E; \theta) dEd\phi \\ &= \sum_j \left[ \left( \sum_{\ell} \Delta t_{\ell} \right) \frac{\int_{\phi_k \times E_j} F(\phi, E; \theta) dEd\phi}{F_{jk}(\theta)} \right] \\ &\quad \times \sum_{\ell} \frac{\Delta t_{\ell}}{(\sum_{\ell} \Delta t_{\ell})} R_{\ell ij}(\hat{r}). \end{aligned} \quad (7)$$

A matrix of source count numbers may thus be evaluated as the dot-product$^{48}$

$$s(\theta) = \left( \sum_{\ell} \Delta t_{\ell} \right) R(\alpha, \beta, \gamma) \cdot F(\theta), \quad (8)$$

where $F(\theta)$ is a matrix of phase-integrated incident photon fluxes.$^{49}$ Note that the elements of the matrix $F(\theta)$ may be approximated using instantaneous fluxes at points within the finite phase intervals instead of explicitly integrating over those intervals, provided that the intervals are determined to be sufficiently small.

We numerically marginalize the likelihood function given by Equation (5) over the subspace of nuisance background parameters $\mathbf{B}$ in order to improve tractability of the sampling process. The target distribution (the posterior) for sampling is written conditional on model $\mathcal{M} \subset \mathcal{M}$ as

$$\pi(\theta | \mathbf{d}, \mathcal{M}) = \int \pi(\theta, \mathbf{B} | \mathbf{d}, \mathcal{M}) d\mathbf{B} \propto p(\theta | \mathcal{M}) \\ \times \frac{\int L(\theta, \mathbf{B}) p(\mathbf{B} | \mathcal{M}) d\mathbf{B}}{L(\theta)} \quad (9)$$

where $L(\theta) = p(\mathbf{d} | \theta, \mathcal{M})$ is the marginal likelihood function supplied as a callback for a sampling process, and the joint prior density distribution $p(\theta, \mathbf{B} | \mathcal{M})$ is separable with respect to $\theta$ and $\mathbf{B}$. The joint prior distribution $p(\mathbf{B} | \mathcal{M})$ is equivalent for all models: jointly flat and separable. Crucially, such a phenomenological background model exhibits a large prior complexity; for instance, the joint density at a background count rate vector $\mathbf{B}$, where the variation between channels is always small relative to the limiting instrument count rate, is equivalent to the joint density at a vector $\mathbf{B}$ whose elements exhibit vast channel-to-channel variations. It follows that a joint flat, separable prior is not considered representative of our prior belief. However, we consider—without proof—the prior to be weakly informative because (i) the conditional likelihood function (given a fixed source vector) exhibits a large curvature relative to the prior density function; (ii) the source photon flux signal always has few extrema in the joint space of energy and phase; and (iii) a posteriori the conditional likelihood function maxima do not wildly fluctuate as a function of channel because such structure does not exist in the data set. Nevertheless, if one were to compare, based on prior predictive performance, the models that we consider in this work with a model invoking a background component with far lower complexity, one should not be surprised if the former are strongly disfavored.

In general, the support of $p(\mathbf{B} | \mathcal{M})$ is compact and bounds can be specified on a channel-by-channel basis to truncate the marginalization integrals. Lower bounds may be derived, for example, from calibration observations of nearby fields that exclude the PSR J0030+0451 and are otherwise devoid of bright sources. Upper limits may, for example, be based on distinct NICER observations of the field containing PSR J0030+0451.

For this work, however, we define the lower bound as zero for each channel, and we eschew definition of an upper bound in each channel because the posterior is considered integrable: non-diverging on joint compact support, and the conditional likelihood function—$L(\theta, \mathbf{B})$ for fixed $\theta$—asymptotes to zero at large background count rates. If a set of sufficiently high upper bounds were specified (e.g., based on NICER count rate limits), the associated normalizing constant for the joint prior, equal to the reciprocal of the products of those bounds, would not modulate relative probability measures defined on $\mathcal{M}$.$^{50}$ We therefore do not view the improperness of the above prior as a misdemeanor, but it does mean that we should not describe our model as *generative* in the strictest sense.

<footnote>$^{48}$ An approximation in many ways, one being the discrete representation of an instrument that responds in a continuous manner to input.</footnote>

<footnote>$^{49}$ In the X-PSI implementation phase integration is performed using splines after a dot-product operation on a matrix of instantaneous incident photon fluxes.</footnote>

<footnote>$^{50}$ The raw event count rate during the exposures used to curate the data set for this work is known to be far below limiting and thus the conditional likelihood function is always relatively small for near-limiting background count rates. A set of upper bounds defined in this limit therefore truncates the evidence integral in a regime where model-dependent sensitivity is negligible.</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 3.** Schematic diagrams of models with single-temperature regions: ST-S defined by antipodal symmetry of the primary and secondary regions, and ST-U defined by the primary and secondary regions not sharing any parameters. Note that the two-surface shown is spherical: the hot regions are defined in spherical coordinates on a (unit) sphere and projected onto the rotationally distorted (oblate) spheroid representing the two-surface of the star, such that the proper areas of the regions are only implicitly defined. To construct the ST-U diagram from the ST-S diagram we displaced the secondary region (defined by having a center colatitude $Θ_s ≥ Θ_p$) away from the antipode of the center of the primary region (with center colatitude $Θ_p$), decreased $ζ_s$ such that $ζ_s < ζ_p$, and assigned the secondary region a distinct temperature parameter. Note that no numeric relationship is implied between the effective temperatures of material within each member—indeed, the prior is separable with respect to the temperatures ($T_p$ and $T_s$), and the support includes combinations $T_s ≤ T_p$ and $T_s > T_p$. For clarity we display for ST-S a projection showing the southern rotational hemisphere; in subsequent diagrams we omit such a projection when antipodal symmetry applies.

The numerical marginalization operation implemented is described in Appendix B of the X-PSI documentation (Riley & Watts 2019).

## 2.5. *Model-specific Definitions*

The surface heating distribution by realistic magnetospheric (return) currents remains uncertain. The global magnetic field may be more complex than a simple dipole (at least in the near vicinity of the surface), while the mapping between currents and surface temperature field is not well determined by existing theoretical models (Harding & Muslimov 2001, 2011; Timokhin & Arons 2013; Philippov et al. 2015a; Gralla et al. 2017;

Lockhart et al. 2019). We thus consider a set of simplified models that are representative of the various theoretical possibilities, albeit restricting our analysis to models with two distinct hot regions. We allow for the possibility of the hot regions being non-antipodal and non-identical (Pavlov & Zavlin 1997; Bogdanov et al. 2007, 2008; Bogdanov 2013); we also consider various hot-region shapes, including circles, rings, and crescents filled with material of uniform local comoving temperature.

While our choices are physically motivated, it is important to emphasize that our inferences are conditional upon these choices. However, posterior computation is computationally intensive and scales with model complexity. The work
presented here is based on usage of a 500,000 core-hour grant on the Dutch national supercomputer Cartesius; we thus find it pragmatic to disseminate information to the community at this point. Further exploration of model variants or execution of higher-resolution calculations requires additional resource allocation on high-performance systems, building upon the information offered here as guidance. In particular, it will be important to explore the sensitivity of marginal posterior estimates of fundamental physical parameters of interest—i.e., exterior and interior spacetime parameters—to expansion of the space of models that we have had the resources to consider.

In this section we detail the properties that distinguish the models in our model space. For each model we give the parameterization details and any remarkable prior details; notes on the *support* of the joint prior distribution are given where appropriate. We also discuss the existence of (continuous linear and nonlinear) degeneracy in posterior modes for a given parameterization, which we interpret as one indicator of unnecessary complexity for a surface radiation field with phenomenological spatial structure. For each model whose associated posterior distribution we compute, we provide a summary table containing a more precise definition of the joint prior distribution; these tables may be found in Section 3 and in Appendix C.

### 2.5.1. Single-temperature Regions with Antipodal Symmetry (ST-S) and with Unshared Parameters (ST-U)

**Parameterization.** The primary hot-region (refer to the leftmost panel of Figure 3) is simply connected and encloses radiating material—a fully ionized hydrogen NSX atmosphere with effective temperature $T$. The boundary of the region is circular: i.e., given a center point on the surface with colatitude $\Theta_p$ and azimuth $\phi_p$ (in a spherical coordinate basis whose polar axis is defined as the stellar rotation axis), the boundary is the locus of points that are equidistant<fnref n="51" /> in angular space from the center point.

Hereafter we use the alias ST-S, parsed as *Single-Temperature-Shared*. The surface radiation field associated with the secondary hot region is derived exactly by applying antipodal symmetry to the primary region: there are no free parameters associated with the secondary region.

Similarly, we use the alias ST-U, parsed as *Single-Temperature-Unshared*. The primary region (refer to the rightmost panel of Figure 3) definition is retained from ST-S as defined above. The secondary region, however, is now endowed with distinct parameters—i.e., the region is *not* derived from the primary region under antipodal symmetry. The parameters of the secondary region have an otherwise equivalent meaning—in terms of surface radiation field specification—to their primary-region counterparts.

**Degeneracy.** We note that a discrete degeneracy—multimodality—can in principle arise for a source such as PSR J0030+0451, but may only be weak when there is detectable asymmetry between the two component pulses over the course of one rotational cycle. There may exist two phase solutions, each corresponding to a distinct mapping between hot regions (distinguished by colatitude) and the pulse

components in the event data. For instance, the primary (lower-colatitude) region could in principle generate either of the component pulses, while the secondary region generates the other. Fortunately, a number of open-source sampling software packages are designed to handle multi-modality efficiently (at least in the absence of nonlinear degeneracy). If the asymmetry between the component pulses is clear, the posterior mass in one mode may be entirely dominant.

**Priors.** For ST-S, we eliminate a region-exchange degeneracy by imposing a constraint $\Theta_p \le \pi/2$ on the prior support. The primary region is uniquely defined as the region whose center subtends the smallest colatitude, $\Theta_p$, to the rotational axis, if the region colatitudes are different.

For ST-U, we eliminate a region-exchange degeneracy by imposing a constraint $\Theta_p \le \Theta_s$ on the support of the joint prior distribution. The primary region is uniquely defined as the region whose center subtends the smallest colatitude, $\Theta_p$, to the rotational axis, if the region colatitudes are different; the regions are distinguishable when $\Theta_p = \Theta_s$ according to the subset of parameters that controls their physical manifestation. The joint prior support is such that the two regions cannot overlap but otherwise are not restricted to be antipodally symmetric.

### 2.5.2. Concentric Single-temperature Regions with Antipodal Symmetry (CST-S) and with Unshared Parameters (CST-U)

**Parameterization.** The primary hot region (see the leftmost panel of Figure 4) is a non-simply connected annulus (or ring) with outer angular radius $\zeta$, which contains material with effective temperature $T$. The non-radiating hole with angular radius $\psi$ is concentric (in angular coordinates) with the radiating annulus. We thus recover the shape defined for the ST-S and ST-U variants in the limit $\psi \to 0$ (which is at the boundary of the prior support).

A generally useful way to distinguish the hole and the annulus—in particular for further increments in complexity—is as follows. Recall the term *member* from Section 2.4.1: consider two simply connected partially overlapping member regions, each wholly filled with radiating material, but impose the logical condition that when evaluating radiating intensities at a spacetime event on the stellar surface, one member—the hole—takes precedence if the event falls within its boundary. In this case (for the models here described), let the temperature of the material in the hole be (effectively) zero so that no signal need be computed for the hole. The statements in Section 2.4.1 pertaining to the proper areas in finite-element representations of radiating regions apply here: the annulus is a subset of a (simply connected circular) ceding member that is *not* *superseeded* by the hole when evaluating local radiation intensities, and its proper area is computed (almost) exactly. These constructions are useful for further extension of the model.

Hereafter we use the alias CST-S, parsed as *Concentric-Single-Temperature-Shared*. For CST-S, the surface radiation field associated with the secondary region is derived exactly by applying antipodal symmetry to the primary region: there are no free parameters associated with the secondary region. The annuli share an outer angular radius $\zeta$, and the holes share a fractional angular radius $f$ such that the hole angular radii are $\psi := f\zeta$.

Similarly, we use the alias CST-U, parsed as *Concentric-Single-Temperature-Unshared*. For CST-U, the primary region (refer to the rightmost panel of Figure 4) definition is retained

<footnote n="51">⁵¹ The ambient spacetime is static, and with respect to a Schwarzschild chart, the points are equidistant in angular coordinates. However, when projected from a spherical two-surface onto that of a rotationally deformed (oblate) spheroid, the spacelike separation—on a Schwarzschild temporal hyperslice— between the center point and boundary points is *not* invariant.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Schematic diagrams of models with concentric single-temperature regions: CST-S defined by antipodal symmetry of the primary and secondary regions, and CST-U defined by the primary and secondary regions not sharing any parameters.</figcaption></figure>

from CST-S as defined above. The secondary region, however, is now endowed with distinct parameters—i.e., it is *not* derived from the primary region under antipodal symmetry. The parameters of the secondary region have an otherwise equivalent meaning—in terms of surface radiation field specification—to their primary-region counterparts.

*Degeneracy*. We now consider the continuous degeneracy labeled I in Figure 5. When the angular extent of a radiating region is small, the signal generated by that region—as registered by a distant detector that does not spatially resolve (image) the star—is insensitive to its shape (detailed spatial structure). Sensitivity is here a measure in terms of the likelihood: i.e., the *total* variation of the parameterized joint sampling distribution of a set of random variables,<fnref n="52" /> in response to motion along a certain set of curves<fnref n="53" /> in parameter

space, is small as summarized by the scalar likelihood. It follows that if the shape is parameterized with more than a single degree of freedom, the likelihood function is degenerate with respect to shapes that satisfy a constraint on the solid angle subtended by the radiating region at the center of the star, and thus which satisfy a constraint on the proper area of the radiating region.

The degeneracy is not in this instance exact, but holds approximately. Consider a ceding member of angular extent $cos ζ ≲ 1$ and hole with angular extent $cos ψ > cos ζ$: the degeneracy is such that the solid angle, $ω(ζ, ψ) = cos ψ - cos ζ$, of the hot region (the annulus), is approximately $ω = 1 - cos ψ̃$, where $ζ = ψ̃$ is the angular radius in the limit $ψ → 0$. In other words, the relative size—and indeed existence—of the hole is at most weakly constrained on small angular scales. This degeneracy in the $(ζ, ψ)$-subspace is nonlinear for $ζ → ψ̃$, but linearizes for increasing $ζ > ψ̃$; however, with increasing $ζ$ the signal generated by the region evolves away from that generated in the limit $ζ → ψ̃$ with $ψ → 0$. Note that if the superseding member is not concentric with the ceding

<footnote n="52">⁵² And whose number usually exceeds the number of parameters.</footnote>

<footnote n="53">⁵³ Where those curves may more generally together generate *m*-dimensional surfaces in an *n*-dimensional space where usually *n* > *m*—if *m* = *n* then no facet of the model is constrainable and the model is arguably not useful unless it can meaningfully be tested in some other manner.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Schematic diagram of (approximate and exact) continuous degeneracies between hot-region structure parameters with respect to the energy-phase-resolved signal incident on a distant telescope. Degeneracy type I is discussed in detail in Section 2.5.2, while types II through IV are discussed in Section 2.5.3; here we focus on defining the components of the diagram. In all cases, let us suppose that there exists a strong posterior constraint on the (coordinate) solid angle $\varpi$ subtended by a component of a hot region (at the stellar origin) with temperature $T$. Type I—we consider a hot region with a single-temperature component, and represent such a region as simply connected with a strongly constrained radius $\zeta = \tilde{\psi}$. If the angular extent of the headed region is sufficiently small that we are insensitive (in the absence of spatial imaging power) to structural details other than solid angle $\varpi$, then we can introduce a hole of angular radius $\psi$ such that the hot region has the topology of a ring and a likelihood function degeneracy $\varpi = \cos\psi - \cos\zeta = \text{const.}$ manifests (on small angular scales). Similar arguments can be made for other parameterized structural modifications to the hot region. Type II—let us consider a hot region constructed from an annulus and a hole (concentric or eccentric), each filled with material of local comoving temperature $T$ and $T$, respectively. Let us assume a strong constraint on the outer angular radius $\zeta = \tilde{\psi}$ of the annulus, where $\tilde{\psi}$ is now not restricted to small angular scales; let us further assume a strong constraint on the structure of the hot region as simply connected circular with uniform temperature $T$, where only one temperature component is useful. It follows that degeneracy manifests for $T \approx T \approx \tilde{T}$ and $\psi + \psi^+ \approx \tilde{\psi}$, where $\psi^+ := \zeta - \psi$. Type III—as type I, but in the limit that $\zeta \gg \psi$, such that the signal generated by the hole, whose temperature is degenerate over a range dependent on $\psi$, is effectively turned off. Type IV—the signal generated by the hole now satisfies the constraints assumed for type II—i.e., $T \approx \tilde{T}$ and $\psi \approx \tilde{\psi}$—while the solid angle, temperature, and coordinates of the annulus exhibit degeneracy with respect to configurations that effectively turn off the signal generated by the annulus.</figcaption></figure>

member, and the overlap is only partial, the form of the degeneracy—the constraint equation satisfied—has additional dependence on the coordinates of the center of the hole relative to the center of the ceding member.

It is important to be aware of such degeneracy for the purpose of efficient posterior computation—inherent to which is accuracy. In our case, if the type of signal that is superior for describing the event data is generated by radiation from a

localized region on the star, the constraint on the solid angle of the region can be viewed as dragging a posterior mode through parameter space along the type I degeneracy direction. Thus, with this parameterization, nonlinearity will exist that will reduce efficiency to some degree.<fnref n="54" />

<footnote n="54">⁵⁴ Which can only be robustly learned during computation.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Schematic diagrams of models with concentric dual-temperature regions: CDT-S defined by antipodal symmetry of the primary and secondary regions, and CDT-U defined by the primary and secondary regions not sharing any parameters. To construct the CDT-U diagram from the CDT-S diagram we displaced the secondary region (defined by having a center colatitude $Θ_s \ge \Theta_p$) away from the antipode of the center of the primary region (with center colatitude $Θ_p$), decreased $ζ_s$ such that $ζ_s < ζ_p$ while leaving $ψ_s = ψ_p$, equivalently meaning $f_s > f_p$, and assigned distinct colors to the hole and annulus of the secondary region. We retained the caveat that no numeric meaning is implied regarding the effective temperature of the material in the members—indeed, the prior is separable with respect to all temperatures ($T_p, T_s, T_s, T_s$), and the support includes combinations $T_s \le T_s$ and $T_s > T_s$.</figcaption></figure>

If the posterior predictive performance is maximal for signals that are generated by localized emission, and degenerate posterior structure is observed, then clearly the most effective manner in which to achieve efficiency increase is to *reduce* the complexity of the structure of the radiating region. One thus inserts a simpler model into the model space, with the caveat that while estimation of ulterior model parameters should be insensitive to this reduction in complexity, the evidence may not be.

Alternatively, working with an integral summary variable such as the solid angle $\varpi$ of the (uniform temperature) radiating region is useful from the perspective of eliminating degeneracy by parameterizing directly in terms of variables to which we are statistically sensitive. On the other hand, the mapping from $\varpi$ to variables that directly control the shape of the radiating regions can behave undesirably, and can thus complicate the action of extending models.

One potential avenue for efficiency improvement by linearizing the degeneracy, is to sample in the space of ($\cos \zeta, \cos \psi$), which will eliminate the emergent small-scale nonlinear degeneracy. The cost is complication of the joint prior definition and implementation: in this case, a singularity<fnref n="55" /> exists in the mapping for $\psi \to 0$, which is the boundary at which the radiating region reduces to being simply connected as for both ST-S and ST-U. If we define finite joint prior density $p(\zeta, \psi)$ at points where $\psi = 0$, the joint density $p(\cos \zeta, \cos \psi)$ is divergent (although integrable, possibly in closed form).

*Priors.* An issue to be aware of is that the mapping from the native sampling space—that of the unit hypercube—would need to avoid introducing nonlinearity, otherwise the effort to

<footnote n="55">⁵⁵ Note that Taylor-expanding on small angular scales to work in the joint space of ($\zeta^2, \psi^2$) by definition cannot bypass the singularity.</footnote>
improve efficiency may be in vain (see the discussion in Section 2.4.1 regarding implementation of the joint prior distribution of *M* and *R*<sub>eq</sub>).

We define the prior as separable on the joint space of *f* and *ζ*. Specifically, we condition on *f* ~ *U(0, 1)* and *ζ* ~ *U(ε<sub>ζ</sub>, π/2 - ε<sub>ζ</sub>)* where ε<sub>ζ</sub> ≪ π/2. For posterior computation, in order to ease prior implementation, we opt to transform the prior onto the joint space of *ψ* and *ζ* and accept the nonlinear degeneracy in the limit *ζ* → *ψ* for small *ψ*. The chosen joint prior is not separable on this chosen space, but remains straightforwardly implementable.

For CDT-U, we eliminate a region-exchange degeneracy by imposing a constraint Θ<sub>p</sub> ≪ Θ<sub>s</sub> on the support of the joint prior distribution. The primary region is uniquely defined as the region whose center subtends the smallest colatitude, Θ<sub>p</sub>, to the rotational axis; the regions are distinguishable when Θ<sub>p</sub> = Θ<sub>s</sub> according to the subset of parameters that controls their physical manifestation. The joint prior support is such that the two regions cannot overlap but otherwise are not restricted by antipodal symmetry. Note that imposing that the regions are non-overlapping modifies the marginal prior density *p*(*ζ*)—and thus *p*(*ψ*)—by redistributing prior mass to lower angular radii.

### 2.5.3. Concentric Dual-temperature Regions with Antipodal Symmetry (CDT-S) and with Unshared Parameters (CDT-U)

**Parameterization.** We extend the CST models by filling the holes (the superseding members) with radiating material at finite temperature and evaluating the signal generated by both the superseding and ceding members. The hot region is simply the union of the members together with an order of precedence, and can be considered to have two heated subregions with distinct temperatures, each of which generates a component of the signal.

Hereafter we use the alias CDT-S, parsed as *Concentric-Dual-Temperature-Shared*. For CDT-S, the surface radiation field associated with the secondary region is derived exactly by applying antipodal symmetry to the primary region: there are no free parameters associated with the secondary region. The annuli share an outer angular radius *ζ*, and the holes share a fractional angular radius *f* such that the hole angular radii are given by *ψ* := *f*ζ.

Similarly, we use the alias CDT-U, parsed as *Concentric-Dual-Temperature-Unshared*. For CDT-U, the primary region (see the rightmost panel of Figure 6) definition is retained from CDT-S. The secondary region, however, is now endowed with distinct parameters—i.e., it is *not* derived from the primary region under antipodal symmetry. The parameters of the secondary region have an otherwise equivalent meaning—in terms of surface radiation field specification—to their primary-region counterparts. It follows that parameter vectors that correspond to, e.g., the primary region having a hotter annulus while the secondary region has a hotter hole, are assigned finite local joint prior density.

**Degeneracy.** Introducing a second temperature component also introduces degenerate structure in the likelihood function—a fundamental problem with such phenomenological descriptions of the surface radiation field. We now consider the continuous degeneracies labeled I, II, III, and IV in Figure 5.

Degeneracy of type I manifests on small angular scales when the temperature *T* of the superseding member is small relative to the temperature *T* of the ceding member. In this case, the

signal from the superseding member is dominated and we effectively recover the structure from the CST models; the reader should imagine the hole in panel I of Figure 5 as colored and labeled *T* ≪ *T*.

Degeneracy of type II manifests because the temperatures of the ceding and superseding members can be approximately and exactly equal. Thus, given a strong constraint on the temperature, the solid angle (π ≈ 1 - cos *ζ*) constraint on the radiating region can be satisfied ∀ψ ≪ *ζ*. Note that unlike degeneracy of type I, degeneracy of type II can be exact (for *T* = *T*) and exists invariantly of angular scale. The smaller the solid angle constraint, the weaker the constraint on topology of the hotter subregion (annulus or hole), and thus degeneracy of type II should be present. On larger angular scales, the signal offers a much stronger constraint on topology: if the hottest subregion is ring-like,<fnref n="56" /> degeneracy of type II cannot be present a posteriori.<fnref n="57" />

Degeneracy of type III manifests because the smaller the superseding member is in relation to the ceding member, the greater the dominance of the signal from the ceding member, and thus the greater the range of values *T* can assume while leaving the signal approximately invariant. Thus, given a strong constraint on the temperature and solid angle (π ≈ 1 - cos *ζ*) of the dominant component, the constraint can be satisfied ∀ψ ≪ *ζ*. Again note that unlike degeneracy of type I, the existence of degeneracy of type II is insensitive to angular scale. However, if there is a strong constraint on the topology of the subregion that dominates the joint signal, and that subregion is ring-like, degeneracy of type II will not be present a posteriori.

We now consider degeneracy IV and its implications for sampling. Suppose that the signal can be dominated by that from a simply connected subregion: if the superseding member satisfies both a temperature constraint and a solid angle constraint (π ≈ 1 - cos *ψ*), the ceding member then forms an annulus whose width *ψ*<sup>+</sup> and whose temperature *T* can together assume a wide range of values. In Figure 5 we illustrate the case of *T* ≪ *T* for varying *ζ* (and, in the last cartoon, for a shift of the center of the ceding member in angular space).<fnref n="58" /> The superseding member need only be described by four parameters, while the non-superseded subset of the ceding member is described by five parameters, three of which are shared with the superseding member. If the ceding member is relatively cool, then the signal it generates is dominated by that of a much hotter superseding member; the signal is then degenerate with respect to *T*. Crucially, if one writes the properties of the superseding member in terms of that of the ceding member, nonlinear degeneracy can arise in what otherwise may appear to be a natural parameterization, as we show below.

For posterior computation (see Appendix A for methodology) we explicitly consider a parameter space in which the superseding member is constructed using the minimal number

<footnote n="56">56 Formed by imposing a relatively large, cooler superseding member (effectively a hole).</footnote>

<footnote n="57">57 Here we assume a unimodal posterior distribution, but more generally, the degeneracy will not be present for a *local* mode in which the topology is strongly constrained as ring-like.</footnote>

<footnote n="58">58 Note that a less important degeneracy occurs, which we do not illustrate: the annulus is arbitrarily hotter (*T* > *T*) but *ψ*<sup>+</sup> → 0 so that the superseding member also dominates the signal. In the *T* ≈ *T* transition zone—which coincides with the type II degeneracy—*ζ* approximately satisfies the solid angle constraint.</footnote>
of parameters necessary, which eliminates needless nonlinear degeneracy. For example, if we work in the joint space of $f$ and $\zeta$, the local direction of approximate signal invariance—which emerges when the superseding member with angular radius $\psi$ is dominant—is not everywhere a basis vector, but is given by the gradient of $f = \psi/\zeta$ for constant $\psi \in \mathbb{R}^+$: $\partial f/\partial \zeta = -\psi/\zeta^2$, so the direction is $(1, -\psi/\zeta^2)$ at point $(\zeta, f)$. This degeneracy is linearized by working in the joint space of $\psi$ and $\zeta$, or alternatively, in the joint space of $\psi$ and the angular annular width $\psi^+ := \zeta - \psi$.

We note that unlike degeneracy of type I, the existence of degeneracy of type IV is insensitive to angular scale. However, if there exists a strong constraint on the topology of the heated subregion that dominates the signal, and that subregion is ring-like, degeneracy of type IV will not be present a posteriori—in the same vein that degeneracy of types II and III are absent. This does not mean that the posterior modes are devoid of degeneracy if the hot regions are both constrained to be ring-like: the temperature of a superseding member can assume a wide range of values (in logarithmic space) lower than that of the ring, because the additional complexity beyond that of a non-radiating hole (i.e., devoid of heated material) is unwarranted. However, the degeneracy is linear and thus handled straightforwardly.

In many cases, degeneracies are characterized by allowing one temperature component (a heated subregion) to dominate the signal, while one temperature—or the solid angle subtended at the stellar origin by a subregion—assumes relatively low values. In these cases the degeneracy is not expected to be particularly problematic because it will be effectively linear, parallel to basis vectors of the parameter space, and truncated by the boundary of the prior support. Thus, provided that the mapping from the native sampling space to the parameter space preserves this behavior, nested sampling efficiency may not be affected as adversely as when posterior modes exhibit strong nonlinear degeneracy.

For instance, consider degeneracy such as type IV: let the angular scale be sufficiently large for the dominant component to be constrained to be simply connected. The degenerate subsets of the $(\psi, \zeta)$- and $(\psi, T)$-subspaces will be linear and orthogonal to the $\hat{\psi}$ basis vector. The notable source of efficiency reduction for type IV will be the boundary of the degenerate subset of the $(\zeta, T)$-subspace: while there will be a hard lower limit on $T$, the boundary is otherwise dependent on combinations of parameters and will not conform trivially to common nested-sampling active-point bounding algorithms.

On the other hand, if degeneracy of type I arises and the hole is filled with material of temperature $T$, we incur a nonlinear degeneracy in the $(\zeta, \psi)$-subspace, and non-trivial boundaries of the degenerate subsets of the $(\zeta, T)$ and $(\psi, T)$-subspaces. In this case the potential for efficiency reduction is greater.

It is interesting to note that degeneracy can also arise when one temperature component does not dominate the signal. Let us suppose that two components with distinct temperatures are favored a posteriori, where those components both contribute non-negligibly to the total signal. We implement a numerical geometrically thin atmosphere where the effective temperature and effective gravity control the local comoving specific intensity as a function of photon energy and direction (i.e., control the spectrum and beaming); it follows that if a hot region exhibits small angular extent, two components (subregions) can in principle generate signals that are

commensurate in total count rate, provided that the solid angles subtended by the subregions at the stellar origin have appropriate relative sizes. For the CDT models, a discrete degeneracy—multi-modality—can then arise on small angular scales: pairs of configurations of the components generate approximately the same total signal. In one configuration the hole hosts a particular component, while in the alternate configuration that component is hosted by the annulus; the relative solid angles subtended by the hole and annulus depend on which component is hosted. However, for a more complex model where the superseding and ceding members are not defined as concentric (see Figure 5, and Section 2.5.5 and beyond), it is clear that (approximate) continuous degeneracies would also arise and thus complicate matters.

We should also be aware of the relative prior masses associated with subsets of parameter space over which the signal is approximately invariant, governed roughly by the dimensionality of the degeneracy. Degeneracy of types II and III generally occupy a subset of parameter space with smaller prior mass than that of type I and IV, the former two being effectively to one-dimensional,<fnref n="59" /> and the latter two being effectively two-dimensional. Note that the type I degeneracy is effectively two-dimensional *because* we fill the hole from the CST models with material whose temperature is finite and a model parameter.

*Priors.* We retain almost all prior definitions from the CST models (Section 2.5.2). We again transform the prior onto the joint space of $\psi$ and $\zeta$ in order ensure a more optimal parameterization for our sampling algorithm of choice, as highlighted above. Crucially, we need to transform from the native space to the parameter space in a manner that preserves linearity of degeneracy IV. We achieve this by inverse sampling the marginal prior density $p(\psi)$ and the condition prior density $p(\zeta | \psi)$; we give these transforms in Appendix B. In some cases we opt to transform to the joint space of $\psi$ and $\psi^+$ space in post-processing for optimality of kernel density estimation when $\psi^+$ approaches a boundary—e.g., due to degeneracy of type IV.

### 2.5.4. Interlude: On Degeneracy and Complexity

*General.* The continuous degeneracies discussed in Sections 2.5.2 and 2.5.3, when in mutual existence, generally form a connected structure in parameter space. These degeneracies could be characterized in more detail, and also with reference to ulterior parameters (including those of the other hot region). We do not study the degeneracies further here for lack of a clear way to simultaneously linearize each of them via transformation.

We instead observe a fundamental aspect of the modeling process *on a source-by-source basis*: it is judicious to evaluate the utility of additional complexity *on a region-by-region basis*, if a posteriori one infers localized emission and/or that a single-temperature component per hot region suffices in describing a particular pulse component in the phase-folded event data. In practice, this could include the following: (i) design a model space with increments in complexity, moving toward more sophisticated region topologies and/or boundaries, and/or an additional temperature component as warranted by the data; (ii) insert models into one's model space

<footnote n="59">⁵⁹ An unillustrated degeneracy, characterized by a hotter ring in the limit $\psi^+ \to 0$, is also effectively one-dimensional.</footnote>
<figure><img src="image_10.png" /></figure>

**Figure 7.** Conditional on CDT-U, we display points (nested samples) {$\theta_k$} that reported (background-marginalized) likelihood function values $L(\theta_k) > L'$, where $\max_k \ln L(\theta_k) - \ln L' \lesssim 13$; the choice of threshold $L'$ is somewhat arbitrary beyond needing a large enough number of points to clearly resolve the form of the degeneracies. We display these points in a selection of two-dimensional spaces to focus on the degeneracies that emerge upon application of CDT-U. These degeneracies pertain to the structure of the primary (lower-colatitude) hot region, which generates a signal describing a particular pulse component visible in Figure 1; the hot region exhibits a small angular extent on the stellar surface, and we are statistically insensitive to the additional complexity offered by a CDT cap over a ST cap. In the left panel we show the points in the joint space ($\zeta_p, \psi_p$)—i.e., in the joint space of the angular radii of the ceding and superseding members. We indicate which type of degeneracy illustrated in Figure 5 corresponds to which structure in parameter space; the curves (blue and red) each correspond to a constraint equation. Types II and III occupy far smaller prior masses than types I and IV, so we choose to indicate the structures associated with the former types with ellipses; in a finite-sample context these structures are allocated points more sparsely, in proportion to the associated prior mass in the full $n$-dimensional parameter space. Note that the region $\psi_p > \zeta_p$ is not a subset of the CDT-U prior support. In the center panel we display the points in the joint space of ($\zeta_p, f_p$), where $f_p := \psi_p/\zeta_p$; our sampling processes (refer to Appendix A for details) were executed in a space in which likelihood function isosurfaces are structured nonlinearly as suggested in this center panel, leading to reduced posterior computation efficiency and higher nested-sampling error. In Sections 2.5.2 and 2.5.4 we discuss the interpretation and treatment of degeneracy: in hindsight, higher efficiency would have been achieved if we had worked in a joint space that linearized degeneracy of type IV. In the right panel we display the corresponding structures in a joint space trivially related to ($T_p, T_s$)—i.e., in the joint space of the effective temperatures of the material enclosed by the ceding and superseding members. The points displayed here are a composite set from two sampling processes, where due to reduced efficiency and increased nested-sampling error, neither process individually resolved all of the degenerate structure present.

that reduce region complexity, if it is apparent a posteriori that a model with unhelpful complexity was first applied; or (iii) reduce region complexity as in point (ii), while increasing the complexity of another facet of the model, based for example on graphical posterior predictive checking.

To provide an example, suppose that a posteriori, given CDT-U, hot region $\mathcal{A}$ is degenerate with respect to multiple parameters controlling its surface radiation field, while hot-region $\mathcal{B}$ is only degenerate with respect to temperature of one relatively cool component. Specifically, suppose $\mathcal{A}$ is constrained to have a small angular extent and a single-temperature component, while offering satisfactory posterior predictive performance in comparison to a particular pulse component. Continuous degeneracies of type I and IV are then clear in the posterior distribution. Additionally, suppose that the secondary region has a strong constraint on its topology, with its dominant component being ring-like with moderate angular extent. On the basis of this inference, one might argue that it is justified to modify the model space: one proceeds to define a model wherein the complexities of $\mathcal{A}$ and $\mathcal{B}$ are unequal. In doing so one can reduce or redistribute complexity while improving computational efficiency. For the aforementioned example, one may define a model ST+CST, parsed as *Single-Temperature + Concentric-Single-Temperature*, which is intermediary in relation to ST-U and CST-U: $\mathcal{A}$ has one simply connected component, while $\mathcal{B}$ has one non-simply connected component (which can reduce to a simply connected component). The dimensionality of the parameterization is then $n=9$, reduced from the CDT-U value of $n=12$. One could also define additional models that increase the complexity of $\mathcal{B}$ (see Sections 2.5.5 and 2.5.6).

As a brief second example: if, given CDT-U, $\mathcal{B}$ was constrained to have two temperature components, one could instead define a model ST+CDT with dimensionality $n=10$. One could then proceed to increment the complexity of $\mathcal{B}$.

An important remark here is that the performance of $\mathcal{A}$ cannot strictly be decoupled from the performance of $\mathcal{B}$ in relation to their respective target pulse components: the regions exist on the same rotating two-surface and share an ambient spacetime solution. It follows that the complexity of one region can in principle affect our conclusions about the level of complexity of the other. Such conclusions cannot be based solely on the existence of degenerate structure in the joint posterior distribution of the region parameters. If $\mathcal{A}$ exhibits degenerate structure while $\mathcal{B}$ does not, the structure may manifest—at least in part—*because* the complexity of $\mathcal{B}$ is insufficient.

A more robust basis for breaking parity between the hot-region complexities would be the existence of degenerate structure, and satisfactory posterior predictive performance of both regions in comparison to their respective pulse components. A practical basis for breaking parity, given the existence of efficiency-reducing degenerate structure a posteriori, may be resource availability and management.

*Modeling PSR J0030+0451.* For the source we focus on in this work, degeneracy emerged upon application of CDT-U. We display the corresponding structures, in parameter space, in Figures 7 and 8. In summary, the primary region was constrained to have a small angular extent and one dominant temperature component, and could generate data structurally consistent with a pulse component visible in the phase-folded event data (Figure 1); the primary region thus here corresponds
<figure><img src="image_11.png" /><figcaption>**Figure 8.** As for Figure 7, conditional on CDT-U, we display points (nested samples) {$\theta_k$} that reported (background-marginalized) likelihood function values $L(\theta_k) > L'$, where $\max_k \ln L(\theta_k) - \ln L' \lesssim 13$. We display these points in a selection of two-dimensional spaces to focus on the secondary (higher-colatitude) hot region. The secondary region exhibits a larger angular extent on the stellar surface than the primary region and, in contrast, we are statistically sensitive to the topology of the secondary region while being statistically insensitive to the additional complexity offered by a CDT cap over a CST cap. In other words, we do not strongly require the material in the hole to radiate when considering the NICER waveband.</figcaption></figure>

to region $\mathcal{A}$ above. Degeneracy of types I and IV arose and dominated in prior mass (and by extension, in posterior mass). The nonlinear degeneracies in the native nested-sampling space also suppressed the efficiency of the sampling processes. We therefore now have grounds to design and apply models that break the parity in complexity between hot regions. We redistribute complexity from the primary region to the secondary region, such that the dimensionality of the parameter space does not exceed that of CDT-U. In Appendix D we define models that preserve parity between regions—a brief continuation of the above scheme for model extension.

In Sections 2.5.5 and 2.5.6 we explain how the complexity of a solitary hot region is extended. In Section 2.5.7 we briefly summarize and provide diagrams of the relevant models formed from two disjoint hot regions with unequal complexities. In Section 2.5.7 we then provide diagrams summarizing the relationships between all models considered in the scope of this work.

### 2.5.5. Eccentric Single- and Dual-temperature Regions (EST and EDT)

**Parameterization.** We extend the CST and CDT hot-region types by not requiring the ceding and superseding members to be concentric. For a CST region, the superseding member was a hole in a ceding member, devoid of radiating material; for a CDT region, the superseding member was similarly a hole, but filled with radiating material. For both CST and CDT regions, the *non-superseded* subset of the ceding member was a radiating annulus with concentric inner- and outer-boundaries; now the radiating annulus is *eccentric*. Further we require—via the joint prior support—that the ceding member is strictly a superset of the superseding member; in other words, the superseding member is a hole in the ceding member, again forming an annulus or ring.<fnref n="60" />

We illustrate these hot-region structures in the topmost panels of Figure 9; we also alluded to such configurations in Figure 5 and in Sections 2.5.2 and 2.5.3. Hereafter we use the aliases EST and EDT, respectively parsed as *Eccentric-Single-Temperature* and *Eccentric-Dual-Temperature*.

An additional two parameters must be defined to specify the eccentricity as the great-circle segment separating the centers of the members: the magnitude and direction of the offset. To do so, we consider a spherical coordinate basis such that: (i) the polar axis subtends an angle $\Theta$ to the stellar rotation axis; (ii) the southern rotational pole subtends angle $\pi - \Theta$ to the polar axis, and has azimuth zero (such that the northern rotational pole has azimuth $\pi$ radians). The azimuth of the annulus center, $\phi$, in the rotated basis uniquely selects the great circle and direction of the offset (depicted in Figure 9). We choose define the magnitude of the angular separation as a fraction $\epsilon$ of the difference $\psi^+ := \zeta - \psi \equiv \zeta(1 - f)$ between the angular radii of the hole and annulus (see Section 2.5.3). The magnitude of the angular offset between the centers—equal to the colatitude of the annulus center in the rotated spherical coordinate basis defined above—is thus $\vartheta := \epsilon\psi^+$.<fnref n="61" />

**Degeneracy.** In a similar vein to the CDT parameterization, optimality of the parameterization for an EDT region should be considered. By this point in the modeling process, we have a handle on a region-by-region basis whether the additional complexity offered by an EST or EDT region is a justified modeling step. Nevertheless, one can construct the sampling space using the minimum number of parameters need to fully control the superseding member. The parameter vector for the regions is thus of the form $\mathbf{v} = (\psi, \Theta, \phi, \zeta, \varepsilon, \varphi)$, and the mapping from the native sampling space to the parameter space is formed in part using the marginal density $p(\psi)$ and the conditional density $p(\zeta | \psi)$.

<footnote n="60">⁶⁰ The hole reduces to a point at the boundary of the prior support for an EST region, while for a EDT region the hole is not permitted to reduce to a point. Also note that because the ceding and superseding members are circular, their boundaries may intersect at a maximum of one point.</footnote>

<footnote n="61">⁶¹ An alternative is to parameterize in terms of the fraction of the sum of their angular radii such that $\vartheta := \epsilon(1 + f)\zeta$; it is this combination that is the limit of zero partial overlap when the superseding member lies partially exterior to the ceding member (see Section 2.5.6). In order to impose that the superseding member is a hole in the ceding member, we would then require the support constraint that $\epsilon(1 + f)\zeta + f\zeta \le \zeta$ and thus that $\epsilon \le (1 - f)/(1 + f)$.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Schematic diagrams of solitary hot regions beyond the CST and CDT complexity levels, as described in Sections 2.5.5 and 2.5.6. We delineate the boundary of the superseding member for the PST region.</figcaption></figure>

*Priors.* As for the CDT type of region (Sections 2.5.3), joint prior support is such that the material within the annulus may exhibit an effective temperature $T < T_\text{eff}$ or $T \ge T_\text{eff}$. Moreover, when an EST/EDT region shares the stellar surface with another hot region C, the joint prior support does not include configurations wherein C overlaps with the ceding (annular) region of the EST/EDT region. Note that a continuous set of coordinate singularities exist in the mapping $(\epsilon, \phi) \mapsto (\Theta', \phi')$, where the coordinates $(\Theta', \phi')$ are the colatitude and azimuth of the center of the annulus in a spherical coordinate system with polar axis defined as the stellar rotation axis. The mapping is singular for all points $(\epsilon, \phi) = (0, \phi)$ at the boundary of the prior support, such that by defining a finite joint prior density at such singular points the joint density at $(\Theta', \phi') = (\Theta, \phi)$ diverges. We consider the spherical coordinates $(\vartheta, \phi)$ to be more natural for increasing complexity and thus we display

density functions with respect to $(\epsilon, \phi)$ to avoid difficulty in representing density functions accurately.

### 2.5.6. Protruding Single- and Dual-temperature Regions (PST and PDT)

*Parameterization.* We extend the EST and EDT hot-region types by not requiring the ceding member to be a strict superset of the superseding member. In other words, the superseding member need not be a hole in the ceding member, leading to radiating component with the topology of a ring. If the superseding member is *not* a hole nor a point, and the boundaries of the ceding and superseding members intersect at two points, then the *non-superseded* subset of the ceding member is simply connected *and* has a non-circular boundary: it cannot therefore always be considered an annulus, and can
assume a crescent- or arc-like morphology.<fnref n="62" /> Note that we do not increase complexity here by incrementing the dimensionality of parameter spaces defined for the EST and EDT regions. We nevertheless choose to distinguish between EST/EDT and PST/PDT in view of a symmetry being broken: that of the circularity of the outer boundary of the region. The superseding member can also now subtend the larger (coordinate) solid angle at the stellar origin, and in the limit that the member boundaries touch at a point, the hot region is simply connected, with a boundary that is either circular or union of two circular subregion boundaries.

We illustrate these hot-region structures in the bottommost panels of Figure 9; we also alluded to such configurations in Sections 2.5.2 and 2.5.3. Hereafter we use the aliases PST and PDT, respectively parsed as *Protruding-Single-Temperature* and *Protruding-Dual-Temperature*: the superseding member can protrude from the ceding member, for parameter vectors within the joint prior support.

We modify the EST (and thus EDT) parameterization in order to permit the superseding member to: (i) protrude<fnref n="63" /> from the ceding member; and (ii) subtend a larger (coordinate) solid angle at the stellar origin than the ceding member. Let us denote the angular radius of the largest member as $\xi$. Let us then consider the interval $f \in [\epsilon_f, 2 - \epsilon_f]$, where $\epsilon_f$ is some small number: when $f \le 1$ the angular radii of the superseding and ceding members are, respectively, $\psi = f\xi$ and $\zeta = \xi$, as in the EDT variants. However, when $f > 1$, the angular radii of the superseding and ceding members are, respectively, $\psi = \xi$ and $\zeta = (2 - f)\xi$. The angular radii are thus piecewise in $f$, and are continuous at the transition point $f = 1$. It follows that the interval for the angular radius $\xi$ of the largest member considered in the previous variants can be maintained, while the member that is largest switches. If one varies $f$ through the interval $f \in [\epsilon_f, 2 - \epsilon_f]$ while all other parameters are fixed, the solid angle subtended by the superseding member increases to match that of the ceding member; at the $f = 1$ transition the superseding member stops expanding, and for $f$ increasing beyond unity, the solid angle subtended by the ceding member decreases.

We choose to define the magnitude of the angular separation as a fraction<fnref n="64" /> $\kappa$ of a combination of the angular radii of the superseding and ceding members. For $f \le 1$, the combination is the *sum* of the angular radii of the superseding and ceding members, $\psi = f\xi$ and $\zeta = \xi$, respectively, for $f \le 1$: $\vartheta(\xi, f \le 1, \kappa) := \kappa(1 + f)\xi = \kappa(\psi + \zeta)$. It is this combination that is the limit of zero partial overlap between the members when the superseding member lies partially exterior to the ceding member. We define the prior support for $\kappa$ simply as the unit interval: we permit only configurations wherein the members partially overlap (including at a point), with maximum separation $\vartheta = \psi + \zeta$.

We require the angular separation $\vartheta$ as piecewise in $f$: for $f > 1$, we let $\vartheta(\xi, f > 1, \kappa) := \psi - \zeta + 2\kappa\zeta = (f - 1)\xi + 2\kappa(2 - f)\xi$. Note that the minimum separation here with respect to $\kappa$, for fixed $\xi$ and fixed $f$, is $\vartheta = \psi - \zeta$, equivalent to

the transition from the ceding member being *partially superseded* to being *wholly superseded*; the minimum separation for the subinterval $f \le 1$ is zero. The maximum separation (with respect to $\kappa$ for fixed $\xi$ and fixed $f$) is $\vartheta = \psi + \zeta$, which is equivalent to the expression for the maximum separation for $f \le 1$. Note that at the transition $f = 1$, the piecewise components of $\vartheta(\xi, f, \kappa)$ continuously match at a value of $2\kappa\zeta$.

The coordinates of the center of the ceding member remain written in terms of those of the center of the superseding member. We cannot define a pair of coordinates as being associated with the larger member (in the same vein that $\xi$ is directly associated with the larger constituent) without generally introducing a discontinuous transition with respect to the configuration of the members at $f = 1$.

Note that the above prescription is equivalent in principle to defining a binary discrete parameter (with associated uniform prior probability mass function) that controls which member $\xi$ in turn directly controls, and thus which member $f$ directly controls, where the upper bound of the prior support of $f$ remains as unity. In this equivalent alternative, the piecewise definitions of $\vartheta$ and the angular radii are required, but are recast with respect to the binary parameter. However, extension of the support of a continuous parameter to include these configurations eases posterior computation.

**Degeneracy.** The two clear exact degeneracies beyond those considered above<fnref n="65" /> are: (i) equal temperatures where the superseding member is wholly enclosed by the ceding member, leading to three-dimensional degeneracy in the subspace of $(f, \kappa, \varphi)$; and (ii) the superseding member *wholly supersedes* the ceding member ($\kappa = 0$ for $f \ge 1$), leading to a two-dimensional degeneracy in the $(f, \varphi)$-subspace. There are also discrete member-exchange degeneracies for a given region when the temperatures are equivalent.

**Priors.** The joint prior support is again such that the material within the ceding member may exhibit an effective temperature $T < T_0$ or $T \ge T_0$. Moreover, when a PST/PDT region shares the stellar surface with another hot region $C$, the joint prior support could be defined to exclude configurations wherein $C$ overlaps with the PST/PDT region. We note that due to choice of parameterization and prior support, the prior density function for the angular separation between the ceding and superseding members has changed from the corresponding EST/EDT density function.

### 2.5.7. Overview of the Model Space

In Figure 10 we provide a diagram of the relationships between the models—spanning the discrete model space—that impose equal complexities. A subset of these models are based on the extensions defined in Sections 2.5.5 and 2.5.6; diagrams of EDT-S, EDT-U, PDT-S, and PDT-U are given in Appendix D because we do not compute posterior distributions that are jointly conditional on these models and the event data.

We also define models wherein an ST region shares the stellar surface with a higher-complexity hot region. We provide diagrams of ST+EST and ST+PST in Figure 11. We provide diagrams of ST+EDT and ST+PDT in the online figure set associated with Figure 11. In Figure 12 we illustrate the relationships between these models. Note that these models

<footnote n="62">⁶² It is now particularly clear why we supplant the *hole* and *annulus* descriptors with those in terms of simply connected *members* and evaluation precedence as noted in Section 2.4.1: one member is distinguished from the other simply by which takes precedence when evaluating radiation intensities along a ray (null geodesic) connecting the stellar surface to a distant observer.</footnote>

<footnote n="63">⁶³ Or partially overlap with the ceding member.</footnote>

<footnote n="64">⁶⁴ Equal to the colatitude of the center of the ceding member in the rotated spherical coordinate basis defined above.</footnote>

<footnote n="65">⁶⁵ The degeneracies considered above should at this stage be of little to no concern if computation of the simpler (nested) models has *not* indicated that the complexity thus far introduced is unhelpful.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Diagram of relationship between a subset of models {$M$}, where each $M \subset \mathcal{M}$ and $\mathcal{M}$ is the model space as defined in Section 2.3. The models are ST-S and ST-U (Section 2.5.1), CST-S and CST-U (Section 2.5.2), CDT-S and CDT-U (Section 2.5.3), EDT-S and EDT-U (Appendix D), and PDT-S and PDT-U (Appendix D). The integers in parentheses are the number of continuous parameters controlling the hot regions; to obtain the total number of parameters constituting the sampling space for each model, add 8 to these numbers—corresponding to the vector $(M, R_{\text{eq}}, i, D, N_{\text{H}}, \alpha, \beta, \gamma)$ shared between all models. Prior predictive complexity increases with dimensionality and/or prior support expansion; complexity increments are achieved by breaking symmetries, allowing the morphology and topology of the radiating regions to change, and adding a second component with parameterized temperature. Solid arrows between model nodes delineate a nested relationship between models: the model at the tail is nested within the model at the head, the latter of which has greater complexity. Adding arrows head-to-tail at a node conserves such relationships—i.e., the model at the tail of the resultant arrow is nested within the more complex model at the head. The dashed arrows are in the background (visualize a third dimension of the graph as indicated by ⊗). The dotted arrows between single- and dual-temperature models indicate that the nested relationship is weaker: in the single-temperature models no material (or material with zero temperature) fills the superseding member and thus no signal is generated, whereas for the dual-temperature models we fill the superseding member with material of finite temperature $T \gtrsim O(10^5)$ K, and a signal is thus physically generated and computed. The blue ring indicates the simplest model that can generate data that is visually comparable to the real data set (refer to Section 3 for discussion). The red annotations are to denote where we learn that a given model is not performing adequately or has unhelpful complexity. A strikethrough denotes that a model was not applied: (i) a solid strikethrough denotes that the model is considered as being incapable of generating synthetic event data that resembles the real XTI event data, based on performance of a simpler model with which some degree of symmetry is shared; (ii) a dotted strikethrough denotes that the model is considered to include unhelpful complexity for at least one hot region.</figcaption></figure>

with unequal complexities exist within the scope of those models present in Figure 10.

When the hot regions are distinguishable by complexity level we can label the regions as *primary* and *secondary* according to that complexity: let the higher-complexity region be *secondary*. Our application of CDT-U suggested that evolving the lower-colatitude (primary) region from ST to CDT did not yield any improvement in describing the event data; the lower-colatitude region in this case describes the pulse component in the event data (Figure 1) that peaks at ~0.5 rotational cycles. Moreover, for both ST+CST and ST+CDT, the secondary region is formed (not necessarily by trivial union) from two members—one or both of which radiate—that are concentric and thus share a single colatitude parameter. Therefore, for ST+CST and ST+CDT it is reasonable to retain the constraint on the prior support that $\Theta_p \le \Theta_s$—i.e., that the colatitude of the *center* of the ST region is at most the colatitude of the *center* of the CST/CDT region.<fnref n="66" />

Notably, we cannot consider a model with an EST/EDT or PST/PDT region as nested within CDT-U, and there is thus not an obvious ordering of regions in colatitude that can be applied based on the CDT-U posterior information. We therefore do not impose such a constraint on the prior support, meaning that the order of the regions in model names ST+EST, ST+EDT, ST+PST, and ST+PDT does *not* indicate an order in colatitude.

Based on the CDT-U posterior information we are now interested in models wherein we couple the lower-complexity (primary) region to the pulse component in the event data (Figure 1) that peaks at ~0.5 rotational cycles. We consider a practical reason to identify a particular hot region with a particular pulse component in the phase-folded event data: to ease posterior computation by focusing sampling resolution in the vicinity of such configurations. If we do not opt for this coupling, a local posterior mode (or modes) can absorb sampling resolution, despite contributing relatively low posterior mass and thus being ultimately uninteresting for both parameter estimation and evidence estimation; such a local mode is characterized by the ST region describing the pulse component at approximately zero rotational cycles (Figure 1).

<footnote n="66">⁶⁶ A model wherein the CST or CDT region is at *higher* colatitude can be continuously connected to this model by simply by forgoing the constraint on the prior support that $\Theta_p \le \Theta_s$.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Schematic diagrams of models wherein an ST region shares the stellar surface with a higher-complexity EST or PST region.</figcaption></figure>

Mitigating such an effect may then require additional active points and/or activation of the mode-separation MULTINEST sampling variant (refer to Appendix A).

To intentionally couple the hot regions to pulse components, we restrict the prior support for the ST region to an interval $\phi_p \in [a, b]$ where $(b-a) < 1$ (such that the boundary of the support is not periodic with respect to $\phi_p$). Remarkably, while such restriction may not affect parameter estimation,<fnref n="67" /> the evidence does clearly depend on changes to the prior support, especially for weakly informative priors whose density does *not* fall to negligible values at the changing boundary of the support. It follows that in order to compare models based on

evidences, we should obtain a lower bound on the evidence of, e.g., ST+EST, by accounting for the increased prior mass in the posterior mode due to contraction of the ST phase support.

We reserve the remaining prior implementation details for Appendix B.2.

### 3. Inferences

In this section we provide (posterior) summary information about each model applied and compare them. Tables giving numerical information for all but one model are available in Appendix C. Numerical files associated with the nested sample sets may be found in the persistent repository of Riley (2019). As regards model comparison, there is no clear maximally optimal measure for relative model performance (and certainly not for performance in an absolute sense). We deem ST+PST

<footnote n="67">⁶⁷ Provided that (i) the dominant posterior mode corresponds to a ST primary region, and (ii) that the marginal prior density $p(\phi_p)$ remains weakly informative relative to the likelihood function</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Diagram of relationship between dual-hot-region models with unequal complexities, constructed via extension of ST-U. These models exist within the scope of those displayed in Figure 10.</figcaption></figure>

to be superior considering the following mixture of measures, both quantitative and qualitative.

**Posterior predictive performance.** Our crude graphical posterior predictive checking procedure suggests that even ST-U can generate synthetic event data that is structurally commensurate with the XTI event data—at least in the channel subset [25, 300]. For reference, see Appendix A.2.1, together with Figure 13 and the associated figure set. Absent are obvious systematic differences in the (Poisson) standardized residuals over phase-channel intervals. Nevertheless, one could study the residual differences for modeling background event arrival processes, instrument operation, or noise properties of the event data—e.g., whether the event arrival processes could detectably deviate from being Poissonian (see also Figure 23). For the purpose of constructing a model that can simply generate event data that is similar to PSR J0030+0451 under visual inspection, ST-U is competitive with all higher-complexity models.

**Visualization.** In Figures 14 and 15 (and the associated online figure sets) we display, in posterior-expected form, various signals (derived quantities) generated by the hot regions; the reader can thus get a handle on the source contribution to the model displayed in Figure 13 (center panel). In Figure 16 we display the model-to-model evolution of the posterior-expected count-rate signals generated. To aid visualization, we provide in Figures 17 and 18 schematics of the hot regions on the surface that generate the aforementioned signals.

**Evidence.** We report evidence estimates in Tables 2 through 6. For the family of models that we have considered, the evidence strongly favors (adopting the guidelines of Kass & Raftery 1995) a model wherein antipodally reflection symmetry

is *not* imposed; in other words, ST-S, and by extrapolation any such model with the -S extension (refer to Figure 10), is considered to be strongly disfavored. Moreover, for all computed models that do not impose antipodal reflection symmetry, posterior modes contain configurations characterized by both hot regions being located in the same rotational hemisphere—the opposite rotational hemisphere to the Earth direction. We also conclude that at least one region should be modeled with more complexity than offered by a single-temperature simply connected circular (ST) region; in other words, there is deemed sufficient evidence to disfavor ST-U. However, there is insufficient evidence to resolve between CDT-U, ST+CST, ST+EST, and ST+PST, especially when considering the evidence not as a scalar estimator, but as a random variable with a (simulated) distribution for estimating error intervals. The estimated error intervals (defined to contain 90% of evidence estimates based on nested sampling process realizations) typically overlap for these models, while the expected log-evidences are within ~2 units.

**Kullback-Leibler divergence.** We report global and marginal divergence estimates in Tables 2 through 6, and in Figures 19–21. Information gain from prior to posterior is defined as a divergence integral over some subset of parameter dimensions: it is a non-negative real scalar measure of the difference between normalized density functions (see Appendix A.2.4 for a more detailed description). The larger the number, the larger the information gain, while a minimum divergence of zero indicates that the density functions are identical. The global information gain—the divergence integral over all dimensions—is comparable for all models. Based on comparison of parameter-by-parameter divergence estimates, and visual comparison of the marginal prior and posterior density functions hence summarized, we consider the joint prior distribution to be weakly informative in the context of the likelihood function for most source parameters; the main exception is the distance, which is strongly prior-dominated. The other parameters whose marginal posterior distributions are entirely prior-dominated are the two of the instrument parameters, α and γ (refer to Figures 21 and 22). These informative prior distributions are however shared by all models, as are the improper flat prior density functions described in the background treatment (refer to Section 2.4.3). Note that more information is gained about the instrument parameter β—the weighting factor between two response matrices—and we can conclude that the instrument manifests a posteriori as a mixture weighted appreciably toward the nominal response for all posterior computations; we reserve discussion on why the instrument calibration may not be accurate for Section 4.1.3.

**Likelihood function.** We generally refrain from reporting parameter vector point measures such as the maximum a posteriori vector or the maximum likelihood vector for the purpose of formal quantitative model comparison (see Appendix A.2.1 for reasoning). We note, however, that the background-marginalized likelihood function exhibits the largest values in any posterior typical set (by ~4 natural logarithmic units), across all models, for ST+PST. It follows that a subset of the additional complexity introduced by ST+PST is helpful—and suggestive of avenues for future modeling—while a subset of the complexity is unhelpful, leading to only commensurate evidence relative to ST+EST. In other words, the (background-marginalized) likelihood
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Count data {$d_{ij}$}, posterior-expected count numbers {$\lambda_{ij}$}, and (Poisson) residuals for ST+PST. Note that we split the count numbers in the upper two panels over two rotational cycles, such that the information on phase interval $\phi \in [0, 1]$ is identical to the information on $\phi \in [0, 2]$; our data sampling distribution, however, is defined as the (conditional) joint probability of all event data grouped into phase intervals on $\phi \in [0, 1]$. We display the standardized (Poisson) residuals in the bottom panel: the residuals for the rotational cycle $\phi \in [0, 1]$ were calculated in terms of all event data on that interval (as for likelihood definition), and simply cloned onto the interval $\phi \in (1, 2]$. In Appendix A.2.1 we elaborate on the information displayed here. The complete figure set (6 images) is available in the online journal, for the ST+PST, ST+EST, ST+CST, ST-U, ST-S, and CDT-U models.</figcaption></figure>

(The complete figure set (6 images) is available.)

function, over a subset of the additional prior support, is at most commensurate with the prior expectation of the likelihood function for the other models; globally, the model does not therefore improve as measured solely by the evidence.

*Tractability and complexity.* Posterior computation accuracy and reproducibility are generally higher as the models become less complex. The ST+CST, ST+EST, and ST+PST models redistributed complexity based on the performance of CDT-U in pursuit of *helpful* complexity and thus higher-efficiency resource consumption. We demonstrate in Figure 16 that effectively all of the improvement of CDT-U over ST-U is captured by the intermediary ST+CST. As a consequence,

ST+PST in particular exhibited substantial improvement in regard to the (background-marginalized) likelihood function in the posterior mode, while reducing the number of parameters by limiting the complexity of one hot region; it follows that it should be clearly preferred over CDT-U.

We remark that although the complexity of one hot region was ultimately limited at the ST-level, this does not mean that the signal generated by that region is fixed or performs maximally in some absolute sense. Indeed, the signal generated by one region cannot be decoupled (as stated in Section 2.5.4) from the signal generated by the other for the purpose of statistical inference: (i) the regions are restricted to exist on the
<figure><img src="image_17.png" /></figure>

**Figure 14.** The posterior-expected signal for ST+PST, both incident on the instrument (top and top-middle) and as registered by the instrument (bottom-middle and bottom). The signal in the top panel has been integrated over the linearly spaced instrument energy intervals, and is effectively proportional to the specific photon flux. The black count-rate curves are the posterior-expected signals generated by each hot region separately, and in combination. We also represent the conditional posterior distribution of the incident photon flux (top-middle) and count-rate (bottom) at each phase as a set of one-dimensional highest-density credible intervals, and connect these intervals over phase via the contours; these distributions are denoted by $\pi$ (photons cm$^{-2}$ s$^{-1}$; $\phi$) and $\pi$ (counts s$^{-1}$; $\phi$). Note that the fractional width of the credible interval at each phase is usually higher for $\pi$ (photons cm$^{-2}$ s$^{-1}$; $\phi$) than for $\pi$ (counts s$^{-1}$; $\phi$) because of the variation permitted for the instrument model; in combination, the signal registered by the instrument is more tightly constrained. To generate the conditional posterior bands we apply the X-PSI package, which in turn wraps the fgivenx (Handley 2018) package. The complete figure set (6 images) is available in the online journal.

(The complete figure set (6 images) is available.)

same two-surface within an ambient spacetime; and (ii) for some subset of rotational phases, and for the configurations that report higher (background-marginalized) likelihood values, the images of both regions are simultaneously visible (a posteriori)

and registered by the instrument in combination. As a result, upon examination of Figure 16, we observe that the posterior-expected count-rate signal generated by the ST region evolves along the ST-U, ST+CST, ST+EST, and ST+PST sequence,
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The posterior-expected spectrum for ST+PST, both incident on the instrument (top) and as registered by the instrument (middle and bottom). The black count-rate curves are the posterior-expected spectra generated by each hot region separately, and in combination. In the top panel we display as black curves the incident photon specific flux spectra both with and without interstellar absorption. We represent the conditional posterior distribution $\pi$ (photons keV$^{-1}$ cm$^{-2}$ s$^{-1}$; E) of the absorbed incident photon specific flux at each energy as a set of one-dimensional highest-density credible intervals, and connect these intervals over phase via the contours (top); the energies displayed are those spanning the waveband of channel subset [25, 300). The credible intervals fan-out at the lowest energies because (i) conditional on event data for channel subset [25, 300), we are relatively insensitive to the details of the signal for $E \in [0.1, 0.2]$; and (ii) the interstellar attenuation factor is stronger the lower the photon energy, and thus the incident signal varies strongly as a function of the neutral hydrogen column density—a free parameter that operates as an exponent. In the top panel we overlay the incident absorbed spectrum inferred by Bogdanov & Grindlay (2009) based on a phase-averaged analysis of low-background XMM observations; the blue band denotes the estimated fitting uncertainty on this spectrum, at each energy, as a Gaussian with fractional standard deviation $\sigma/\mu = 0.15$, and the three opacity levels indicate intervals $1\sigma$ through $3\sigma$. There is additional systematic XMM flux calibration uncertainty at the ~10%-level that is not included here. In the middle panel we display the background-marginalized posterior-expectation of the source count-rate signal, plus the background count-rate terms that maximize the conditional likelihood function; the signal is equivalent to that displayed in the middle panel of Figure 13. In the bottom panel we display the posterior-expected count-rate spectra generated by the hot regions in combination and individually; we opt not to render the conditional posterior count-rate distribution for each channel because it is too narrow about the expected spectrum to be useful. Moreover, the topmost black step function is the phase-average of the center panel—it is effectively, but not exactly, the observed count-number spectrum divided by the total exposure time $T_{exp}$. We combine the XMM-derived count-rate spectrum (and its associated uncertainty) with the marginal NICER instrument posterior on parameters $\alpha, \beta$, and $\gamma$ to simulate a conditional probability distribution $\pi$ (counts s$^{-1}$; i) for the count-rate in the ith channel; these conditional distributions are connected via the contours in a manner congruent to the top panel described above. The complete figure set (6 images) is available in the online journal.</figcaption></figure>

(The complete figure set (6 images) is available.)

both at phases where the non-ST region is invisible, and at phases where both regions are simultaneously visible. The reader should not therefore interpret our decision to limit the

complexity of one region to ST as a statement that an ST region is forecasted to be sufficient for the purpose of future modeling efforts nor that exploration of different models will
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Top panel: the posterior-expected channel-by-channel source count rate signal generated as a function of rotational phase by ST+PST surface emission: $s := E_{\pi(\theta)}[g(\theta)]$, where $g(\theta)$ is a map from parameter space to the corresponding count rate signal. Others panels: the difference between signal $s$ (for ST+PST) and the posterior-expected signal conditional on some other model: let $s' := E_{\pi(\theta)}[h(\theta)]$, where $h(\theta)$ is a map from parameter space to a count-rate signal conditional on some model other than ST+PST. The signals $s$ and $s'$ represent Poissonian arrival processes; we therefore opt to display the difference as $(s - s')/\sqrt{s}$, and scale the signals by $T_{\text{exp}}/32$. A grayscale value is then locally representative of the absolute change (shown exclusively on the phase interval $\phi \in [1, 2]$), in units of the Poisson standard deviation, within a phase interval of width 1/32 cycles as defined during event data pre-processing.</figcaption></figure>
Northern rotational hemisphere
 view( Earth inclina tion)

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Top panel: schematic diagram of a surface heating configuration—and Earth inclination—representative of those corresponding to points in the ST+PST posterior mode. The configuration rendered here corresponds to the sample that reported the highest background-marginalized likelihood function value across all models (among the values reported by the set of all nested samples). We project the hot regions onto a (unit) sphere and view from the Earth inclination with no ambient gravitational field. The regions are constrained to exist in the same hemisphere, but with remarkably different morphologies. The hot regions are approximately equal in effective temperature and thus we define a new temperature symbol *T*<sub>eff</sub> that is common to both. We also display the channel-summed count rate pulse generated by the source emission and indicate which region generates which component; we refer the reader to Figure 14, where this signal is also displayed, for more information. Bottom panel: note that we impose (via the prior support) that the Earth inclination lies within the northern rotational hemisphere, but an identical configuration (in terms of the physics that we consider and thus signal generation) is given via an equatorial reflection of both the Earth direction and radiating regions; we render this alternative configuration as viewed from the equatorial plane, and display the Earth inclination (but not azimuth) as the shaded angular interval bounded by the 16% and 84% quantiles in marginal posterior mass (see Figure 21 and Table 2 for the numerical interval).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Surface heating configurations conditional on ST+ST (equivalent to ST-U), ST+CST, and ST+EST. For each model the hot regions correspond to the sample assigned the greatest posterior weight and is thus representative, being a draw from the posterior typical set. For each model we display the Earth inclination (but not azimuth) as an angular interval bounded by the 16% and 84% quantiles in marginal posterior mass.</figcaption></figure>

fail to yield improvement; indeed, one could even consider
models such as PST-U given that our modeling route for one
region was based on CDT-U.

*Telescope cross-checking.* In Figure 15 we display the posterior information for the phase-integrated spectrum generated by the ST +PST hot regions; corresponding figures for the other models may
<table><caption>**Table 2**
Summary Table for ST+PST</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>C&#773;I<sub>68%</sub></th><th>D&#773;<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87,<sup>a</sup> fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass<sup>b</sup></td><td>M ~ U(1, 3)<sup>c</sup></td><td>1.34<sup>+0.15</sup><sub>-0.16</sub></td><td>1.26<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius<sup>d</sup></td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1), 16]<sup>e</sup></td><td>12.71<sup>+1.14</sup><sub>-1.19</sub></td><td>0.78<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td></td><td>Compactness condition</td><td>R<sub>eq</sub>/r<sub>g</sub>(M) > 3</td><td></td><td></td></tr><tr><td></td><td>Compactness condition<sup>f</sup></td><td>R<sub>polar</sub>(M, R<sub>eq</sub>, Ω) ≥ 3r<sub>g</sub>(M)</td><td></td><td></td></tr><tr><td></td><td>Enforce elliptical two-surface cross-section</td><td>function of (M, R<sub>eq</sub>, Ω)</td><td></td><td></td></tr><tr><td>Θ<sub>p</sub> (radians)</td><td>ST region center colatitude<sup>g</sup></td><td>Θ<sub>p</sub> ~ U(0, π)</td><td>2.23<sup>+0.10</sup><sub>-0.10</sub></td><td>2.98<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>φ<sub>p</sub> (cycles)</td><td>ST region initial phase (from Earth)</td><td>φ<sub>p</sub> ~ U(a, a + 0.2)<sup>h</sup></td><td>0.46<sup>+0.00</sup><sub>-0.00</sub></td><td>4.44<sup>+0.03</sup><sub>-0.04</sub></td></tr><tr><td>ζ<sub>p</sub> (radians)</td><td>ST region angular radius</td><td>ζ<sub>p</sub> ~ U(0, π/2)</td><td>0.09<sup>+0.01</sup><sub>-0.01</sub></td><td>4.62<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>PST region superseding center colatitude</td><td>Θ<sub>s</sub> ~ U(0, π)</td><td>2.91<sup>+0.05</sup><sub>-0.05</sub></td><td>3.85<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>PST region initial phase (from Earth antipode)</td><td>φ<sub>s</sub> ~ U(-0.5, 0.5), periodic<sup>i</sup></td><td>-0.59<sup>+0.04</sup><sub>-0.04</sub></td><td>2.43<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>ψ<sub>s</sub><sup>+</sup> (radians)</td><td>PST region angular radii difference<sup>j</sup></td><td>ξ<sub>s</sub> ~ U(0, π/2), ψ<sub>s</sub><sup>+</sup> = ψ<sub>s</sub><sup>+</sup>(ξ<sub>s</sub>, f<sub>s</sub>)<sup>k</sup></td><td>-0.01<sup>+0.02</sup><sub>-0.02</sub></td><td>3.16<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>ψ<sub>s</sub><sup>-</sup> (radians)</td><td>PST region superseding angular radius</td><td>f<sub>s</sub> ~ U(0, 2), ψ<sub>s</sub> = ψ<sub>s</sub>(ξ<sub>s</sub>, f<sub>s</sub>)</td><td>0.47<sup>+0.13</sup><sub>-0.13</sub></td><td>1.52<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>χ<sub>s</sub></td><td>PST region ceding fractional angular offset</td><td>χ<sub>s</sub> ~ U(0, 1)</td><td>0.10<sup>+0.05</sup><sub>-0.03</sub></td><td>2.67<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>φ<sub>s</sub> (radians)</td><td>PST region ceding azimuthal offset</td><td>φ<sub>s</sub> ~ U(0, 2π)</td><td>0.71<sup>+0.10</sup><sub>-0.10</sub></td><td>2.28<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>log<sub>10</sub>(T<sub>p</sub> (K))</td><td>Non-overlapping hot regions<sup>l</sup></td><td>function of Θ<sub>p</sub> through φ<sub>s</sub> above</td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>ST region NSX effective temperature</td><td>log<sub>10</sub>(T<sub>p</sub> ~ U(5.1, 6.8), NSX limits)</td><td>6.11<sup>+0.01</sup><sub>-0.01</sub></td><td>5.29<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>i (radians)</td><td>PST region NSX effective temperature</td><td>log<sub>10</sub>(T<sub>s</sub> ~ U(5.1, 6.8), NSX limits)</td><td>6.11<sup>+0.01</sup><sub>-0.01</sub></td><td>5.75<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>D (kpc)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>0.94<sup>+0.11</sup><sub>-0.10</sub></td><td>1.91<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>N<sub>H</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)<sup>m</sup></td><td>0.33<sup>+0.01</sup><sub>-0.01</sub></td><td>0.07<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>NICER α</td><td>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>0.78<sup>+0.17</sup><sub>-0.16</sub></td><td>2.86<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>NICER β</td><td>Calibrated matrix scaling</td><td>α ~ N(1, 0.1), α ∈ [0.5, 1.5]</td><td>0.99<sup>+0.08</sup><sub>-0.09</sub></td><td>0.04<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>NICER γ</td><td>Reference-to-calibrated matrix weighting</td><td>β ~ U(0, 1)</td><td>0.19<sup>+0.17</sup><sub>-0.12</sub></td><td>0.96<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td></td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>0.98<sup>+0.08</sup><sub>-0.08</sub></td><td>0.08<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td colspan="5"><strong>Sampling process information<sup>n</sup></strong></td></tr><tr><td colspan="5">Number of free parameters:<sup>o</sup> 19</td></tr><tr><td colspan="5">Number of runs:<sup>o</sup> 1</td></tr><tr><td colspan="5">Number of live points: 10<sup>3</sup></td></tr><tr><td colspan="5">Inverse hypervolume expansion factor:<sup>q</sup> 0.8</td></tr><tr><td colspan="5">Termination condition: 10<sup>-1</sup></td></tr><tr><td colspan="5">Evidence:<sup>f</sup> ln Z&#773; = -36368.28<sup>+0.49</sup><sub>-0.46</sub></td></tr><tr><td colspan="5">Global KL divergence: D&#773;<sub>KL</sub> = 68.9<sup>+0.9</sup><sub>-0.8</sub> bits</td></tr><tr><td colspan="5">Number of core<sup>g</sup> hours: 42453</td></tr><tr><td colspan="5">Likelihood evaluations: 78343018</td></tr><tr><td colspan="5">Nested replacements: 57972</td></tr><tr><td colspan="5">Weighted posterior samples:<sup>i</sup> 20177</td></tr></tbody></table>

**Notes.** We provide (i) the parameters that constitute the sampling space, with symbols, units, and short descriptions; (ii) any notable derived or fixed parameters; (iii) the joint prior distribution, including hard truncation bounds and constraint equations that define the hyperboundary of the support; (iv) one-dimensional (marginal) 68.3% credible interval estimates symmetric in posterior mass about the median (C̄I<sub>68%</sub>); and (v) KL divergence estimates in bits (D̄<sub>KL</sub>) representing prior-to-posterior information gain (a scalar with an associated error calculated as the 68.3% distributional interval about the median with respect to sampling process realizations; see Appendix A.2.4 for high-level description of the divergence). Constraint equations in terms of two or more parameters result in marginal distributions that are not equivalent to those inverse-sampled.

<footnote><sup>a</sup> Arzoumanian et al. (2018).</footnote>

<footnote><sup>b</sup> Interpreted as a rotationally perturbed mass monopole moment (e.g., Hartle 1967), but the perturbation is small for the spin frequency of PSR J0030+0451 (see Section 4.2).</footnote>

<footnote><sup>c</sup> Hard lower bound based loosely on plausible astrophysical formation channels (see, e.g., Strobel et al. 1999).</footnote>

<footnote><sup>d</sup> An alternative two-dimensional space of g(M, R<sub>eq</sub>) and h(M, R<sub>eq</sub>)—where g ≠ h—on which we could choose to specify a joint flat prior density distribution is that of g := M and h := R<sub>eq</sub>/r<sub>g</sub>(M).</footnote>

<footnote><sup>e</sup> The function r<sub>g</sub>(M) denotes the gravitational radius explicitly in dimensions of length.</footnote>

<footnote><sup>f</sup> The coordinate polar radius of the source two-surface, R<sub>polar</sub>(M, R<sub>eq</sub>, Ω), is a quasi-universal function adopted from AlGendy & Morsink (2014), where Ω := 2π/P is the coordinate angular rotation frequency. This compactness condition, together with the elliptical surface requirement below are unimportant for the P ~ 5 ms spin period of PSR J0030+0451, and therefore we ignore these constraints in Tables 3 through 7.</footnote>

<footnote><sup>g</sup> Note that for parameters where a lower bound of zero would correspond to absence of pulsations, we use some small finite number as a lower bound.</footnote>

<footnote><sup>h</sup> Where φ<sub>p</sub> = a is an arbitrary phase dependent on event data pre-processing. We set a = 0.35.</footnote>

<footnote><sup>i</sup> The periodic boundary is admitted and handled by MULTINEST. However, this is an unnecessary measure because we straightforwardly define the mapping from the native sampling space to the space of φ<sub>s</sub> such that the likelihood function maxima are not in the vicinity of this boundary.</footnote>

<footnote><sup>j</sup> The difference is defined as ψ<sub>s</sub><sup>+</sup> := ζ<sub>s</sub> - ψ<sub>s</sub><sup>-</sup>, where ζ<sub>s</sub> is the angular radius of the ceding member that radiates where it is not superseded.</footnote>

<footnote><sup>k</sup> See Section 2.5.6.</footnote>

<footnote><sup>l</sup> Refer to Appendix B.2.4.</footnote>

<footnote><sup>m</sup> Constructed to approximate the information in the measurement (statistical information) reported by Arzoumanian et al. (2018). The support D ∈ [0.235, 0.415] is equivalent to the interval μ ± 10σ.</footnote>

<footnote><sup>n</sup> Refer to Appendix A for definitions.</footnote>

<footnote><sup>o</sup> In the sampling space; the number of background count-rate parameters is equal to the number of channels defined by the data set.</footnote>

<footnote><sup>p</sup> The mode-separation MULTINEST variant was deactivated, meaning that isolated modes are not evolved independently and nested sampling threads contact multiple modes.</footnote>

<footnote><sup>q</sup> For this sampling process, nor any such process reported in this work, we did not activate constant-efficiency MULTINEST active-point bounding variant.</footnote>

<footnote><sup>r</sup> Defined as the prior predictive probability p(d | ST + PST). We report the interval about the median containing ±45% of 10<sup>3</sup> joint bootstrap-weight replications for the combined run. Note, however, that in order to complete the reported evidence for comparison to models other than those defined in this work, upper bounds for the background parameters need to be specified as described in Section 2.4.3.</footnote>

<footnote><sup>s</sup> Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) processors on the SURFsara Cartesian supercomputer. Note that these are physical cores—i.e., hyper-threading technology is not invoked.</footnote>

<footnote><sup>t</sup> Excludes samples with important weight smaller than 10<sup>-6</sup> times the largest such weight among samples.</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** One- and two-dimensional marginal posterior density distributions for the MSP spacetime parameters conditional on the each in the sequence ST-U, ST+CST, ST+EST, and ST+PST. From leftmost to rightmost in each panel, the parameters are equatorial radius, (equatorial) compactness, and gravitational mass. We display the marginal prior density distributions for each parameter as the dashed-dotted functions. For the less expensive ST-U and ST+CST models we executed two production runs and combined them, while for each of ST+EST and ST+PST we executed a single production run. We also display an ST-U run in which the mode-separation MULTINEST variant was activated, but because neither the theory nor the software exists for combining such a sampling process with runs in default mode (where sampling threads can migrate between posterior modes), it is not included in the combined run. The mode-separation run allocated ~one-third of the sampling resolution (i.e., sampling threads, or active points) to a second posterior mode with negligible local mass; this mode corresponds to a distinct phase configuration, with the hot regions—which are distinguished by their order in colatitude—transposed in their coupling to the pulse components visible in the phase-folded event data. We report the KL divergence, $D_{KL}$, from prior to posterior in bits for each parameter, together with an error interval containing 68.3% of $\hat{D}_{KL}$ estimates based on simulated nested sampling process realizations. The shaded credible intervals CI$_{68\%}$ for each parameter are symmetric in marginal posterior mass about the median, containing 68.3% of the mass; the (barely discernible) darker intervals at the CI$_{68\%}$ boundaries contain 68.3%, respectively, of the 15.85% and 84.15% quantiles in posterior mass, again based on simulated nested sampling process realizations. The credible regions in the off-diagonal panels, on the other hand, are uniquely the highest-density—and thus the smallest possible—credible regions, containing 68.3%, 95.4%, and 99.7% of the posterior mass. In Appendix A we provide additional information regarding posterior kernel density estimation, error analysis, and the estimators displayed here. The complete figure set (7 images) is available in the online journal.</figcaption></figure>

(The complete figure set (7 images) is available.)
be found in the online figure set. We overlay a model PSR J0030+0451 spectrum, derived via analysis of low-background low-signal-to-noise ratio phase-integrated *XMM* observations (see Bogdanov & Grindlay 2009); this spectrum is considered as a guiding upper bound for all emission in the *NICER* waveband from the PSR J0030+0451 system. The ST-U, CDT-U, ST+CST, ST+EST, and ST+PST models do not violate this condition in channels [25, 300] spanned by the event data considered, nor at incident photon energies that couple strongly to this channel subset; the ST-S model clearly violates this condition. The *XMM*-derived model is more uncertain than shown at low energies in Figure 15, due to unresolved likelihood function degeneracy between $N_H$ and the thermal components. The *XMM* spectral analysis also included a power-law component, which is subsumed into the phase-invariant (background) terms in our analysis of *NICER* data (see the topmost black count-rate step function in the bottom panel of Figure 15); the discrepancy at higher energies is thus accounted for. Nevertheless, a pertinent question remains regarding how much of the signal generated by surface emission is captured by the phase-invariant count-rate terms; these terms combine linearly with the signals from the hot regions and are intended to parameterize the background contribution. We refer the reader to Section 4.1 for more detailed discussion on this topic.

**Marginalization.** For parameters that are shared between discrete models {$M_m$}, where $M_m \subset \mathcal{M}$, we could in principle marginalize over the discrete parameter *m* (see Appendix C of Riley & Watts 2019, for a formulation consistent with Section 2.3), provided we accept a marginal prior mass distribution of *m*. If two or more models are competitive, the marginal joint posterior distribution of the shared parameters is not dominated by the information from a single model. As stated above, we cannot distinguish between four of the highest-evidence models (ST +PST, ST+EST, ST+CST, and CDT-U). Moreover, for ST+CST, ST+EST, and ST+PST, the marginal joint posterior distribution of the shared spacetime parameters of interest is only mildly sensitive to model choice—i.e., the compactness constraints are commensurate, while the mass and radius are only weakly constrained individually, and the joint credible regions exhibit high partial overlap. For ST+PST, the marginal joint posterior distribution of the shared spacetime parameters discernibly evolves, but the mass and radius remain weakly constrained individually. Given that *m* labels models that are often (approximately) nested and which differ only in phenomenological complexity, we opt not to marginalize shared parameters over those models; instead, we report headline parameter estimates for ST+PST, which exhibits the largest background-marginalized likelihood function values in any posterior typical set.

# 4. Discussion

In this section we highlight how our inferences may be sensitive to the modeling assumptions made, and discuss the implications of our inferences for both dense matter physics and NS astrophysics.

## 4.1. Modeling Assumptions

The inferences that we report are conditional upon a number of modeling assumptions. These assumptions were physically motivated, and reassuringly there are no obvious large discrepancies or structures in the pulse-profile residuals to indicate a major problem. Nonetheless, the sensitivity of our results to these assumptions should be explored in future work, given additional

computational resource allocations. With regard to the MSP, the biggest model-dependencies are (i) the atmosphere (see Section 4.1.1); (ii) the treatment of phase-invariant components of the total signal (see Section 4.1.2); and (iii) the assumption there exist two disjoint hot regions, each of which is radiatively contiguous and has a temperature field that is adequately represented by one of the models described in Section 2 (see 4.1.2). The instrument model (see Section 4.1.3) is of less concern.

### 4.1.1. Atmosphere

Two properties that could affect the atmosphere models used in our analyses are chemical composition (hydrogen as opposed to helium) and ionization state (fully versus partially ionized). In this Letter we have considered only a hydrogen composition: hydrogen would dominate the composition of matter accreted from the interstellar medium (Blaes et al. 1992), while matter accreted from a binary companion star would be predominantly hydrogen or helium. There are other processes that may drive changes in composition. For instance: hydrogen would result from spallation (Bildsten et al. 1992), diffusive nuclear burning could convert hydrogen to helium (Chang & Bildsten 2003, 2004), and significant pulsar wind excavation could make visible an underlying heavy element layer (Chang & Bildsten 2004). If the atmosphere were in fact dominated by helium, we could expect changes because helium atmospheres radiate differently from those of hydrogen. For example, hydrogen and helium model specific intensities, at an atmosphere effective temperature of $10^6$ K (the approximate inferred temperature for the hot regions for this source), have fractional differences of at most 2%–5% at 0.5–1 keV in the (maximal) forward direction (Bogdanov et al. 2019c).

Regarding ionization state, our atmosphere models are constructed assuming the atmospheric plasma is fully ionized, such that the dominant opacity in regimes of interest is that due to electron free-free absorption (Ho & Lai 2001). While opacity tables for partially ionized matter exist (Iglesias & Rogers 1996; Badnell et al. 2005; Colgan et al. 2016), they do not cover the full range of energies and temperatures needed for our analysis. However, the hydrogen neutral fraction in the atmosphere at $10^6$ K is low, and a comparison of our fully ionized hydrogen atmosphere model with that constructed using the OP (Badnell et al. 2005) opacity table yields specific intensity fractional differences of at most 1%–2% at 0.5–1 keV in the (maximal) forward direction at an atmosphere effective temperature of $10^6$ K (Bogdanov et al. 2019c). The importance of including partial ionization, in comparison to developing other aspects of the model (e.g., Section 4.1.2), is not clear; but partially ionized models will be part of future re-analysis as updated opacity tables become available.

### 4.1.2. Surface Heating Estimation and Treatment of Phase-invariant Components

Predicting the MSP surface temperature field from ab initio$^{68}$ calculations of energy deposition by magnetospheric currents is, as described in Section 2, notoriously challenging. In this work we have assumed that there are only two distinct hot regions, motivated by the fact that there appear to be two pulsed components in the pulse profile. If the (surface) field

<footnote>$^{68}$ Terminology adopted from the series of studies by Philippov et al. (2015a, 2015b) and Philippov & Spitkovsky (2018) on pulsar magnetospheric simulations.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** One- and two-dimensional marginal posterior density distributions for the MSP parameters conditional on ST+PST. From leftmost to rightmost: equatorial radius $R_{eq}$; gravitational radius $M$; ST center colatitude $\Theta_p$; ST angular radius $\zeta_p$; ST NSX effective temperature $\log_{10} T_p$; PST superseding member center colatitude $\Theta_s$; PST superseding angular radius $\psi_s$; PST superseding member angular radius difference $\psi_s^\dagger = \zeta_s - \psi_s$; where $\zeta_s$ is the angular radius of the ceding member; PST ceding member fractional angular offset $\kappa_s$ (labeled as $\kappa_s$); PST ceding member azimuthal offset $\varphi_s$; and PST ceding NSX effective temperature $\log_{10} T_s$. For descriptions of the information displayed, refer to Figure 19; note that here we display the marginal posterior density distribution for each parameter as a single solid function due to the number of panels. We choose not to display joint posterior distributions for all pairs of model parameters because the number of panels is prohibitive; moreover, the posterior azimuthal separation of the ST and PST regions is displayed in Figure 21. The complete figure set (6 images) is available in the online journal.</figcaption></figure>

(The complete figure set (6 images) is available.)

structure involves higher-order multipoles,<fnref n="69" /> additional polar caps—and thus additional disjoint surface heating—may be

possible, and/or more complicated polar cap topologies (such as ring-like) may be possible. We have also assumed specific forms for each hot region; the true temperature field is more complex, and both the physical complexity and our statistical sensitivity to such complexity should be investigated further. In lieu of a physical emission model for the stellar surface exterior

<footnote n="69">⁶⁹ With respect to some coordinate system that simplifies vector spherical harmonic field expansion; this coordinate system will generally be rotated and displaced from the system with stellar spin axis defined as the polar axis.</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure 21.** One- and two-dimensional marginal posterior density distributions for parameters pertaining mostly to the MSP observation, conditional on ST+PST. From leftmost to rightmost: distance *D*; Earth inclination *i*; ST center azimuth (phase) relative to Earth direction *φ*<sub>*p*</sub>; PST superseding member center azimuth (phase) relative to Earth antipode *φ*<sub>*s*</sub>; interstellar neutral hydrogen column density *N*<sub>H</sub>; and NICER instrument parameters α, β, and γ. For descriptions of the information displayed, refer to Figure 19; note that here we display the marginal posterior density distribution for each parameter as a single solid function due to the number of panels. The complete figure set (6 images) is available in the online journal.</figcaption></figure>

(The complete figure set (6 images) is available.)

of the hot regions, we subsumed the non-pulsed component of
any such emission (which is expected a priori to be dominant if
hot regions with smaller angular extent are favored) within
phase-invariant count-rate terms; however, this would be a

minor concern on the premise that outside of the footpoints of
open magnetic field lines at the polar caps, there is no
energetically comparable heating to which we are sensitive
when observing with NICER.
<figure><img src="image_25.png" /><figcaption>**Figure 22.** We update the instrument prior (displayed in Figure 2) conditional on ST+PST. The conditional posterior distributions are represented by the orange bands; the conditional prior distributions are displayed in blue to indicate information gain about the instrument, and the black lines and points are equivalent to those in Figure 2. The complete figure set (6 images) is available in the online journal.</figcaption></figure>

(The complete figure set (6 images) is available.)

All emission from sources other than the hot regions—a combination of astrophysical and instrumental—was left free in our models in the form of a set of phase-invariant count-rate terms (background parameters), one per channel, which we collected under the envelope of *background contribution* (refer to Section 2.4.3). Moreover, an improper joint flat prior was implemented that was separable with respect to these background parameters; no upper bounds (or lower bounds) were defined for the prior support, in lieu of a physical (generative) model for the total contribution from the hot regions (and thus surface if emission exterior of the regions is considered unimportant). Such a model would need to account for the combined (phase-invariant) signal attributed to off-surface emission and any astrophysical backgrounds in the field of view, in the NICER waveband, based on previous observations of PSR J0030+0451 and/or theoretical modeling.

Remarkably, by neglecting any physical (generative) model for the total counts attributed to *all* surface emission, the phase-invariant terms can even capture emission from the hot regions. One can reason that by permitting the phase-invariant terms to capture all or most of the phase-invariant signal components, there may then exist *background-marginalized* likelihood function maxima corresponding to signals that are (i) dominated by phase-invariant terms over a fraction of a rotational cycle, and (ii) elsewhere found to describe the pulse-profile adequately (in combination with the phase-invariant terms). Near pulse minima, the hot region contribution can be entirely dominated in linear combination with the phase-invariant terms, if the regions are both close to the visible limb of the star, or even partially or wholly non-visible. Such a heating configuration may exhibit systematic bias in the sense that (i) it is not considered an adequate approximation of the configuration inferred when physical limits are modeled for the contribution from the stellar surface (or specifically the hot regions), or (ii) that a configuration alluded to in point (i) is not

encompassed by a given posterior credible region boundary. On the other hand, it could be viewed that our treatment of the phase-invariant component of the total signal is in some aspects conservative: posterior credible regions may be appreciably larger than if such physical limits on surface contribution are imposed. In the *X-PSI* documentation (Riley & Watts 2019) a simple parameter estimation workflow is demonstrated using the same default background treatment implemented in this present work; while there is no evidence for systematic bias nor credible region inflation in that specific case (where the true data-generating process is known), guarantees cannot be made universally.

Upon examination of Figure 14 we see that the combined signal from the hot regions falls to near zero at its minimum—such that the fractional amplitude of the signal is near unity—and is thus *not inconsistent* with a phase-invariant component of the combined signal being subsumed in the background. Further, examination of Figure 15, in which we cross-check the ST+PST region spectrum with a PSR J0030+0451 spectrum inferred by Bogdanov & Grindlay (2009), further suggests that a fraction of the contribution from the hot regions is subsumed in the background. Interestingly, Bogdanov & Grindlay (2009) inferred a spectrum with two thermal components of different temperatures.<fnref n="70" />

We attribute the discrepancy to a number of factors. First, the *XMM* photon event set had size $O(10^4)$, nearly two orders of magnitude smaller than the number of NICER events used in this work. Second, it is known that Bogdanov & Grindlay (2009) did not fully resolve degeneracy between the thermal components and the neutral hydrogen column density $N_H$. Lastly, *XMM* is an imaging telescope: the background signal

<footnote n="70">⁷⁰ Note that Bogdanov & Grindlay (2009) explicitly calculate the phase-resolved signal generated by rotating single-temperature circular hot spots, and then phase-average; it is thus not the case that the discrepancy is explained by *our* rotating circular hot spots (refer to ST-U), whose effective temperatures are commensurate, mimicking a dual-temperature *incident* spectrum due to relativistic rotation.</footnote>
<figure><img src="image_26.png" /></figure>

**Figure 23.** Summary of the two-dimensional residual panels in the online figure set associated with Figure 13. For each model we identify the sample (parameter vector) $\hat{\theta}_{\text{BMML}}$ that reported the highest background-marginalized likelihood function value among identified posterior modes, together with the background count-rate vector $\hat{\mathbf{B}}_{\text{CML}}$ that maximizes the conditional likelihood function given $\hat{\theta}_{\text{BMML}}$. Given the expectations $E[x_{ij}] = \lambda_{ij}(\hat{\theta}_{\text{BMML}}, \hat{\mathbf{B}}_{\text{CML}})$, we evaluate the residuals $(\lambda - d)/\sqrt{\lambda}$ for each phase-channel interval, where $\{d_{ij}\}$ are variates of independent Poisson-random variables $\{x_{ij}\}$. We display in each panel the distribution of these random variates (the real count numbers). The minimum number of expected counts in any phase-channel interval is $\sim 30$, and thus the $\{x_{ij}\}$ drawn from the sampling distribution are, approximately, identically and normally distributed if transformed as $y := (\lambda - x)/\sqrt{\lambda}$. The smooth black distribution in each panel is a normal distribution. The data sampling-distributions for the ST-U model and models higher in complexity do not exhibit systematic structural inaccuracies in the context of the variates: there are no clear signs of noise-model inaccuracy or over-fitting, nor are there clear signs of residual correlations in, e.g., Figure 13. These models can thus perform adequately as point-measure predictors of structure in the real count-number data. The residuals conditional on the ST-S model are distributed with larger variance than unity (normal): more weight is visible in the wings, and less in the near vicinity of zero.

(The complete figure set (3 images) is available.)

was well determined by imaging nearby source-free regions of the sky, and was subsequently used to impose that the surface hot spots generate the remaining signal from the imaged (point-source) MSP. Such a model is distinctly different from those we consider here for NICER, a non-imaging telescope; we do not impose (e.g., via some informative background prior) the signal to be generated by the surface hot regions. As a result, we may inaccurately subsume a cooler contribution from the hot regions into the phase-invariant likelihood terms as suggested above; Bogdanov & Grindlay (2009), on the other hand, may miss some non-diffuse radiative component(s) in the near vicinity of the surface of PSR J0030+0451, and thus require their hot spots to generate additional cooler emission than they do in physical reality.

With these considerations in hand, we conclude that without further work, it is unclear which inferred signal is more physically accurate. Indeed, the existence and treatment of additional X-ray emission by PSR J0030+0451—or from its circumstellar vicinity—in the NICER waveband is considered an open question for future modeling. Moreover, the sensitivity of parameter estimation to requiring a certain phase-invariant contribution specifically from the hot regions could be investigated with newly allocated computing resources.

A (superficially) straightforward alternative to a full physical generative model of the non-pulsed surface and

off-surface emission would be to define bounds—upper in particular—on the prior support of the channel-by-channel background count-rate parameters (refer to related discussion on limits in Section 2.4.3). One must then address the question of how. One option is to move the difficulty of defining a prior density for the background to another level in a Bayesian hierarchy via a hyperparameter (e.g., the upper bound of a flat density function) and a corresponding hyperprior for each channel; each hyperparameter and background parameter pair must then be jointly numerically marginalized over to ensure sampling is tractable in a lower-dimensional parameter space. Another option would be to define a conditional prior distribution for each background parameter, where the upper limit of the support is a function of the source parameters—e.g., some fraction of the phase-average source counts in a channel.

We did not consider at the outset of this work the possibility of jointly modeling NICER and XMM event data—a strategy both tractable and arguably more rigorous given that we do not have likelihood function information nor posterior information from Bogdanov & Grindlay (2009) that is compatible as a prior information to be updated conditional on the NICER data. We could thus consider constructing a joint likelihood function over observations with both telescopes (requiring definition of at least one more nuisance parameter). On the other hand, the