NICER event data comprise far more photons than do the XMM data, and so inferences may be dominated by the (background-marginalized) NICER likelihood function. There are, however, good prospects for improving our understanding of the NICER background (particle radiation and diffuse sky terms; see Section 2.4.3 and Bogdanov et al. 2019a), even without new input from XMM observations or other imaging capabilities.

Observations with a future high-time resolution soft X-ray spectroscopic imaging telescope, such as the *Athena X-ray Observatory* (Nandra et al. 2013), would be synergistic with archival NICER data—and/or observations with a future large-area soft spectroscopic timing telescope such as the *enhanced X-ray Timing and Polarimetry mission (eXTP)*; Zhang et al. 2019; Watts et al. 2019) or the *Spectroscopic Time-Resolving Observatory for Broadband Energy X-rays (STROBE-X)*; Ray et al. 2019). The high-sensitivity (read large-area) imaging capabilities of *Athena* (with ~5"–10" angular resolution<sup>71</sup>) and its two instruments, the *X-ray Integral Field Unit (X-IFU)*; 10 µs time resolution<sup>72</sup>) and the *Wide Field Imager* (80 µs time resolution), will permit modeling of the background emission not originating from the line of sight of the PSR J0030+0451<sup>73</sup> system (as permitted by the point-spread function) in a similar manner to Bogdanov & Grindlay (2009) for XMM (and also Bogdanov 2013, for PSR J0437−4715). Furthermore, the microcalorimeter spectral resolution of *X-IFU* (~2.5 eV resolution for photons energies <7 keV) will permit precise measurements of various edges of the interstellar medium (e.g., oxygen K edge at 0.55 keV, iron L edge at 0.72 keV), thereby enabling derivation of independent tight constraints on the hydrogen column density $N_H$ for a given MSP. Archival NICER MSP data, on the other hand, will remain valuable far into the future because the principal NICER mission science objectives focus purely on MSPs; NICER will have compiled $O(10^6)$ s integrated exposures on the primary MSP targets, meaning that despite the smaller total effective area, the size of the event data set will be synergistic with an advanced multifaceted mission such as *Athena* whose science objectives are broader. Consequently, observations with the *Athena* instruments will enable us to jointly model data sets, and make definitive progress on disentangling the signal generated by surface hot regions from the complicated phase-invariant emission detected by NICER (or a future non-imaging telescope).

### 4.1.3. Instrument

We implemented a specific instrument response model (Section 2.4.2): a simple ad hoc parameterization designed to combine several available calibration products and thus account for both energy-independent (absolute flux calibration) and energy-dependent uncertainty. Sensitivity to our choice of astrophysical calibration source, the Crab, also needs further study; it is known that operation in response to an incident radiation field is a function of its properties (see Section 2.4.2).

Future study could explore whether using different calibration sources affects inferences. Another interesting question may be posed as to how the posterior information about this instrument model evolves as we update our knowledge via analysis of other sources (MSPs and otherwise).

As our understanding of the NICER instrument improves, it should be possible to improve on this model. One option might be to adopt a more sophisticated approach based on Lee et al. (2011) and Xu et al. (2014) to estimating the uncertainty in instrumental response for Chandra. At present such a sophisticated approach may be unjustified because computational expense would be amplified by increasing the complexity of the instrument model instead of—or in addition to—the complexity of model astrophysical sources whose nature is more uncertain.

### 4.2. Mass and Radius Constraints in Context

For ST+PST, the inferred mass *M* and equatorial radius $R_{eq}$ are $1.34_{-0.16}^{+0.15} M_{\odot}$ and $12.71_{-1.19}^{+1.14}$ km, where the credible interval bounds are approximately the 16% and 84% quantiles in marginal posterior mass, given relative to the median. The marginal credible intervals may thus be considered as $1\sigma$ intervals (containing 68.3% of the posterior mass), where $\sigma_M/M \sim 11\%$ and $\sigma_R/R \sim 9.2\%$. For completeness, and to assist comparison with results derived using other methods, we also give the values for the 90% credible interval ($M = 1.34 \pm 0.24 M_{\odot}$, $R_{eq} = 12.71_{-1.85}^{+1.83}$ km) and the 95% credible interval ($M = 1.34_{-0.27}^{+0.28} M_{\odot}$, $R_{eq} = 12.71_{-2.14}^{+2.15}$ km). The (equatorial) compactness $GM/R_{eq}c^2 = 0.156_{-0.010}^{+0.008}$ is more tightly constrained than both *M* and $R_{eq}$ individually, at the ~6% level; the 90% credible interval is $GM/R_{eq}c^2 = 0.156_{-0.017}^{+0.013}$, and the 95% credible interval is $GM/R_{eq}c^2 = 0.156_{-0.021}^{+0.015}$.

The effect of a rotation rate of 200 Hz on NSs is a small deformation of the star into an oblate spheroid. The deformation enters in two ways. A star rotating at 200 Hz will have an equatorial radius that is at most 2% larger than a nonrotating star with the same mass for the stiffest equations of state, with smaller increases in radius for soft equations of state (Cook et al. 1994).<sup>74</sup> The polar radius will be smaller than the equatorial radius, an effect that depends on the equatorial radius, mass, and spin, with very little dependence on the EOS (AlGendy & Morsink 2014, and references therein); for the a posteriori most probable exterior spacetimes inferred in this present work, the polar radius is ~1% smaller than the equatorial radius. Rotation also increases the mass of a star compared to a nonrotating star with the same number of baryons, but this is at most a 0.2% effect for the rotation-powered pulsars observed by NICER. The effect of rotation on the location of the innermost stable circular orbit (ISCO) is larger (van Doesburgh et al. 2018): the (prograde) ISCO typically ranges from $5.6r_g-5.8r_g$ for stars spinning at 200 Hz for most theoretical equations of state, compared to the Schwarzschild ISCO radius of $6r_g$. The 95% compactness

<footnote n="71"><sup>71</sup> Refer to Bogdanov et al. (2019a) for an XMM image of the PSR J0030+0451 field; *Athena* promises to improve on the XMM point-spread function half-energy width by a factor of ~1.5–3.</footnote>

<footnote n="72"><sup>72</sup> Similar to that expected for the *eXTP Spectroscopic Focusing Array* (~10 µs), but not as good as spectro-timing dedicated missions, like NICER or STROBE-X with ~0.1 µs time resolution. Note that event time-tagging resolution offered by X-IFU will still achieve a level of ~1/500 of the spin period of PSR J0030+0451.</footnote>

<footnote n="73"><sup>73</sup> And other MSPs targeted by a mission such as NICER.</footnote>

<footnote n="74"><sup>74</sup> Note that a one-parameter sequence of nonrotating stars deforms into a sequence of stars rotating at some rate, but a unique deterministic map between the central densities that parameterize those sequences does not exist. This reflects the absence of a unique physical mode for stars to evolve in rotation rate. For comparative purposes, one often considers sequences—parameterized by rotation rate—that conserve a quantity such as total mass, total baryon number, or central density. In the accompanying Letter of Raaijmakers et al. (2019), it is the central density that is explicitly defined as a model parameter (to be marginalized out), and which is held constant for comparison of nonrotating stars to rotating stars in the context of likelihood function evaluation.</footnote>
credible interval corresponds to equatorial radii in the range of $5.84r_g-7.41r_g$, meaning that the most probable stellar surfaces<fnref n="75" /> have radii close to or larger than the radius of the ISCO, and thus that the innermost permitted stable orbit is typically at (or just exterior to) the surface.

In this subsection we proceed to discuss how the constraints derived in this work compare to existing constraints on mass, radius, and compactness derived using other methods and independent observations. We also consider prospects for improving NICER constraints on PSR J0030+0451. When comparing radii, it should be remembered that most of the published radius determinations using other methods have assumed that the star is spherical. For luminosity radius determinations of rapidly rotating NSs (spins of a few hundred Hz), systematic errors of 5% could be introduced (Bauböck et al. 2015), although other systematic errors may well dominate over an inaccurate treatment of rotational surface and metric deformation.

The present constraints on the radius are consistent with the previous radius lower limit for this pulsar, $R_{eq} > 10.7$ km (95% confidence, assuming a $1.4 M_\odot$ NS), obtained from early pulse-profile modeling of XMM data (Bogdanov & Grindlay 2009). It is also in agreement with the radii inferred from XMM observations of other MSPs: $R_{eq} > 11$ km (3$\sigma$ confidence) for PSR J0437-4715 (Bogdanov 2013); and $R_{eq} > 7.8$ km (68% confidence) for PSR J2124-3358 (Bogdanov et al. 2008). However, these early light-curve models only considered uniform-temperature circular hot spots, neglected stellar oblateness, and we consider the statistical computation described in this Letter as more advanced.

Early constraints on the NS radius were obtained from the X-ray spectroscopic modeling of the thermal emission originating from the entire surface of isolated NS RX J1856.5-3754 (Burwitz et al. 2001; Drake et al. 2002; Pons et al. 2002). However, uncertainties due to the distance and due to the calculations of radiative transfer in the magnetized atmosphere hindered an accurate radius constraint for this NS (Ho et al. 2007). More robust constraints, on the other hand, can be extracted from the X-ray spectra of quiescent low-mass X-ray binaries (qLMXBs) hosted in globular clusters. Not only are their distances known to better than ~10% precision, but their purely thermal emission is thought to emerge from non-magnetic NS atmosphere models (similar to the NSX model used in the present work) in order to estimate the apparent radius $R_\infty = R_{eq}(1+z) = R_{eq}(1 - 2GM/R_{eq}c^2)^{-1/2}$. However, the degeneracy between $R_{eq}$ and $M$ in the estimation of $R_\infty$ precluded obtaining useful constraints on the EOS (e.g., Webb & Barret 2007; Heinke et al. 2006, 2014; Guillot et al. 2011), due to the typical shapes of the $M-R_{eq}$ confidence contours that made them compatible with many families of EOS. These results prompted the simultaneous analysis of a set of sources.

In these combined analyses, the degeneracy between $M$ and $R_{eq}$ was lifted by assuming a parameterized shape for the EOS, either a toy-model (constant-radius EOS, as a simplistic representation of nucleonic EOS) or an analytical representation using polytropes. A handful of qLMXBs in globular clusters (up to seven) and of Type-I X-ray bursters (four or five) have been combined to produce constraints on the EOS and/or the radii of NSs (given the

assumed EOS shape). Early works using only qLMXBs produced rather small NS radii, $R_{eq} \approx 9-10$ km with ~10%–15% uncertainties (90% credible interval; Guillot et al. 2013; Guillot & Rutledge 2014; Guillot 2016), but the addition of new data and the use of more recent globular cluster distance measurements resulted in higher values, $R_{eq}$ in the 9.9–11.2 km range for a $1.5 M_\odot$ NS (2$\sigma$ credible interval; Bogdanov et al. 2016). Combining qLMXBs and Type-I X-ray bursts, other works have found radii in a wide range of values: $R_{eq} \approx 10.4-12.9$ km (95% credible interval; Steiner et al. 2010, 2013), $R_{eq} \approx 10.5-12.7$ km (90% credible interval; Lattimer & Steiner 2014), and $R_{eq} \approx 9.8-11.4$ km (95% credible interval, Özel et al. 2016), for a $1.4 M_\odot$ NS. However, these analyses may have been affected by systematics, such as those due to the modeling of piled-up X-ray photons in the Chandra data<fnref n="76" />—which was not considered for all qLMXBs in these early analyses—or those due to the choice of atmospheric composition for qLMXBs (generally H versus He). The mass-radius constraints extracted from Type-I X-ray bursts are dependent on the color-correction factors used (between the measured blackbody temperatures and the modeled effective temperatures of the burning atmospheres), which have been debated in the literature (e.g., Suleimanov et al. 2011a; Güver et al. 2012; Güver & Özel 2013; Kajava et al. 2014). These issues cast doubt on the robustness of the error intervals reported in these early works.

Nättilä et al. (2017) recently analyzed Type-I X-ray bursts from 4U 1702-429 by fitting bursting atmosphere models directly to spectra during the cooling tail of the bursts, hence avoiding the use of color-correction factors. These authors found $R_{eq} = 12.4 \pm 0.4$ km and $M = 1.9 \pm 0.3 M_\odot$ (68% credible interval), although the posterior distributions also allow smaller radii ~10 km for higher masses ~2.1 $M_\odot$.

Steiner et al. (2018) considered the effects of pile-up on qLMXB spectra, the possibility of pure helium atmospheres (instead of pure hydrogen atmospheres), as well as non-uniform surface temperature distributions. They obtained constraints on polytropic EOS via Bayesian inference of mass-radius probability distributions of seven NS qLMXBs. These constraints translate to a NS radius in the 10.0–14.4 km range (95% credible interval, assuming a $1.4 M_\odot$ NS), when considering all of the models tested.

More recently, a physically justified parameterization of the EOS was proposed as an alternative to polytropes (Margueron et al. 2018a, 2018b). In that work, the EOS is a meta-model expressed as a Taylor expansion of nuclear physics parameters, and was applied to a combined spectral analysis of seven qLMXBs to directly extract values of nuclear physics parameters $L_{\text{sym}}$, $K_{\text{sym}}$, and $Q_{\text{sat}}$ (Baillot d'Etivaux et al. 2019). Radius estimates were also derived: $R_{eq} = 12.35 \pm 0.37$ km (2$\sigma$ credible interval) assuming a $1.45 M_\odot$ NS.

The cold emission from the MSP PSR J0437-4715 is detectable in the far ultraviolet (Durant et al. 2012) and in the soft X-ray band (0.1–0.3 keV; Guillot et al. 2016), and its mass and distance are known precisely from radio timing (Reardon et al. 2016). Applying NS atmosphere models<fnref n="77" /> has permitted estimation of this pulsar's radius: $R_{eq} = 13.1 \pm 0.8$ km (68%

<footnote n="75">⁷⁵ Which in this present work are only embedded in the ambient spacetime while neglecting rotational metric deformation (see Section 2.4.1) and are thus not self-consistently computed with global numerical solutions to the field equations given interior conditions (including an EOS).</footnote>

<footnote n="76">⁷⁶ It was shown that an unmodeled pile-up fraction as low as ~1% could affect the radius estimated via modeling of spectral data by as much as 10% (Bogdanov et al. 2016).</footnote>

<footnote n="77">⁷⁷ Similar but not identical to the atmosphere models used in our NICER analysis—the models used in Gonzalez-Caniulef et al. (2019) incorporated various effects that are important at lower temperatures such as partial ionization and plasma frequency effects.</footnote>
credible interval; Gonzalez-Caniulef et al. 2019). We note that
this was a phase-averaged spectral analysis, wherein the
emission is assumed to originate from the ~10⁵ K stellar
surface exterior to the (heated) polar cap regions.

An indirect method to constrain NS radii is to use emission
features from an inner accretion disk in accreting LMXBs.
Narrow emission lines, such as Fe K, arising from this rotating
material are asymmetrically broadened—to which there is a
strong relativistic contribution based upon proximity of the
inner disk to the NS surface (Fabian et al. 2000). The accretion
disk must truncate at the stellar surface or at a larger radius:
spectral modeling of these emission lines enables derivation of
a statistical constraint on the inner radius of the (prograde) disk
*in units of* the gravitational radius or the spin-dependent ISCO
radius (Cackett et al. 2008, 2010; Miller et al. 2013; Degenaar
et al. 2015; Ludlam et al. 2017). To do so, an approximative
ambient spacetime solution (e.g., Schwarzschild or Kerr) is
typically invoked—as in this present work—but without the
embedding of a NS surface.<fnref n="78" /> The constrained inner radius can
then be translated into a lower limit on the stellar compactness
for that particular NS.

