# A NICER View of PSR J1231-1411: A Complex Case

Tuomo Salmi¹,² ID, Julia S. Deneva³ ID, Paul S. Ray⁴ ID, Anna L. Watts¹ ID, Devarshi Choudhury¹ ID, Yves Kini¹ ID, Serena Vinciguerra¹ ID, H. Thankful Cromartie⁵,<fnref n="14" /> ID, Michael T. Wolff⁴ ID, Zaven Arzoumanian⁶ ID, Slavko Bogdanov⁷ ID, Keith Gendreau⁶ ID, Sebastien Guillot⁸ ID, Wynn C. G. Ho⁹ ID, Sharon M. Morsink¹⁰ ID, Ismaël Cognard¹¹,¹² ID, Lucas Guillemot¹¹,¹² ID, Gilles Theureau¹¹,¹²,¹³ ID, and Matthew Kerr⁴ ID

¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands

² Department of Physics, University of Helsinki, P.O. Box 64, FI-00014, Finland; tuomo.salmi@helsinki.fi

³ George Mason University, Resident at the U.S. Naval Research Laboratory, Washington, DC 20375, USA

⁴ Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375, USA

⁵ National Research Council Research Associate, National Academy of Sciences, Washington, DC 20001, USA

⁶ X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA

⁷ Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA

⁸ Institut de Recherche en Astrophysique et Planétologie, UPS-OMP, CNRS, CNES, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France

⁹ Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA

¹⁰ Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB, T6G 2E1, Canada

¹¹ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, 45071 Orléans Cedex 02, France

¹² Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, 18330 Nançay, France

¹³ Laboratoire Univers et Théories LUTH, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

Received 2024 August 26; revised 2024 September 20; accepted 2024 September 21; published 2024 November 13

## Abstract

Recent constraints on neutron star mass and radius have advanced our understanding of the equation of state (EOS) of cold dense matter. Some of them have been obtained by modeling the pulses of three millisecond X-ray pulsars observed by the Neutron Star Interior Composition Explorer (NICER). Here, we present a Bayesian parameter inference for a fourth pulsar, PSR J1231-1411, using the same technique with NICER and XMM-Newton data. When applying a broad mass-inclination prior from radio timing measurements and the emission region geometry model that can best explain the data, we find likely converged results only when using a limited radius prior. If limiting the radius to be consistent with the previous observational constraints and EOS analyses, we infer the radius to be $12.6 \pm 0.3$ km and the mass to be $1.04_{-0.03}^{+0.05} M_{\odot}$, each reported as the posterior credible interval bounded by the 16% and 84% quantiles. If using an uninformative prior but limited between 10 and 14 km, we find otherwise similar results, but $R_{eq} = 13.5_{-0.5}^{+0.3}$ km for the radius. In both cases, we find a nonantipodal hot region geometry where one emitting spot is at the equator or slightly above, surrounded by a large colder region, and where a noncircular hot region lies close to southern rotational pole. If using a wider radius prior, we only find solutions that fit the data significantly worse. We discuss the challenges in finding the better fitting solutions, possibly related to the weak interpulse feature in the pulse profile.

Unified Astronomy Thesaurus concepts: X-ray astronomy (1810); Neutron stars (1108); Millisecond pulsars (1062)

Materials only available in the online version of record: figure sets

## 1. Introduction

Modeling the X-ray pulses from a neutron star (NS) can be used to infer the NS's mass, radius, and the geometric properties of the X-ray emitting regions (see e.g., G. G. Pavlov & V. E. Zavlin 1997; A. L. Watts et al. 2016; S. Bogdanov et al. 2019b, and the references therein). The inferred mass and radius can then be used to constrain the equation of state (EOS) of the supranuclear dense matter in NS cores (see e.g., J. M. Lattimer & M. Prakash 2007). On the other hand, the emitting regions are associated with the NS's polar caps, and their properties can be used to probe the NS's magnetic field structure (e.g., A. V. Bilous et al. 2019; F. Carrasco et al. 2023).

Previously, X-ray pulses of three NSs have been modeled using Neutron Star Interior Composition Explorer (NICER; K. C. Gendreau et al. 2016) observations: PSR J0030+0451

(M. C. Miller et al. 2019; T. E. Riley et al. 2019; T. Salmi et al. 2023; S. Vinciguerra et al. 2024), the massive pulsar PSR J0740+6620 (M. C. Miller et al. 2021; T. E. Riley et al. 2021; T. Salmi et al. 2022, 2023; A. J. Dittmann et al. 2024; T. Salmi et al. 2024), and the bright pulsar PSR J0437-4715 (D. Choudhury et al. 2024a). These have already provided useful constraints on dense matter models (see, e.g., M. C. Miller et al. 2021; G. Raaijmakers et al. 2021a; B. Biswas 2022; E. Annala et al. 2023; M.-Z. Han et al. 2023; J. Takátsy et al. 2023; N. Rutherford et al. 2024). However, measuring masses and radii for a larger set of NSs and over a broad range of masses is expected to allow both tighter and more robust constraints.

PSR J1231-1411 is a rotation-powered millisecond pulsar (RMP; discovered by S. M. Ransom et al. 2011), like many NSs previously analyzed with NICER. Its thermally emitting surface regions are expected to be heated by the bombardment of charged particles from a magnetospheric return current (see, e.g., M. A. Ruderman & P. G. Sutherland 1975; J. Arons 1981; A. K. Harding & A. G. Muslimov 2001). Soft X-ray pulsations at the known spin frequency of 271 Hz were found by P. S. Ray et al. (2019) and further analyzed in S. Bogdanov et al. (2019a).

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

<footnote n="14">¹⁴ Resident at Naval Research Laboratory, Washington, DC 20375, USA.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The phase-folded PSR J1231−1411 event data for two rotational cycles (for clarity). The top panel shows the pulse profile summed over the channels. The total number of counts is given by the sum over all phase-channel pairs (over both cycles). For the modeling, all the event data are grouped into a single rotational cycle instead, and thus each phase-channel bin has twice the number of counts shown here.</figcaption></figure>

The pulsed emission shows a broad and slightly asymmetric main pulse and a much fainter secondary interpulse (see Figure 1). This feature differs from those seen in NICER RMPs analyzed so far (see S. Bogdanov et al. 2019a for PSR J0030+0451 and PSR J0437−4715, and M. T. Wolff et al. 2021 for PSR J0740+6620). Similar to PSR J0740+6620 and PSR J0437−4715, PSR J1231−1411 is located in a binary system, allowing radio timing measurements to constrain the NS mass and observer inclination based on the relativistic Shapiro delay (E. Fonseca et al. 2021; D. J. Reardon et al. 2024). However, the obtained constraints in the preliminary radio timing results (available at the time of this work) are far less restrictive for PSR J1231−1411 (T. Cromartie et al. 2024, in preparation) compared to PSR J0740+6620 and PSR J0437−4715. In this work, we use these measurements as prior distributions and explore the robustness of NS parameter inference for PSR J1231−1411.

The remainder of this paper is structured as follows. In Section 2, we introduce the NICER and XMM-Newton data sets used for PSR J1231−1411. In Section 3, we summarize the modeling procedure, and in Section 4 we present the results. We discuss the implications of the results in Section 5 and conclude in Section 6.

## 2. X-Ray Event Data

### 2.1. NICER

The NICER data set for this work was produced using HEASoft 6.32.1¹⁵ (NASA High Energy Astrophysics Science Archive Research Center (Heasarc), 2014) and CALDB xti20221001. We selected data from the full NICER mission up through 2023 May 14 (ObsID 6060060730), which is before NICER was affected by the light leak.<fnref n="16" /> After standard L2 processing with the task nicer12, but before any of our filtering, the initial data set contained about 3.53 Ms of

exposure time. We applied further cuts to exclude data with magnetic cutoff rigidity <1.5 GeV c⁻¹; planetary K-index Kₚ > 5; overshoot rate<fnref n="17" /> > 1.5 counts per second (cps) per Focal Plane Module (FPM); or a median undershoot rate<fnref n="18" /> >200 cps/FPM. These cuts reduce the exposure to 81% of the L2 exposure. Subsequently, we generated a background light curve using only 2–10 keV photons (where the pulsar flux is small so the count rate is background dominated) with 16 s bins and discarded times where that rate exceeded 1.25 cps. This removes the times polluted by high levels of trapped electron (TREL, in the nomenclature of the SCORPEON background model, see below) and low-energy electron (LEEL) backgrounds, reducing the exposure to 75% of the L2 exposure. The final exposure was 2651393.0413 s.

The timing model to assign pulse phases to photon measured times was developed by first fitting a timing model to the radio time of arrivals from the Nançay radio telescope (NRT). Details on the NRT and on pulsar observations with this instrument can be found in L. Guillemot et al. (2023). This determined the astrometric and binary parameters (including Shapiro delay) and a constant dispersion measure (DM). These were then held fixed while the frequency and frequency first derivative were fit to the Fermi-LAT gamma-ray data, as done in FERMI-LAT Collaboration et al. (2022). Finally, the reference epoch was set to the values used in the Third Fermi-LAT Pulsar Catalog (3PC; D. A. Smith et al. 2023) so the absolute pulse phases will be consistent with that. This model provides pulse phases accurate to <0.1% of a pulse period over the entire duration of the NICER data set.

The NICER science tool nicer13-spect was used to compute a pulse height spectrum, response matrix, and background estimate using the SCORPEON background model.<fnref n="19" /> We used XSPEC (K. A. Arnaud 1996) to fit the SCORPEON model parameters simultaneously with different source spectral models, including two blackbody components, one blackbody and one power-law component (with a power-law index fixed at 1.8 as in S. M. Ransom et al. 2011), just a power-law component, one NS atmosphere component (nsatmos), and two NS atmosphere components. These were all modified by interstellar absorption (e.g., TBabs (BB + powerlaw)). For these data we found that, in addition to the standard free parameters in the SCORPEON model, we needed to free the normalization of the O VII line at 0.57 keV (presumably originating from solar wind charge exchange). In prior work before the development of SCORPEON (P. S. Ray et al. 2019), this line was added as an ad hoc Gaussian line. We found that the obtained SCORPEON background estimates were notably dependent on the chosen source model, and therefore we predicted the background to range from the smallest estimate to the highest at each channel (after including also one standard deviation uncertainty for each model). This range was used only in the end to compare our inferred background to it (see Section 4.2).

In the pulse profile analysis, we used the pulse invariant (PI) channel subset [30, 150], corresponding to the nominal photon

<footnote n="15">¹⁵ https://heasarc.gsfc.nasa.gov/ftools</footnote>

<footnote n="16">¹⁶ https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/light-leak-overview/</footnote>

<footnote n="17">¹⁷ https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/overshoot-intro/</footnote>

<footnote n="18">¹⁸ https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/undershoot-intro/</footnote>

<footnote n="19">¹⁹ https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/scorpeon-overview/</footnote>
energy range [0.3, 1.5] keV, and 32 rotational phase bins. The data split over two rotational cycles are visualized in Figure 1.

## 2.2. XMM-Newton

The XMM-Newton EPIC-MOS (M. J. L. Turner et al. 2001; 1 and 2) and EPIC-pn (L. Strüder et al. 2001) data (ObsID 0605470201) used in this study were obtained on 2009 July 15, as part of a Fermi/LAT unidentified source observation campaign (S. M. Ransom et al. 2011). We utilize the Scientific Analysis System software (SAS)<fnref n="20" /> to analyze and build our data products, which in turn are input into the inference software. The fields of the new unidentified LAT sources were observed using the XMM-Newton imaging mode, thus with insufficient time resolution to look for pulsations. Using the original raw data products from the EPIC-MOS and EPIC-pn observations, we extracted X-ray events from the now-known pulsar position. We set the PATTERN ≤ 12 for the EPIC-MOS1 and EPIC-MOS2 instruments and PATTERN ≤ 4 for the EPIC-pn along with FLAG = 0 for all three instruments. This resulted in 28.8 ks of exposure from the EPIC-MOS instruments and 17.7 ks of exposure for the EPIC-pn instrument. For the EPIC-MOS1 and EPIC-MOS2 event data, we chose source+background extraction regions of 32" and source-free background extraction regions of 90". The extraction regions for source+background and background events for the EPIC-pn were 38" and 60", respectively. We filtered out times of higher-than-average particle noise from each event stream. The resulting final list of events are stored in FITS files that are then converted into ASCII files containing just the event times and PI value. However, no pulse phases are included in the XMM-Newton event files because the observations were made with a detector frame time that was too long to detect pulsations. Once the event filtering is selected, we build X-ray spectra using the SAS tool “evselect” and then use the SAS tools “rmfgen” and “arfgen” to build response matrices. The resulting response matrices are then turned into ASCII text files for ingesting into the inference code that is described below. We also implicitly make the assumption that the average pulsar X-ray properties have not changed from the epoch of the XMM-Newton observation to the epoch of the NICER observations. The assumption is justified since no evidence for long-term flux variability of PSR J1231-1411 was found in S. Bogdanov et al. (2019a, see their Section 7.3).

In the analysis, we used the channel subset [50, 200) for the EPIC-pn instrument, and the channel subset [20, 100) for the both EPIC-MOS1 and EPIC-MOS2 instruments. These correspond to the nominal photon energy range [0.27–1.02 keV] for EPIC-pn and [0.3–1.5 keV] for EPIC-MOS1 and EPIC-MOS2. The data are visualized in Figure 2.

# 3. Modeling Procedure

We use the X-ray Pulse Simulation and Inference<fnref n="21" /> (X-PSI) code, with version v2.1.1 for the inference runs and v2.2.7 for producing the figures (T. E. Riley et al. 2023). Complete information of each run, data products, posterior sample files, and all of the analysis files can be found in the Zenodo repository: doi:10.5281/zenodo.13358349 (T. Salmi et al. 2024). In the next sections, we summarize the modeling

procedure and focus on how it differs from that used in previous work.

## 3.1. Effective-area Scaling Models

As in previous works (T. E. Riley et al. 2021; D. Choudhury et al. 2024a; T. Salmi et al. 2024; S. Vinciguerra et al. 2024), we capture the uncertainty in the effective area of both NICER and all XMM-Newton detectors by defining absolute energy-independent effective-area scaling factors. Each factor consists of a shared and telescope specific component (see Section 2.2. in T. Salmi et al. 2024), and, as in most recent works, we select the uncertainties in those to be 10% and 3%, correspondingly, leading to a 10.4% uncertainty in the full scaling factor. As before, we commonly assume all scaling factors to be identical between EPIC-pn, EPIC-MOS1, and EPIC-MOS2 (leaving then only the scaling factors $\alpha_{NICER}$ and $\alpha_{XMM}$ as free parameters). Since this may not be true, we have also performed a test run treating $\alpha_{XMM}$ with $\alpha_{pn}$, $\alpha_{MOS1}$, and $\alpha_{MOS2}$ as free independent parameters, although with all of them defined to have the same uncertainty (10.4% as mentioned above) and the same correlation with $\alpha_{NICER}$ (see Section 4.1).

## 3.2. Pulse Profile Modeling Using X-PSI

As in previous NICER analyses, we use the “Oblate Schwarzschild” approximation for the spacetime and NS shape model to calculate the energy-resolved X-ray pulses (see e.g., M. C. Miller & F. K. Lamb 1998; N. R. Nath et al. 2002; J. Poutanen & M. Gierlinski 2003; C. Cadeau et al. 2007; S. M. Morsink et al. 2007; K. H. Lo et al. 2013; M. AlGendy & S. M. Morsink 2014; S. Bogdanov et al. 2019b; A. L. Watts 2019). The spectrum and the beaming pattern of the surface radiation are in most cases assumed to follow those from a fully ionized hydrogen atmosphere model NSX (W. C. G. Ho & D. Lai 2001), but we also performed a run assuming fully ionized helium NSX atmosphere (see Section 4.2). In addition, for one case we also modeled the NS surface outside the hot regions, assuming it emits blackbody radiation (for computational efficiency). This assumption is more realistic than having no surface emission, although more accurate still would be to use a numerical (and thus more expensive) partially ionized atmosphere model for the likely relatively cold surface. To account for the interstellar attenuation, we use the TBabs model (J. Wilms et al. 2000, updated in 2016).

To characterize the geometry (and number) of the hot regions, we apply two different models: Single Temperature-Unshared (ST-U) and Protruding Dual Temperature-Unshared (PDT-U). In the former model two hot regions are described by two single-temperature disjoint spherical caps with unshared parameters. The latter model is the most complex model applied with X-PSI so far, where each of the two hot regions is described by two overlapping spherical caps with different temperatures (see T. E. Riley et al. 2019 and Figure 1 of S. Vinciguerra et al. 2023 for more details). Where the two overlap, emission of only one of them is accounted for, and that component is referred to as being “superseding” (the other being “ceding”).

For the pulse profile modeling resolution with X-PSI, we mostly follow the settings<fnref n="22" /> used in the previous PSR J0740+6620 analyses (T. E. Riley et al. 2021; T. Salmi et al. 2022, 2024).

<footnote n="20">²⁰ https://www.cosmos.esa.int/web/xmm-newton/sas</footnote>

<footnote n="21">²¹ https://github.com/xpsi-group/xpsi</footnote>

<footnote n="22">²² For more details, see Section 2.3.1 in S. Vinciguerra et al. (2024) or X-PSI documentation at https://xpsi-group.github.io/xpsi/hotregion.html.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** XMM-Newton EPIC-pn camera PSR J1231-1411 count number spectrum as a function of nominal photon energies (blue-stepped curve). The measured rescaled background (red-stepped curve) and the background limits used in the analysis (purple dashed and dotted-dashed stepped curves) are also shown. The choice of the limits is explained in Section 3.4. The complete figure set (3 images), including also the spectra and backgrounds for EPIC-MOS1 and EPIC-MOS2 instruments is available in the online journal (HTML version).</figcaption></figure>

(The complete figure set (3 images) is available in the online article.)

For example, we calculate multiple images (photons reaching the observer from the same point at the NS surface with many deflection angles due to high gravity) up to third order in case they are visible. The initial deflection angle grid is again calculated for 512 different emission angles at each colatitude. However, to speed-up the calculation we reduced the number of model energies (from 128 to 64) and phases (from 64 to 32) for the more expensive PDT-U model. This change was tested and found to have no effects on the results for a ST-U run. In contrast, the resolution of the hot region cell mesh (sqrt_num32, min10, maxsqrt_num�64) was not reduced to avoid loss of accuracy in the computed waveforms for extreme configurations, as found by D. Choudhury et al. (2024b).

3.3. Priors

3.3.1. Mass and Inclination Priors

We use mass *M* and inclination *i* priors that are based on the preliminary Shapiro delay measurement results from T. Cromartie et al. (2024, in preparation), which were obtained through a Bayesian timing model fit to data from the Robert C. Byrd Green Bank Telescope and from the NRT. To approximate the 2D probability density function, we first fitted a Gaussian function to the marginalized probability distribution of *M* using a subset of the Markov chain Monte Carlo samples from the Shapiro delay analysis (those having *M* between 1 and 3 *M*<sub>☉</sub>). We obtained a mean of ~1.002 *M*<sub>☉</sub> and a standard deviation of ~0.930 *M*<sub>☉</sub>. We inverse sampled the mass from this distribution (truncating it between 1 and 3 *M*<sub>☉</sub>).<sup>23</sup> After that, we inverse sampled the cos *i*

<figure><img src="image_4.png" /><figcaption>**Figure 3.** A 2D cos *i*–*M* prior distribution based on the radio timing measurements from T. Cromartie et al. (2024, in preparation). The color indicates the probability of a cos *i*–*M* cell.</figcaption></figure>

from another Gaussian, whose mean and standard deviations
were set to depend on the value of the already drawn mass;
μ<sub>cos i</sub>(M), σ<sub>cos i</sub>(M). The values were found by fitting Gaussians
to the cos i probability distributions for a grid of masses and by
fitting their dependence on mass with a polynomial function.
This resulted in the following:

$$
\mu_{\cos i}(M) = a_{\mu} \left( \frac{M}{M_{\odot}} \right)^2 + b_{\mu} \frac{M}{M_{\odot}} + c_{\mu} \quad (1)
$$

and

$$
\sigma_{\cos i}(M) = a_{\sigma} \left( \frac{M}{M_{\odot}} \right)^2 + b_{\sigma} \frac{M}{M_{\odot}} + c_{\sigma}, \quad (2)
$$