Ludlam et al. (2017) inferred, for two LMXB systems, inner
radii that are consistent with the disks extending down to the
ISCOs of their respective ambient spacetimes. In the absence of
a constraint on the NS mass that is independent<fnref n="79" /> of the inner
radius of the disk, a mass merely has to be assumed to obtain an
upper limit on the stellar radius. For example, suppose that the
4U 1636-53 system contains a nonrotating<fnref n="80" /> minimally
compact NS (i.e., disk truncation by surface): for a mass of
1.4 M<sub>☉</sub>, the upper limit on the stellar radius would lie close to
the Schwarzschild ISCO, at 12.4-13.1 km (1σ confidence
interval; Ludlam et al. 2017), which is not inconsistent with the
constraint derived conditional on the NICER data.<fnref n="81" /> Lastly, if
the surface of an NS in the 4U 1636-53 system does not
approximately truncate the disk—and is thus more compact
than the ISCO—and/or the (prograde) ISCO is more compact
due to NS rotation, the NS could only be viewed as inconsistent
with our compactness estimate for PSR J0030+0451 if the
following are true: (i) the NS masses are both tightly
constrained and happen to be highly commensurate, and (ii)
the NSs do not occupy a segment of the mass-radius sequence
along which the radius is highly sensitive to increasing mass
due to EOS softening or phase transitions (e.g., Drago et al.
2014; Alford & Han 2016; Alford & Sedrakian 2017).
However, one mass is unconstrained (disk modeling) while
the other is constrained at the ~10%-level (NICER pulse-
profile modeling).

Overall, these recent publications have estimated NS radii in
the range of 12–14 km, which is compatible with our
PSR J0030+0451 radius estimate. We note that equating
accurately measured radii of distinct NSs should elicit
agreement if (i) the EOS is shared from core to crust, meaning
that perturbative effects attributed to crust composition,
temperature, and magnetic field strength are sufficiently small
in the context of measurement precision; (ii) the EOS is of a

nucleonic composition that supports NSs with similar radii over
a wide range in mass (0.8 to ~2.0 M☉); and (iii) differences due
to spin-dependent rotational deformations are accounted for or
are small enough to justify neglecting (see the discussion in the
second and third paragraphs of Section 4.2, and in Raaijmakers
et al. 2019). For other families of EOS, such as those involving
quarks or hyperons in hybrid stars (e.g., Zdunik &
Haensel 2013) or baryon resonances (Drago et al. 2014), the
radius may be (highly) sensitive to increasing mass. In such
cases, we would expect to find NSs whose radii differ by
several km, reinforcing the importance of jointly estimating
both the radius and the mass of each member of a population
of NSs.

Constraints on NS masses and tidal deformabilities are now
also being reported based on the first binary NS merger
gravitational wave event, GW170817 (Abbott et al.
2018, 2019). These can be translated into constraints on mass
and radius. The inferred values vary somewhat depending on
modeling and prior assumptions. Abbott et al. (2018), for
example, employed two methods: the first, which did not
assume that both stars had the same EOS, yielded radii of
10.8<sup>+2.0</sup><sub>-1.7</sub> and 10.7<sup>+2.1</sup><sub>-1.5</sub> km for the two stars, with masses in the
range 1.16–1.62 M<sub>☉</sub>; the second, which assumed a common
EOS, yielded a radius of 11.9 ± 1.4 km for both stars and
masses in the range 1.18–1.58 M<sub>☉</sub> (all results 90% credible
intervals). De et al. (2018), who also assume a common EOS,
report radii of 10.7<sup>+2.1</sup><sub>-1.6</sub> ± 0.2 km and masses in the range
1.12–1.67 M<sub>☉</sub> (90% credible interval). Constraints taking into
account additional information tend to support slightly larger
mean values for the radius, and smaller uncertainties, e.g.,
those derived from the electromagnetic counterpart
(12.4<sup>+1.1</sup><sub>-0.4</sub> km, 2σ confidence interval; Most et al. 2018); and
those invoking a theoretical minimum (Tews et al. 2017) for
neutron matter pressure (11.4<sup>+1.9</sup><sub>-0.8</sub> km, 90% credible interval
marginalizing over mass; T. Zhao & J. M. Lattimer 2019, in
preparation). The results of Abbott et al. (2018) and De et al.
(2018), which employ the assumption that both stars share the
same EOS, suggest that the radii of the two stars are nearly
equal despite the fact that the mass ratio of the stars could lie
between 0.7 and 1.0 with almost uniform probability. Let us
make the assumption that PSR J0030+0415 also has the same
radius as the binary members to well within the posterior
uncertainty on each star: the radius reported in this Letter,
inferred from NICER data, is more consistent with values in the
upper ranges emerging from the gravitational wave analysis of
GW170817. Larger radii would only be consistent with mass
ratios closer to unity. It follows that if the three stars have
nearly the same radius, and if the mass ratio of GW170817 was
near its lower limit of 0.7, the common radius should be at the
lower end of the NICER range. Electromagnetic observations
might suggest relatively large amounts of dynamical ejecta,
which would favor mass ratios considerably less than unity
(Radice et al. 2018).

*What are the prospects for improving constraints on mass and radius for PSR J0030+0451?* Unfortunately we cannot obtain an independent constraint on the mass for PSR J0030+0451 because it is not in a binary, unlike some of NICER’s other MSP targets. Our model of the NICER background (particle radiation and diffuse sky terms) is, however, expected to improve without needing to wait for input from other telescope missions; understanding this background accurately will prove crucial, as discussed in Section 4.1.2, and may well

<footnote n="78">⁷⁸ The real surface could thus in principle enclose the ISCO associated with the ambient spacetime solution.</footnote>

<footnote n="79">⁷⁹ For example, via a binary mass function together with classification of the companion star.</footnote>

<footnote n="80">⁸⁰ Or at least a small dimensionless spin, despite the 581.0 Hz spin frequency (Ludlam et al. 2017, and references therein).</footnote>

<footnote n="81">⁸¹ Where as discussed above the rotational deformation of the surface is small for PSR J0030+0451.</footnote>
impact mass–radius estimation. A longer total exposure time could also certainly be accumulated (and indeed the NICER team anticipate doing this). Previous studies that have examined how posterior estimation of mass and radius is sensitive to factors such as geometry, spin rate, and the number of source counts in the event data, indicate that constraining power increases as the square root of the number of counts (Lo et al. 2013; Psaltis et al. 2014). However, those studies all assumed a single circular single-temperature hot spot, not a more complex hot region configuration such as those we have considered and inferred here. While it is likely that gathering more data will improve the joint constraint on mass and radius (without reference to an EOS model), the precise observing time required to achieve a given level of precision cannot be estimated robustly without further study.

Our report here is encouraging in terms of prospects for other NICER targets such as PSR J0437−4715: the pulsar mass is constrained independently to within a few percent via radio pulsar timing because it is in a binary system (Reardon et al. 2016). If such a constraint had been available for PSR J0030+0451, it is clear that we could have obtained a comparable posterior uncertainty on the inferred radius of a few percent.

### 4.3. EOS Implications

One of the primary goals of mass–radius inference is to use posterior information<fnref n="82" /> to infer the properties of the dense matter EOS, if such information is deemed sufficiently likelihood-dominated to warrant the study. Studies have utilized joint mass–radius posteriors inferred from X-ray spectral modeling of bursting and quiescent NSs (Steiner et al. 2010, 2013; Özel et al. 2016; Raithel et al. 2017). Moreover, studies have utilized mass and tidal deformability constraints derived from analysis of the NS binary merger event GW170817 (Abbott et al. 2018; Annala et al. 2018; Lim & Holt 2018; Malik et al. 2018; Most et al. 2018; Tews et al. 2018b; Carson et al. 2019; Li & Sedrakian 2019; Montaña et al. 2019), and consideration is already being given to combining constraints from electromagnetic and gravitational wave analysis (Kumar & Landry 2019; McNeil Forbes et al. 2019; Weih et al. 2019).

Given a suitable model for the EOS (see, e.g., Read et al. 2009; Raithel et al. 2016; Lindblom 2018; Tews et al. 2018a) there are two approaches to EOS inference: one is to jointly infer the EOS parameters (and central densities) directly from the data (e.g., pulse-profile data); the other is to jointly infer EOS parameters from per-source nuisance-marginalized likelihood functions of exterior-spacetime parameters (e.g., gravitational mass and equatorial radius). The former approach is *at least* as computationally intensive as the direct mass–radius inference reported in this Letter; the latter approach is less computationally intensive *given* archival likelihood function information about exterior spacetime parameters (Riley et al. 2018). In any case, care is required in both overall approach and the selection of (interior source matter and exterior spacetime) model parameterization and priors (Carney et al. 2018; Raaijmakers et al. 2018; Riley et al. 2018; Greif et al. 2019; Landry & Essick 2019). For our analysis here we deliberately defined a joint flat prior density function for *M* and *R*<sub>eq</sub>, with the intention that the posterior density function can be invoked as a likelihood function marginalized over all nuisance

parameters. We explore the dense matter EOS implications of the inferred mass, radius, and compactness for PSR J0030+0451 in an accompanying Letter (Raaijmakers et al. 2019), following the approach to EOS inference outlined in Greif et al. (2019).

### 4.4. Implications of the Surface Heating Configuration

We constructed a sequence of simple models for the properties of the two hot regions, nevertheless motivated by (numerical) pulsar theory. We considered models in which the regions were related via antipodal reflection symmetry with respect to the stellar origin, and models that do not impose such symmetry, meaning that their properties and location were described with distinct parameters (with the restriction that the regions cannot overlap). The models included simply connected circular and crescent regions, and rings (whose hole and annulus are concentric or eccentric), each filled with single-temperature material. The models also included annular (ring) regions whose concentric hole is filled with material of finite temperature distinct from that of the material in the annulus.

We were able to rule out the hot regions being antipodal and identical based on clear systematic structure in the residuals between data and model a posteriori; moreover, a model wherein the regions are both assumed to be simply connected circular single-temperature spots was strongly disfavored. We inferred that the regions are configured to exist in the same rotational hemisphere: one region subtends an angular extent of only a few degrees (in spherical coordinates with origin at the stellar center) but whose other structural details we are insensitive to; the other region is far more azimuthally extended,<fnref n="83" /> in the form of a narrow hot crescent. The inferred effective temperature of the NSX atmosphere was remarkably consistent across all models considered—for both regions—at ~1.3 × 10⁶ K. The ST+PST model exhibited the largest background-marginalized likelihood function values in the typical set of a posterior mode. Figure 17 renders a representative configuration from the posterior mode, and Table 2 reports the marginal credible intervals for the hot-region parameters. Note that ST+PST includes within prior support, configurations wherein the regions are similar or even congruent in shape, and a priori favors (albeit weakly) smaller angular extents—the heating asymmetry is emergent in spite of this.

One of the principal astrophysical questions arising is how such a heating configuration can occur. It appears to be incompatible with magnetospheric current heating at the footpoints of a simple near-centered dipole magnetic field, and is likely to require some higher-order multipole structure (Barnard & Arons 1982; Gralla et al. 2017; Lockhart et al. 2019). We now need to determine the type of field configuration required, the magnitude of the different moments, and whether this is feasible on physical grounds. Consideration will also need to be given to how magnetospheric currents actually map to temperature fields on the stellar surface. We note that there are clear similarities between the inferred ST+PST configuration and the current heating distribution contemporaneously derived by Lockhart et al. (2019, see their Figure 6 in particular) via quadrupolar extension of the magnetic field, considering that the heating ring is asymmetric

<footnote n="82">⁸² Strictly, likelihood information (Riley et al. 2018).</footnote>

<footnote n="83">⁸³ In spherical coordinates with polar axis coincident with the stellar rotation axis.</footnote>
with respect to the dipole axis, and closely resembles a large-scale arc-like hot region. There are, however, also some differences: Lockhart et al. (2019) restricted their study to configurations where the center of the heated ring is antipodal to the heated spot (both of which emit as approximate blackbodies), an assumption that would need to be relaxed to recover our preferred configuration.

There are also implications for pulsar emission in wavebands other than the X-ray (see the reviews by Grenier & Harding 2015; Cerutti & Beloborodov 2017). If a multipolar field structure is required to explain the surface temperature field, how does this affect radio and gamma-ray emission generated further out in the magnetosphere? Could the multipole structures persist out to the point where emission in these wavebands is thought to happen? Quadrupole fields fall off faster than dipole fields, as the inverse fifth power of the radius, but the radius beyond which the field is predominantly dipolar would depend on the ratio of quadrupole to dipole components. Most current models of radio and gamma-ray emission assume that the field structure is a centered dipole (e.g., Radhakrishnan & Cooke 1969; Gil et al. 1984; Kijak & Gil 2003; Dyks & Harding 2004; Johnson et al. 2014), and this would need to be revisited.

There are also questions pertaining to stellar evolution. NSs are born with a field structure that could be quite complex as a result of the supernova process (Ardeljan et al. 2005; Obergaulinger & Aloy 2017) but various diffusive evolutionary processes can subsequently modify field structure even for isolated NSs (Reisenegger 2009; Viganò et al. 2013; Mitchell et al. 2015; Gourgouliatos & Hollerbach 2018). In addition, rotation-powered MSPs are thought to go through an extended period of accretion-induced spin-up to reach the observed spin rates. The accretion process may also act to modify the field structure (Romani 1990; Melatos & Phinney 2001; Payne & Melatos 2004). It remains to be determined whether a complex multipolar field structure could emerge and survive from birth, or be generated during the accretion process. If such a field structure is present or evolves during the accretion phase, there will also be implications for the spin-up process and for X-ray emission during that phase of NS evolution. If the magnetic field were to channel accreting material onto two magnetic polar caps on the same hemisphere, for example, this would certainly affect the emission from accreting MSPs (Long et al. 2007, 2008; Patruno & Watts 2012). Whether the star is even visible as an accreting pulsar will depend not only on the geometry of the hot regions where accreting material impacts the star, but also on whether the observer views the hemisphere containing the polar caps or the other one. The flow at the inner-edge of the accretion disk, a strong source of potentially variable X-ray emission, would also be affected by a multipolar field structure. Finally, extremely off-center dipoles or strong non-centered multipole fields, for instance, will produce asymmetries in the Poynting flux of low-frequency radiation parallel to the spin axis, with consequences that could include a large space velocity (Harrison & Tademaru 1975). Lommen et al. (2006) found that PSR J0030+0451 has a relatively low transverse space velocity—a property that is potentially in contention with a field far from that of a centered dipole—although since PSR J0030+0451 is isolated its space velocity would also depend on how the binary was disrupted after the spin-up phase. Discussion on some of these issues is reserved for an accompanying Letter (Bilous et al. 2019). However,

further work on the implications of the inferred configuration for pulsar field structure, emission mechanisms, and stellar evolution is certainly required. It is clear that the mass and radius inferred for PSR J0030+0451 depend strongly on the surface radiation field models (including prior support) that we have explored. Further study may show that our models are either not general enough or too general (for example with respect to the prior support). For example, we may find that temperature gradients in the hot regions cannot be neglected; or that field stability considerations impose a minimum angular separation between the polar caps—and possibly, by extension, the hot regions—larger than found in our analysis; or that no magnetospheric model can generate a heating configuration in which one region is a small-scale spot while the other an azimuthally extended crescent.

There are a number of computational aspects for pulsar theorists to consider when developing surface heating models suitable for statistical inference. We were only able to consider the configurations offered uniquely by ST+PST as we approached the limit of our computational resource allocation; we can therefore provide stronger guarantees about the accuracy of the posterior computation for the lower-complexity models with more than one run (ST-U and ST+CST). In the future, we suggest that more resources be devoted to models at the ST+PST-level of complexity, in particular for work on (i) parameterization of hot regions with more complex topologies and/or boundaries, (ii) their efficient numerical resolution for likelihood function evaluation, and (iii) exploration of the associated parameter space via sampling or other methods.

We also suggest that additional resources be devoted to research avenues such as self-consistent theory and computation of surface heating by magnetospheric currents, for the purpose of statistical computation. A question may be posed as to the generation of highly non-dipolar surface temperature fields—e.g., arc- or ring-like heating distributions—which can be (i) parameterized such that approximate representations find compromise between accuracy and complexity, capturing the facets considered most crucial to signal generation, perhaps as a sequence of models increasing in complexity; and (ii) built into efficient software implementations. Progress on such fronts should encourage a bridge to form between phenomenological efforts and more self-consistent theory for the purpose of efficient statistical computation, and may also offer a way in which to connect distinct theoretical models (in an approximative manner) on a continuous space.

### 4.5. Pulse-profile Modeling for Other Types of NS

NICER is the first mission designed specifically to use the pulse-profile modeling technique to infer the mass and radius of NSs. As a soft X-ray telescope with an effective area of less than a square meter, it is optimized for applying the technique to MSPs, which have soft, stable pulse profiles meaning we can use multiple exposures taken over a long baseline to accumulate a sufficient number of events to statistically probe MSP physics. However, the technique can also be applied to other NSs with emission modulated by rapid rotation: accretion-powered pulsars and thermonuclear burst oscillation sources.

In accretion-powered pulsars (see Patruno & Watts 2012, for a review), accreting material is channeled by the magnetic field onto the magnetic polar caps and the pulsed emission has two main components: thermal emission from the heated vicinity of
the accretion impact zones, and nonthermal emission from the shock in the accretion funnel (Poutanen & Gierliński 2003). A third pulsed component may arise due to reflection from the accretion disk (Wilkinson et al. 2011). Thermonuclear burst oscillations (see Watts 2012, for a review) are generated by rotational modulation of global asymmetries that form in a surface radiation field during thermonuclear (Type-I X-ray) bursts; such a burst occurs in the ocean of an accreting NS, driven by unstable burning of accreted hydrogen, helium, or carbon (see, e.g., Galloway et al. 2008). The precise mechanism driving the *detectable* asymmetry (oscillations) is not clear: possibilities include disrupted flame spread (Spitkovsky et al. 2002; Cavecchi et al. 2013), large-scale waves in the burning ocean (Heyl 2004; Piro & Bildsten 2005; Chambers et al. 2019), or patterns triggered by convection (Garcia et al. 2018).

Accretion-powered pulsations and thermonuclear burst oscillations are radiatively harder (~1–30 keV) than the pulsations of MSPs. Accumulating the requisite number of photons for tight constraints, in a realistic observing time, also requires a telescope with an effective area of several square meters (Watts et al. 2016; Watts 2019). Several mission concepts are currently being developed for large-area broadband X-ray timing telescopes that would access a larger, fainter population of MSPs than we can observe with NICER, and fuel pulse-profile modeling for accretion-powered pulsars and thermonuclear burst oscillators: these include the eXTP (Watts et al. 2019; Zhang et al. 2019), and the STROBE-X (Ray et al. 2019). For an idea of the constraints that can be delivered by pulse-profile modeling and inference using existing data, see the Salmi et al. (2018) analysis of Rossi X-ray Timing Explorer observations of the accretion-powered MSP SAX J1808.4–3658.

Our uncertainty toward surface heating physics affects pulse-profile modeling not only for rotation-powered MSPs, but also accretion-powered pulsars and thermonuclear burst oscillators. Although the general mechanism that gives rise to the pulsed components in accretion-powered pulsars is clear, the surface and off-surface temperature field and local comoving beaming function—particularly from the accretion-funnel shock—are a priori highly uncertain for any given source. Thermonuclear burst emission has a well-understood (local comoving) beaming function due to the sub-surface thermal origin (Suleimanov et al. 2011b), but the mechanism for generating asymmetries in the global surface radiation field, which in turn generate such rotational oscillations, remains highly uncertain. Reducing the remaining theoretical uncertainties, and developing physically motivated parameterized models of the surface temperature field, will be important. However, our analysis of NICER data provides an important real-world demonstration that pulse-profile modeling is a viable technique for constraining masses and radii of NSs, and that the analysis machinery can operate on somewhat flexible models, with weakly informative priors, for both source and background emission.

## 5. Conclusion