where a_μ = -0.00835942, b_μ = 0.10890304, c_μ = 0.01118702,
a_σ = -0.00097777, b_σ = 0.01013241, and c_σ = 0.01509739.

The resulting 2D prior for *M* and cos *i* is shown in Figure 3. We
see that highest probabilities are at the smallest allowed masses but
that the distribution fades very slowly toward higher masses.
Inclination is slightly better constrained, peaking above 80°
(cos *i* ≲ 0.18). We note that the final prior adopted in X-PSI is
additionally modified by the compactness and surface gravity
conditions (as e.g., in T. E. Riley et al. 2021) to prevent the polar
radius from being smaller than the photon sphere (R<sub>polar</sub>/r<sub>g</sub>(*M*) > 3,
where r<sub>g</sub> is the Schwarzschild radius), and to not extrapolate
intensities beyond the precomputed atmosphere grid (e.g.,
13.7 ≲ log<sub>10</sub> g ≲ 15.0, for the fully ionized hydrogen grid, where
g is the surface gravity). This makes the prior probability at smaller
masses even higher because fewer NS samples fail the compact-
ness condition there, hence more are retained.

3.3.2. Radius Prior

For the NS (equatorial) radius, X-PSI analyses have usually
assumed a prior that is initially flat between ~4.4 km (i.e., 3r<sub>g</sub> for a
1 M<sub>☉</sub> star) and 16 km but then modified by the compactness and
surface gravity conditions, mentioned also in Section 3.3.1. In this
work, we apply also more informative radius priors. In particular,
we consider four cases when choosing a radius prior. In the first
scenario, we use the standard X-PSI radius prior. In the second
scenario, we set the lower limit of radius to be 8 km to avoid
extremely compact solutions. In the third scenario, we restrict the
radius to be between 10 and 14 km. In the fourth scenario (called

<footnote n="23"><sup>23</sup> We note that the resulting prior is still quite close to uniform.</footnote>
"R21"), we require the radius to be consistent with previous observational and theoretical EOS constraints by using the posterior samples from G. Raaijmakers et al. (2021a; on Zenodo, G. Raaijmakers et al. 2021b) to numerically build up a 1D prior distribution for radius. We selected the samples from the "combined constraints" case obtained with the "PP" parameterization of the EOS corresponding roughly a Gaussian centered around 12 km (see the lower left-hand panel in Figure 5 of G. Raaijmakers et al. 2021a). We converted the samples to a cumulative probability distribution using the ` numpy.cumsum` function and interpolated from that distribution with `scipy.interpolate.Akima1DInterpolator` when sampling.

### 3.3.3. Other Priors

For all of the other model parameters, we use fairly uninformative priors. For the geometry parameters describing the hot regions, we follow the choices of T. E. Riley et al. (2021) and S. Vinciguerra et al. (2024), i.e., we have initially uniform priors in the spherical cap angular radii (limited to be below $\pi/2$), in their phase coordinates and in cosines of their colatitudes. These priors are then modified by the requirements that primary and secondary hot regions cannot overlap (but in case of `PDT-U` the two caps forming a hot region must overlap), and that the primary hot region is always the one with smaller colatitude. The hot region temperatures are uniform in $log_{10} T_{\text{eff}}$ and bound based on the limits of the precomputed atmosphere model grid (e.g., $log_{10} T_{\text{eff}} \in [5.1, 6.8]$ for fully ionized hydrogen). When including blackbody emission from the rest of the star's surface, we sample its temperature $log_{10} T_{\text{else}}$ uniformly between 5.0 and 6.5, as in S. Vinciguerra et al. (2024).

Unlike the previous NSs analyzed with NICER, the distance D to PSR J1231-1411 is known rather poorly. J. M. Yao et al. (2017) estimated the D to be 420 pc using the DM of PSR J1231-1411 and an electron density model of the Galaxy. However, they also reported a distance range between 350 and 510 pc based on independent estimates from optical measurements of the white dwarf companion. In this work, we use slightly more conservative limits and sample D uniformly between 100 and 700 pc. We note, though, that preliminary parallax distance measurements indicate a distance of roughly 600 ± 100 pc (T. Cromartie et al. 2024, in preparation).

Since the interstellar hydrogen column density $N_H$ is also not well known, we sample that uniformly between 0 and $10^{21} \text{cm}^{-2}$, as in T. E. Riley et al. (2021). In addition, the effective-area scaling factor priors are those described already in Section 3.1.

### 3.4. Likelihood and Background

We construct a joint NICER and XMM-Newton likelihood function and marginalize the likelihoods over background variables (one for each detector channel) in the same way as described in Section 2.5 of T. E. Riley et al. (2021). We do not directly impose any informative limits on the NICER background, but for XMM-Newton instruments we bound the background marginalization to be done within $B_X \pm 4\sqrt{B_X}$ (but always above 0), where $B_X$ is the measured number of background counts for a given channel (see the dashed lines in Figure 2 after being rescaled based on the BACKSCAL factor, i.e., the difference in the extraction area between the PSR J1231-1411 and background observation). The prior in the background is assumed to be uniform between the lower and upper limits. Constraining the background using XMM-Newton also offers an indirect constraint on the NICER

background by restricting the number of source photons that can satisfy both NICER and XMM-Newton data simultaneously. In the end, we cross-check the inferred NICER background with that predicted by the SCORPEON background model for NICER (see Sections 2.1 and 4.2).

### 3.5. Posterior Computation

We compute the posterior samples using PYMULTINEST (J. Buchner et al. 2014) and MULTINEST (F. Feroz & M. P. Hobson 2008; F. Feroz et al. 2009, 2019), as in the previous X-PSI analyses. In all runs, we set both the sampling efficiency (defined as described in Appendix B of T. Salmi et al. 2024) and the evidence tolerance to 0.1. We also keep the multimode/mode-separation setting of MULTINEST turned on to get more information about the multiple modes found in the posterior volume of PSR J1231-1411. We checked for one `ST-U` test run that turning the multimode setting off had no significant effects on the results or on the computation time. For `ST-U` runs, we used 4000 live points, which was deemed sufficient based on preliminary test runs with older data and mass-inclination prior versions. For `PDT-U` runs, we increased the number of live points to 10,000 and to 20,000. However, these were found to provide likely converged results only in the case of the two most restricted radius priors (the third and fourth cases described in Section 3.3.2). Due to the already high computational cost (above 200,000 core hours for a 20,000 live point run), better convergence by increasing the live points further or by dropping the sampling efficiency was not attempted.

## 4. Inferences

### 4.1. ST-U

We start with the inference for PSR J1231-1411 using the `ST-U` model. We applied the model using three different versions of it: (1) Basic `ST-U`, (2) `ST-U` where the rest of the star's surface is also emitting (as mentioned in Section 3.2), (3) `ST-U` where all three XMM-Newton instruments have their effective-area scaling factors as free parameters (as mentioned in Section 3.1). The corresponding posteriors for the spacetime parameters are shown in the Appendix for the standard X-PSI radius prior (posteriors for the other parameters are available in the Zenodo, T. Salmi et al. 2024 (doi:10.5281/zenodo.13358349)). We found that in all three cases the inferred radius is very small, around 5 km, and the NS is extremely compact, even though the mass is only slightly above $1 M_\odot$. In fact, all these solutions strongly hit our prior upper bound on surface gravity set by the hydrogen atmosphere table limits ($log_{10} g = 15.0$).

While no very obvious patterns were observed in the residuals between the inferred phase-energy resolved signals and the PSR J1231-1411 data (which would indicate a poor fit to the data), we found that the fit to the phase-averaged (or bolometric) NICER pulse profile is rather bad for all the `ST-U` models (see the left-hand panel of Figure 4 for an example). The bolometric $\chi^2$ values for maximum-likelihood solutions are around 80 for 32 phase bins.<fnref n="24" /> In addition, if forcing the radius to be above 8 km

<footnote n="24">²⁴ The degree of freedom (d.o.f.) is not obvious, but it should be a few smaller than 32 depending on how many of the model parameters influence the bolometric pulse profile (A. J. Dittmann et al. 2024). We also note that better bolometric $\chi^2$ values are found when looking at the 1000 highest likelihood samples in each `ST-U` run because the best bolometric NICER $\chi^2$ value does not exactly coincide with the maximum joint (2D) NICER and XMM-Newton likelihood sample. However, these $\chi^2$ values are still always above 70.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Bolometric pulse profiles (counts summed over energy channels) for the maximum-likelihood parameters (red-stepped curves) of the ST−U H T<sub>else</sub> (upper left-hand panel) and PDT−U H R<sub>eq</sub> ∈ [10, 14] km (upper right-hand panel) km models. See Table 1 and Section 4.2 for model explanations. The bolometric NICER data are shown with black-stepped curves. The differences (data − model) are shown in the bottom panels also with black-stepped curves. The expected Poisson fluctuation levels (for NICER data) are shown there with blue-dashed curves. We note that the model curves are rather similar between all the ST−U versions and between all the PDT−U versions.</figcaption></figure>

(making a full run with the second radius prior case mentioned
in Section 3.3.2), both the maximum-likelihood and the
corresponding bolometric χ² values became significantly
worse. The former dropped by more than 20 in ln-units and
the latter increased to ~160 (results of this run are available in
the Zenodo, T. Salmi et al. 2024). Most notably, the best-fit
model was no longer able to produce the interpulse around
phase 0.5 that is seen in the data (see Figure 1). This was easier
for the high compactness star, due to multiple images of the hot
spots contributing to the interpulse. The highest Bayesian
evidence for a ST−U model (the one including rest of the star
emission) was also found to be significantly smaller than for
any of the PDT−U models (by 8–33 in ln-units), so we therefore
focus on the latter next.

4.2. PDT-U

We present next the inference results for PSR J1231−1411 using the PDT−U model. In Figure 5, we show the radius, mass, compactness, and cos i posteriors for three radius priors: (1) using a lower limit of 8 km (with the upper limit being 16 km as usual); (2) limiting the radius between 10 and 14 km, and (3) using an informative radius prior based on the results of G. Raaijmakers et al. (2021a; the R21 case from Section 3.3.2). The last was performed using both hydrogen (H) and helium (He) atmospheres (He R21 with 10,000 live points instead of the 20,000 used otherwise, but the H R21 results were checked and found to be close to identical with both settings). In all cases, we report the 68.3% credible intervals around the median. For the $R_{\text{eq}} \in [8, 16]$ km case, we obtain $R_{\text{eq}} = 14.74_{-0.69}^{+0.62}$ km and $M = 1.07_{-0.05}^{+0.06} M_{\odot}$. For the $R_{\text{eq}} \in [10, 14]$ km case, we obtain