We reported on pulse-profile modeling efforts for the rotation-powered millisecond X-ray pulsar PSR J0030+0451, conditional on NICER data. We focused on PSR J0030+0451, a challenging source due to the absence of an independent constraint on the gravitational mass (compared to NICER's other primary target PSR J0437−4715). Nevertheless, PSR J0030+0451 was selected as the optimal source to

demonstrate simultaneous inference of gravitational mass and equatorial radius given weakly informative priors, and to develop our analysis procedures.

The mass and radius each have marginal posterior 68% credible interval half-widths at the ~9%–11% level, conditional on the NICER XTI event data. These constraints are consistent with those emerging from both gravitational wave analysis and X-ray spectral modeling, and are expected to improve with further exposure. The compactness is constrained more tightly, at the ~6% level. Prospects for NICER delivering tight constraints for rotation-powered MSPs where the mass is known independently to uncertainties of a few percent are clearly excellent.

In addition to inferring properties of the spacetime (mass, radius, and compactness), we were also able to infer the properties of the thermally emitting hot regions that we assume generate the pulsations. For the specific set of models that we considered, the inferred configuration has both hot regions in the same rotational hemisphere, with one hot region being a small spot and the other an azimuthally extended narrow crescent. Models wherein the hot regions are antipodal are strongly disfavored, implying a complex offset dipolar and multipolar field structure that, if accurate, has major implications for both pulsar emission and stellar evolution.

We thank the anonymous referees for their suggested improvements to this work. This work was supported in part by NASA through the NICER mission and the Astrophysics Explorers Program. T.E.R., A.L.W., and A.V.B. acknowledge support from ERC Starting grant No. 639217 CSINEUTRON-STAR (PI: Watts). A.L.W. would also like to thank Ralph Wijers for environmental support. This work was sponsored by NWO Exact and Natural Sciences for the use of supercomputer facilities, and was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This research has made extensive use of NASA's Astrophysics Data System Bibliographic Services (ADS) and the arXiv. R.M.L. acknowledges the support of NASA through Hubble Fellowship Program grant HST-HF2-51440.001. S.G. acknowledges the support of the Centre National d'Études Spatiales (CNES). W.C.G.H. appreciates use of computer facilities at the Kavli Institute for Particle Astrophysics and Cosmology. S.M.M. thanks NSERC for support. J.M.L. acknowledges support from NASA through Grant 80NSSC17K0554 and the U.S. DOE from grant DE-FG02-87ER40317.

**Facility:** NICER (Gendreau et al. 2016).

**Software:** Python/C language (Oliphant 2007), GNU Scientific Library (GSL; Gough 2009), NumPy (van der Walt et al. 2011), Cython (Behnel et al. 2011), SciPy (Jones et al. 2001), OpenMP (Dagum & Menon 1998), MPI (Forum 1994), MPI for Python (Dalcín et al. 2008), Matplotlib (Hunter 2007; Droettboom et al. 2018), IPython (Perez & Granger 2007), Jupyter (Kluyver et al. 2016), TEMPO2 (photons; Hobbs et al. 2006), PINT (photonphase; https://github.com/nanograv/PINT), MULTINEST (Feroz et al. 2009), PYMUL-TINEST (Buchner et al. 2014), GetDist (https://github.com/cmbant/getdist), nestcheck (Higson 2018; Higson et al. 2018, 2019), fgivenx (Handley 2018), X-PSI (v0.1; https://github.com/ThomasEdwardRiley/xpsi; Riley & Watts 2019).
# Appendix A
## Posterior Computation

Here we describe how we derived posterior inferences conditional on the models defined in Section 2.

### A.1. Nested Sampling

We implemented nested sampling using the open-source software MULTINEST<sup>84</sup> (Feroz & Hobson 2008; Feroz et al. 2009, 2013; Buchner et al. 2014). For a subset of model nodes in Figures 10 and 12 we executed two production runs—i.e., we generated two realizations of a particular stochastic sampling process—with the following resolution settings: the number of active (or live) points was $10^3$; the bounded-hypervolume expansion factor was (at least<sup>85</sup>) $0.3^{-1}$; the termination condition as a function of iteration number $i$ was, schematically,

$$ \ln(\hat{Z}_i + \Delta\hat{Z}_i) - \ln(\hat{Z}_i) < x \Rightarrow \frac{\hat{Z}_i}{\hat{Z}_i + \Delta\hat{Z}_i} > e^{-x}, \quad (10) $$

where $\hat{Z}_i$ is an estimator for the evidence integral up to iteration $i$ over estimated prior mass ($1 - \hat{X}_i$), and $\Delta\hat{Z}_i$ is an estimator for the maximum remaining contribution to the evidence over the complementary estimated prior mass $\hat{X}_i$. As $x \to 0$, $e^{-x} \to 1^{-}$ and the estimated amassed evidence entirely dominates the estimated remaining evidence; we generally set $x = 10^{-1}$, which is five times smaller than suggested by Feroz et al. (2009).<sup>86</sup> The number of active points is chosen to be a number that is larger than will be typically reported in the literature for similar dimensional problems (perhaps with simpler distributional structure), or recommended by the authors of MULTINEST, for a compromise between resource consumption and accuracy. The hypervolume expansion factor is, for all but one run,<sup>87</sup> greater than the number recommended by the authors of MULTINEST for accurate evidence estimation. The combination of a number of active points—which is between 40 and 80 times larger than the dimensionality of the sampling space—and the recommended expansion factor, targets posterior computation with an implementation-specific error<sup>88</sup> that is smaller than the error due to the inherent stochasticity of Monte Carlo sampling.

The constant-efficiency MULTINEST bounding variant should in general be avoided where tractability is not compromised to reduce risk of under-sampling when a hyper-ellipsoidal decomposition does not conform well to likelihood level hypersurfaces; we did not use this bounding variant.

However, in the higher-dimensional contexts that we are approaching in this work, it could be useful to activate this bounding variant if integration is forecasted to consume too many computing resources (e.g., the acceptance fraction drops too low below $10^{-3}$ when using a high-performance system and an expensive likelihood function). Imposing a target acceptance fraction may be useful for initial exploratory runs in order to probe for configurations with high (marginal) likelihood—a weak indicator of potential model performance if future efforts achieve greater computational efficiency. When constant-efficiency mode is activated the evidence estimates should be assumed to be positively biased unless importance nested sampling is also activated (which can consume an appreciable fraction of the memory on a typical supercomputer node for ~$10^8$ likelihood function evaluations). The acceptance fraction decaying to such levels is indicative that the minimum-bounding hyper-ellipsoidal decomposition does not conform sufficiently well to the nested likelihood function level hypersurfaces, resulting in (i) a large excess of hypervolume exterior to the likelihood surface but within the bounding union of ellipsoids, and/or (ii) the union of bounding ellipsoids exhibits a large fractional overlap. In such cases one could construct an alternative parameterization of the problem that is more tractable with the MULTINEST bounding algorithm (in the native sampling space), or use an alternative (nested) sampling algorithm.

The mode-isolation sampling variant was not activated unless stated otherwise. This variant isolates the evolution of local modes whose hyper-ellipsoidal clustering decompositions are mutually non-overlapping; upon isolation, the active points of sampling threads constituting each mode cannot migrate between modes at subsequent iterations—they are locked in, which can alleviate premature mode deactivation. Sampling resolution, however, can be absorbed by unimportant modes, which is undesirable. The mode separation sampling variant generates distinct statistics for each mode and labels actives points according to mode association.

With the above settings, among ulterior numerical settings for marginal likelihood function evaluation, X-PSI typically executes $O(10^7)$ marginal likelihood function evaluations in $O(10^4)$ core hours via MULTINEST.

### A.2. Summary of Information Presentation

We now summarize how we opt to present probabilistic information in the form of figures and tables. We condition on a set of model variants that are treated in effectively the same manner, and the information presentation is consistent for these models; we thus here describe the information once. Figures and tables are displayed in the main body—see Section 3—for the model we deem to perform best considering both prior and posterior predictive measures and checks; we refer the reader to Table 2 and Figures 13 through 21 for reference.

#### A.2.1. Graphical Posterior-predictive Checking

In a figure such as Figure 13 we display salient information for assessing performance of a model in isolation. Our generative modeling process is fundamentally built on the statement that we do not believe the true data-generating process exists within the model space considered (refer to Section 2.3); however, the models may, for the purpose of generating data, be deemed adequate approximations. We aim

<footnote n="84"><sup>84</sup> MULTINEST v3.11 can be located at https://github.com/farhanferoz/MultiNest with SHA1-hash 4b3709c. PYMULTINEST v2.6 can be located at https://github.com/JohannesBuchner/PyMultiNest SHA1-hash 5d8c103.</footnote>

<footnote n="85"><sup>85</sup> A lower bound for reasons pertaining to the implementation of the joint prior distribution; the lower bound itself is numerically transformed into an appropriate MULTINEST setting in order to achieve at least this desired expansion factor. Further detail is beyond the scope of this description (refer to Riley & Watts 2019).</footnote>

<footnote n="86"><sup>86</sup> The minimum tolerance used for any process was $x = 10^{-3}$ to confirm that termination was not premature.</footnote>

<footnote n="87"><sup>87</sup> The ST+PST model was implemented as we were exhausting our computational resources and thus we lowered resolution—increased the nested-sample acceptance fraction—to ensure completion of an exploratory run—see Table 2.</footnote>

<footnote n="88"><sup>88</sup> The error due to not sampling from the entire prior hypervolume subject to a given likelihood function constraint.</footnote>
to graphically approximate an answer to the question: *does the model generate*, a posteriori, synthetic data that emulates structure in the real event data?

Graphical posterior predictive checking here relies on the power of human identification of systematic structural differences, which if physically characterized can drive future model development toward better-performing approximations. Structural differences include spectro-temporal correlations between random variables assumed to be statistically independent (e.g., the ST-S version of Figure 13 in the online figure set), and inaccurate noise modeling leading to under- or over-estimation of the variance of random variables; the latter also may manifest due to excessive, non-physical predictive complexity (over-fitting). In the top panel of such a figure we display the data set for convenience. In the middle panel we choose to display the posterior-mean Poisson expected count numbers:

$$ \lambda_{ij} := \mathbb{E}_{\pi(\theta)}[c_{ij}(\theta)] = \int_S c_{ij}(\theta, B_i) \pi(\theta, B_i | d) d\theta \\ \approx \sum_k w_k c_{ij}(\theta_k), \qquad (11) $$

where $\theta_k \sim \pi(\theta | d)$ are samples with normalized importance weights $w_k$ drawn from the background-marginalized posterior density $\pi(\theta | d)$, and $S$ is the prior support. Crucially, we do not sample the joint posterior distribution $\pi(\theta, B | d)$ because of prohibitive scaling of expense with dimensionality;<fnref n="89" /> therefore, in order to compute a data-space posterior-mean quantity, we opt to maximize the conditional likelihood function $L(B_i; \theta_k)$ with respect to each background count-rate parameter $B_i$, generating an estimator $\hat{B}_i$, such that $c_{ij}(\theta_k) := c_{ij}(\theta_k, \hat{B}_i)$. In the bottom panel we display standardized residuals between data count numbers and the quantities $\lambda_{ij}$, where we define Poisson-random variables $x_{ij} \sim p(x_{ij} | \lambda_{ij})$. We consider the set of figures associated with Figure 13 as graphical posterior-checking plots. If no clear systematic structure manifests during posterior predictive checking and the posterior-predictive distribution is well approximated by a sampling distribution conditional on a parameter vector, we consider the model to be appropriate for predicting future observations against which the model may be falsified.

Bayesian global performance measures—such as the evidence—are useful for determining the utility of increments in the complexity of a generative model. However, while being a target of our posterior computation, prior predictive probabilities are known to not be universally robust and can be sensitive to prior choices and parameterization, especially in phenomenological contexts. Thus in order to assess performance and identify model features that conform well to data structure, it is

necessary—and natural—to also visually inspect predictions a posteriori. We consider other measures such as the expected posterior utility (i.e., power for future statistical falsification) beyond the scope of this work (Vehtari & Ojanen 2012).

In this work we refrain as much as possible from invoking point measures—based on parameter vector point estimates—to quantitatively summarize and compare models. The reason we do this is twofold: (i) a point estimate is usually the parameter vector that is estimated to globally optimize some quantity, and subsequent to estimation all information encoded by the posterior distribution is not explicitly regarded; and (ii) the target of our statistical computation (nested sampling) is not any particular point measure, but instead to draw samples from the posterior typical set for the purpose of estimating posterior integrals. It follows that stronger guarantees can be made about the accuracy of statistical estimators that are posterior integrals than can be made about those which are point estimates; we therefore view integral estimators as generally more helpful and robust than point estimators.

As an example, consider the canonical estimation of the parameter vector that globally maximizes the likelihood function conditional on some model. In our case, our posterior computation does not target accurate maximization, and we are forced to marginalize over the phase-invariant (background) count-rate parameters in order to define a nested-sampling space whose dimensionality is not prohibitive. We thus do not guarantee that the sample (parameter vector) $\hat{\theta}_{BMML}$ that reported the highest background-marginalized likelihood function value among identified posterior modes, together with the background count-rate vector $\hat{B}_{CML}$ that maximizes the conditional likelihood function given $\hat{\theta}_{BMML}$, is an adequate estimator of the joint vector $(\theta, \hat{B})_{ML}$ that globally maximizes the likelihood function. In Figure 23 we consider the joint vector $(\hat{\theta}_{BMML}, \hat{B}_{CML})$: empirically we see that all but one model (ST-S) has, within prior support, parameter vectors that perform adequately as point-measure predictors of structure in the real count-number data.

### A.2.2. Parameter Kernel Density Estimation and Credible Regions

We applied the post-processing module of the X-PSI package (v0.1; Riley & Watts 2019). X-PSI wraps—combines and adds functionality—to several other packages for statistical computation: relevant here for Gaussian kernel density estimation (KDE) is GetDist.<fnref n="90" /> We do not use GetDist to compute numerical one-dimensional credible intervals on each parameter, but GetDist is used whenever KDE is necessary for post-processing (posterior) samples into estimators: for calculation of the plotted one- and two-dimensional density functions, for calculation of the plotted two-dimensional (joint) credible regions, and for one-dimensional KL divergence estimation where both the marginal posterior and marginal prior density functions need to be estimated from samples.

GetDist can execute smoothing using either a manual Gaussian kernel bandwidth or an automatic optimized bandwidth. For all parameters and models, for simplicity and consistency, we invoke a manual bandwidth of 0.4 times the

$$ p(x_{ij} | d) = \int p(x_{ij} | \theta, B_i) \pi(\theta, B | d) dB d\theta \\ \approx \sum_k w_k p(x_{ij} | \theta_k, B_{i,k}), $$

<footnote n="89">⁸⁹ If our set of samples was drawn from the joint posterior distribution $\pi(\theta, B | d)$ we could better approximate the expectation integral in Equation (11), or at greater cost, approximate the posterior-predictive probability mass distribution in data space:</footnote>

<footnote>whose expectation $\mathbb{E}[x_{ij}]$ and variance $\mathbb{V}[x_{ij}]$ are indicators of residual structure in the data relative to the model a posteriori. The posterior predictive distribution may also be constructed by jointly generating samples from the Bayesian joint distribution: $x_k \sim p(x | \theta_k, B_k)$ and $(\theta_k, B_k) \sim \pi(\theta, B | d)$.</footnote>

<footnote n="90">⁹⁰ https://github.com/cmbant/getdist. Specifically, tag v0.3.1, with some minor plotting customization locatable at https://github.com/Thomas EdwardRiley/getdist/tree/customization, with SHA1-hash 61f69d0. The technical KDE notes for GetDist are located at https://cosmologist.info/notes/GetDist.pdf.</footnote>
estimated parameter standard deviations, based on the *GetDist* guidelines for choosing analysis settings. Moreover, *GetDist* performs various smoothing corrections to avoid smoothing bias as one attempts to mitigate finite-sample noise: (i) smoothed-density correction near hard one-dimensional parameter bounds defined through the prior support, with capability to estimate a finite local density function gradient; and (ii) iterative multiplicative bias-correction to the estimated density function to nullify over-smoothing. Density estimation near non-trivial prior support boundaries in two parameters does not account for the local boundary, but this is not problematic for our work here. We roll with the default bias-correction settings: a linear boundary kernel (in both one- and two-dimensional spaces) and zeroth-order multiplicative bias-correction (a single application with no iterations).

### A.2.3. Error Analysis for Statistical Estimators

Two frameworks now require distinction: (i) the parametric probabilistic framework of the generative model (a Bayesian context) for the data; and (ii) a non-parametric probabilistic framework that operates with realizations of a stochastic sampling process, which in turn operates on a deterministic target probability distribution defined in the parametric framework given data and a generative model (e.g., Skilling 2006; Higson et al. 2018). The sampling process has a mixture of known and approximated properties, including fixed non-physical settings. Estimators derived from the process output are stochastic and we are interested in their distributions.

The purpose of executing some number of computationally expensive *repeats* is for error analysis—specifically the estimation of implementation-specific error (Higson et al. 2019) pertaining to sampling from the joint prior subject to a likelihood function constraint (a more thorough review of error analysis techniques may be found in Riley & Watts 2019). One way to approach this problem is to compute posterior marginal density functions for many bootstrapped realizations of each stochastic run and then graphically probe for the manifestation of implementation-specific error (Higson et al. 2019).

*X-PSI* wraps the package *nestcheck*<fnref n="91" /> (Higson 2018; Higson et al. 2018, 2019) to access existing error analysis routines. As an example of application in this work, consider Figure 20: in the associated online figure set, where we supply a higher-resolution version of each panel, we display the marginal posterior density distributions for each parameter as a set of shaded error bands. The shaded bands represent the distribution of posterior density, at each parameter value, based on simulated nested sampling process realizations; the posterior density function is estimated for each realization with *GetDist* as described above in Appendix A.2.2. The colorbar denotes the percentage of realizations spanned by a band with a given shade, where each band connects intervals (at each parameter value) containing the highest *realization* density of *posterior density*. Note that the colorbar is not associated with the shaded joint density distributions in the off-diagonal plots. The contours in the on-diagonal panels thus encode information on the variation of the parameter kernel density estimator due to the inherent stochasticity of each of sampling processes; the estimator distributions are connected as a function of each

respective parameter to delineate the behavior of the probability mass. If the member processes are deemed to exhibit consistency under visual inspection, the combined process may be invoked to estimate distributions of estimators.

We apply *nestcheck* routines to bootstrap re-sample threads and simulate weights for the following: the one-dimensional quantiles in posterior mass for each parameter, which are in turn used to report the credible intervals; the global and parameter-by-parameter KL divergences; and the evidence. The numerical values we report are in some cases (ST-S, ST-U, and ST+CST) derived by combining two realizations into a single realization with $O(10^3)$ active points (sampling threads), if the runs are considered sufficiently consistent and exhibit sufficient resolution for our purposes here. However, due to computational expense we could only afford *at most* two runs, and for both ST+EST and ST+PST we were limited to a single run with $O(10^3)$ active points. Moreover, note that processes executed by the mode-separation algorithm are incompatible with the notion of process combination: the theory nor software implementation exists at the time of writing, and therefore we only display such a run—when relevant and available—in the posterior figures, but do not use it to calculate numerical estimators.

### A.2.4. Estimating Posterior Information Gain

In a Bayesian context, the KL divergence (Kullback & Leibler 1951) can be applied as a non-negative real scalar<fnref n="92" /> measure of posterior information *gain* about a parameter—or jointly about parameters—of a generative model, conditioned on the data set; it is also known via information-theoretic interpretation as *relative entropy*. Equivalently it is the *posterior-expected additional* number of bits necessary to encode the value of a parameter sample for lossless communication between agents, *if* the (marginal) prior distribution is invoked to design an optimal encoding. Alternate interpretations exist to satiate a variety of readers—e.g., Shlens (2014).

KL-divergence maximization is central to an information-theoretic—but often in practice intractable—definition of a minimally informative prior via reference to the generative model, but without reference to the data (via data-space marginalization). While we cannot feasibly determine the maximal KL divergence with respect to the space of all proper prior density functions, the number of bits of information gain is a useful indicator of the degree to which the likelihood function dominates the information encoded in the posterior. The KL divergence has, as an example, recently been applied by The LIGO Scientific Collaboration et al. (2019) to probe posterior information gain and sensitivity to prior assumptions.

For continuous symbols, the KL divergence is defined in the limit that a discrete symbol becomes a continuous subset of $\mathbb{R}^n$; in practice, of course, computer representation of the reals is discrete. Given the existence of a known optimal prior encoding, the number of additional bits for lossless communication of a posterior sample at a given precision is, in the context of our inference problem, far less than the number of bits required to store the sample. Mathematically, in units of

<footnote n="91">⁹¹ Specifically, tag v0.2.0, with customization to support use of *GetDist* KDE, locatable at https://github.com/ThomasEdwardRiley/nestcheck/tree/feature/getdist_KDE with SHA1-hash 4555df0.</footnote>

<footnote n="92">⁹² And thus parameterization invariant.</footnote>
bits:

$$D_{\text{KL}}(\pi \| p) = \int_S \pi(\theta | \mathbf{d}) \log_2 \left[ \frac{\pi(\theta | \mathbf{d})}{p(\theta)} \right] d\theta \\ \approx \sum_k w_k \log_2 \left[ \frac{\pi(\theta_k | \mathbf{d})}{p(\theta_k)} \right], \quad (12)$$

where $\theta_k \sim \pi(\theta | \mathbf{d})$ are samples with normalized importance weights $w_k$ drawn from the background-marginalized density $\pi(\theta | \mathbf{d})$. Computation of a single scalar divergence for the n-dimensional joint posterior is straightforwardly given by

$$D_{\text{KL}}(\pi \| p) = \int_S \pi(\theta | \mathbf{d}) \log_2 \left[ \frac{\pi(\theta | \mathbf{d})}{p(\theta)} \right] d\theta \\ = \int_S \pi(\theta | \mathbf{d}) \log_2 \left[ \frac{1}{\mathcal{Z}} \frac{p(\mathbf{d} | \theta) p(\theta)}{p(\theta)} \right] d\theta \\ = -\log_2 \mathcal{Z} + \int_S \pi(\theta | \mathbf{d}) \log_2 L(\theta) d\theta \\ \approx -\log_2 \mathcal{Z} + \sum_k w_k \log_2 L(\theta_k). \quad (13)$$

It is more useful, however, to compute a marginal divergence for each parameter, yielding a handle on which parameters the (marginal) likelihood function is most insensitive to in the context of the (marginal) prior, and thus to which prior assumptions global posterior inferences may be most sensitive. The parameters that exhibit the lowest marginal posterior information gain are those for which prior assumptions are generally more important to be aware of and should thus be an accurate representation of prior belief. Divergence estimation for each parameter (or jointly for $m < n$ parameters) requires more involved post-processing because kernel density estimation is performed for evaluation of the quotient of marginal densities appearing in the integrand.

Note that if the marginal KL divergence for some parameter $\theta$ is small (relative to the divergences of other parameters) but the divergence is close to the theoretical maximum expected divergence, the prior exists in the minimally informative limit while being relatively informative. For example, if the Fisher information for $\theta$ is everywhere relatively small—meaning the experiment is at most a weak probe of the parameter—the marginal posterior for $\theta$ should be dominated even by a minimally informative prior, but posterior inferences about other parameters should be insensitive to all information about $\theta$. If, on the other hand, the likelihood function is a useful probe of $\theta$ while the marginal divergence is relatively small, the global posterior may be sensitive to prior information about $\theta$. In practice, we determine that the parameters that typically exhibit the smallest divergences—pulsar distance $D$ and NICER instrument parameters $\alpha$ and $\gamma$—are assigned marginal priors that are not weakly informative in the context of the likelihood function; moreover, we do not need to calculate Fisher information to understand that the likelihood function itself is a useful probe of these parameters, which exhibit degeneracies with other pulsar parameters. Therefore, we conclude that small divergences in these cases do indicate that our global posterior inferences are strongly conditional on these prior assumptions.

## Appendix B
### Prior Transforms

For nested sampling we aim to transform from a native sampling space—a unit hypercube $\mathcal{H} = [0, 1]^n$—to a (physical) parameter space $\mathbb{R}^n$ according to an inverse transformation of the joint prior density distribution $p(\theta)$ defined on $\mathbb{R}^n$ (typically with some compact support). We thus require implementation of a mapping $\mathbf{x} \mapsto \theta$ where $\mathbf{x} \in \mathcal{H}$ and $\theta \in \mathbb{R}^n$. In this appendix we give the prior transforms implemented in order to facilitate reproduction of the sampling processes, and also to provide a demonstration of some of the necessary architectural work for parameter estimation via nested sampling.

In this work we provide a summary table for each model; as an example, refer to Table 2. Within each table the joint prior density function $q(\theta)$ and its support $S$ are reported. Here we give a prescription for drawing a sample from a given prior: (i) draw a sample for each parameter according to the listed sampling distribution $p(\theta)$ with support $S^\dagger \subset \mathbb{R}^n$—using one-dimensional inverse sampling—to generate a candidate vector $\tilde{\theta} \in \mathbb{R}^n$; (ii) systematically evaluate the constraint equations to determine whether $\tilde{\theta} \in S \subset \mathbb{R}^n$ or whether $S \not\supset \tilde{\theta} \in \mathbb{R}^n$; and (iii) accept the candidate sample if $\tilde{\theta} \in S$, otherwise reject the sample. The form—i.e., the relative marginal density at two values of parameter $\theta$—is, for a subset of parameters, not given by the one-dimensional distribution $\theta \sim q(\theta)$ (whose definition in some cases, as the distance $D$, requires explicitly stated truncation bounds) explicitly written in the prior column: the set of constraint equations defining the joint compact support $S$ often non-trivially modulate the density distribution $q(\theta)$ if $\theta$ appears in constraint equations jointly with some subset of the other parameters. If the joint prior is separable with respect $\theta$ to then by definition $p(\theta) \equiv q(\theta)$.

### B.1. Gravitational Mass and Equatorial Radius

We defined our joint prior distribution of gravitational mass and equatorial radius in Section 2.4.1. We apply a technique from Appendix E of Riley & Watts (2019): inverse sampling of a joint flat density function $q(M, R_{eq})$ with a trivial rectangular boundary $M \in [M_a, M_b]$ and $R_{eq} \in [R_a, R_b]$, and subsequent rejection only if $R_{eq} \notin [3r_g, 16]$ km.<fnref n="93" /> In this case, let $\mathcal{H} = [0, 1] \times [0, 1]$ and let the support of $q(\theta) = q(M, R_{eq})$ be $S^\dagger \subset \mathbb{R}^2$. The mapping is then $\mathcal{H} \to S^\dagger$, $\mathbf{x} \mapsto \theta$, where $S^\dagger \supset S$. For $(M, R_{eq}) \in S$, the joint prior density $p(M, R_{eq}) \propto q(M, R_{eq})$ because the constraint equation in compactness is dependent only on $M$ and $R_{eq}$, and not on any ulterior source parameters. Such a procedure is also summarized in the preamble of Appendix B above.

A standard transform (in the context of the nested sampling software) on the other hand would take the form $\mathcal{H} \to S$, $\mathbf{x} \mapsto \theta$. To construct such a map, one may write $p(M, R_{eq}) = p(R_{eq} | M)p(M)$ where

$$p(M) = \int_{3r_g(M)}^{R_b} p(M, R_{eq}) dR_{eq} \propto R_b - 3r_g(M). \quad (14)$$

<footnote n="93">⁹³ Strictly speaking, we also impose several ulterior constraint equations. However, for the spin of PSR J0030+0451 these constraint equations are unimportant for defining the prior support on the joint space of $M$ and $R_{eq}$.</footnote>
Then define $\mathbf{x} := (x_M, x_R)$ where

$$x_M(M) = \int_{M_a}^{M} p(M')dM' \quad (15)$$

and

$$\begin{aligned} x_R(R_{eq}; M) &= \int_{3r_g(M)}^{R_{eq}} p(R'_{eq} | M) dR'_{eq} \\ &= \int_{3r_g(M)}^{R_{eq}} \frac{p(M, R'_{eq})}{p(M)} dR'_{eq} = \frac{R_{eq} - 3r_g(M)}{R_b - 3r_g(M)}. \end{aligned} \quad (16)$$

Inverting, one has

$$R_{eq}(x_R; M) = x_R R_b + 3(1 - x_R) r_g(M), \quad (17)$$

and similarly one obtains a nonlinear function $M(x_M)$.

The problem with such a standard transformation is clear from inspection of $x_R(R_{eq}; M)$. The common $(M, R_{eq})$ degeneracy in pulse-profile modeling is linear due to sensitivity to compactness $M/R_{eq}$. Therefore, requiring a constant compactness $r_g(M)/R_{eq} = \text{const.}$ implies

$$x_R(R_{eq}; M) = \frac{\text{const.} - 3}{R_b/r_g(M) - 3}, \quad (18)$$

where $M = M(x_M)$, meaning that $x_R = x_R(x_M; M/R_{eq})$ is generally a nonlinear function. An optimal mapping would preserve the linearity of the degeneracy, and thus we do not opt for a standard transformation where $\mathcal{H} \to \mathcal{S}$, $\mathbf{x} \mapsto \mathbf{\theta}$; instead we inverse sample with rejection as described above.

## B.2. Hot Regions

The joint prior distribution for the parameters of the members comprising both hot regions is non-trivial to implement, requiring a number of considerations. Moreover, the difficulty scales with the complexity of the hot regions. Here we break down the implementation into a series of steps.

### B.2.1. Parameter Space

The joint parameter space for the members may be denoted $\mathbf{v} = (\Theta_p, \phi_p, \zeta_p, \Theta_s, \phi_s, \zeta_s, \dots)$, where the six parameters explicitly written are inherent to every model wherein antipodal reflection symmetry is not imposed, and any ulterior parameters depend on the model. For ST-U, these six parameters are sufficient.

### B.2.2. Region-exchange Degeneracy

If the two hot regions sharing the stellar surface are of equivalent complexity, the prior support for the coordinates of the regions can be defined so as to avoid degeneracy of the likelihood function under exchange of the region positions.

For regions related via antipodal reflection symmetry, the prior support for the colatitude of the regions can be defined such that only at (or near to) the support boundary are there configurations wherein a pair of regions mutually map onto one another via a rotation about the stellar rotation axis. The same condition applies if the regions are not related via antipodal reflection symmetry, but have equivalent complexities—e.g., ST-U meaning two ST regions. For instance, for ST-U, we can impose that $\Theta_p \le \Theta_s$, where $\Theta_p$ and $\Theta_s$ are, respectively,

the colatitudes of the centers of the primary and secondary ST regions.

If the regions are not related via antipodal reflection symmetry and do not have equivalent complexities—e.g., ST +CST or ST+EST—then by definition it is not true that an arbitrary point in parameter space yields exactly the same system configuration as distinct point in parameter space due to region exchange. Indeed, region exchange degeneracy may only be exist for a subset of parameter space, and thus one need not impose a joint constraint on the prior support for the colatitudes of the regions.

### B.2.3. Ceding- and Superseding-member Radii

As highlighted in Section 2.5.3, it is advisable to define one's native nested-sampling space in order to linearize certain continuous degeneracies where possible.

The angular extent of the regions are remarkable in this respect for models constituted by at least one region at the CST complexity level or beyond—i.e., when a hot region is constituted by a superseding member and a ceding member, each with an angular radius. While it is not clear how to fully linearize degeneracy of type I illustrated in Figure 5 (see Section 2.5.2), the important<fnref n="94" /> degeneracy of type IV can be linearized by working with the joint space of $\psi$ and $\zeta$, the radii of the superseding and ceding members, respectively (see Section 2.5.3). Although such a choice may be appear obvious in isolation, it may not be the space on which one chooses to intuitively define a joint prior density distribution.

Transforms for a CST or an EST region. If the superseding member subtends smaller angular extent than the ceding member, it is useful to consider $\psi := f\zeta$ where $f \in [\epsilon_f, 1]$ for small (or zero) $\epsilon_f$. Moreover, it is common to invoke uniform prior density distributions for parameters with the intention of choosing a weakly informative prior but without rigorous proof. We therefore consider a flat separable density for $f$ and $\zeta$, such that $q(f, \zeta) = q(f)q(\zeta)$ where $f \sim U(\epsilon_f, 1)$ and $\zeta \sim U(\epsilon_\zeta, \pi/2 - \epsilon_\zeta)$ for small (or zero) $\epsilon_\zeta$. More generally, one might choose $\zeta \in [\epsilon_\zeta, b_\zeta]$; in our case $b_\zeta := \pi/2 - \epsilon_\zeta$. We then require the marginal density function $q(\psi)$ and the conditional density function $q(\zeta | \psi)$ in order to define a map $\mathcal{H} \to \mathcal{S}^\dagger$, where $\mathcal{H} = [0, 1] \times [0, 1]$ and where $\mathcal{S}^\dagger$ indicates that the support $\mathcal{S}$ of $p(\psi, \zeta)$, after all considerations in this appendix (see Appendix B.2.4), will be such that $\mathcal{S}^\dagger \supset \mathcal{S}$.

We must now consider the size of $\epsilon_f$ and $\epsilon_\zeta$: these limits determine the boundary of the prior support in the joint space of $\psi$ and $\zeta$. For a single-temperature hot region, we simply choose $\epsilon_f = 0$, such that either the superseding or ceding member can subtend zero angular extent at the boundary of the support; for a dual-temperature region one might choose a small finite value for $\epsilon_f$ given that it is filled with material of finite temperature. Given the choice $\epsilon_f = 0$ the choice of $\epsilon_\zeta$ is unimportant for deriving the prior distributions of interest.

The joint density $q(\psi, \zeta)$ is given by

$$q(\psi, \zeta) = q(f, \zeta) \left| \frac{\partial f}{\partial \psi} \right| = \zeta^{-1} q(\psi/\zeta) q(\zeta). \quad (19)$$

<footnote n="94">⁹⁴ In terms of prior mass in comparison to that associated with other degenerate structures.</footnote>
$$x_{\zeta}(\zeta; \psi) := \int_{a(\psi)}^{\zeta} q(\zeta' | \psi) d\zeta' \\ = \begin{cases} \frac{\ln(\zeta/\epsilon_{\zeta})}{\ln(b_{\zeta}/\epsilon_{\zeta})} & \text{if } (\psi < \epsilon_{\zeta}) \land (\epsilon_{\zeta} \le \zeta \le b_{\zeta}) \\ \frac{\zeta}{\psi}[1 + \ln(b_{\zeta}/\psi)]^{-1} & \text{if } (\epsilon_{\zeta} \le \psi \le b_{\zeta}) \land (0 \le \zeta \le \psi) \\ [1 + \ln(\zeta/\psi)][1 + \ln(b_{\zeta}/\psi)]^{-1} & \text{if } (\epsilon_{\zeta} \le \psi \le b_{\zeta}) \land (\psi < \zeta \le b_{\zeta}). \end{cases} \quad (32)$$

The function $\zeta(x_\zeta; \psi)$ is written in closed form as

$$\zeta(\mathbf{x}) = \begin{cases} \epsilon_{\zeta} \exp(x_{\zeta} \ln(b_{\zeta}/\epsilon_{\zeta})) & \text{if } (\psi < \epsilon_{\zeta}) \land (\epsilon_{\zeta} \le \zeta \le b_{\zeta}) \\ x_{\zeta} \psi(x_{\psi})[1 + \ln(b_{\zeta}/\psi)] & \text{if } (\epsilon_{\zeta} \le \psi \le b_{\zeta}) \land (0 \le \zeta \le \psi) \\ \psi(x_{\psi}) \exp(x_{\zeta}[1 + \ln(b_{\zeta}/\psi)] - 1) & \text{if } (\epsilon_{\zeta} \le \psi \le b_{\zeta}) \land (\psi < \zeta \le b_{\zeta}). \end{cases} \quad (33)$$

### B.2.4. Non-overlapping Hot Regions

We implicitly define the support $\mathcal{S}$ of the joint density $p(\mathbf{v})$, where $\mathbf{v} = (\Theta_p, \phi_p, \zeta_p, \Theta_s, \phi_s, \zeta_s, \dots)$ is the vector of parameters controlling both hot regions sharing the stellar surface, by imposing a constraint equation in terms of $\mathbf{v}$: we require that the regions are *non-overlapping*. More explicitly: two radiating regions associated with distinct regions cannot overlap—if they were to overlap, additional logical conditions would be required to specify an order of precedence for intensity evaluation at spacetime events at the stellar surface.

For single-temperature regions with a single member—i.e., simply connected and circular—one need only determine whether the ceding members, with their simple (circular) boundaries, overlap; the same condition is true if there exists a superseding member that is a *hole* in a ceding member—i.e., ST, CST, or EST. For dual-temperature hot regions, whose boundary is always constituted by a maximum of two simple (circular) boundaries, one need only determine whether any pair of members—from two distinct regions—overlap. However, for a PST region the boundary of the radiating region is more unwieldy for evaluating whether or not the radiating region of the PST region overlaps with another region sharing the stellar surface; in this case we simply define an overlap condition only in terms of the simple boundary of the ceding member, a subset of which is superseded by non-radiating surface.

To derive the joint prior density $p(\mathbf{v})$ and its support $\mathcal{S}$, we (i) define $q(\mathbf{v})$ as a product of density functions, (ii) inverse sample as $\mathcal{H} \to \mathcal{S}^\dagger$, $\mathbf{x} \mapsto \mathbf{v}$, and (iii) accept the sample $\tilde{\mathbf{v}}$ if the regions are determined to be non-overlapping. Algorithmically, we identify a set of pairs of members with simple boundaries whose center coordinates and angular radii are defined by vector $\tilde{\mathbf{v}}$, and evaluate via a spherical coordinate transformation whether or not the angular separation of the member centers is at least equal to the sum of the angular radii. If any of the pairs of members overlap, then $\tilde{\mathbf{v}} \notin \mathcal{S}$.

As an example, for ST+EST, $\mathbf{v} = (\Theta_p, \phi_p, \zeta_p, \Theta_s, \phi_s, \zeta_s, \epsilon_s, \phi_s)$ where the coordinates of the center of the ceding member of the EST region is derived from the vector $(\Theta_s, \phi_s, \psi_s, \zeta_s, \epsilon_s, \phi_s)$. Overlap is then evaluated for this ceding member in relation to the ST spot whose boundary is derived from the vector $(\Theta_p, \phi_p, \zeta_p)$. The marginal density function for every parameter constituting (the space of) vector $\mathbf{v}$ is thus modulated by excluding overlaps.

## Appendix C
### Model Summary Tables

In this appendix we provide posterior summary tables for all models applied to the PSR J0030+0451 event data, other than ST+PST (Table 2). For ST+EST see Table 3. For ST+CST see Table 4. For ST-U see Table 5. For ST-S see Table 6. For CDT-U see Table 7.
<table><caption>**Table 3**
Summary Table for ST+EST, Introduced in Section 2.5.7 and Illustrated in Figure 11</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87, fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass</td><td>M ~ U(1, 3)</td><td>1.46<sup>+0.17</sup><sub>-0.18</sub></td><td>1.15<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius</td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1), 16]</td><td>13.89<sup>+1.14</sup><sub>-1.30</sub></td><td>0.69<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>Θ<sub>p</sub> (radians)</td><td>ST region center colatitude</td><td>Θ<sub>p</sub> ~ U(0, π)</td><td>2.22<sup>+0.09</sup><sub>-0.10</sub></td><td>3.01<sup>+0.03</sup><sub>-0.03</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>a</sup></td><td>0.45<sup>+0.00</sup><sub>-0.00</sub></td><td>6.59<sup>+0.02</sup><sub>-0.03</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.07<sup>+0.01</sup><sub>-0.01</sub></td><td>4.78<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>EST region hole center colatitude</td><td>Θ<sub>s</sub> ~ U(0, π)</td><td>2.66<sup>+0.07</sup><sub>-0.09</sub></td><td>3.36<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>EST region initial phase (from Earth antipode)</td><td>φ<sub>s</sub> ~ U(-0.5, 0.5), periodic</td><td>-0.51<sup>+0.01</sup><sub>-0.01</sub></td><td>4.28<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>ψ<sub>s</sub><sup>+</sup> (radians)</td><td>EST region angular radii difference</td><td>ζ<sub>s</sub> ~ U(0, π/2), ψ<sub>s</sub><sup>+</sup> := ζ<sub>s</sub> - ψ<sub>s</sub></td><td>0.03<sup>+0.01</sup><sub>-0.01</sub></td><td>2.93<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>ψ<sub>s</sub> (radians)</td><td>EST region hole angular radius</td><td>f<sub>s</sub> ~ U(0, 1), ψ<sub>s</sub> := f<sub>s</sub>ζ<sub>s</sub></td><td>0.25<sup>+0.04</sup><sub>-0.04</sub></td><td>2.29<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>ε<sub>s</sub></td><td>EST region annulus fractional angular offset</td><td>ε<sub>s</sub> ~ U(0, 1)<sup>b</sup></td><td>0.45<sup>+0.45</sup><sub>-0.33</sub></td><td>0.08<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>φ<sub>s</sub> (radians)</td><td>EST region annulus azimuthal offset</td><td>φ<sub>s</sub> ~ U(0, 2π)</td><td>0.79<sup>+0.35</sup><sub>-0.26</sub></td><td>0.64<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td></td><td>Non-overlapping hot-region annuli</td><td>function of (Θ<sub>p</sub>, Θ<sub>s</sub>, φ<sub>p</sub>, φ<sub>s</sub>, ζ<sub>p</sub>, ζ<sub>s</sub>)</td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</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.03<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>EST 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.60<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>i (radians)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>1.01<sup>+0.07</sup><sub>-0.07</sub></td><td>2.40<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>D (kpc)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)</td><td>0.33<sup>+0.01</sup><sub>-0.01</sub></td><td>0.02<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>N<sub>H</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>0.61<sup>+0.18</sup><sub>-0.16</sub></td><td>2.84<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.09</sup><sub>-0.09</sub></td><td>0.02<sup>+0.01</sup><sub>-0.00</sub></td></tr><tr><td>NICER β</td><td>Reference-to-calibrated matrix weighting</td><td>β ~ U(0, 1)</td><td>0.16<sup>+0.17</sup><sub>-0.11</sub></td><td>1.07<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>NICER γ</td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>0.99<sup>+0.09</sup><sub>-0.09</sub></td><td>0.04<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td colspan="5"><strong>Sampling process information</strong></td></tr><tr><td colspan="5">Number of free parameters: 19</td></tr><tr><td colspan="5">Number of runs: 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: 0.3</td></tr><tr><td colspan="5">Termination condition: 10<sup>-1</sup></td></tr><tr><td colspan="5">Evidence:<sup>c</sup> ln Z̄ = -36367.81<sup>+0.48</sup><sub>-0.43</sub></td></tr><tr><td colspan="5">Global KL divergence: D̂<sub>KL</sub> = 62.1<sup>+0.8</sup><sub>-0.8</sub> bits</td></tr><tr><td colspan="5">Number of core<sup>d</sup> hours: 61210</td></tr><tr><td colspan="5">Likelihood evaluations: 88965106</td></tr><tr><td colspan="5">Nested replacements: 53149</td></tr><tr><td colspan="5">Weighted posterior samples: 17671</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Where $\phi_p = a$ is an arbitrary phase dependent on event data pre-processing. We set $a = 0.35$.</footnote>

<footnote><sup>b</sup> If we were to parameterize the eccentricity in terms of the sum of angular radii of the superseding (hole) and ceding regions, $\vartheta := \varepsilon(1 + f)\zeta$, a conditional prior such as $\varepsilon|f \sim U(0, (1-f)/(1+f))$ would be necessary, where the upper bound imposes that the radiating region is not simply connected (i.e., is an annulus): $\varepsilon(1+f)\zeta + f\zeta \le \zeta \Rightarrow \varepsilon \le (1-f)/(1+f)$.</footnote>

<footnote><sup>c</sup> Defined as the prior predictive probability $p(d | ST + EST)$. We report the interval about the median containing ±45% of 10³ joint bootstrap-weight replications for the combined run.</footnote>

<footnote><sup>d</sup> Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) processors on the SURFsara Cartesius supercomputer.</footnote>
<table><caption>**Table 4**
Summary Table for ST+CST, Introduced in Section 2.5.7</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87, fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass</td><td>M ~ U(1, 3)</td><td>1.44<sup>+0.18</sup><sub>-0.19</sub></td><td>1.07<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius</td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1, 16)]</td><td>13.88<sup>+1.23</sup><sub>-1.38</sub></td><td>0.61<sup>+0.01</sup><sub>-0.00</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>Θ<sub>p</sub> (radians)</td><td>ST region center colatitude</td><td>Θ<sub>p</sub> ~ U(0, π)</td><td>2.24<sup>+0.09</sup><sub>-0.09</sub></td><td>5.14<sup>+0.03</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(-0.5, 0.5), periodic</td><td>0.46<sup>+0.00</sup><sub>-0.00</sub></td><td>6.66<sup>+0.03</sup><sub>-0.02</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.07<sup>+0.01</sup><sub>-0.01</sub></td><td>4.59<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>CST region center colatitude</td><td>Θ<sub>s</sub> ~ U(0, π)</td><td>2.60<sup>+0.05</sup><sub>-0.06</sub></td><td>2.79<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>CST region initial phase (from Earth antipode)</td><td>φ<sub>s</sub> ~ U(-0.5, 0.5), periodic</td><td>-0.50<sup>+0.00</sup><sub>-0.00</sub></td><td>7.33<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>ψ<sub>s</sub><sup>+</sup> (radians)</td><td>CST region annulus angular width</td><td>ζ<sub>s</sub> ~ U(0, π/2), ψ<sub>s</sub><sup>+</sup> := ζ<sub>s</sub> - ψ<sub>s</sub></td><td>0.04<sup>+0.01</sup><sub>-0.01</sub></td><td>2.71<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>ψ<sub>s</sub> (radians)</td><td>CST region hole angular radius</td><td>f<sub>s</sub> ~ U(0, 1), ψ<sub>s</sub> := f<sub>s</sub>ζ<sub>s</sub></td><td>0.23<sup>+0.03</sup><sub>-0.03</sub></td><td>2.52<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td></td><td>Enforce ST and CST colatitude order<sup>a</sup></td><td>Θ<sub>s</sub> ≥ Θ<sub>p</sub></td><td></td><td></td></tr><tr><td></td><td>Non-overlapping hot regions</td><td>function of (Θ<sub>p</sub>, Θ<sub>s</sub>, φ<sub>p</sub>, φ<sub>s</sub>, ζ<sub>p</sub>, ζ<sub>s</sub>)</td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</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.10<sup>+0.01</sup><sub>-0.01</sub></td><td>4.99<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>CST 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.58<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>i (radians)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>1.02<sup>+0.07</sup><sub>-0.08</sub></td><td>2.32<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>D (kpc)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)</td><td>0.33<sup>+0.01</sup><sub>-0.01</sub></td><td>0.01<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td>N<sub>H</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>0.62<sup>+0.19</sup><sub>-0.18</sub></td><td>2.75<sup>+0.02</sup><sub>-0.02</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.10</sup><sub>-0.10</sub></td><td>0.01<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td>NICER β</td><td>Reference-to-calibrated matrix weighting</td><td>β ~ U(0, 1)</td><td>0.16<sup>+0.18</sup><sub>-0.11</sub></td><td>0.98<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>NICER γ</td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>0.99<sup>+0.09</sup><sub>-0.09</sub></td><td>0.01<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td colspan="5"><strong>Sampling process information</strong></td></tr><tr><td colspan="5">Number of free parameters: 17</td></tr><tr><td colspan="5">Number of runs: 2</td></tr><tr><td colspan="5">Number of live points per run: 10<sup>3</sup></td></tr><tr><td colspan="5">Inverse hypervolume expansion factor: 0.3</td></tr><tr><td colspan="5">Termination condition: 10<sup>-1</sup></td></tr><tr><td colspan="5">Evidence:<sup>b</sup> ln Z̄ = -36368.00<sup>+0.34</sup><sub>-0.33</sub></td></tr><tr><td colspan="5">Global KL divergence: D̂<sub>KL</sub> = 62.9<sup>+0.6</sup><sub>-0.5</sub> bits</td></tr><tr><td colspan="5">Combined number of core<sup>c</sup> hours: 23010</td></tr><tr><td colspan="5">Combined likelihood evaluations: 39501475</td></tr><tr><td colspan="5">Combined nested replacements: 105264</td></tr><tr><td colspan="5">Combined weighted posterior samples: 32839</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Based on learning that additional complexity (of the form we consider in our model space) beyond ST is not warranted for the one region, but is warranted for the other region.</footnote>

<footnote><sup>b</sup> Defined as the prior predictive probability p(d | ST + CST). We report the interval about the median containing ±45% of 10³ joint bootstrap-weight replications for the combined run.</footnote>

<footnote><sup>c</sup> Approximate equal-partition between Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) and E5-2690v3 (2.60 GHz; Haswell) processors on the SURFsara Cartesius supercomputer.</footnote>
<table><caption>**Table 5**
Summary Table for ST-U, Introduced in Section 2.5.1 and Illustrated in Figure 3</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87, fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass</td><td>M ~ U(1, 3)</td><td>1.09<sup>+0.11</sup><sub>-0.07</sub></td><td>2.09<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius</td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1, 16)]</td><td>10.44<sup>+1.10</sup><sub>-0.86</sub></td><td>1.44<sup>+0.03</sup><sub>-0.00</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>Θ<sub>p</sub> (radians)</td><td>p region center colatitude</td><td>Θ<sub>p</sub> ~ U(0, π)</td><td>2.48<sup>+0.06</sup><sub>-0.06</sub></td><td>6.92<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>s region center colatitude</td><td>Θ<sub>s</sub> ~ U(0, π)</td><td>2.78<sup>+0.02</sup><sub>-0.02</sub></td><td>4.08<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>φ<sub>p</sub> (cycles)</td><td>p region initial phase (from Earth)</td><td>φ<sub>p</sub> ~ U(-0.5, 0.5), periodic</td><td>0.46<sup>+0.00</sup><sub>-0.00</sub></td><td>7.51<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>s region initial phase (from Earth antipode)</td><td>φ<sub>s</sub> ~ U(-0.5, 0.5), periodic</td><td>-0.50<sup>+0.00</sup><sub>-0.00</sub></td><td>8.05<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>ζ<sub>p</sub> (radians)</td><td>p region angular radius</td><td>ζ<sub>p</sub> ~ U(0, π/2)</td><td>0.14<sup>+0.02</sup><sub>-0.02</sub></td><td>3.95<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>ζ<sub>s</sub> (radians)</td><td>s region angular radius</td><td>ζ<sub>s</sub> ~ U(0, π/2)</td><td>0.29<sup>+0.04</sup><sub>-0.03</sub></td><td>2.98<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td></td><td>Eliminate region-exchange degeneracy</td><td>Θ<sub>s</sub> ≥ Θ<sub>p</sub></td><td></td><td></td></tr><tr><td></td><td>Non-overlapping hot regions</td><td>function of (Θ<sub>p</sub>, Θ<sub>s</sub>, φ<sub>p</sub>, φ<sub>s</sub>, ζ<sub>p</sub>, ζ<sub>s</sub>)</td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</sub>(K))</td><td>p 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.53<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub>(K))</td><td>s region NSX effective temperature</td><td>log<sub>10</sub>T<sub>s</sub> ~ U(5.1, 6.8), NSX limits</td><td>6.10<sup>+0.01</sup><sub>-0.01</sub></td><td>5.71<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>i(radians)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>1.04<sup>+0.07</sup><sub>-0.08</sub></td><td>2.38<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>D (kpc)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)</td><td>0.33<sup>+0.01</sup><sub>-0.01</sub></td><td>0.35<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>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>1.23<sup>+0.17</sup><sub>-0.17</sub></td><td>2.82<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>NICER α</td><td>Calibrated matrix scaling</td><td>α ~ N(1, 0.1), α ∈ [0.5, 1.5]</td><td>0.96<sup>+0.10</sup><sub>-0.10</sub></td><td>0.11<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.23<sup>+0.23</sup><sub>-0.16</sub></td><td>0.57<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>NICER γ</td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>0.91<sup>+0.10</sup><sub>-0.09</sub></td><td>0.56<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td colspan="5"><strong>Sampling process information</strong></td></tr><tr><td colspan="5">Number of free parameters: 16</td></tr><tr><td colspan="5">Number of runs:<sup>a</sup> 2</td></tr><tr><td colspan="5">Number of live points per run: 10<sup>3</sup></td></tr><tr><td colspan="5">Inverse hypervolume expansion factor: 0.3</td></tr><tr><td colspan="5">Termination condition: 10<sup>-3</sup></td></tr><tr><td colspan="5">Evidence:<sup>b</sup> ln Z̄ = -36377.60<sup>+0.36</sup><sub>-0.35</sub></td></tr><tr><td colspan="5">Global KL divergence: D̂<sub>KL</sub> = 63.7<sup>+0.5</sup><sub>-0.5</sub> bits</td></tr><tr><td colspan="5">Combined number of core<sup>c</sup> hours: 9588</td></tr><tr><td colspan="5">Combined likelihood evaluations: 25346841</td></tr><tr><td colspan="5">Combined nested replacements: 121617</td></tr><tr><td colspan="5">Combined weighted posterior samples: 49481</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> The mode-separation MULTINEST variant was deactivated for these two runs that were combined to compute estimators. Mode separation means that modes are not evolved independently and nested sampling threads contact multiple modes; a mode-separation run was executed and is displayed in Figure 19 and in the figure sets available in the online corresponding to Figures 20 and 21. The theory nor software implementation exists for combining this run with the two reported in the table.</footnote>