$R_{\text{eq}} = 13.52_{-0.50}^{+0.32}$ km and $M = 1.06_{-0.04}^{+0.06} M_{\odot}$. For the H R21 case, we obtain $R_{\text{eq}} = 12.60_{-0.33}^{+0.31}$ km and $M = 1.04_{-0.03}^{+0.05} M_{\odot}$. Finally, for the He R21 case, we obtain $R_{\text{eq}} = 12.48_{-0.34}^{+0.35}$ km and $M = 1.09_{-0.06}^{+0.13} M_{\odot}$.

As mentioned also in Section 3.5, only the R21 and
$R_{\text{eq}} \in [10, 14]$ km results are likely converged in terms of
improving sampler settings. Both of them infer practically the
same geometry solution (see the remaining posteriors in the
Appendix) and yield almost equally good fit quality and
evidence (see Table 1). On the contrary, the maximum-
likelihood and evidence values are significantly worse for the
$R_{\text{eq}} \in [8, 16]$ km prior run: by 7 and 17, respectively (both in
ln-units) when comparing to the corresponding R21 run. This
indicates that the wider radius prior run did not explore the
parameter space sufficiently, and thus inferred a completely
different (wrong) geometry solution. This was also the case in
our early stage 4000 live point run using the standard
(broadest) radius prior for PDT−U (with some older settings).
In that case, both 5 and 15 km radius modes were visibly
detected in the posteriors,<fnref n="25" /> but both of them are significantly
worse than the R21 and $R_{\text{eq}} \in [10, 14]$ km results in terms of
likelihood and evidence. From Table 1, we also see that the
evidence and maximum likelihood are highest for the He R21
run: by 8 and 3 in ln-units, respectively, compared to the H R21
run. All of the credible intervals for the two atmosphere options
are, however, very consistent with each other (with He giving
often broader intervals).

<footnote n="25"><sup>25</sup> The former was better in the likelihood but the latter in the evidence.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Radius, mass, compactness, and inclination posterior distributions using the NICER and XMM-Newton data sets for the PDT-U model. The results with three different radius priors are shown, but the run with widest radius prior ($R_{eq} \in [8, 16]$ km) seems not to have converged because it finds only significantly worse fits to the data than the others. Dashed-dotted curves represent the marginal prior probability density functions (PDFs). The vertical dashed orange line shows $R_{eq} = 14$ km to guide the eye. The shaded vertical bands show the 68.3% credible intervals (for the posteriors corresponding to the red curves), and the contours in the off-diagonal panels show the 68.3%, 95.4%, and 99.7% credible regions. See the captions of Figure 5 of T. E. Riley et al. (2021) for additional details about the figure elements.</figcaption></figure>

The performance of the PDT-U model was found to be good, both in terms of the 2D NICER residuals (see Figure 6) and the bolometric $\chi^2$ value of the maximum-likelihood parameter

vector (see the right-hand panel of Figure 4). The $\chi^2$ value was around 18 for the He R21 run, 19–20 for the H R21 and $R_{eq} \in [10, 14]$ km runs, and 21 for the nonconverged
<table><caption>**Table 1**
Model Performance Measures</caption><thead><tr><th>Model</th><th>Δln Z̃</th><th>Δln ML̃</th><th>χ<sup>2</sup><sub>ML,bol</sub></th></tr></thead><tbody><tr><td>PDT-U He R<sub>eq</sub> R21</td><td>0</td><td>0</td><td>18.2</td></tr><tr><td>PDT-U H R<sub>eq</sub> ∈ [10, 14] km</td><td>-7.35</td><td>-2.39</td><td>19.8</td></tr><tr><td>PDT-U H R<sub>eq</sub> R21</td><td>-8.13</td><td>-2.72</td><td>19.5</td></tr><tr><td>PDT-U H R<sub>eq</sub> ∈ [8, 16] km</td><td>-24.73</td><td>-10.16</td><td>21.7</td></tr><tr><td>ST-U H T<sub>else</sub></td><td>-32.56</td><td>-48.05</td><td>77.4</td></tr><tr><td>ST-U H 3α</td><td>-35.98</td><td>-58.78</td><td>80.7</td></tr><tr><td>ST-U H</td><td>-36.56</td><td>-58.83</td><td>85.5</td></tr></tbody></table>

**Note.** The evidence (ln Z̃) and maximum-likelihood (ln ML̃) differences are shown compared to the PDT-U He R<sub>eq</sub> R21 model. We also show the bolometric NICER χ² values (χ²<sub>ML,bol</sub>) corresponding the maximum-likelihood sample. Note that the exact d.o.f. is not known, but it is something smaller than 32 (see Footnote 24). For ST-U, only runs with the broadest radius prior are considered here; besides the standard case (ST-U H), values are also shown when including emission from rest of the star surface (ST-U H T<sub>else</sub>) and when having three different XMM-Newton effective-area scaling factors (ST-U H 3α).

H R<sub>eq</sub> ∈ [8, 16] km run (all again with 32 phase bins). We also checked if the inferred NICER background (based mostly on simultaneous XMM-Newton fits) is consistent with the predicted SCORPEON background in Figure 7. We see that it generally is, except at the lowest energies (close to 0.3 keV) where the inferred background is somewhat below, at 0.4–0.6 keV where it is somewhat above, and at 1.0–1.4 keV where it is slightly below the SCORPEON range. Inferring a slightly higher background is not surprising because the SCORPEON estimate does not include any unpulsed and unmodeled emission from the NS nor any emission from nearby sources (see the XMM-Newton image in Figure 8 of S. Bogdanov et al. 2019a). In addition, as noted in Section 2.1, the SCORPEON background estimate was found to depend on the chosen source spectral model (especially at the lowest energies), and none of the tried XSPEC models matches exactly the model used in X-PSI and its inferred parameters. Therefore, we consider the agreement reasonable. We also note that in our early NICER-only runs (without XMM-Newton) the inferred background was found to be substantially too large (covering all the unpulsed signal).

The contribution from the two different hot regions to the total pulse can be viewed from Figure 8, both for the incident signal and for that registered by NICER (without background). We see that the pulse from the primary region has a smoother shape than that from the secondary region, which remains mostly flat between the pulse peaks. The corresponding inferred hot region geometry is shown in Figure 9 for the maximum-likelihood sample of the PDT-U H R<sub>eq</sub> ∈ [10, 14] km run. The primary region is located on the northern hemisphere and consists of a rather small hot component surrounded by a very large cold component. The secondary region is located close to the southern pole and consists of almost equally hot but barely connected components with similar sizes. We note that the maximum-likelihood geometry is very similar for the PDT-U H R<sub>eq</sub> R21 and PDT-U He R<sub>eq</sub> R21 runs, except that the small hot component of the primary region is located at the equator in the latter case (which is still consistent with the H cases considering the width of the credible intervals).

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The NICER count data (top panel), posterior-expected count numbers (middle panel, averaged from 200 equally weighted posterior samples), and (Poisson) residuals (bottom panel) for the PDT-U H R<sub>eq</sub> ∈ [10, 14] km model. See Figure 6 of T. E. Riley et al. (2021) for additional details about the figure elements. The complete figure set (seven images), which includes the residuals for all of the models listed in Table 1, is available in the online journal (HTML version).</figcaption></figure>

(The complete figure set (seven images) is available in the online article.)

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Inferred NICER background for the joint NICER and XMM-Newton analysis using the PDT-U H R<sub>eq</sub> ∈ [10, 14] km model. The blue solid stepped curve shows the total NICER count-rate spectrum. The orange stepped curves show the backgrounds that maximize the likelihood for 1000 equally weighted posterior samples. The red curve shows the corresponding background for the maximum-likelihood sample. The black-dashed (dotted-dashed) curve shows the highest (lowest) background when considering the model uncertainty, as explained in Section 2.1. The inferred backgrounds of the other presented models are very similar to those shown here.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Posterior-expected pulse profiles for the joint NICER and XMM-Newton analysis using the *PDT-U* H $R_{\text{eq}} \in [10, 14]$ km model. The top and top-center panels show the signal incident on the telescopes, and the bottom-center and bottom panels show the signal registered by NICER after accounting for the instrument response. The blue-dashed curves show the posterior-expected signals generated by the two hot regions separately and the black-solid curves in combination. The conditional posterior distribution of the incident photon flux and the NICER count rate are also presented (purple contours in the top-center and bottom panels, respectively). See Figure 11 of T. E. Riley et al. (2021) for additional details about the figure elements.</figcaption></figure>

The inferred values for the remaining parameters, e.g., for *D* and *N*<sub>H</sub>, are found to be well consistent with expectations (see Section 5.1). For all the models, we find the distance to be around *D* = 0.65 ± 0.05 kpc. Inferred column densities show slightly more scatter, but all of the converged *PDT-U* runs provide roughly *N*<sub>H</sub> = 1.8 ± 0.5 × 10<sup>20</sup> cm<sup>-2</sup>.

## 5. Discussion

As shown in Section 4, our results are rather inconclusive due to the sampling challenges encountered when exploring the *PDT-U* model parameters with a wide prior for NS radius. However, at least using a radius prior based on the existing observations and theory (G. Raaijmakers et al. 2021a), or limiting the radius between 10 and 14 km, seems to lead to converged results and provides the best fits we have found to the PSR J1231−1411 data. Next, we discuss what these results imply (in Section 5.1), what factors may make analyzing

PSR J1231−1411 harder than the previous NICER sources (in Section 5.2), and some caveats (in Section 5.3).

### 5.1. Implications

One of the more robust findings of this work is that the preferred hot region model is *PDT-U* rather than *ST-U*. This means that assuming two circular uniformly emitting hot regions is not sufficient and may lead to wrong estimates of the NS parameters. Assuming a NS radius consistent with the previous EOS analysis and using the *PDT-U* model, the inferred hot region geometry is also constrained to be rather nonantipodal (see Figure 9), which implies that purely centered dipolar magnetic field is unlikely for PSR J1231−1411 (as for other NICER RMPs analyzed so far). The large cold emission component surrounding the primary hot spot may also be associated with emission from the rest of the star’s surface, hinting that using a *ST+PDT* model (see S. Vinciguerra et al. 2024 for a definition) and allowing the rest of the star surface to emit might lead to similar results. However, this is expected to speed-up the computation only mildly, and was thus not tried here.