<footnote><sup>b</sup> Defined as the prior predictive probability p(d | ST - U). We report the interval about the median containing ±45% of 10³ joint bootstrap-weight replications for the combined run. See the footnote in Table 2.</footnote>

<footnote><sup>c</sup> Approximate equal-partition between Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) and E5-2690v3 (2.60 GHz; Haswell) processors on the SURFsara Cartesius supercomputer.</footnote>
<table><caption>**Table 6**
Summary Table for ST-S, Introduced in Section 2.5.1 and Illustrated in Figure 3</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87, fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass</td><td>M ~ U(1, 3)</td><td>2.93<sup>+0.01</sup><sub>-0.01</sub></td><td>6.75<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius</td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1), 16]</td><td>15.97<sup>+0.02</sup><sub>-0.04</sub></td><td>6.11<sup>+0.04</sup><sub>-0.04</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>Θ<sub>p</sub> (radians)</td><td>p region center colatitude</td><td>Θ<sub>p</sub> ~ U(0, π/2)</td><td>1.26<sup>+0.02</sup><sub>-0.02</sub></td><td>4.07<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>φ<sub>p</sub> (cycles)</td><td>p region initial phase (from Earth)</td><td>φ<sub>p</sub> ~ U(-0.5, 0.5), periodic</td><td>-0.09<sup>+0.00</sup><sub>-0.00</sub></td><td>8.57<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>ζ (radians)</td><td>ST region angular radius (shared)</td><td>ζ ~ U(0, π/2)</td><td>0.09<sup>+0.01</sup><sub>-0.00</sub></td><td>6.55<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>log<sub>10</sub>(T(K))</td><td>ST region NSX effective temperature (shared)</td><td>log<sub>10</sub>T ~ U(5.1, 6.8), NSX limits</td><td>6.08<sup>+0.00</sup><sub>-0.00</sub></td><td>6.97<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>s region center colatitude</td><td>Θ<sub>s</sub> = π - Θ<sub>p</sub>, derived</td><td></td><td></td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>s region initial phase</td><td>φ<sub>s</sub> = φ<sub>p</sub> + 0.5, derived</td><td></td><td></td></tr><tr><td>i (radians)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>1.23<sup>+0.02</sup><sub>-0.03</sub></td><td>3.98<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>D (kpc)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)</td><td>0.32<sup>+0.01</sup><sub>-0.01</sub></td><td>0.02<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td>N<sub>H</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>0.02<sup>+0.03</sup><sub>-0.02</sub></td><td>5.93<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>NICER α</td><td>Calibrated matrix scaling</td><td>α ~ N(1, 0.1), α ∈ [0.5, 1.5]</td><td>1.00<sup>+0.09</sup><sub>-0.09</sub></td><td>0.01<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td>NICER β</td><td>Reference-to-calibrated matrix weighting</td><td>β ~ U(0, 1)</td><td>0.36<sup>+0.22</sup><sub>-0.19</sub></td><td>0.35<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>NICER γ</td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>1.01<sup>+0.09</sup><sub>-0.09</sub></td><td>0.03<sup>+0.00</sup><sub>-0.00</sub></td></tr><tr><td colspan="5">Sampling process information</td></tr><tr><td colspan="5">Number of free parameters: 12</td></tr><tr><td colspan="5">Number of runs: 2</td></tr><tr><td colspan="5">Number of live points per run: 10<sup>3</sup></td></tr><tr><td colspan="5">Inverse hypervolume expansion factor: 0.3</td></tr><tr><td colspan="5">Termination condition: 10<sup>-3</sup></td></tr><tr><td colspan="5">Evidence:<sup>a</sup> ln Z̄ = -37211.71<sup>+0.29</sup><sub>-0.31</sub></td></tr><tr><td colspan="5">Global KL divergence: D̂<sub>KL</sub> = 51.3<sup>+0.5</sup><sub>-0.5</sub> bits</td></tr><tr><td colspan="5">Combined number of core<sup>b</sup> hours: 1494</td></tr><tr><td colspan="5">Combined likelihood evaluations: 1666483</td></tr><tr><td colspan="5">Combined nested replacements: 87488</td></tr><tr><td colspan="5">Combined weighted posterior samples: 30202</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Defined as the prior predictive probability $p(d | ST-S)$. We report the interval about the median containing ±45% of 10³ joint bootstrap-weight replications for the combined run. See the footnote in Table 2.</footnote>

<footnote><sup>b</sup> Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) processors on the SURFsara Cartesius supercomputer.</footnote>
<table><caption>**Table 7**
Summary Table for CDT-U, Introduced in Section 2.5.3 and Illustrated in Figure 6</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P (ms)</td><td>Coordinate spin period</td><td>P = 4.87,fixed</td><td>...</td><td>...</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>Gravitational mass</td><td>M ~ U(1, 3)</td><td>1.44<sup>+0.17</sup><sub>-0.18</sub></td><td>1.15</td></tr><tr><td>R<sub>eq</sub> (km)</td><td>Coordinate equatorial radius</td><td>R<sub>eq</sub> ~ U[3r<sub>g</sub>(1, 16)]</td><td>13.86<sup>+1.16</sup><sub>-1.26</sub></td><td>0.68</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>Θ<sub>p</sub> (radians)</td><td>p region center colatitude</td><td>Θ<sub>p</sub> ~ U(0, π)</td><td>2.24<sup>+0.08</sup><sub>-0.08</sub></td><td>4.68</td></tr><tr><td>φ<sub>p</sub> (cycles)</td><td>p region initial phase (from Earth)</td><td>φ<sub>p</sub> ~ U(-0.5, 0.5), periodic</td><td>0.46<sup>+0.00</sup><sub>-0.00</sub></td><td>6.75</td></tr><tr><td>ψ<sub>p</sub><sup>+</sup> (radians)</td><td>p region annulus angular width</td><td>ζ<sub>p</sub> ~ U(0, π/2), ψ<sub>p</sub><sup>+</sup> := ζ<sub>p</sub> - ψ<sub>p</sub></td><td>0.08<sup>+0.16</sup><sub>-0.06</sub></td><td>2.59</td></tr><tr><td>ψ<sub>p</sub><sup>-</sup> (radians)</td><td>p region hole angu- lar radius</td><td>f<sub>p</sub> ~ U(0, 1), ψ<sub>p</sub><sup>-</sup> := f<sub>p</sub> ζ<sub>p</sub></td><td>0.07<sup>+0.01</sup><sub>-0.01</sub></td><td>5.28</td></tr><tr><td>Θ<sub>s</sub> (radians)</td><td>s region center colatitude</td><td>Θ<sub>s</sub> ~ U(0, π)</td><td>2.61<sup>+0.05</sup><sub>-0.06</sub></td><td>2.98</td></tr><tr><td>φ<sub>s</sub> (cycles)</td><td>s region initial phase (from Earth antipode)</td><td>φ<sub>s</sub> ~ U(-0.5, 0.5), periodic</td><td>-0.50<sup>+0.00</sup><sub>-0.00</sub></td><td>7.38</td></tr><tr><td>ψ<sub>s</sub><sup>+</sup> (radians)</td><td>s region annulus angular width</td><td>ζ<sub>s</sub> ~ U(0, π/2), ψ<sub>s</sub><sup>+</sup> := ζ<sub>s</sub> - ψ<sub>s</sub><sup>+</sup></td><td>0.04<sup>+0.01</sup><sub>-0.01</sub></td><td>5.23</td></tr><tr><td>ψ<sub>s</sub><sup>-</sup> (radians)</td><td>s region hole angular radius</td><td>f<sub>s</sub> ~ U(0, 1), ψ<sub>s</sub><sup>-</sup> := f<sub>s</sub> ζ<sub>s</sub></td><td>0.23<sup>+0.03</sup><sub>-0.03</sub></td><td>3.66</td></tr><tr><td></td><td>Eliminate region-exchange degeneracy</td><td>Θ<sub>s</sub> ≥ Θ<sub>p</sub></td><td></td><td></td></tr><tr><td></td><td>Non-overlapping hot-region annuli</td><td>function of (Θ<sub>p</sub>, Θ<sub>s</sub>, φ<sub>p</sub>, φ<sub>s</sub>, ζ<sub>p</sub>, ζ<sub>s</sub>)</td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</sub> (K))</td><td>p region annulus NSX effective temperature</td><td>log<sub>10</sub>T<sub>p</sub> ~ U(5.1, 6.8), NSX limits</td><td>5.43<sup>+0.23</sup><sub>-0.20</sub></td><td>1.11</td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>p region hole 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.06</td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>s region annulus 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.56</td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> (K))</td><td>s region hole NSX effective temperature</td><td>log<sub>10</sub>T<sub>s</sub> ~ U(5.1, 6.8), NSX limits</td><td>5.47<sup>+0.19</sup><sub>-0.23</sub></td><td>1.30</td></tr><tr><td>i (radians)</td><td>Earth inclination to rotation axis</td><td>i ~ U(0, π/2)</td><td>1.02<sup>+0.07</sup><sub>-0.07</sub></td><td>2.37</td></tr><tr><td>D (kpc)</td><td>Earth distance</td><td>D ~ N(0.325, 0.009)</td><td>0.33<sup>+0.01</sup><sub>-0.01</sub></td><td>0.02</td></tr><tr><td>N<sub>H</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>Interstellar neutral H column density</td><td>N<sub>H</sub> ~ U(0, 5)</td><td>0.70<sup>+0.19</sup><sub>-0.18</sub></td><td>2.70</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.09</sup><sub>-0.09</sub></td><td>0.02</td></tr><tr><td>NICER β</td><td>Reference-to-calibrated matrix weighting</td><td>β ~ U(0, 1)</td><td>0.16<sup>+0.17</sup><sub>-0.11</sub></td><td>1.05</td></tr><tr><td>NICER γ</td><td>Reference matrix scaling</td><td>γ ~ N(1, 0.1), γ ∈ [0.5, 1.5]</td><td>0.99<sup>+0.09</sup><sub>-0.09</sub></td><td>0.03</td></tr><tr><td></td><td>Sampling process information</td><td></td><td></td><td></td></tr><tr><td colspan="5">Number of free parameters: 20</td></tr><tr><td colspan="5">Number of runs: 1</td></tr><tr><td colspan="5">Number of live points:<sup>a</sup> 10<sup>3</sup></td></tr><tr><td colspan="5">Inverse hypervolume expansion factor:<sup>b</sup> 0.3</td></tr><tr><td colspan="5">Termination condition: 10<sup>-1</sup></td></tr><tr><td colspan="5">Evidence:<sup>c</sup> ln Z̄ = -36366.76 ± 0.21</td></tr><tr><td colspan="5">Global KL divergence: D̂<sub>KL</sub> = 64.0 bits</td></tr><tr><td colspan="5">Combined number of core<sup>d</sup> hours: 101917</td></tr><tr><td colspan="5">Likelihood evaluations: 156707329</td></tr><tr><td colspan="5">Nested replacements: 54610</td></tr><tr><td colspan="5">Weighted posterior samples: 17503</td></tr></tbody></table>

**Notes.** In this table we only give the numerical details for one completed run. However, we attempted to perform a higher-efficiency higher-resolution second run, without activation of the mode-separation MULTINEST sampling algorithm. Our attempt to improve the sampling efficiency by linearizing the degeneracy (type IV in Figures 5 and 7) observed in the first run; however, we failed to fully apply the necessary transformations (later applied to models ST+CST and beyond, as described in Section 2.5.3 and Appendix B.2), and the mapping from the native sampling space to the physical parameter space inadvertently preserved the nonlinearity of the degeneracy. We thus did not attain higher efficiency, which coupled with higher-resolution calculation, meant that this run was nearing—but did not reach—termination according to the standard criterion used for the other sampling processes reported in this work. Due to the low sampling efficiencies being reported, we ceased computation after consumption of ~160540 core hours in order to preserve resources and redesign our modeling route as described in Sections 2.5.4 and 2.5.7. The highest-likelihood points from the posterior mode in this second run were however utilized, in combination with those from the first run, to map out all of the degenerate posterior structure in Figures 7 and 8. The inefficiency suffered during CDT-U posterior computation served as a stark reminder that we should design our problems as carefully as possible in order to avoid resource wastage, and was the motivation behind considering in detail the question of how much hot-region complexity is helpful. While we could in principle conclude this adjourned run with additional computing resources, we have argued in this work that would not be fruitful do so.

<sup>a</sup> The mode-separation MULTINEST variant was activated, meaning that isolated modes are evolved independently and nested sampling threads migrate between multiple modes. A local posterior mode was identified, corresponding to a weaker phase solution in which the primary and secondary hot regions transpose—relative to the global posterior mode—in their coupling to pulse components visible in the phase-folded event data. The number of live points locked into the mode with dominant posterior mass was 637, a number assigned according to the prior mass distribution upon mode separation and under the influence of Monte Carlo noise. The posterior mass ratio (or ratio of local evidences) is estimated to be ~1100. As stated in Appendix A.2.3, because of activation of mode-separation sampling, we cannot perform error analysis via process bootstrapping, and we would not be able to combine with another run, supposing that one was available. We nevertheless have an error on the log-evidence reported natively by MULTINEST. Another consequence of the activation of mode-separation is that sampling resolution was absorbed by the local posterior mode with much lower mass—in other words, the active points were sparser in the dominant mode. In combination with sampling error due to likelihood isosurface nonlinearity, the consequence was that the dominant and degenerate posterior mode—which forms a large connected structure in parameter space as discussed in Section 2.5.4—was not fully resolved. In particular, the type I degeneracy branch (refer to Figures 5 and 7) was not fully resolved, with sampling threads (active points) migrating to and densely populating the type IV branch (again refer to Figure 7) in an unbalanced manner. The second, albeit incomplete, run described in the caption of this table exhibited much improved resolution of the dominant mode (see the points in Figure 7); the resolution remained incomplete, however, which we attribute to sampling error due to the clear nonlinear degeneracy present in the mode.

<sup>b</sup> We decreased the expansion factor to 0.8<sup>-1</sup> when the acceptance rate decayed to below 5 × 10<sup>-4</sup>, which slightly decreased the rate of decay. At this point the process was sampling from the typical set and nearing termination (2 × 10<sup>4</sup> core hours remaining at low acceptance fraction).

<sup>c</sup> Defined as the prior predictive probability p(d | CDT - U).

<sup>d</sup> Intel® Xeon E5-2697Av4 (2.60 GHz; Broadwell) processors on the SURFsara Cartesius supercomputer.
<figure><img src="image_27.png" /><figcaption>**Figure 24.** Schematic diagrams of models wherein an ST region shares the stellar surface with a higher-complexity EDT or PDT region. The complete figure set (3 images) is available in the online journal.</figcaption></figure>

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

# Appendix D
## Supplementary Ideas for Model Extension

We now make note of models within the scope of this work that were either clearly (without need for explicit posterior computation due to the posterior properties of simpler models) not competitive for PSR J0030+0451 or unhelpfully complex in phenomenologically describing the structure of a hot region. Excess complexity does not mean that the modeling has been optimized, but indicates that a particular extension to a model is not warranted because we are insensitive to a subset of parameters (or combinations of parameters), *and* that the

signals that maximize the likelihood function are signals that are effectively generated by a simpler (nested) model.

*ST+EDT* and *ST+PDT*. Obtained via simple extension of *ST+EST* and *ST+PST*: let the superseding member (the hole for *EST*) contain radiating material. Refer to Figure 24.

*EDT-S* and *EDT-U*. Refer to the online figure set associated with Figure 24. For *EDT-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. Conversely, for *EDT-U*, the secondary region is endowed with distinct parameters—i.e., it is *not* derived from the primary
region under antipodal symmetry. However, the parameters of the secondary region have an otherwise equivalent meaning—in terms of surface radiation field specification—to their primary-region counterparts. As an example, the azimuth of the center of the secondary ceding member is defined relative to the meridian passing through the center of the secondary hole and through the rotational poles.

**PDT-Sand PDT-U.** Refer to the online figure set associated with Figure 24. For PDT-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. Conversely, for PDT-U, the secondary region is endowed with distinct parameters—i.e., it is *not* derived from the primary region under antipodal symmetry. However, the parameters of the secondary region have an otherwise equivalent meaning—in terms of surface radiation field specification—to their primary-region counterparts.

## ORCID iDs

T. E. Riley https://orcid.org/0000-0001-9313-0493

A. L. Watts https://orcid.org/0000-0002-1009-2354

S. Bogdanov https://orcid.org/0000-0002-9870-2742

P. S. Ray https://orcid.org/0000-0002-5297-5278

R. M. Ludlam https://orcid.org/0000-0002-8961-939X

S. Guillot https://orcid.org/0000-0002-6449-106X

A. V. Bilous https://orcid.org/0000-0002-7177-6987

D. Chakrabarty https://orcid.org/0000-0001-8804-8946

A. K. Harding https://orcid.org/0000-0001-6119-859X

W. C. G. Ho https://orcid.org/0000-0002-6089-6836

S. M. Morsink https://orcid.org/0000-0003-4357-0575

T. E. Strohmayer https://orcid.org/0000-0001-7681-5845

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, PhRvL, 121, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, PhRvX, 9, 011001

Alford, M., & Sedrakian, A. 2017, PhRvL, 119, 161104

Alford, M. G., & Han, S. 2016, EPJA, 52, 62

AlGendy, M., & Morsink, S. M. 2014, ApJ, 791, 78

Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, PhRvL, 120, 172703

Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Sci, 340, 448

Ardeljan, N. V., Bisnovatyi-Kogan, G. S., & Moiseenko, S. G. 2005, MNRAS, 359, 333

Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37

Arzoumanian, Z., Gendreau, K. C., Baker, C. L., et al. 2014, Proc. SPIE, 9144, 914420

Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458

Baillot d'Etivaux, N., Guillot, S., Margueron, J., et al. 2019, arXiv:1905.01081

Barnard, J. J., & Arons, J. 1982, ApJ, 254, 713

Bauböck, M., Berti, E., Psaltis, D., & Özel, F. 2013, ApJ, 777, 68

Bauböck, M., Özel, F., Psaltis, D., & Morsink, S. M. 2015, ApJ, 799, 22

Bauböck, M., Psaltis, D., & Özel, F. 2019, ApJ, 872, 162

Baym, G., Hatsuda, T., Kojo, T., et al. 2018, RPPh, 81, 056902

Becker, W., Trümper, J., Lommen, A. N., & Backer, D. C. 2000, ApJ, 545, 1015

Behnel, S., Bradshaw, R., Citro, C., et al. 2011, CSE, 13, 31

Bildsten, L., Salpeter, E. E., & Wasserman, I. 1992, ApJ, 384, 143

Bilous, A. V., Watts, A. L., Harding, A. K., et al. 2019, ApJL, 887, L23

Blaes, O. M., Blandford, R. D., Madau, P., & Yan, L. 1992, ApJ, 399, 634

Bogdanov, S. 2013, ApJ, 762, 96

Bogdanov, S., & Grindlay, J. E. 2009, ApJ, 703, 1557