We can also note a few interesting details about the inferred mass and radius. We see from Figure 5 that all of the inferred masses are quite small and close to our $1 M_{\odot}$ lower limit (this limit is motivated by supernova theory, see B. Müller et al. 2024, and by the description of the early evolution of a NS, see K. Strobel et al. 1999). These mass constraints are notably tighter than the mass prior which we used from radio timing measurements but they are still in the region where the prior is highest (see Figure 3). The radius constraints, on the other hand, do not provide as much information gain over the prior, especially in the R21 case. There the Kullback-Leibler (KL) divergence<sup>26</sup> is 0.55, implying a nonzero information given by the data. In the H $R_{\text{eq}} \in [10, 14]$ km case, the KL-divergence of radius is notably better (1.40) but still worse than that of mass (2.55).

Additionally, our results offer some insights into the NS atmosphere composition by finding better evidence and likelihoods for helium rather than hydrogen when using the *PDT-U* model and R21 radius prior (see Section 4.2, the same was, however, also detected in preliminary nonrestricted *ST-U* analyses). However, the difference is not large enough to be conclusive considering that hydrogen composition is a priori more probable based on the rapid sinking of heavier elements via diffusive gravitational separation (C. Alcock & A. Illarionov 1980; J. M. Hameury et al. 1983; E. F. Brown et al. 2002; V. E. Zavlin & G. G. Pavlov 2002). A helium atmosphere could occur if all the hydrogen was converted to helium via diffusive nuclear burning (M. J. P. Wijngaarden et al. 2019) or if the star had only accreted from a completely hydrogen-depleted companion star (see S. Bogdanov et al. 2021; T. Salmi et al. 2023, and the references therein). Nevertheless, there is evidence that the companion of PSR J1231−1411 is a white dwarf with a hydrogen envelope (C. G. Bassa et al. 2016).

The $N_{\text{H}}$ column density can be estimated via independent means and compared to the inferred values. Neutral hydrogen maps predict values ~3.4 × 10<sup>20</sup> cm<sup>-2</sup> (HI4PI Collaboration et al. 2016), but this should be considered an approximate upper limit because this represents the column density integrated to the edge of the Galaxy in that direction.<sup>27</sup> Other

<footnote n="26"><sup>26</sup> Representing prior-to-posterior information gain.</footnote>

<footnote n="27"><sup>27</sup> For this, we used the HEASARC $N_{\text{H}}$ tool at https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The maximum-likelihood geometry configuration for the joint NICER and XMM-Newton analysis using the PDT-U H $R_{\text{eq}} \in [10, 14]$ km model. The viewing angle represents the Earth's inclination to the spin axis. The maximum-likelihood geometry is very similar for the PDT-U H $R_{\text{eq}}$ R21 and PDT-U He $R_{\text{eq}}$ R21 models, except that the primary superseding component is located at the equator for the latter.</figcaption></figure>

estimates of $N_{\text{H}}$, at the distance of the pulsar, can be obtained from 3D maps of infrared or optical extinction $A_{\text{V}}$ or reddening $E(B-V)$. At the position and distance of PSR J1231−1411, ~0.6 kpc, we find ~3.5 × 10²⁰ cm⁻² based on the 3D dust maps of G. M. Green et al. (2019)<fnref n="28" /> this 3D map shows little dependence of $E(B-V)$ on the distance beyond 0.2 kpc). Alternatively, with the 3D dust maps of J. L. Vergely et al. (2022) and R. Lallement et al. (2022)<fnref n="29" /> we find ~3.4 × 10²⁰ cm⁻². Both results make use of the relation of D. R. Foight et al. (2016) between $N_{\text{H}}$ and $A_{\text{V}}$. Another estimate can exploit the relation between the DM (number density of free electrons in the line of sight) and $N_{\text{H}}$ (C. He et al. 2013). For PSR J1231−1411, $\text{DM} = 8.09 \text{ pc cm}^{-3}$, which gives ~2.4 × 10²⁰ cm⁻².

The $N_{\text{H}}$ values inferred in this work are globally consistent with these independent measurements. It is important to keep in mind that the measurements deduced above should be considered carefully because the employed relations may be uncertain, especially at low-absorption values because the sampling of objects is much smaller for low values of $N_{\text{H}}$. In other words, these independent measurements mostly bring support that the chosen prior on $N_{\text{H}}$ is broad enough and does not bias the inferred parameters.

Our inferred distance is also found to be well consistent with the preliminary parallax distance measurement (T. Cromartie et al. 2024, in preparation), both being around or above 600 pc. We note though that our distance upper limit (700 pc) might be too restrictive because our posteriors peak only slightly below it and the parallax distance uncertainties can also allow higher distances, especially when accounting for the Lutz-Kelker bias (T. E. Lutz & D. H. Kelker 1973; J. P. W. Verbiest et al. 2012). Using a more sophisticated distance prior is, however, left to future work.

In one of the ST-U runs we also explored the effect of separate effective-area scaling factors for all the XMM-Newton instruments. We found $\alpha_{\text{pn}} = 0.93 \pm 0.09$, $\alpha_{\text{MOS1}} = 0.98 \pm 0.09$, and $\alpha_{\text{MOS2}} = 1.00 \pm 0.09$. These are similar enough with each other

that using one single $\alpha_{\text{XMM}}$ seems justified for the final PDT-U models.

## 5.2. Comparison to Other Pulsars

Compared to the other NICER RMPs analyzed so far, constraining the NS parameters for PSR J1231−1411 appears more challenging. There are several possible reasons for this. First of all, the available priors on mass and observer inclination are much less restrictive for PSR J1231−1411 than for PSR J0740+6620 and PSR J0437−4715. For the isolated pulsar PSR J0030+0451, there is even less prior information about mass and inclination, but the distance is known much more precisely based on radio timing parallax measurements (Z. Arzoumanian et al. 2018). These factors, together with the relatively large number of free hot region geometry parameters needed for the PDT-U model, likely make the sampling problem more difficult. In the previous X-PSI analyses, the simple ST-U model was found to be adequate for PSR J0740+6620 (T. E. Riley et al. 2021) and the CST+PDT model, with three parameters less than for PDT-U, for PSR J0437−4715 (D. Choudhury et al. 2024a). In S. Vinciguerra et al. (2024) PDT-U was one of the models applied to PSR J0030+0451, however its convergence was not investigated.

PSR J1231−1411 also differs from the other pulsars in terms of the pulse shape. As seen from Figures 1 and 8, there is only a weak interpulse appearing in the bolometric pulse profile. This feature is especially weak in the 2D data, leading possibly to a multitude of solutions that can more or less fit it but that are too distinct in the parameter space to be explored thoroughly if one does not shrink the prior space. For PSR J0030+0451 and PSR J0740+6620, the secondary pulse is much stronger, and for PSR J0437−4715 there is no peaky interpulse (see S. Bogdanov et al. 2019a and M. T. Wolff et al. 2021). However, it could also be that the source properties of PSR J1231−1411 simply differ from the others in such a way that they are more difficult to recover.

It is nevertheless encouraging that despite the uncertainty, we are able to obtain solutions for this source that fit the data well and are also consistent with the results from other NICER RMPs and multimessenger analyses. There are clearly various steps that might improve the situation: tighter priors, additional

<footnote n="28">²⁸ Query tool is available at http://argonaut.skymaps.info/.</footnote>