Bogdanov, S., Grindlay, J. E., & Rybicki, G. B. 2008, ApJ, 689, 407

Bogdanov, S., Guillot, S., Ray, P., et al. 2019a, ApJL, 887, L25

Bogdanov, S., Guillot, S., Ray, P., et al. 2019c, Constraining the Neutron Star Mass-Radius Relation and Dense Matter Equation of State with NICER. I. The Millisecond Pulsar X-Ray Data Set, Zenodo, doi:10.5281/zenodo. 35244457

Bogdanov, S., Heinke, C. O., Özel, F., & Güver, T. 2016, ApJ, 831, 184

Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019b, ApJL, 887, L26

Bogdanov, S., Rybicki, G. B., & Grindlay, J. E. 2007, ApJ, 670, 668

Bchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Burwitz, V., Zavlin, V. E., Neuhäuser, R., et al. 2001, A&A, 379, L35

Cackett, E. M., Miller, J. M., Ballantyne, D. R., et al. 2010, ApJ, 720, 205

Cackett, E. M., Miller, J. M., Bhattacharyya, S., et al. 2008, ApJ, 674, 415

Cadeau, C., Morsink, S. M., Leahy, D., & Campbell, S. S. 2007, ApJ, 654, 458

Carney, M. F., Wade, L. E., & Irwin, B. S. 2018, PhRvD, 98, 063004

Carson, Z., Steiner, A. W., & Yagi, K. 2019, PhRvD, 99, 043010

Cavecchi, Y., Watts, A. L., Braithwaite, J., & Levin, Y. 2013, MNRAS, 434, 3526

Cerutti, B., & Beloborodov, A. M. 2017, SSRv, 207, 111

Chambers, F. R. N., Watts, A. L., Keek, L., Cavecci, Y., & Garcia, F. 2019, ApJ, 871, 61

Chang, P., & Bildsten, L. 2003, ApJ, 585, 464

Chang, P., & Bildsten, L. 2004, ApJ, 605, 830

Colgan, J., Kilcrease, D. P., Magee, N. H., et al. 2016, ApJ, 817, 116

Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1994, ApJ, 424, 823

Cromartie, H., Fonseca, E., Ransom, S. M., et al. 2019, NatAs, in press (arXiv:1904.06759)

Dagum, L., & Menon, R. 1998, CSE, 5, 46

Dalcín, L., Paz, R., Storti, M., & D'Elía, J. 2008, JPDC, 68, 655

De, S., Finstad, D., Lattimer, J. M., et al. 2018, PhRvL, 121, 091102

Degenaar, N., Miller, J. M., Chakrabarty, D., et al. 2015, MNRAS, 451, L85

Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Natur, 467, 1081

Drago, A., Lavagno, A., & Pagliara, G. 2014, PhRvD, 89, 043014

Drake, J. J., Marshall, H. L., Dreizler, S., et al. 2002, ApJ, 572, 996

Droettboom, M., Caswell, T. A., Hunter, J., et al. 2018, matplotlib/matplotlib v2.2.2, Zenodo, doi:10.5281/zenodo.1202077

Durant, M., Kargaltsev, O., Pavlov, G. G., et al. 2012, ApJ, 746, 6

Dyks, J., & Harding, A. K. 2004, ApJ, 614, 869

Fabian, A. C., Iwasawa, K., Reynolds, C. S., & Young, A. J. 2000, PASP, 112, 1145

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2013, arXiv:1306.2144

Forum, M. P. 1994, MPI: A Message-Passing Interface Standard, Tech. Rep., MPIF, https://www.mpi-forum.org/docs/mpi-1.0/mpi-10.ps

Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360

Garcia, F., Chambers, F. R. N., & Watts, A. L. 2018, PhRvF, 3, 123501

Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Gil, J., Gronkowski, P., & Rudnicki, W. 1984, A&A, 132, 312

Gonzalez-Caniulef, D., Guillot, S., & Reisenegger, A. 2019, MNRAS, 490, 5848

Gough, B. 2009, GNU Scientific Library Reference Manual (3rd ed.; Network Theory Ltd.)

Gourgouliatos, K. N., & Hollerbach, R. 2018, ApJ, 852, 21

Gralla, S. E., Lupsasca, A., & Philippov, A. 2017, ApJ, 851, 137

Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., & Watts, A. L. 2019, MNRAS, 485, 5363

Grenier, I. A., & Harding, A. K. 2015, CRPhy, 16, 641

Guillot, S. 2016, MmSAI, 87, 521

Guillot, S., Kaspi, V. M., Archibald, R. F., et al. 2016, MNRAS, 463, 2612

Guillot, S., & Rutledge, R. E. 2014, ApJL, 796, L3

Guillot, S., Rutledge, R. E., & Brown, E. F. 2011, ApJ, 732, 88

Guillot, S., Servillat, M., Webb, N. A., & Rutledge, R. E. 2013, ApJ, 772, 7

Güver, T., & Özel, F. 2013, ApJL, 765, L1

Güver, T., Psaltis, D., & Özel, F. 2012, ApJ, 747, 76

Hamaguchi, K. & NICER Calibration Team 2019, in 14th Workshop of Int. Astronomical Consortium for High Energy Calibration, https://indico2. riken.jp/event/2910/contributions/12756/attachments/8511/10339/ NICER-IACHEC-2019.pdf

Handley, W. 2018, JOSS, 3, 849

Harding, A. K., & Muslimov, A. G. 2001, ApJ, 556, 987

Harding, A. K., & Muslimov, A. G. 2011, ApJ, 743, 181

Harrison, E. R., & Tademaru, E. 1975, ApJ, 201, 447

Hartle, J. B. 1967, ApJ, 150, 1005

Hartle, J. B., & Thorne, K. S. 1968, ApJ, 153, 807

Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2010, PhRvL, 105, 161102

Heinke, C. O., Cohn, H. N., Lugger, P. M., et al. 2014, MNRAS, 444, 443

Heinke, C. O., Rybicki, G. B., Narayan, R., & Grindlay, J. E. 2006, ApJ, 644, 1090

Heyl, J. S. 2004, ApJ, 600, 939
Higson, E. 2018, JOSS, 3, 916

Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2018, BayAn, 13, 873

Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2019, MNRAS, 483, 2044

Ho, W. C. G., & Heinke, C. O. 2009, Natur, 462, 71

Ho, W. C. G., Kaplan, D. L., Chang, P., van Adelsberg, M., & Potekhin, A. Y. 2007, MNRAS, 375, 821

Ho, W. C. G., & Lai, D. 2001, MNRAS, 327, 1081

Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

Hunter, J. D. 2007, CSE, 9, 90

Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943

Johnson, T. J., Venter, C., Harding, A. K., et al. 2014, ApJS, 213, 6

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python (http://www.scipy.org/) /

Kajava, J. J. E., Nättilä, J., Latvala, O.-M., et al. 2014, MNRAS, 445, 4218

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kijak, J., & Gil, J. 2003, A&A, 397, 969

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (Amsterdam: IOS Press), 87

Kramer, A., & Wex, N. 2009, CQGra, 26, 073001

Kullback, S., & Leibler, R. A. 1951, Ann. Math. Statist., 22, 79

Kumar, B., & Landry, P. 2019, PhRvD, 99, 123026

Kurkela, A., Fraga, E. S., Schaffner-Bielich, J., & Vuorinen, A. 2014, ApJ, 789, 127

Landry, P., & Essick, R. 2019, PhRvD, 99, 084049

Lattimer, J. M., & Prakash, M. 2016, PhR, 621, 127

Lattimer, J. M., & Steiner, A. W. 2014, ApJ, 784, 123

Lee, H., Kashyap, V. L., van Dyk, D. A., et al. 2011, ApJ, 731, 126

Li, J. J., & Sedrakian, A. 2019, ApJL, 874, L22

Lim, Y., & Holt, J. W. 2018, PhRvL, 121, 062701

Lindblom, L. 2018, PhRvD, 97, 123019

Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19

Lockhart, W., Gralla, S. E., Özel, F., & Psaltis, D. 2019, MNRAS, 490, 1774

Lommen, A. N., Kipphorn, R. A., Nice, D. J., et al. 2006, ApJ, 642, 1012

Lommen, A. N., Zepka, A., Backer, D. C., et al. 2000, ApJ, 545, 1007

Long, M., Romanova, M. M., & Lovelace, R. V. E. 2007, MNRAS, 374, 436

Long, M., Romanova, M. M., & Lovelace, R. V. E. 2008, MNRAS, 386, 1274

Ludlam, R. M., Miller, J. M., Arzoumanian, Z., et al. 2018, ApJL, 858, L5

Ludlam, R. M., Miller, J. M., Bachetti, M., et al. 2017, ApJ, 836, 140

Malik, T., Alam, N., Fortin, M., et al. 2018, PhRvC, 98, 035804

Margueron, J., Hoffmann Casali, R., & Gulminelli, F. 2018a, PhRvC, 97, 025805

Margueron, J., Hoffmann Casali, R., & Gulminelli, F. 2018b, PhRvC, 97, 025806

McNeil Forbes, M., Bose, S., Reddy, S., et al. 2019, PhRvD, 100, 083010

Melatos, A., & Phinney, E. S. 2001, PASA, 18, 421

Miller, J. M., Parker, M. L., Fuerst, F., et al. 2013, ApJL, 779, L2

Miller, M. C. 2013, arXiv:1312.0029

Miller, M. C., & Lamb, F. K. 1998, ApJL, 499, L37

Miller, M. C., & Lamb, F. K. 2015, ApJ, 808, 31

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24

Mitchell, J. P., Braithwaite, J., Reisenegger, A., et al. 2015, MNRAS, 447, 1213

Montaña, G., Tolós, L., Hanauske, M., & Rezzolla, L. 2019, PhRvD, 99, 103009

Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244

Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018, PhRvL, 120, 261103

Nandra, K., Barret, D., Barcons, X., et al. 2013, arXiv:1306.2307

Nättilä, J., Miller, M. C., Steiner, A. W., et al. 2017, A&A, 608, A31

Nättilä, J., & Pihajoki, P. 2018, A&A, 615, A50

Nättilä, J., Steiner, A. W., Kajava, J. J. E., Suleimanov, V. F., & Poutanen, J. 2016, A&A, 591, A25

Obergaulinger, M., & Aloy, M. Á. 2017, JPhCS, 932, 012043

Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, RvMP, 89, 015007

Oliphant, T. E. 2007, CSE, 9, 10

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Özel, F., & Psaltis, D. 2009, PhRvD, 80, 103003

Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28

Patruno, A., & Watts, A. L. 2012, arXiv:1206.2727

Pavlov, G. G., & Zavlin, V. E. 1997, ApJL, 490, L91

Payne, D. J. B., & Melatos, A. 2004, MNRAS, 351, 569

Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJL, 815, L19

Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94

Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015b, ApJL, 801, L19

Piro, A. L., & Bildsten, L. 2005, ApJ, 629, 438

Pons, J. A., Walter, F. M., Lattimer, J. M., et al. 2002, ApJ, 564, 981

Poutanen, J., & Beloborodov, A. M. 2006, MNRAS, 373, 836

Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301

Psaltis, D., & Özel, F. 2014, ApJ, 792, 87

Psaltis, D., Özel, F., & Chakrabarty, D. 2014, ApJ, 787, 136

Raaijmakers, G., Riley, T. E., & Watts, A. L. 2018, MNRAS, 478, 2177

Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, ApJL, 887, L22

Radhakrishnan, V., & Cooke, D. J. 1969, ApL, 3, 225

Radice, D., Perego, A., Hotokezaka, K., et al. 2018, ApJL, 869, L35

Raithel, C. A., Özel, F., & Psaltis, D. 2016, ApJ, 831, 44

Raithel, C. A., Özel, F., & Psaltis, D. 2017, ApJ, 844, 156

Rangelov, B., Pavlov, G. G., Kargaltsev, O., et al. 2017, ApJ, 835, 264

Ray, P. S., Arzoumanian, Z., Ballantyne, D., et al. 2019, arXiv:1903.03035

Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, PhRvD, 79, 124032

Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751

Reisenegger, A. 2009, A&A, 499, 557

Riley, T. E. 2019, A NICER View of PSR J0030+0451: Nested Samples for Millisecond Pulsar Parameter Estimation, v1.0, Zenodo, doi:10.5281/zenodo.3386449

Riley, T. E., Raaijmakers, G., & Watts, A. L. 2018, MNRAS, 478, 1093

Riley, T. E., & Watts, A. L. 2019, ApJS, submitted

Romani, R. W. 1990, Natur, 347, 741

Salmi, T., Nättilä, J., & Poutanen, J. 2018, A&A, 618, A161

Shaw, A. W., Heinke, C. O., Steiner, A. W., et al. 2018, MNRAS, 476, 4713

Shlens, J. 2014, CoRR, arXiv:1404.2000

Skilling, J. 2006, BayAn, 1, 833

Spitkovsky, A., Levin, Y., & Ushomirsky, G. 2002, ApJ, 566, 1018

Steiner, A. W., Heinke, C. O., Bogdanov, S., et al. 2018, MNRAS, 476, 421

Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33

Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, ApJL, 765, L5

Stevens, A. L., Fiege, J. D., Leahy, D. A., & Morsink, S. M. 2016, ApJ, 833, 244

Strobel, K., Schaab, C., & Weigel, M. K. 1999, A&A, 350, 497

Suleimanov, V., Poutanen, J., Revnivtsev, M., & Werner, K. 2011a, ApJ, 742, 122

Suleimanov, V., Poutanen, J., & Werner, K. 2011b, A&A, 527, A139

Tews, I., Carlson, J., Gandolfi, S., & Reddy, S. 2018a, ApJ, 860, 149

Tews, I., Lattimer, J. M., Ohnishi, A., & Kolomeitsev, E. E. 2017, ApJ, 848, 107

Tews, I., Margueron, J., & Reddy, S. 2018b, PhRvC, 98, 045804

The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, B. P., et al. 2019, PhRvX, 9, 031040

Timokhin, A. N., & Arons, J. 2013, MNRAS, 429, 20

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van Doesburgh, M., van der Klis, M., & Morsink, S. M. 2018, MNRAS, 479, 426

Vehtari, A., & Ojanen, J. 2012, Statist. Surv., 6, 142

Viganò, D., Rea, N., Pons, J. A., et al. 2013, MNRAS, 434, 123

Vincent, F. H., Bejger, M., Różańska, A., et al. 2018, ApJ, 855, 116

Watts, A. L. 2012, ARA&A, 50, 609

Watts, A. L. 2019, in AIP Conf. Ser. 2127, Xiamen-CUSTIPEN Workshop on the EOS of Dense Neutron-Rich Matter in the Era of Gravitational Wave Astronomy, ed. A. Li, B.-A. Li, & F. Xu (Melville, NY: AIP) 020008

Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, RvMP, 88, 021001

Watts, A. L., Yu, W., Poutanen, J., et al. 2019, SCPMA, 62, 29503

Webb, N. A., & Barret, D. 2007, ApJ, 671, 727

Weih, L. R., Most, E. R., & Rezzolla, L. 2019, ApJ, 881, 73

Wilkinson, T., Patruno, A., Watts, A., &uttley, P. 2011, MNRAS, 410, 1513

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Xu, J., van Dyk, D. A., Kashyap, V. L., et al. 2014, ApJ, 794, 97

Zavlin, V. E., Pavlov, G. G., & Shibanov, Y. A. 1996, A&A, 315, 141

Zdunik, J. L., & Haensel, P. 2013, A&A, 551, A61

Zhang, S., Santangelo, A., Feroci, M., et al. 2019, SCPMA, 62, 29502