<footnote n="29">²⁹ Data are available on Vizier at (https://vizier.cds.unistra.fr) or via an online tool at https://explore-platform.eu/sdas/about/gtomo.</footnote>
computational resources, or a longer data set with more
photons. However, before devoting additional resources (in the
form of observing or computational time) to this source, we
would recommend a program of simulations using synthetic
data to determine whether this is likely to lead to improvements
or whether there are fundamental issues with this type of pulse
profile that are likely to prevent efforts to obtain tighter and
more robust constraints.

5.3. Caveats

The caveats of this work mostly follow those present in the
previous NICER analyses as well (see e.g., S. Bogdanov et al.
2021, for discussion). Additionally, we are now dependent on
more restrictive radius priors, e.g., based on the results of
G. Raaijmakers et al. (2021a), because the radius could not be
constrained fully independently in this work. Therefore, the
existence of better solutions far from the R21 prior peak, which
is around 12 km, or at least outside from the [10, 14] km range,
cannot be ruled out (even though the hard limits for the R21
case were still kept at ~4.4 and 16 km).

We also note that our results rely on $M - \cos i$ priors from a
preliminary radio timing inference (T. Cromartie et al. 2024, in
preparation). In the early analyses without these priors, we still
found mostly similar (small) masses but the inclination angles were
found to be much smaller, with $\cos i$ peaking between 0.4 and 0.9
(depending on the model) instead of at around 0.15 as predicted by
the radio timing. However, these analyses did not yet employ any
restricted radius prior for the PDT-U model. We also note that the
inclination angle from radio timing refers to the angle between the
line of sight and the orbital axis, whereas we use it as the angle
between line of sight and NS spin axis. However, previous
accretion of matter is expected to have gradually aligned the spin
and orbital axes of the system (see, e.g., D. Bhattacharya &
E. P. J. van den Heuvel 1991), and any remaining tilt between them
is likely to be insignificant compared to the width of the prior.

We have not tested if setting the lower limit of mass below $1 M_{\odot}$ would affect our results. Recent spectral modeling of HESS J1731-347 has suggested that some NSs could have such low masses (V. Doroshenko et al. 2022). However, this interpretation relies on several assumptions, including the distance associated with the star and the use of a uniform-temperature carbon atmosphere model, as discussed by J. A. J. Alford & J. P. Halpern (2023). More importantly, this model is a poor fit to longer, better quality XMM-Newton data from 2014 of the same source, while a two-blackbody model with a mass of, e.g., $1.4 M_{\odot}$, provides a good fit to the same data (J. A. J. Alford & J. P. Halpern 2023). Therefore, current data do not support a mass less than $1 M_{\odot}$ for HESS J1731-347.

Finally, we note that the inferred temperatures of some of the
hot regions are rather cold (notably below $10^6$ K), making the
effects or partial ionization more relevant in the NS atmosphere,
whereas we assume it to be fully ionized both for the hydrogen
and helium case. However, the effect is still expected to be minor
because the signal is dominated by the hottest regions. In addition,
switching between hydrogen and helium was not found to
significantly alter the results (for the R21 case), implying only
minor dependence on the atmosphere choice.

**6. Conclusions**

We have jointly analyzed NICER and XMM-Newton data
for PSR J1231−1411 using X-ray pulse profile modeling. We

found that the simple ST-U model with two uniformly emitting
spherical caps is not sufficient to explain the data, especially
the bolometric pulse profile. Instead, good fits both in
bolometric and energy-resolved data and significantly larger
Bayesian evidence were found using the PDT-U model, where
both hot regions consist of two circular components that are
allowed to have different temperatures. Unfortunately, we did
not manage to get converged results using a wide NS radius
prior with this computationally expensive model. However,
forcing the radius to be consistent with previous multi-
messenger observations and nuclear theory, we found
$R_{\text{eq}} = 12.6 \pm 0.3$ km and $M = 1.04_{-0.03}^{+0.05} M_{\odot}$ and a substantial
improvement in evidence and likelihoods. When assuming
helium (instead of hydrogen) composition for the NS
atmosphere, we were also able reproduce the data well and
obtained very similar parameter constraints. Equally good fits,
and very similar results, were also found for the hydrogen
atmosphere when applying an uninformative prior limited
between 10 and 14 km, except that the radius was inferred to be
$R_{\text{eq}} = 13.5_{-0.5}^{+0.3}$ km. In all the best-fitting models we infer a
quite similar nonantipodal hot region geometry, implying that a
pure dipolar magnetic field is not likely for PSR J1231−1411.

Acknowledgments

This work was supported in part by NASA through the
NICER mission and the Astrophysics Explorers Program. T.S.,
A.L.W., D.C., Y.K., and S.V. acknowledge support from ERC
Consolidator grant No. 865768 AEONS (PI: Watts). The use of
the national computer facilities in this research was subsidized
by NWO Domain Science. Astrophysics research at the Naval
Research Laboratory is supported by the NASA Astrophysics
Explorer Program. S.G. acknowledges the support of the
CNES. W.C.G.H. acknowledges support through grant
80NSSC23K0078 from NASA. This research has made use
of data products and software provided by the High Energy
Astrophysics Science Archive Research Center (HEASARC),
which is a service of the Astrophysics Science Division at
NASA/GSFC and the High Energy Astrophysics Division of
the Smithsonian Astrophysical Observatory. We thank Thomas
Riley for help in the early stages of this project and Jason
Farquhar for discussions about statistical fits to weak features.

Facilities: NICER, XMM.

Software: Cython (S. Behnel et al. 2011), fgivenx (W. Handley 2018), GetDist (A. Lewis 2019), GNU Scientific Library (GSL; M. Galassi et al. 2009), HEASoft (NASA High Energy Astrophysics Science Archive Research Center (Heasarc), 2014), Matplotlib (J. D. Hunter 2007), MPI for Python (L. Dalcin et al. 2008), MULTINEST (F. Feroz et al. 2009), nestcheck (E. Higson 2018), NumPy (S. van der Walt et al. 2011), PYMULTINEST (J. Buchner et al. 2014), Python/C language (T. E. Oliphant 2007), SciPy (E. Jones et al. 2001), X-PSI (T. E. Riley et al. 2023).

**Appendix**
**Additional Tables and Figures**

We summarize the results for the best-fitting hydrogen models (PDT-U H $R_{\text{eq}} \in [10, 14]$ km and PDT-U H $R_{\text{eq}}$ R21) in Tables 2 and 3. Additional PDT-U posteriors, including the parameters not shown in Figure 5, are presented in Figures 10, 11, and 12. The main ST-U posterior distributions for radius, compactness, and mass are shown in Figure 13.
<table><caption>**Table 2**
Summary Table for the PDT-U H $R_{\text{eq}} \in [10, 14]$ km Run Described in Section 4.2</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior PDF (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th><th>M̂<sub>KL</sub></th></tr></thead><tbody><tr><td>P [ms]</td><td>Coordinate spin period</td><td>P = 3.6839, fixed</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M [M<sub>☉</sub>]</td><td>Gravitational mass</td><td>M ~ N(1.00233782, 0.93010502)<sup>a</sup></td><td>1.06<sup>+0.06</sup><sub>-0.04</sub></td><td>2.55</td><td>1.069</td></tr><tr><td>cos(i)</td><td>Cosine Earth inclination to spin axis</td><td>cos(i) ~ N(&mu;<sub>cos i</sub>(M), &sigma;<sub>cos i</sub>(M))<sup>b</sup></td><td>0.15<sup>+0.01</sup><sub>-0.02</sub></td><td>0.95</td><td>0.172</td></tr><tr><td>R<sub>eq</sub> [km]</td><td>Coordinate equatorial radius<br>With compactness condition<br>With surface gravity condition</td><td>R<sub>eq</sub> ~ U(10, 14)<br>R<sub>polar</sub>/r<sub>g</sub>(M) > 3<br>13.7 &le; log<sub>10</sub> g(&theta;) &le; 15.0, &forall;&theta;</td><td>13.52<sup>+0.32</sup><sub>-0.50</sub></td><td>1.40</td><td>13.007</td></tr><tr><td>&Theta;<sub>p</sub> [radians]</td><td>p superseding component center colatitude</td><td>cos(&Theta;<sub>p</sub>) ~ U(-1, 1)</td><td>1.29<sup>+0.21</sup><sub>-0.21</sub></td><td>0.93</td><td>1.239</td></tr><tr><td>&Theta;<sub>c,p</sub> [radians]</td><td>p ceding component center colatitude</td><td>cos(&Theta;<sub>c,p</sub>) ~ U(-1, 1)</td><td>0.23<sup>+0.13</sup><sub>-0.09</sub></td><td>3.32</td><td>0.123</td></tr><tr><td>&Theta;<sub>s</sub> [radians]</td><td>s superseding component center colatitude</td><td>cos(&Theta;<sub>s</sub>) ~ U(-1, 1)</td><td>2.78<sup>+0.03</sup><sub>-0.03</sub></td><td>4.72</td><td>2.761</td></tr><tr><td>&Theta;<sub>c,s</sub> [radians]</td><td>s ceding component center colatitude</td><td>cos(&Theta;<sub>c,s</sub>) ~ U(-1, 1)</td><td>2.81<sup>+0.04</sup><sub>-0.04</sub></td><td>4.46</td><td>2.844</td></tr><tr><td>&phi;<sub>p</sub> [cycles]</td><td>p superseding component initial phase</td><td>&phi;<sub>p</sub> ~ U(-0.25, 0.75), wrapped</td><td>0.38<sup>+0.01</sup><sub>-0.01</sub></td><td>5.25</td><td>0.375</td></tr><tr><td>&phi;<sub>s</sub> [cycles]</td><td>s superseding component initial phase</td><td>&phi;<sub>s</sub> ~ U(-0.25, 0.75), wrapped</td><td>0.53<sup>+0.01</sup><sub>-0.01</sub></td><td>5.29</td><td>0.536</td></tr><tr><td>&chi;<sub>p</sub> [radians]</td><td>Azimuthal offset between the p components</td><td>&chi;<sub>p</sub> ~ U(-&pi;, &pi;)</td><td>-0.09<sup>+0.26</sup><sub>-0.24</sub></td><td>0.41</td><td>0.010</td></tr><tr><td>&chi;<sub>s</sub> [radians]</td><td>Azimuthal offset between the s components</td><td>&chi;<sub>s</sub> ~ U(-&pi;, &pi;)</td><td>-1.63<sup>+0.05</sup><sub>-0.05</sub></td><td>6.89</td><td>-1.622</td></tr><tr><td>&zeta;<sub>p</sub> [radians]</td><td>p superseding component angular radius</td><td>&zeta;<sub>p</sub> ~ U(0, &pi;/2)</td><td>0.10<sup>+0.02</sup><sub>-0.02</sub></td><td>2.85</td><td>0.094</td></tr><tr><td>&zeta;<sub>c,p</sub> [radians]</td><td>p ceding component angular radius</td><td>&zeta;<sub>c,p</sub> ~ U(0, &pi;/2)</td><td>1.44<sup>+0.09</sup><sub>-0.15</sub></td><td>4.28</td><td>1.548</td></tr><tr><td>&zeta;<sub>s</sub> [radians]</td><td>s superseding region angular radius</td><td>&zeta;<sub>s</sub> ~ U(0, &pi;/2)</td><td>0.26<sup>+0.03</sup><sub>-0.02</sub></td><td>3.15</td><td>0.237</td></tr><tr><td>&zeta;<sub>c,s</sub> [radians]</td><td>s ceding component angular radius</td><td>&zeta;<sub>c,s</sub> ~ U(0, &pi;/2)</td><td>0.28<sup>+0.02</sup><sub>-0.02</sub></td><td>3.16</td><td>0.266</td></tr><tr><td></td><td>No region-exchange degeneracy</td><td>&Theta;<sub>s</sub> &ge; &Theta;<sub>p</sub></td><td></td><td></td><td></td></tr><tr><td></td><td>Nonoverlapping hot regions</td><td>function of all &Theta;, &phi;, &chi;, and &zeta;</td><td></td><td></td><td></td></tr><tr><td></td><td>Overlapping hot region components<sup>c</sup></td><td>function of all &Theta;, &phi;, &chi;, and &zeta;</td><td></td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</sub> [K])</td><td>p superseding component effective temperature</td><td>log<sub>10</sub>(T<sub>p</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.90<sup>+0.02</sup><sub>-0.03</sub></td><td>4.05</td><td>5.918</td></tr><tr><td>log<sub>10</sub>(T<sub>c,p</sub> [K])</td><td>p ceding component effective temperature</td><td>log<sub>10</sub>(T<sub>c,p</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.53<sup>+0.03</sup><sub>-0.03</sub></td><td>3.82</td><td>5.531</td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> [K])</td><td>s superseding component effective temperature</td><td>log<sub>10</sub>(T<sub>s</sub>) ~ U(5.1, 6.8), NSX limits</td><td>6.02<sup>+0.01</sup><sub>-0.01</sub></td><td>5.50</td><td>6.027</td></tr><tr><td>log<sub>10</sub>(T<sub>c,s</sub> [K])</td><td>s ceding component effective temperature</td><td>log<sub>10</sub>(T<sub>c,s</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.96<sup>+0.02</sup><sub>-0.02</sub></td><td>4.54</td><td>5.987</td></tr><tr><td>D [kpc]</td><td>Earth distance</td><td>D ~ U(0.1, 0.7)</td><td>0.65<sup>+0.04</sup><sub>-0.05</sub></td><td>2.19</td><td>0.620</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, 10)</td><td>1.77<sup>+0.50</sup><sub>-0.50</sub></td><td>2.29</td><td>2.101</td></tr><tr><td>&alpha;<sub>NICER</sub></td><td>NICER effective-area scaling</td><td>&alpha;<sub>NICER</sub>, &alpha;<sub>XMM</sub> ~ N(&mu;, &Sigma;)</td><td>0.98<sup>+0.09</sup><sub>-0.07</sub></td><td>0.14</td><td>1.071</td></tr><tr><td>&alpha;<sub>XMM</sub></td><td>XMM-Newton effective-area scaling<br>With joint prior PDF N(&mu;, &Sigma;)</td><td>&alpha;<sub>NICER</sub>, &alpha;<sub>XMM</sub> ~ N(&mu;, &Sigma;)<br>&mu; = [1.0, 1.0]<sup>T</sup><br>&Sigma; = [0.104<sup>2</sup> 0.100<sup>2</sup><br>0.100<sup>2</sup> 0.104<sup>2</sup>]</td><td>0.96<sup>+0.08</sup><sub>-0.08</sub></td><td>0.22</td><td>0.970</td></tr></tbody></table><p><strong>Sampling Process Information</strong></p><p>Number of free parameters: 23</p><p>Number of processes (multimodes): 4</p><p>Number of live points: 2 &times; 10<sup>4</sup></p><p>Sampling efficiency (SE): 0.1</p><p>Termination condition: 0.1</p><p>Evidence: ln Z = -21352.42 &plusmn; 0.56</p><p>Number of core hours: 5.34533 &times; 10<sup>5</sup></p><p>Likelihood evaluations: 4.83928845 &times; 10<sup>8</sup></p>

**Notes.** We show the prior PDFs, 68.3% credible intervals around the median $\hat{\text{CI}}_{68\%}$, KL-divergence $\hat{D}_{\text{KL}}$ in *bits*, and the maximum-likelihood nested sample $\hat{\text{ML}}$ for all the parameters. The subscripts *p* and *s* denote for primary and secondary hot region parameters, respectively. Note that $\phi_p$ ($\phi_s$) is measured with respect to the meridian on which Earth (Earth antipode) lies.

a Truncated between 1 and 3 $M_{\odot}$. However, note that for each $\cos i$ the *M* prior is different, as seen in Figure 3.

b See the definitions of $\mu_{\cos i}$ and $\sigma_{\cos i}$ in Equations (1) and (2).

c Ensuring that superseding components never engulf the corresponding ceding components.
<table><caption>**Table 3**
Summary Table for the PDT-U H R<sub>eq</sub> R21 Run Described in Section 4.2</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior PDF (Density and Support)</th><th>CI<sub>68%</sub></th><th>D̂<sub>KL</sub></th><th>M̂<sub>L</sub></th></tr></thead><tbody><tr><td>P [ms]</td><td>Coordinate spin period</td><td>P = 3.6839, fixed</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M [M<sub>☉</sub>]</td><td>Gravitational mass</td><td>M ~ N(1.00233782, 0.93010502)<sup>a</sup></td><td>1.04<sup>+0.05</sup><sub>-0.03</sub></td><td>3.00</td><td>1.002</td></tr><tr><td>cos(i)</td><td>Cosine Earth inclination to spin axis</td><td>cos(i) ~ N(&mu;<sub>cos i</sub>(M), &sigma;<sub>cos i</sub>(M))<sup>b</sup></td><td>0.15<sup>+0.01</sup><sub>-0.02</sub></td><td>1.02</td><td>0.165</td></tr><tr><td>R<sub>eq</sub> [km]</td><td>Coordinate equatorial radius<br>With compactness condition<br>With surface gravity condition</td><td>G. Raaijmakers et al. (2021a)<sup>c</sup><br>R<sub>polar</sub>/r<sub>g</sub>(M) > 3<br>13.7 &le; log<sub>10</sub> g(&theta;) &le; 15.0, &forall;&theta;</td><td>12.60<sup>+0.31</sup><sub>-0.33</sub></td><td>0.55</td><td>12.859</td></tr><tr><td>&Theta;<sub>p</sub> [radians]</td><td>p superseding component center colatitude</td><td>cos(&Theta;<sub>p</sub>) ~ U(-1, 1)</td><td>1.34<sup>+0.19</sup><sub>-0.21</sub></td><td>1.01</td><td>1.172</td></tr><tr><td>&Theta;<sub>c,p</sub> [radians]</td><td>p ceding component center colatitude</td><td>cos(&Theta;<sub>c,p</sub>) ~ U(-1, 1)</td><td>0.26<sup>+0.16</sup><sub>-0.10</sub></td><td>2.90</td><td>0.155</td></tr><tr><td>&Theta;<sub>s</sub> [radians]</td><td>s superseding component center colatitude</td><td>cos(&Theta;<sub>s</sub>) ~ U(-1, 1)</td><td>2.77<sup>+0.03</sup><sub>-0.03</sub></td><td>4.67</td><td>2.716</td></tr><tr><td>&Theta;<sub>c,s</sub> [radians]</td><td>s ceding component center colatitude</td><td>cos(&Theta;<sub>c,s</sub>) ~ U(-1, 1)</td><td>2.79<sup>+0.04</sup><sub>-0.04</sub></td><td>4.35</td><td>2.827</td></tr><tr><td>&phi;<sub>p</sub> [cycles]</td><td>p superseding component initial phase</td><td>&phi;<sub>p</sub> ~ U(-0.25, 0.75), wrapped</td><td>0.38<sup>+0.01</sup><sub>-0.01</sub></td><td>5.30</td><td>0.375</td></tr><tr><td>&phi;<sub>s</sub> [cycles]</td><td>s superseding component initial phase</td><td>&phi;<sub>s</sub> ~ U(-0.25, 0.75), wrapped</td><td>0.53<sup>+0.01</sup><sub>-0.01</sub></td><td>5.31</td><td>0.541</td></tr><tr><td>&chi;<sub>p</sub> [radians]</td><td>Azimuthal offset between the p components</td><td>&chi;<sub>p</sub> ~ U(-&pi;, &pi;)</td><td>-0.04<sup>+0.25</sup><sub>-0.26</sub></td><td>0.37</td><td>0.219</td></tr><tr><td>&chi;<sub>s</sub> [radians]</td><td>Azimuthal offset between the s components</td><td>&chi;<sub>s</sub> ~ U(-&pi;, &pi;)</td><td>-1.65<sup>+0.05</sup><sub>-0.05</sub></td><td>6.95</td><td>-1.550</td></tr><tr><td>&zeta;<sub>p</sub> [radians]</td><td>p superseding component angular radius</td><td>&zeta;<sub>p</sub> ~ U(0, &pi;/2)</td><td>0.12<sup>+0.02</sup><sub>-0.02</sub></td><td>2.82</td><td>0.093</td></tr><tr><td>&zeta;<sub>c,p</sub> [radians]</td><td>p ceding component angular radius</td><td>&zeta;<sub>c,p</sub> ~ U(0, &pi;/2)</td><td>1.44<sup>+0.09</sup><sub>-0.14</sub></td><td>4.40</td><td>1.412</td></tr><tr><td>&zeta;<sub>s</sub> [radians]</td><td>s superseding region angular radius</td><td>&zeta;<sub>s</sub> ~ U(0, &pi;/2)</td><td>0.27<sup>+0.03</sup><sub>-0.02</sub></td><td>3.20</td><td>0.234</td></tr><tr><td>&zeta;<sub>c,s</sub> [radians]</td><td>s ceding component angular radius</td><td>&zeta;<sub>c,s</sub> ~ U(0, &pi;/2)</td><td>0.29<sup>+0.02</sup><sub>-0.02</sub></td><td>3.21</td><td>0.285</td></tr><tr><td></td><td>No region-exchange degeneracy</td><td>&Theta;<sub>s</sub> &ge; &Theta;<sub>p</sub></td><td></td><td></td><td></td></tr><tr><td></td><td>Nonoverlapping hot regions</td><td>function of all &Theta;, &phi;, &chi;, and &zeta;</td><td></td><td></td><td></td></tr><tr><td></td><td>Overlapping hot region components<sup>d</sup></td><td>function of all &Theta;, &phi;, &chi;, and &zeta;</td><td></td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>p</sub> [K])</td><td>p superseding component effective temperature</td><td>log<sub>10</sub>(T<sub>p</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.90<sup>+0.02</sup><sub>-0.02</sub></td><td>4.09</td><td>5.920</td></tr><tr><td>log<sub>10</sub>(T<sub>c,p</sub> [K])</td><td>p ceding component effective temperature</td><td>log<sub>10</sub>(T<sub>c,p</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.53<sup>+0.02</sup><sub>-0.03</sub></td><td>3.77</td><td>5.556</td></tr><tr><td>log<sub>10</sub>(T<sub>s</sub> [K])</td><td>s superseding component effective temperature</td><td>log<sub>10</sub>(T<sub>s</sub>) ~ U(5.1, 6.8), NSX limits</td><td>6.02<sup>+0.01</sup><sub>-0.01</sub></td><td>5.55</td><td>6.018</td></tr><tr><td>log<sub>10</sub>(T<sub>c,s</sub> [K])</td><td>s ceding component effective temperature</td><td>log<sub>10</sub>(T<sub>c,s</sub>) ~ U(5.1, 6.8), NSX limits</td><td>5.96<sup>+0.02</sup><sub>-0.02</sub></td><td>4.50</td><td>5.980</td></tr><tr><td>D [kpc]</td><td>Earth distance</td><td>D ~ U(0.1, 0.7)</td><td>0.65<sup>+0.03</sup><sub>-0.05</sub></td><td>2.41</td><td>0.633</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, 10)</td><td>1.62<sup>+0.47</sup><sub>-0.48</sub></td><td>2.35</td><td>1.905</td></tr><tr><td>&alpha;<sub>NICER</sub></td><td>NICER effective-area scaling</td><td>&alpha;<sub>NICER</sub>, &alpha;<sub>XMM</sub> ~ N(&mu;, &Sigma;)</td><td>0.98<sup>+0.09</sup><sub>-0.08</sub></td><td>0.13</td><td>1.040</td></tr><tr><td>&alpha;<sub>XMM</sub></td><td>XMM-Newton effective-area scaling<br>With joint prior PDF N(&mu;, &Sigma;)</td><td>&alpha;<sub>NICER</sub>, &alpha;<sub>XMM</sub> ~ N(&mu;, &Sigma;)<br>&mu; = [1.0, 1.0]<sup>T</sup><br>&Sigma; = [ 0.104<sup>2</sup> 0.100<sup>2</sup> 0.100<sup>2</sup> 0.104<sup>2</sup> ]</td><td>0.95<sup>+0.08</sup><sub>-0.08</sub></td><td>0.26</td><td>0.950</td></tr><tr><td colspan="6"><strong>Sampling Process Information</strong></td></tr><tr><td colspan="6">Number of free parameters: 23</td></tr><tr><td colspan="6">Number of processes (multimodes): 4</td></tr><tr><td colspan="6">Number of live points: 2 &times; 10<sup>4</sup></td></tr><tr><td colspan="6">Sampling efficiency (SE): 0.1</td></tr><tr><td colspan="6">Termination condition: 0.1</td></tr><tr><td colspan="6">Evidence: ln Z&#773; = -21353.20 &plusmn; 0.56</td></tr><tr><td colspan="6">Number of core hours: 2.27230 &times; 10<sup>5</sup></td></tr><tr><td colspan="6">Likelihood evaluations: 5.47847609 &times; 10<sup>8</sup></td></tr></tbody></table>

**Notes.** We show the prior PDFs, 68.3% credible intervals around the median CI<sub>68%</sub>, KL-divergence D̂<sub>KL</sub> in bits, and the maximum-likelihood nested sample M̂<sub>L</sub> for all the parameters. The subscripts *p* and *s* denote for primary and secondary hot region parameters, respectively. Note that φ<sub>p</sub> (φ<sub>s</sub>) is measured with respect to the meridian on which Earth (Earth antipode) lies.

<sup>a</sup> Truncated between 1 and 3 M<sub>☉</sub>. However, note that for each cos i the *M* prior is different, as seen in Figure 3.

<sup>b</sup> See the definitions of μ<sub>cos i</sub> and σ<sub>cos i</sub> in Equations (1) and (2).

<sup>c</sup> Truncated between 3r<sub>g</sub>(1) ≈ 4.4 and 16 km. The prior density corresponds to the R21 case that was explained in Section 3.3.2.

<sup>d</sup> Ensuring that superseding components never engulf the corresponding ceding components.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** PDT-U posterior distributions for radius and the primary hot region parameters. See the caption of Figure 5 for additional details.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11. PDT-U posterior distributions for radius and the secondary hot region parameters. See the caption of Figure 5 for additional details.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** PDT-U posterior distributions for radius and the other parameters. See the caption of Figure 5 for additional details.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Posterior distributions for radius, mass, compactness, the logarithm of the surface gravity, and the cosine of the inclination for the ST-U runs. See the caption of Figure 5 for additional details.</figcaption></figure>

ORCID iDs

Tuomo Salmi https://orcid.org/0000-0001-6356-125X

Julia S. Deneva https://orcid.org/0000-0003-1226-0793

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

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

Devarshi Choudhury https://orcid.org/0000-0002-2651-5286

Yves Kini https://orcid.org/0000-0002-0428-8430

Serena Vinciguerra https://orcid.org/0000-0003-3068-6974

H. Thankful Cromartie https://orcid.org/0000-0002-6039-692X

Michael T. Wolff https://orcid.org/0000-0002-4013-5650

Zaven Arzoumanian https://orcid.org/0009-0008-6187-8753

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

Keith Gendreau https://orcid.org/0000-0001-7115-2819
Sebastien Guillot https://orcid.org/0000-0002-6449-106X

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

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

Ismaël Cognard https://orcid.org/0000-0002-1775-9692

Lucas Guillemot https://orcid.org/0000-0002-9049-8716

Gilles Theureau https://orcid.org/0000-0002-3649-276X

Matthew Kerr https://orcid.org/0000-0002-0893-4073

## References

Alcock, C., & Illarionov, A. 1980, ApJ, 235, 534

Alford, J. A. J., & Halpern, J. P. 2023, ApJ, 944, 36

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

Annala, E., Gorda, T., Hirvonen, J., et al. 2023, NatCo, 14, 8451

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Arons, J. 1981, ApJ, 248, 1099

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

Bassa, C. G., Antoniadis, J., Camilo, F., et al. 2016, MNRAS, 455, 3806

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

Bhattacharya, D., & van den Heuvel, E. P. J. 1991, PhR, 203, 1

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

Biswas, B. 2022, ApJ, 926, 75

Bogdanov, S., Dittmann, A. J., Ho, W. C. G., et al. 2021, ApJL, 914, L15

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

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

Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920

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

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

Carrasco, F., Pelle, J., Reula, O., Vigano, D., & Palenzuela, C. 2023, MNRAS, 520, 3151

Choudhury, D., Salmi, T., Vinciguerra, S., et al. 2024a, ApJL, 971, L20

Choudhury, D., Watts, A. L., Dittmann, A. J., et al. 2024b, ApJ, 975, 202

Dalcin, L., Paz, R., Storti, M., & D'Elia, J. 2008, JPDC, 68, 655

Dittmann, A. J., Miller, M. C., Lamb, F. K., et al. 2024, ApJ, 974, 295

Doroshenko, V., Suleimanov, V., Pühlhofer, G., & Santangelo, A. 2022, NatAs, 6, 1444

FERMI-LAT Collaboration, Ajello, M., Atwood, W. B., et al. 2022, Sci, 376, 521

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. 2019, OJAp, 2, 10

Foight, D. R., Güver, T., Özel, F., & Slane, P. O. 2016, ApJ, 826, 66

Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJL, 915, L12

Galassi, M., Davies, J., Theiler, J., et al. 2009, GNU Scientific Library Reference Manual (3rd ed.; Godalming: Network Theory Ltd.)

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

Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93

Guillemot, L., Cognard, I., van Straten, W., Theureau, G., & Gérard, E. 2023, A&A, 678, A79

Hameury, J. M., Heyvaerts, J., & Bonazzola, S. 1983, A&A, 121, 259

Han, M.-Z., Huang, Y.-J., Tang, S.-P., & Fan, Y.-Z. 2023, SciBu, 68, 913

Handley, W. 2018, JOSS, 3, 849

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

He, C., Ng, C. Y., & Kaspi, V. M. 2013, ApJ, 768, 64

HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116

Higson, E. 2018, JOSS, 3, 916

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

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

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

Lallement, R., Vergely, J. L., Babusiaux, C., & Cox, N. L. J. 2022, A&A, 661, A147

Lattimer, J. M., & Prakash, M. 2007, PhR, 442, 109

Lewis, A. 2019, arXiv:1910.13970

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

Lutz, T. E., & Kelker, D. H. 1973, PASP, 85, 573

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

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

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJL, 918, L28

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

Müller, B., Heger, A., & Powell, J. 2024, arXiv:2407.08407

NASA High Energy Astrophysics Science Archive Research Center (Heasarc), 2014 HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

Nath, N. R., Strohmayer, T. E., & Swank, J. H. 2002, ApJ, 564, 353

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

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

Poutanen, J., & Gierlinski, M. 2003, MNRAS, 343, 1301

Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021a, ApJL, 918, L29

Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021b, Constraints on the Dense Matter Equation of State and Neutron Star Properties from NICER's Mass-Radius Estimate of PSR J0740+6620 and Multimessenger Observations: Posterior Samples and Scripts for Generating Plots, v1, Zenodo, doi:10.5281/zenodo.4696232

Ransom, S. M., Ray, P. S., Camilo, F., et al. 2011, ApJL, 727, L16

Ray, P. S., Guillot, S., Ransom, S. M., et al. 2019, ApJL, 878, L22

Reardon, D. J., Bailes, M., Shannon, R. M., et al. 2024, ApJL, 971, L18

Riley, T. E., Choudhury, D., Salmi, T., et al. 2023, JOSS, 8, 4977

Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21

Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27

Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51

Rutherford, N., Mendes, M., Svensson, I., et al. 2024, ApJL, 971, L19

Salmi, T., Choudhury, D., Kini, Y., et al. 2024, ApJ, 974, 294

Salmi, T., Deneva, J. S., Ray, P. S., et al. 2024, Data and Software for: 'A NICER View of PSR J1231-1411: A Complex Case', v1.0.0, Zenodo, doi:10.5281/zenodo.13358349

Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2022, ApJ, 941, 150

Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2023, ApJ, 956, 138

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

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

Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18

Takátsy, J., Kovács, P., Wolf, G., & Schaffner-Bielich, J. 2023, PhRvD, 108, 043002

Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27

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

Verbiest, J. P. W., Weisberg, J. M., Chael, A. A., Lee, K. J., & Lorimer, D. R. 2012, ApJ, 755, 39

Vergely, J. L., Lallement, R., & Cox, N. L. J. 2022, A&A, 664, A174

Vinciguerra, S., Salmi, T., Watts, A. L., et al. 2023, ApJ, 959, 55

Vinciguerra, S., Salmi, T., Watts, A. L., et al. 2024, ApJ, 961, 62

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

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

Wijngaarden, M. J. P., Ho, W. C. G., Chang, P., et al. 2019, MNRAS, 484, 974

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

Wolff, M. T., Guillot, S., Bogdanov, S., et al. 2021, ApJL, 918, L26

Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29

Zavlin, V. E., & Pavlov, G. G. 2002, in Proc 270 WE-Heraeus Seminar on Neutron Stars, Pulsars, and Supernova Remnants, ed. W. Becker, H. Lesch, & J. Trümper (Garching bei München: Max-Plank-Institut für extraterrestrische Physik), 263