# A NICER view of PSR J0030+0451: Implications for the dense matter equation of state

G. Raaijmakers, T.E. Riley, A.L. Watts, S.K. Greif, S.M. Morsink, K.

Hebeler, A. Schwenk, T. Hinderer, S. Nissanke, S. Guillot, et al.

► To cite this version:

G. Raaijmakers, T.E. Riley, A.L. Watts, S.K. Greif, S.M. Morsink, et al., A NICER view of PSR J0030+0451: Implications for the dense matter equation of state. The Astrophysical Journal Letters, 2019, 887 (1), pp.L22. 10.3847/2041-8213/ab451a. [hal-02440009]

HAL Id: hal-02440009

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

Submitted on 29 Aug 2025

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

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

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

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

# A NICER View of PSR J0030+0451: Implications for the Dense Matter Equation of State

G. Raaijmakers¹ ID, T. E. Riley¹ ID, A. L. Watts¹ ID, S. K. Greif²,³, S. M. Morsink⁴ ID, K. Hebeler²,³, A. Schwenk²,³,⁵,

T. Hinderer¹,⁶ ID, S. Nissanke¹,⁶ ID, S. Guillot⁷,⁸ ID, Z. Arzoumanian⁹, S. Bogdanov¹⁰ ID, D. Chakrabarty¹¹ ID, K. C. Gendreau⁹,

W. C. G. Ho¹²,¹³ ID, J. M. Lattimer¹⁴, R. M. Ludlam¹⁵,<fnref n="17" /> ID, and M. T. Wolff¹⁶ ID

¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090GE Amsterdam, The Netherlands; G.Raaijmakers@uva.nl

² Institut für Kernphysik, Technische Universität Darmstadt, D-64289 Darmstadt, Germany

³ ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, D-64291 Darmstadt, Germany

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

⁵ Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg, Germany

⁶ GRAPPA Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

⁷ IRAP, CNRS, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France

⁸ Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France

⁹ 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

¹¹ MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

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

¹³ Mathematical Sciences, Physics and Astronomy, and STAG Research Centre, University of Southampton, Southampton, SO17 1BJ, UK

¹⁴ Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA

¹⁵ Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA

¹⁶ Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375-5352, USA

Received 2019 August 7; revised 2019 September 6; accepted 2019 September 16; published 2019 December 12

## Abstract

Both the mass and radius of the millisecond pulsar PSR J0030+0451 have been inferred via pulse-profile modeling of X-ray data obtained by NASA's *Neutron Star Interior Composition Explorer (NICER)* mission. In this Letter we study the implications of the mass-radius inference reported for this source by Riley et al. for the dense matter equation of state (EoS), in the context of prior information from nuclear physics at low densities. Using a Bayesian framework we infer central densities and EoS properties for two choices of high-density extensions: a piecewise-polytropic model and a model based on assumptions of the speed of sound in dense matter. Around nuclear saturation density these extensions are matched to an EoS uncertainty band obtained from calculations based on chiral effective field theory interactions, which provide a realistic description of atomic nuclei as well as empirical nuclear matter properties within uncertainties. We further constrain EoS expectations with input from the current highest measured pulsar mass; together, these constraints offer a narrow Bayesian prior informed by theory as well as laboratory and astrophysical measurements. The NICER mass-radius likelihood function derived by Riley et al. using pulse-profile modeling is consistent with the highest-density region of this prior. The present relatively large uncertainties on mass and radius for PSR J0030+0451 offer, however, only a weak posterior information gain over the prior. We explore the sensitivity to the inferred geometry of the heated regions that give rise to the pulsed emission, and find a small increase in posterior gain for an alternative (but less preferred) model. Lastly, we investigate the hypothetical scenario of increasing the NICER exposure time for PSR J0030+0451.

*Unified Astronomy Thesaurus concepts: Neutron stars (1108); Nuclear astrophysics (1129); Neutron star cores (1107); Pulsars (1306); Bayesian statistics (1900); X-ray stars (1823); Millisecond pulsars (1062)*

## 1. Introduction

The cores of neutron stars (NSs) provide a unique environment for exploring matter at densities above nuclear saturation density ($\rho_s = 2.7 \times 10^{14} \text{ g cm}^{-3}$). Theoretical predictions in this regime are diverse, ranging from nucleonic matter under extreme neutron-rich conditions, to stable states of strange matter such as hyperons or deconfined quarks, color superconducting phases, and Bose-Einstein condensates (for recent reviews see Hebeler et al. 2015; Lattimer & Prakash 2016; Oertel et al. 2017; Baym et al. 2018). Our uncertainty about the nature of cold supranuclear-density matter is often encoded in the equation of state (EoS) through general parametric extensions to high densities with an associated prior distribution. Each EoS maps via the stellar structure equations to sequences of stable spacetime solutions given interior boundary conditions (see the review by

Paschalidis & Stergioulas 2017). Properties such as total (or gravitational) mass *M* and equatorial radius *R*<sub>eq</sub> of the NS surface feature strongly in the exterior spacetime solution.<fnref n="17" /> Observational phenomena that are sensitive to the structure of the exterior spacetime, such as the propagation of radiation from the stellar surface to a distant observer, can thus be used to probe the EoS and hence the microphysics of dense matter.

NASA's *Neutron Star Interior Composition Explorer* (NICER; Gendreau et al. 2016), a soft X-ray telescope installed on the *International Space Station* in 2017, was developed to estimate masses and radii of NSs using pulse-profile modeling of nearby rotation-powered millisecond pulsars (MSPs). The magnetic polar caps of MSPs, thought to be heated by (return) currents in the pulsar magnetosphere, produce thermal emission in the soft X-ray band (Harding & Muslimov 2002). As the

<footnote n="26">²⁶ See Appendix A.2.4 of Riley et al. (2019) for supplementary detail about the KL divergence.</footnote>

<footnote n="17">¹⁷ Both in terms of metric functions, and the spatial domain of those functions.</footnote>

<footnote n="17">¹⁷ Einstein Fellow.</footnote>
MSP rotates, this emission gives rise to perceived pulsations, and relativistic effects encode information about the spacetime into the phase-energy resolved pulse-profile.<fnref n="18" /> Pulse-profile modeling employs relativistic ray-tracing and Bayesian inference software to jointly infer mass and radius (see Bogdanov 2016; Watts et al. 2016; Watts 2019, for an overview of the technique).

Riley et al. (2019) jointly estimated the mass *M* and radius *R*<sub>eq</sub> of the MSP PSR J0030+0451 conditional on NICER X-ray Timing Instrument (XTI) photon event data curated by Bogdanov et al. (2019). The results derived are also conditional upon the modeling choices made in the analysis, e.g.: the assumption of two disjoint surface hot regions, each with some local comoving effective temperature field but no magnetic field physics; a fully ionized hydrogen atmosphere; and a specific parameterization of the uncertainty in the NICERXTI instrument response.<fnref n="19" /> The restriction to two disjoint hot regions was motivated by the presence of two distinct pulses in the observed pulse profile. Riley et al. (2019) allowed for the possibility of the hot regions being non-antipodal and non-identical, and considered various shapes for the hot regions including circles, rings (with the centers both concentric and offset), and crescents, filled with material of a single local comoving temperature. Model comparison enabled the identification of a favored configuration, using a combination of performance measures including the evidence (the prior predictive probability of the data) and graphical posterior predictive checking (to verify whether or not an updated model generates synthetic data<fnref n="20" /> without obvious residual systematic structure in comparison to the real data).

For the family of models considered, the favored configuration is one in which the hot regions consist of a small hot spot with angular extent of only a few degrees, and a more extended hot crescent, both in the same rotational hemisphere (referred to in Riley et al. 2019 as ST+PST). For this configuration, the inferred mass and equatorial radius<fnref n="21" /> are $M = 1.34^{+0.15}_{-0.16} M_{\odot}$ and $R_{\text{eq}} = 12.71^{+1.14}_{-1.19}$ km. The compactness $GM/R_{\text{eq}}c^2 = 0.156^{+0.008}_{-0.010}$ is more tightly constrained.

The credible bounds reported here are approximately the 16% and 84% quantiles in marginal posterior mass. The spin frequency of PSR J0030+0451 is only 205 Hz: *M* and *R*<sub>eq</sub> can, due to the size of the credible intervals, therefore be identified as those of a non-rotating star with an equivalent number of baryons. The effects of rotation are discussed in more detail in Section 4.

If the extended hot region is restricted to have ring-like topology rather than that of a simply connected crescent (ST+CST in Riley et al. 2019), the inferred mass and equatorial radius are $M = 1.44^{+0.18}_{-0.19} M_{\odot}$ and $R_{\text{eq}} = 13.89^{+1.22}_{-1.39}$ km. The compactness $GM/R_{\text{eq}}c^2 = 0.16 \pm 0.01$ is however shared with ST+PST—at the quoted precision. Although ST+CST was not the favored configuration a posteriori, it provides a useful illustration of the sensitivity of dense matter inferences to the nuisance parameters controlling the surface radiation field. As pulsar theory develops, dense matter inferences therefore need to be re-examined in step. Fortunately, such calculations are less expensive to execute given posterior samples because nuisance-parameter marginalization is thereby approximated.

In this Letter, we examine how the constraints on NS mass and radius translate into constraints on the dense matter EoS. Ultimately, we intend to carry out a population-level analysis conditional on all NICER MSP targets, in order to report a joint summary for NICER. We propose to inject as little information as is reasonable from other statistical constraints derived from astronomical data sets—the exception being information from a radio pulsar mass measurement. Eventually, we aim to combine the joint NICER constraints with those derived using other missions, where appropriate. However, for now we have information for a single source, and in this Letter we address how the joint mass–radius information derived by Riley et al. (2019) for PSR J0030+0451 maps to constraints on the dense matter EoS. The second principal aim of this Letter is therefore to formalize a plan for post-processing posterior information derived via pulse-profile modeling, into posterior information about dense matter. The post-processing phase for dense matter study is far less computationally expensive than the preceding X-ray analysis in which the likelihood information relevant for dense matter study is computed. We can therefore effectively update our posterior information on-the-fly as new information becomes available, by jointly compiling NICER source-by-source nuisance-marginalized likelihood functions into posterior constraints about a common EoS.

## 2. NICER EoS Constraints

### 2.1. EoS Parameterizations

Following the methods described in Greif et al. (2019), we model the interior of PSR J0030+0451 using two distinct EoS parameterizations: the piecewise-polytropic (PP) model from Hebeler et al. (2013) and a speed of sound (CS) model introduced in Greif et al. (2019)—see also Tews et al. (2018a). These parameterizations were matched at 1.1 $\rho_s$ to either the upper limit or the lower limit of a calculated EoS range based on chiral effective field theory (cEFT) interactions including theoretical uncertainties (for details, see Hebeler & Schwenk 2010; Hebeler et al. 2013). This discrete matching leads to a bimodality in the prior of the EoS, which we mitigate here by introducing an additional parameterization of the EoS inside the chiral EFT band. For simplicity we assume a single polytrope, i.e.,

$$P(\rho) = K \left( \frac{\rho}{\rho_s} \right)^{\Gamma}, \qquad (1)$$

where *P* is the pressure, $\rho$ is the baryon mass density, *K* (in units of MeV fm⁻³) is a free parameter, and $\Gamma$ is the adiabatic

<footnote n="18">¹⁸ A pulse-profile consists of X-ray counts per rotational phase bin per instrument detector channel, curated by phase-folding X-ray events according to a pulsar timing ephemeris.</footnote>

<footnote n="19">¹⁹ For an independent analysis of the same data set using different modeling choices and methodology see Miller et al. (2019b), which follows the approach outlined in (Miller et al. 2019a).</footnote>

<footnote n="20">²⁰ For illustration, a pulse-profile count-number data set can be simulated given specific instances of the following components: a spacetime solution; a surface hot-region configuration (effective temperature, geometry); an atmospheric beaming function (composition, ionization); background contribution (astrophysical, instrumental); an instrument response function; and a noise model (Poissonian). Source emission is propagated via relativistic ray-tracing through the spacetime toward a distant observer inclined to the stellar spin axis, and is subsequently operated on by the instrument response function; the product is a joint sampling distribution for photon count numbers, which is intrinsic to the definition of a likelihood function. The notion of synthetic data generation is a vital part of the Bayesian inference framework. See Riley et al. (2019) and Watts (2019) for more discussion.</footnote>

<footnote n="21">²¹ With respect to a Schwarzschild coordinate chart, see Section 2.3.1 of Riley et al. (2019) for more details.</footnote>
exponent. We determine Γ and the bounds of K by fitting Equation (1) to the lower and upper limits of the cEFT band for densities between 0.5ρs and 1.1 ρs and find that these limits are well approximated by Γ = 2.5, Kmin = 1.70, and Kmax = 2.76. At densities below 0.5 ρs we match to a single crust EoS (BPS; see?). Comparing the full range of masses and radii permitted by EoS under both parameterizations with a continuous matching to the the cEFT band against the upper/lower limit case, we find that they are consistent, although a slightly larger range is obtained for the PP model for the continuous case. This is due to the polytropic fit to the cEFT band allowing for a small set of additional EoSs that are soft enough in the low-density regime to result in small NS radii but stiff enough at larger densities to comply with the pulsar mass constraint (see Section 2.2.1 and Greif et al. 2019).

## 2.2. Bayesian Framework

To derive constraints on the EoS from a single-star mass-radius posterior density distribution, we use the Bayesian framework outlined in Greif et al. (2019) and Riley et al. (2018). Let us combine the EoS parameters and the central density εc of PSR J0030+0451 into a vector θ. The posterior distribution of θ is proportional to the product of the prior distribution of θ and the nuisance-marginalized likelihood function of θ (Bayes' theorem):

$$ p(\theta|\mathcal{d}, \mathcal{M}) \propto p(\theta|\mathcal{M}) p(\mathcal{d}|\theta, \mathcal{M}) \\ \propto p(\theta|\mathcal{M}) p(\mathcal{M}, R|\mathcal{d}, \mathcal{M}), \quad (2) $$

where d denotes the NICER PSR J0030+0451 data set, and M denotes the model. The model includes the physics of processes both interior and exterior to the star: the EoS and central conditions (present work); and X-ray emission, propagation, and detection (Riley et al. 2019, and references therein). Implicitly, the model includes all Bayesian prior information. The parameters θ map deterministically to the mass M = M(θ; Ω) and radius R = R(θ; Ω), where the coordinate angular rotation frequency Ω = 0 (see Section 4.3). All parameters apart from M and R are, for the purposes of the discussion that follows, termed nuisance parameters and are marginalized out.<fnref n="22" />

To obtain the second line of Equation (2) we equated the nuisance-marginalized likelihood function of M and R to the nuisance-marginalized joint posterior density distribution of M and R reported by Riley et al. (2019). This proportionality holds exactly because the marginal joint prior distribution of M and R chosen by Riley et al. (2019) is jointly flat. Our numerical nuisance-marginalized likelihood function is an approximation to the exact nuisance-marginalized likelihood function because we are post-processing posterior samples, and because post-processing involves kernel density estimation (KDE) of the posterior density function. We then sample from the posterior density p(θ|d, M) in Equation (2) using the nested sampling software MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2013; Buchner et al. 2014).

## 2.2.1. Priors

The prior density, p(θ|M), in Equation (2) is identical to the prior described in Section 3.1.1 of Greif et al. (2019), with the exception of the implementation of the continuous range within the cEFT band. Here we introduce an additional uniform prior on the parameter K in Equation (1) with support K ∈ [1.70, 2.76]. The prior for all parameters in the PP and CS model is summarized in Table 2. Moreover, we require that every EoS assigned a finite local prior density can support a stable 1.97 M⊙ NS, equal to the lower 1σ limit on the mass of the most massive NS measured to date (PSR J0348+0432; Antoniadis et al. 2013).<fnref n="23" />

The joint prior density p(θ|M) is written (for both parameterizations) as

$$ p(\theta|\mathcal{M}) = p(\varepsilon_c|\text{EoS}, \mathcal{M})p(\text{EoS}|\mathcal{M}), \quad (3) $$

where the conditional prior density of the central density εc, p(εc|EoS, M), is numerically evaluated to impose global spacetime stability. Given an EoS, we inverse sample the conditional density with rejection: we reject θ if θ ↦ (M, R) does not yield a stable spacetime solution that exists within the support of the posterior density p(M, R|d, M). Outside of this support, the nuisance-marginalized likelihood function has not been estimated in the preceding X-ray analysis. Riley et al. (2019) imposed prior support with hard bounds of M ∈ [1.0, 3.0] M⊙ and R ∈ [3r<sub>g</sub>, 16] km, where r<sub>g</sub> = GM/c². Note that, except for the lower bound on the mass, these bounds already have zero support from the prior on the EoS model. The lower bound on the mass is implemented in MULTINEST by assigning a likelihood value below a certain threshold to any mass-radius pair outside this bound. Any mass-radius pair with a likelihood lower than this threshold will then be ignored in the nested sampling process by MULTINEST (see T. E. Riley & A. L. Watts 2019, in preparation, for a discussion on prior density implementation options for use with MULTINEST). Besides the matching to the cEFT band, the PP model is constrained by causality and the requirement to support a 1.97 M⊙ NS. Thus, the PP model has prior support for the central density up to the maximum mass (which is required to be reached before the speed of light c<sub>s</sub> = c). For the CS model we also consider the EoS for central densities up to the maximum mass and impose the constraints that

(i) the speed of sound for all energy densities is less than the speed of light;

(ii) the speed of sound of each EoS converges to (c<sub>s</sub>/c)² = 1/3 from below at ~50ρ<sub>s</sub>, following the calculations of the speed of sound for asymptotically high densities by perturbative quantum chromodynamics (pQCD; Fraga et al. 2014);

(iii) the bulk properties of matter at densities ρ ≤ 1.5ρ<sub>s</sub> can be described as a normal Fermi liquid, which restricts the speed of sound at these densities to be (c<sub>s</sub>/c)² ≤ 0.163 (see Greif et al. 2019, and references therein, for more detail).

While both our prior for the CS model and the PP model allow for phase transitions at certain density ranges, they do not cover

<footnote n="22">²² Note that they do, however, describe important and interesting physics on the surface of the star and exterior to it, to which our pulse-profile modeling is extremely sensitive. See for example Bilous et al. (2019) on the implications of some of these inferred “nuisance parameters” for our understanding of pulsar magnetospheres.</footnote>

<footnote n="23">²³ Note that this is an ad hoc interpretation of the information encoded in the pulsar mass measurement. We discuss this in more detail in Section 4.2.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Prior probability distributions for PSR J0030+0451 transformed to the joint space of central pressure and central energy density (upper panels), and the space of mass and radius (lower panels), for the PP (left panels) and CS parameterization (right panels). The dotted and dashed contours bound the highest-density two-dimensional regions respectively containing 68% and 95% of the prior mass. The peak in the CS model just above 1.97 $M_{\odot}$ is due to how the model is constructed: all EoSs are forced to soften at high densities to comply with causality and at asymptotic densities with the constraint from pQCD, causing all mass-radius sequences to have $\partial M/\partial R \lesssim 0$ near the maximum mass, thereby overlapping each other (see Greif et al. 2019). Note that these priors are specifically for PSR J0030+0451 because adjustments are made to match the priors in the analysis in Riley et al. (2019); i.e., $M \in [1.0, 3.0] M_{\odot}$ and $R \in [3r_g, 16]$ km where $r_g(M) = GM/c^2$. For comparison the three representative EoSs from Hebeler et al. (2013) are shown as solid curves: HLPS Soft, Intermediate, and Stiff. Note that the discernible small-scale structure is due to: (i) the behavior of the (numerical) transformation from interior matter parameters to exterior spacetime parameters; and (ii) finite sampling noise.</figcaption></figure>

all possibilities of transitions to other forms of matter. However, EoSs that mimic hybrid stars for which the transition is smooth do exist within our prior bounds. We discard any EoS model that allows for two disconnected stable branches on the mass-radius sequence (similar to Alford & Han 2016), which occurred in our sampling for certain parameter sets of the PP model but is also possible in the CS model.

In order to understand our prior choices and the effect of the continuous matching to the cEFT band in more detail, we randomly sample ~$10^5$ points from the prior distributions and calculate for each point: (i) a central energy density and central pressure pair; and (ii) a mass-radius pair. The resulting distributions for both the PP and CS model are shown in Figure 1. We note that the bimodality observed in the prior in Greif et al. (2019) has been smoothed out by the addition of the continuous matching. The darker region at low density, or equivalently at low masses, comes as a result of all EoSs being matched to the cEFT band. However, in the analysis of PSR J0030+0451 this region is outside of the prior support ($M \ge 1.0 M_{\odot}$) from the analysis in Riley et al. (2019). Less intuitive is the darker region for the prior of the CS model just

above 1.97 $M_{\odot}$, which can be explained by investigating the individual mass-radius sequences that contribute to this clustering: the CS model is constructed to be causal at all densities and is only truncated when $dM/d\epsilon_c \le 0$. This causes the mass-radius sequences to bend over on the $M-R$ plot after the 1.97 $M_{\odot}$ constraint is fulfilled and extend horizontally toward smaller radii, overlapping with each other. The PP model differs in this feature because individual EoSs are truncated when the EoS reaches $c_s = c$, so that such an EoS does not need to bend over. This allows mass-radius sequences with steep slopes to exist up to the density where $c_s = c$.

## 2.3. EoS Constraints Based on PSR J0030+0451

We consider two distinct mass-radius posterior distributions supplied by Riley et al. (2019), each conditional on assumptions about the thermally emitting hot regions on the surface. These assumption sets are identified as the ST+PST and ST+CST models, which yielded 68% credible intervals on the equatorial radius of $R = 12.71_{-1.19}^{+1.14}$ km and $R = 13.89_{-1.38}^{+1.23}$ km, respectively. We stress that the favored configuration identified
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Marginal posterior distributions of the pressure *P* conditional on energy density *ε*, *p*(*P*|*ε*), *d*, *M*), for the PP model (left) and the CS model (right), and for both the preferred ST+PST model (upper panels) and the alternative ST+CST model (lower panels). At each value of *ε*, there exist 68% and 95% posterior credible intervals for the pressure *P*; we connect these intervals to form the shaded bands. The black dotted and dashed lines, respectively, indicate the joined 68% and 95% credible interval bands, but for the conditional and marginal prior distribution, *p*(*P*|*ε*, *M*). The red contours in each panel indicate the 68% and 95% highest-density posterior credible regions of central energy density and central pressure. Constraints on the EoS for densities higher than these contours are only determined by our choice of parameterization and are not directly informed by the mass–radius likelihood function (and thus in turn, the data). The lower-right inset panels illustrate the evolution of the Kullback–Leibler (KL) divergence with respect to the energy density, showing that most posterior information is gained for densities below 10¹⁵ g cm⁻³, not coincidentally the highest possible central density reached in PSR J0030+0451. Note that due to finite sampling noise the precise features of the evolution of the KL divergence might be disputed, but the global trend of the curve is unaffected.</figcaption></figure>

by Riley et al. (2019) is ST+PST; we explore the other configuration only to illustrate the potential sensitivity of our conclusions to developments in pulsar theory or additional observations.

We show the 68% and 95% credible regions of the resulting posterior distributions transformed to the EoS space in Figure 2. Sensitivity to compactness manifests strongly as a constraint on central conditions—the density and pressure. Note that the marginal posterior credible interval on the central density of PSR J0030+0451 is dependent (or conditional) on the assumed EoS model. The inset panels in Figure 2 show that most information gain about matter pressure, measured through the Kullback–Leibler (KL) divergence (Kullback & Leibler 1951, see below), is in the vicinity of densities found in the core, but in absolute terms the gain is negligible. The EoS at higher densities than the central density (see also Table 1) is not directly informed by the data because a posteriori, matter at such densities does not exist in the star; any information gain is due to dependence on our choice of EoS parameterization which couples low- and high-density regimes via simple functional forms. Also rendered is a comparison of the

<table><caption>**Table 1**
Central Energy Density and Central Pressure of PSR J0030+0451</caption><thead><tr><th>Parameterization</th><th>Hot Region Model</th><th>log<sub>10</sub>(ε<sub>c</sub>)</th><th>log<sub>10</sub>(P<sub>c</sub>)</th></tr></thead><tbody><tr><td rowspan="2">PP model</td><td>ST+PST</td><td>14.80<sup>+0.04</sup><sub>-0.07</sub></td><td>34.87<sup>+0.06</sup><sub>-0.08</sub></td></tr><tr><td>ST+CST</td><td>14.78<sup>+0.04</sup><sub>-0.06</sub></td><td>34.84<sup>+0.06</sup><sub>-0.09</sub></td></tr><tr><td rowspan="2">CS model</td><td>ST+PST</td><td>14.88<sup>+0.05</sup><sub>-0.06</sub></td><td>34.94<sup>+0.07</sup><sub>-0.09</sub></td></tr><tr><td>ST+CST</td><td>14.86<sup>+0.05</sup><sub>-0.06</sub></td><td>34.90<sup>+0.07</sup><sub>-0.1</sub></td></tr></tbody></table>

posterior distributions with the prior distributions (compare the green shaded bands and the black dashed and dotted bands), which suggests no remarkable reduction in the degree of uncertainty—the prior distribution of the EoS functions is dominant (just as in Greif et al. 2019). In addition, the sensitivity of the posterior distributions to the chosen EoS parameterization (PP or CS) is stronger than to the model chosen for the hot regions (ST+PST or ST+CST). This can also be concluded from the inferred values of central energy densities and corresponding pressures in Table 1; larger differences occur between chosen parameterizations than
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 68% and 95% highest-density credible regions of the PSR J0030+0451 likelihood function (normalized by a flat joint density) are bounded by the blue contours: the ST+PST model is featured in the upper panels, and the ST+CST model is featured in the lower panels. The black contours are identical to those shown in Figure 1, being associated with a population-level prior transformed to mass–radius space; here we show how this prior is updated as a result of our analysis of PSR J0030+0451. To provide context, we are considering the implications for a hypothetical future analysis of data from a different pulsar *S* (i.e., not PSR J0030+0451). We assume that the EoS of all matter (core to crust) is shared between *S* and PSR J0030+0451; as a joint prior for the EoS of matter in *S*, we thus invoke the joint posterior distribution of EoS parameters conditional on NICER observations of PSR J0030+0451. We assume the central density of *S* is drawn from the same population-level density *p*(*ε*<sub>c</sub>|EoS) as PSR J0030+0451. We then transform the joint prior of the EoS parameters and central density of *S* to the joint space of the mass and radius of *S*; we render the two-dimensional regions enclosing 68% and 95% of the updated prior mass in green. We note that the updated prior distributions are still mostly dominated by the prior for PSR J0300+0451, with slightly more support for higher radii in the ST+CST model.</figcaption></figure>

models for hot regions. As a consequence it is difficult to
investigate the sensitivity of our inference to the assumed hot
region model.

In Figure 3 we show the nuisance-marginalized likelihood
function of the mass and radius of PSR J0030+0451. We also
show how the analysis of the dense matter in this source
modifies our population-level prior, when transformed to
mass–radius space (an alternative representation of the poster-
ior information on the EoS). Consider the hypothetical future
analysis of some other observed NS, *S*, which shares an EoS—
from core to crust—with PSR J0030+0451, and whose central
density ε<sub>c</sub> is drawn from the same population-level density
p(ε<sub>c</sub>|EoS) as that of PSR J0030+0451.<sup>24</sup> The joint prior
for analysis of *S* is then based on the PSR J0030+0451 posterior
(see the caption for details). The figure clearly illustrates that
*after* learning about the EoS from PSR J0030+0451, the prior
distribution of mass and radius of *S* remains dominated by the

original prior information invoked for analysis of PSR J0030
+0451. Note that considering the ST+CST model for PSR J0030
+0451 shifts the prior mass–radius sequences toward slightly
higher radii.<sup>25</sup>

The posterior distributions on the speed of sound for the CS
model (not shown in this Letter) are similar to the distributions
shown in Figure 8 of Greif et al. (2019), again showing no
evidence of the speed of sound reaching the asymptotic limit
(c_s/c)² = 1/3 within the range of energy densities relevant
for NSs.

To quantify the information gain (in bits) of the posterior
over the prior, we compute for each model the KL divergence,
an asymmetric measure of how different one probability
distribution is from another (Kullback & Leibler 1951).<sup>26</sup> The
errors on the divergences are obtained by repeated calculation
for each in a set of posterior realizations; each realization is

<footnote n="24"><sup>24</sup> Note that this a Bayesian hierarchical model, where the shared EoS parameters appear in the likelihood function but also effectively appear as hyperparameters of the central density prior distribution.</footnote>

<footnote n="25"><sup>25</sup> To compute the 68% and 95% highest-density credible regions in Figure 3 we have applied a Gaussian KDE with the bandwidth parameter according to Scott's rule, and determined that the results and conclusion presented here are not affected by this parameter.</footnote>
<table><caption>**Table 2**
Summary Table of the Prior, KL Divergences, and Evidences for all Inferred Posterior Distributions</caption><thead><tr><th>Parameter</th><th>Prior Density</th><th>ST+PST (T)</th><th>ST+CST (T)</th><th>ST+PST (2T)</th><th>ST+CST (2T)</th></tr><tr><td></td><td></td><td>D̂<sub>KL</sub> (10<sup>-3</sup> bits)</td><td>D̂<sub>KL</sub> (10<sup>-3</sup> bits)</td><td>D̂<sub>KL</sub> (10<sup>-3</sup> bits)</td><td>D̂<sub>KL</sub> (10<sup>-3</sup> bits)</td></tr></thead><tbody><tr><td rowspan="6">PP model</td><td>K</td><td>U (1.7, 2.76)</td><td>2.47 ± 0.27</td><td>26.24 ± 1.13</td><td>3.72 ± 0.37</td><td>74.80 ± 7.94</td></tr><tr><td>Γ<sub>1</sub></td><td>U (1, 4.5)</td><td>16.32 ± 1.00</td><td>72.77 ± 1.89</td><td>36.83 ± 1.31</td><td>129.98 ± 11.62</td></tr><tr><td>Γ<sub>2</sub></td><td>U (0, 8)</td><td>3.09 ± 0.40</td><td>7.57 ± 0.61</td><td>5.27 ± 0.52</td><td>6.72 ± 1.83</td></tr><tr><td>Γ<sub>3</sub></td><td>U (0.5, 8)</td><td>6.64 ± 0.24</td><td>8.02 ± 0.37</td><td>7.41 ± 0.31</td><td>20.51 ± 1.56</td></tr><tr><td>ρ<sub>12</sub></td><td>U (1.5, 8.3) n<sub>0</sub></td><td>8.13 ± 0.57</td><td>8.71 ± 0.61</td><td>8.45 ± 0.64</td><td>8.38 ± 1.64</td></tr><tr><td>ρ<sub>23</sub></td><td>U (ρ<sub>12</sub>, 8.3) n<sub>0</sub></td><td>3.53 ± 0.43</td><td>6.08 ± 0.45</td><td>3.86 ± 0.42</td><td>7.64 ± 1.82</td></tr><tr><td></td><td>log(ε<sub>c</sub>)</td><td>U (14.6, max(ε<sub>c</sub>))</td><td>982.44 ± 4.20</td><td>1074.48 ± 4.60</td><td>1089.29 ± 4.20</td><td>1213.18 ± 18.97</td></tr><tr><td colspan="2">Log-evidence ln Ẑ</td><td></td><td>-1.43 ± 0.02</td><td>-1.96 ± 0.02</td><td>-0.92 ± 0.02</td><td>-1.70 ± 0.02</td></tr><tr><td rowspan="6">CS model</td><td>K</td><td>U (1.7, 2.76)</td><td>7.37 ± 0.67</td><td>16.40 ± 0.92</td><td>11.32 ± 0.71</td><td>43.83 ± 1.16</td></tr><tr><td>a<sub>1</sub></td><td>U (0.1, 1.5)</td><td>7.74 ± 0.64</td><td>6.80 ± 0.52</td><td>6.92 ± 0.51</td><td>7.63 ± 0.54</td></tr><tr><td>a<sub>2</sub></td><td>U (1.5, 12.0)</td><td>5.94 ± 0.64</td><td>11.5 ± 0.88</td><td>10.18 ± 0.69</td><td>46.53 ± 1.23</td></tr><tr><td>a<sub>3</sub>/a<sub>2</sub></td><td>U (0.05, 2.0)</td><td>42.48 ± 1.63</td><td>40.77 ± 1.28</td><td>43.54 ± 1.30</td><td>52.06 ± 1.46</td></tr><tr><td>a<sub>4</sub></td><td>U (1.5, 37.0)</td><td>20.10 ± 1.13</td><td>16.48 ± 0.92</td><td>17.43 ± 0.98</td><td>23.35 ± 1.00</td></tr><tr><td>a<sub>5</sub></td><td>U (0.1, 1.0)</td><td>2.02 ± 0.31</td><td>2.16 ± 0.35</td><td>1.12 ± 0.23</td><td>1.09 ± 0.16</td></tr><tr><td></td><td>log(ε<sub>c</sub>)</td><td>U (14.6, max(ε<sub>c</sub>))</td><td>1331.8 ± 3.77</td><td>1376.7 ± 3.34</td><td>1440.8 ± 3.5</td><td>1508.7 ± 3.12</td></tr><tr><td colspan="2">Log-evidence ln Ẑ</td><td></td><td>-1.83 ± 0.02</td><td>-2.65 ± 0.02</td><td>-1.11 ± 0.02</td><td>-2.46 ± 0.02</td></tr><tr><td colspan="2">Bayes' factor</td><td></td><td>1.49</td><td>1.99</td><td>1.21</td><td>2.14</td></tr></tbody></table>

**Note.** The parameters in the PP and CS model with their corresponding priors, where $U(a, b)$ denotes uniformly sampled between $a$ and $b$. The upper prior bound on the central energy densities is a variable that depends on the other EoS parameters, see Section 2.2.1. Also shown are the estimated log-evidences $\ln Z$ for the four models considered in this Letter, and parameter-by-parameter KL divergences $\hat{D}_{KL}$ (mean values and standard deviation), for each individual parameter. We report the parameter information gain for: (i) likelihood functions $L(M, R)$ associated with the Bogdanov et al. (2019) $T = 1.94$ Ms NICER data set; and (ii) crudely simulated likelihood functions $L^*(M, R)$ for an exposure time of $2T$. Errors are estimated by calculating these quantities for a set of equally weighted realizations of the nested sampling process, where each realization has an associated posterior distribution. The low absolute value of the KL divergence in bits is indicative that not much information is gained over the prior in each model considered here. The Bayes' factors show that for the ST+CST model there is substantially more support for the PP model, caused by the stricter prior constraints on higher radii for the CS model. One should be careful in comparing evidences in combinations other than the reported Bayes' factors: the nuisance-marginalized mass-radius likelihood functions need to be normalized appropriately (by another evidence), and therefore only if a likelihood function is shared between models does this normalization cancel exactly. Evidence ratios between models with different surface hot regions depend on evidences estimated by Riley et al. (2019); the error intervals for the ST+PST and ST+CST models overlap substantially, however, and can be safely equated. The normalization for a likelihood function $L^*(M, R)$ defined by simulating extension of observing time is unknown.

simply simulated by bootstrap resampling with replacement from the samples according to their importance weights. Note that this is intended as a fast and simple approximation to the nested sampling error theory treated generally in the literature (namely, Skilling 2006; Higson et al. 2018, 2019; Higson 2018), in order to get a handle on the magnitude of the noise.

Inspecting individual KL divergences for each parameter (see Table 2) reveals that most information is gained in the distribution of central densities, while other parameters have KL divergences closer to zero. This can be visualized by comparing the posterior on central energy densities and pressures in Figure 2 with the prior distribution in Figure 1: the most information is gained along direction of the central energy density.

We further quantify the posterior distributions by performing a model comparison between the PP and CS model using Bayes' factors. The Bayes' factor is the ratio between the model evidences;<fnref n="27" /> if we were to accept a uniform prior mass distribution over the discrete models, the posterior odds ratios are equal to the Bayes' factors. In Table 2 we report the evidences as well as the Bayes' factors, computed here as the ratio of the PP model evidence to the CS model evidence for

some exterior-physics (PST or CST) likelihood function. Following the interpretation of Bayes' factors of Kass & Raftery (1995), we observe that there is no preference for either of the two parameterizations when using both the ST+PST and the ST+CST model. The Bayes' factors show however slightly increased support for the PP model when using the ST+CST model, a consequence of the tighter constraint on large radii for the CS model. As a result there is a tension between the inferred radius $R = 13.89^{+1.23}_{-1.38}$ km for the ST+CST model and the allowed range of radii in the CS model, which is less evident in the PP model.

### 3. Effect of Increased Exposure Time

The dense matter information yield, conditional on NICER observations of PSR J0030+0451, is weak in the context of existing knowledge, both theoretical and observational. An important aspect of any telescope mission is resource management: we therefore now investigate a scenario wherein the integrated observing time is increased. The analysis presented by Riley et al. (2019) used a data set with an integrated exposure time of 1.94 Ms, curated by Bogdanov et al. (2019).

Previous studies that have examined how posterior estimation of mass and radius is sensitive to factors including the number of source counts in the event data indicate that

<footnote n="27">²⁷ The model evidence being the principal computational target of the MULTINEST algorithm.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Similar to Figure 3 but for the hypothetical scenario where the observing time of PSR J0030+0451 is increased. We crudely simulate the evolution in the nuisance-marginalized likelihood function by contracting the *p*(M, R|d) credible regions by a factor of two, while retaining the coordinates of the point which reports the highest nuisance-marginalized likelihood value.</figcaption></figure>

constraining power increases as the square root of the number of counts (Lo et al. 2013; Psaltis et al. 2014) and thus observing time *T*. It is unclear, however, whether the credible region areas should be expected to scale as ~1/√*T* or ~1/*T*. The mass and radius are correlated in Lo et al. (2013) and Miller & Lamb (2015), but the areal reduction is given for approximate uncertainty in *M* and *R* separately. Let us simply suppose that the credible regions contract along some dimension *M*/*R* ≈ constant by a factor of ~1/√*T*, and also enjoy a ~1/√*T* scaling along the local (orthogonal) compactness direction; the overall scaling of area is then ~1/*T*. Considering this scaling as optimistic, and a scaling of ~1/√*T* as conservative, we would require ~2–4 times the exposure on PSR J0030+0451 in order to halve the credible region area.

Let us make the assumption that the credible regions for mass and radius *halve* in area. In Figure 4 we show these speculative posterior distributions conditional on the ST+PST and ST+CST models, corresponding to some unknown extension to the observing time. We contract the credible regions, simply assuming that the posterior distributional mean vector is insensitive to continued observation (refer to the Appendix). The effect is to artificially increase the absolute curvature of the mass–radius nuisance-marginalized likelihood function.

For the ST+PST model we note that the distributions remain very similar to the distributions in the upper panels of Figure 3: the 68% highest-density credible region of the updated mass–radius posterior<fnref n="28" /> spans almost the entire region containing 95% of the prior mass. For the ST+CST model the distributions do show slightly more support for higher radii, while the range of the 68% and 95% credible intervals decreases. Examination of the one-dimensional KL divergences for individual parameters exhibit, in most cases, a small increase. In some cases, however, the KL divergence for individual parameters decreases, but only when the KL divergence is close to zero where the error intervals due to sampling noise substantially overlap. The Bayes’ factors indicate that there is still only substantial posterior support for the PP model when considering the ST+CST model.

# 4. Discussion and Conclusions

## 4.1. Prospects for PSR J0030+0451 as an EoS Probe

Following the reported mass–radius posterior distribution from data obtained with NICER on PSR J0030+0451 by (Riley et al. 2019) we have explored the implications for the dense

<footnote n="28">²⁸ Normalized, nuisance-marginalized likelihood function.</footnote>
matter EoS. Two distinct hot-region models were considered, the superior ST+PST and also (for illustrative purposes) ST +CST, both yielding a different constraint on the NS radius. We have inferred the EoS using two high-density parameterizations, the PP model and CS model, which are matched to the cEFT band just above nuclear saturation density. The posterior distributions and corresponding KL divergences have shown that not much information is gained over the, already narrow, prior for either hot-region model. This can be attributed to the relatively large uncertainty in both mass–radius likelihood functions compared to the highest-density region of our prior, and the substantial overlap between the two. From the distributions shown in Figure 2 we observe that the changes from the prior to the posterior are also insignificant at nuclear densities, so that from the present analysis we cannot draw conclusions about further constraints on dense matter interactions within cEFT.

The Bayes' factors indicate that for this particular source neither parameterization is preferred, although when using the ST+CST model there is slightly more support for the PP model, a result of the CS model not being able to produce large enough radii to be somewhat consistent with the likelihood function. When more information on NS radii will become available in the future we would hopefully be able to discard certain parameterizations or construct new functional EoS forms based on the information encoded in the Bayes' factors.

In all these cases the posterior distributions are mostly informed by our prior choices, the two most stringent constraints being: (i) the lower limit of each EoS supporting a 1.97 $M_\odot$ NS, which is a simple way to include the information encoded in the pulsar radio-timing (and companion modeling) mass measurement by Antoniadis et al. (2013); and (ii) the computed cEFT band of possible EoS around nuclear saturation density. The first constraint causes a lower limit on the NS radius, while the latter has the opposite effect, resulting in a strongly peaked prior between ~11–13 km. With the tentative measurement of a $2.17_{-0.10}^{+0.11} M_\odot$ pulsar (PSR J0740+6620; Cromartie et al. 2019) the prior might even be updated to further exclude EoSs with small radii. We discuss the implementation of these priors in detail in Section 4.2.

As discussed in Riley et al. (2019) there are good prospects for advancing our understanding of the NICER background (from sources other than the MSP), and therefore tightening the PSR J0030+0451 mass–radius constraint via re-analysis. Any such re-analysis may also involve more sophisticated modeling of the MSP (and its near vicinity). It is not possible—at least at present—to robustly forecast how the mass–radius nuisance-marginalized function would change in response to such modeling efforts.

Substantial extension of the PSR J0030+0451 observing time is feasible given that the NICER mission has recently been extended for three more years. In Section 3 we crudely simulated the evolution of the nuisance-marginalized likelihood function with observing time. Several remarks must clearly be made in regards to this: (i) we neglected the notion of re-analysis of the NICER data curated by Bogdanov et al. (2019); and (ii) the studies by Lo et al. (2013) and Psaltis et al. (2014) on credible interval scaling assumed a single circular or infinitesimal single-temperature hot spot, not a more complex hot-region geometry such as that which emerged in Riley et al. (2019). While constraints should certainly improve with increased exposure time, we cannot confirm the precise

scaling for these more complex hot-region geometries without further study, and re-analysis may yield more constraining power.

Our modification of the nuisance-marginalized likelihood function—which may require substantial future observing time to formally realize—does not promise to enhance the dense matter information yield in the context of the radio-timing pulsar mass likelihood function and of the calculations of nuclear interactions based on cEFT around the nuclear saturation density. Ultimately, in order to improve synergy with the radio-timing probe of dense matter, an independent and tight mass constraint for PSR J0030+0451 would need to be combined with our compactness constraint conditional on NICER observations. However, the pulsar is not in a binary and thus there is no known prospect of an independent tight mass constraint—based on radio observations or otherwise.

In conclusion, we cannot yet robustly justify further allocation of observing time to PSR J0030+0451 for the purpose of dense matter study—at least without first exhausting modeling avenues for the 1.94 Ms Bogdanov et al. (2019) data set, and without similar modeling of the other primary NICER targets, principally PSR J0437−4715, which has a tight mass constraint derived via radio timing. Such a statement would arguably be valid even if the projected information gain was deemed substantial because it would be based on a crude forecast of the response to extended observing time, and the current state of knowledge may yet evolve via remodeling.

## 4.2. Prior Robustness and Implementation Approximations

cEFT allows for a systematic expansion of nuclear forces between neutrons and protons at low energies in terms of long-range pion-exchange contributions and short-range interactions. Within cEFT it is possible to determine contributions to nucleon–nucleon and many-body forces at different orders in the low-energy expansion and to provide estimates of theoretical uncertainties due to neglected higher-order terms. While theoretical predictions for systems with a significant proton fraction, such as atomic nuclei and isospin-symmetric matter, generally depend more sensitively on properties of the employed interactions, the scheme dependence of results for pure neutron systems or very neutron-rich systems exhibit a remarkable insensitivity to such details; see Hebeler et al. (2015) for a review, and also Lynn et al. (2016) and Drischler et al. (2019). In addition, the theoretical results for the nuclear symmetry energy are in good agreement with experimental constraints (Hebeler et al. 2013). These findings suggest that the predictions for the EoS of neutron-rich matter up to about nuclear saturation density are robust and rather well constrained. Current efforts aim at determining more systematically the upper density limit for such calculations, which is suspected to be closely related to the breakdown scale of cEFT. These studies may allow us to extend the calculations in a reliable way to higher densities and by this reduce also the EoS range of the PP and CS extensions at high densities. We note that in the analyses in this Letter we have approximated the EoS around nuclear densities with a polytropic fit to the cEFT band calculated by Hebeler et al. (2013). The true EoS within the cEFT band is however unknown, causing additional uncertainty in the presented posteriors that remains to be quantified.

One other assumption that merits comment is the prior requirement that each EoS support a 1.97 $M_\odot$ NS (Section 2.2).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left panel: effect of spin on the mass vs. equatorial radius curves for six representative EoSs allowing for a wide range of stiffness. Two curves with dark and light lines are shown for each EoS. For each pair, the lighter curve with smaller radii corresponds to the zero-spin mass–radius curve, while the darker curve with larger radii is the mass–radius curve for stars spinning at 205 Hz. The EoS in order of increasing stiffness (i.e., in order of increasing radius for a 1.4 $M_\odot$ star) are HLPS Soft (Hebeler et al. 2013), one of the softest EoSs allowed by nuclear physics; BBB2 (Baldo et al. 1997) is a soft EoS just marginally ruled out by the observation of a 1.97 $M_\odot$ pulsar; APR (Akmal et al. 1998) includes boost corrections; HLPS Intermediate and Stiff are representative EoSs from Hebeler et al. (2013); and L (Pandharipande & Smith 1975), a very stiff EoS, is most likely ruled out by the Laser Interferometer Gravitational-Wave Observatory (LIGO) observation of GW170817 (Abbott et al. 2017). The HLPS EoSs are also shown in Figures 1 and 2. Right panel: effect of spin on the equatorial compactness ratio $M/R$ vs. central energy density. The order of curves from left to right at a value of $M/R_{eq} = 0.15$ is from stiffest to softest EoS. Two curves, corresponding to 0 and 205 Hz are plotted for each EoS, however the difference between the curves is smaller than the line width so it is difficult to see the difference by eye. The gray horizontal box shows the compactness range of $M/R_{eq} = 0.156_{-0.010}^{+0.008}$ for the ST+PST model reported in Riley et al. (2019).</figcaption></figure>

This hard binary restriction discards the information encoded in the radio pulsar mass measurement: the information in the *shape* of the likelihood function is not included; and information is lost by assigning zero (rather than the correct) posterior density to EoS that do not support a 1.97 $M_\odot$ star; this is not compatible with the notion of future Bayesian updating. A more accurate approach would use the pulsar mass measurement as a nuisance-marginalized likelihood function and would not truncate until far into the tails (Raaijmakers et al. 2018).<fnref n="29" />

We considered the conclusions of this present work to be sufficiently insensitive to the above likelihood implementation detail to warrant discussion instead of recalculation. The approach taken here has the advantage that it enabled a fast, approximate separation of posterior information gain based on PSR J0030+0451 from the mass information gain based on PSR J0348+0432; it also avoided an additional parameter (a central density) in the sampling problem, leading to a minor reduction in computational expense. In future applications, however, we advocate for the default interpretation of (radio) pulsar mass measurements as nuisance-marginalized likelihood functions.

### 4.3. Effect of Rotation on the Accuracy of Likelihood Evaluation

The effect of rotation on an NS is to deform it into an oblate spheroid with a larger equatorial radius and to increase the mass compared to a non-rotating NS. These effects are correctly included in the mass and radius inference reported by Riley et al. (2019). The joint mass–radius inference made

use of an empirical quasi-universal relation for the oblate shape of the surface; the surface was then embedded in the Schwarzschild metric, and higher-order metric and shape corrections were neglected as is acceptable for stars spinning near 200 Hz (AlGendy & Morsink 2014). In particular, the inferred radius in the work of Riley et al. (2019) is the rotating star’s equatorial radius. However, in this Letter we assume that the NS is spherical, with a radius equal to the equatorial radius of the rotating star. It is worthwhile to consider the effect of this simplifying assumption on the accuracy of these results by computing sequences of rotating axisymmetric stars in hydrostatic equilibrium using the RNS code (Stergioulas & Friedman 1995).

The pulsar PSR J0030+0451 rotates with a spin frequency of 205 Hz. Figure 5 demonstrates the effect of rotation at this spin frequency on some representative EoS. In the left panel of Figure 5 we plot mass versus equatorial radius curves for stars spinning at 0 and 205 Hz for a representative set of tabulated EoS that cover a wide range of stiffness (see the figure caption). This shows that the increase in radius at fixed mass is at most 0.2 km for the stiffest EoS. More physically plausible EoS, such as the representative EoS of Hebeler et al. (2013), are deformed by a smaller coordinate distance. The difference of 0.2 km is much smaller than the precision, ±1.2 km, of the marginal radius estimate in Riley et al. (2019). The EoS constructed in this Letter using the PP and CS models will respond to rotation in a similar fashion. This suggests that this analysis is insensitive to the effects of rotation on the properties of this pulsar.

If two stars with the same central energy density and EoS but different spin rates are compared, the mass of the faster rotating star is larger (e.g., Hartle 1967; Hartle & Thorne 1968). Given that the tightest constraints reported by (Riley et al. 2019) are on the equatorial compactness ratio of $M/R$, it is more constructive to compute the changes in the equatorial compactness ratio, instead of the mass or radius, as a star spins. For zero-spin stars, given an EoS, each possible value of

<footnote n="29">²⁹ To clarify with an alternative perspective, suppose that one is simultaneously combining independent likelihood functions (e.g., from different stars) in a population-level analysis. The information encoded in a subset of the other likelihood functions might strongly weight EoSs that would be outside of the prior support according to a binary mass-threshold condition. Consequently, independent likelihood information would be censored.</footnote>
central energy density is mapped by the equations of hydrostatic equilibrium to a unique compactness ratio. The right panel of Figure 5 shows the curves of compactness versus central energy density for the same representative set of six different EoSs.<fnref n="30" /> Note that Figure 5 actually shows 12 curves, corresponding to two different spin rates of 0 and 205 Hz for each EoS. The curves for 205 Hz differ from the zero-spin curves by an amount that is smaller than 0.1%, which is smaller than the thickness of the line. This suggests that if the compactness for a particular central energy density and EoS is computed, it does not matter whether or not spin is included in the calculation. This property appears to extend to higher values of spin and will be investigated in more detail elsewhere.

### 4.4. Consistency with Previous EoS Constraints

There have been several attempts to constrain the parameters of dense matter EoS models using joint posterior information about mass and radius, where that information was derived via phase-averaged X-ray spectral modeling of bursting and quiescent accreting NSs (see Section 4 of Riley et al. 2019). Bogdanov et al. (2016) and Özel et al. (2016) used a piecewise-polytropic EoS model, inferring both EoS parameters and the associated mass–radius bands. Steiner et al. (2010, 2013, 2018) and Lattimer & Steiner (2014) considered a larger range of parameterized EoS models to infer dense matter parameters including those associated with the nuclear symmetry energy and the mass–radius bands. Baillot d’Etivaux et al. (2019) have performed a similar analysis using the parametrized EoS model of Margueron et al. (2018a, 2018b), inferring symmetry energy parameters, speed of sound profiles, and mass–radius bands.

A direct comparison between the EoS constraints derived in these works and ours is difficult: they use different models and priors, and the influence of the priors is not always clear from the published analysis. Due to the form of the models being used, we expect that the models used in these publications should also have a clear peak in the prior distributions in both pressure-energy density and mass–radius space. An exception is Model C of Steiner et al. (2013), which is formulated to give a flat prior distribution in pressure-energy density space (although it may still have a peak in the mass–radius space prior distribution due to the non-linear nature of the Tolman–Oppenheimer–Volkoff equation mapping). That the EoS model and priors have a major influence is clear from Figure 4 of Steiner et al. (2013), which shows the inferred posterior distributions for one of the symmetry energy parameters. There is clear variation between models, and the posterior for Model C is noticeably broader. To perform a consistent comparison with the EoS parameters inferred in our work, we would need to know the prior distributions on the pressure-energy density and mass–radius space associated with the EoS parameterizations and priors used in these previous works (the equivalent, for those models, of our Figure 1). However, the necessary details are not shown in these earlier works, so further study will be required to make a robust comparison.

Recently, the gravitational wave (GW) observation of the NS binary inspiral GW170817 (Abbott et al. 2017), where the progenitor is widely accepted to be a NS–NS system, has

provided an independent method for constraining the EoS by measuring tidal effects of the NS in the evolution of GW phase. The dominant tidal GW imprint depends on the characteristic tidal deformability parameters $\Lambda = (2/3)k_2 (c^2 R_{\text{eq}})^5 / (GM)^5$, where $k_2$ is the EoS-dependent Love number (Flanagan & Hinderer 2008). Several studies, as discussed below, infer $\Lambda$ from the GW data (e.g., Abbott et al. 2017, 2018a, 2019; De et al. 2018). For example, work performed by the LIGO-Virgo Scientific Collaborations (Abbott et al. 2019) quantifies the impact of the choice of spin priors and systematic uncertainties in the waveform models, and find these to be non-negligible yet smaller than the statistical errors. Specializing to a low-spin prior with a dimensionless value of less than 0.05 (as expected from extrapolating the spin-down of observed Galactic binary pulsars that will merge within a Hubble time), a representative GW model, and the case where both binary components are assumed to be NSs and have the same EoS, Abbott et al. (2018a) inferred constraints on the EoS and the radius using two methods: a parameterized EoS and approximate EoS-insensitive relations. Although the inferred masses from the GW data involved in GW170817 (Abbott et al. 2017, 2019) were similar to PSR J0030+0451 within the measurement uncertainties, the compatibility of EoS results is difficult to assess in detail because of the different priors imposed in the analysis. However, the inferred 90% credible intervals for the radii of the two components $R \in [9.1, 12.8]$ km and $R \in [10.5, 13.3]$ km with the two methods are consistent with the results in this Letter. Independent analyses of GW170817-only data (De et al. 2018) and results (e.g., Annala et al. 2018; Lim & Holt 2018; Most et al. 2018; Raithel et al. 2018; Tews et al. 2018b) obtained compatible constraints of the radius in a broad range of $R \in [9, 14]$ km. With the ongoing third observing run of the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo detectors at a higher detector sensitivity and further improvements planned (Abbott et al. 2018b; Shoemaker 2019), GW measurements of a greater number of NS binaries (encompassing both NS–NS and NS–black hole systems) will yield tighter EoS constraints in the coming years.

### 4.5. Final Remarks

We have studied the implications of the available PSR J0030+0451 mass–radius likelihood information for dense matter EoS knowledge. The likelihood function of mass and radius is predominantly sensitive to their combination in the compactness ratio, conditional on the current NICER data set and X-ray pulse-profile modeling. The posterior information gain over our choice of prior knowledge is weak in the joint context of both prior constraints imposed by cEFT interactions at nuclear densities, and all EoSs being able to support a 1.97 $M_{\odot}$ NS. This is a consequence of the substantial overlap between the relatively broad mass–radius likelihood function and our narrowly peaked prior. However, we have shown how our methods can be applied to data obtained through pulse-profile modeling of MSPs. Our understanding of the nature of dense matter is expected to improve in the near future with the constraining power offered by the NICER mission: NICER is concurrently observing rotation-powered MSPs such as PSR J0437–4715 that have an independent mass measurement derived via radio timing. Moreover, joint radio and X-ray information from these MSPs promises synergism with the radio information from high-mass pulsars such as PSR J0348+0432.

<footnote n="30">³⁰ Note that for the three representative EoSs of Hebeler et al. (2013) there are some values of central energy density where the slope of the compactness curve changes. These values of density correspond where the piecewise polytropes are matched in these EoSs.</footnote>
This work was supported in part by NASA through the NICER mission and the Astrophysics Explorers Program. T.E.R. and A.L.W. acknowledge support from ERC Starting Grant No. 639217 CSINEUTRONSTAR (PI: Watts). A.L.W. would like to thank Andrew Steiner for useful discussions on the role of priors in previously published results. The authors would also like to thank Kent Wood for helpful comments. 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. S.K.G., K.H., and A.S. acknowledge support from the DFG through SFB 1245. G.R., T.H., and S.N. are grateful for support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) through the VIDI and Projectruimte grants (PI: Nissanke). S.G. acknowledges the support of the CNES. 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. R.M.L. acknowledges the support of NASA through Hubble Fellowship Program grant HST-HF2-51440.001. This research has made extensive use of NASA's Astrophysics Data System Bibliographic Services (ADS) and the arXiv.

*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), 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), MULTINEST (Feroz et al. 2009), PYMULTINEST (Buchner et al. 2014), RNS (Stergioulas & Friedman 1995).

## Appendix

### Likelihood Function Modification

Here we give the prescription for modifying the likelihood function to crudely simulate the effect of increased exposure time. Let

$$L^{\dagger}(M, R) \propto L(M, R) \frac{p^{\dagger}(M, R|d)}{p(M, R|d)},$$

where $p^{\dagger}(M, R|d)$ is given by isotropic compression of mass with flat $p(M, R)$:

$$X = \int_{C^{\dagger}} p^{\dagger}(M, R | d) dM dR = \int_{C} p(M, R|d) dM dR,$$

such that $\forall X \le 1$, region $C^\dagger$ has Euclidean area $A^\dagger$ and $A^\dagger/A \in [1/\sqrt{n}, 1/n]$, where $A$ is the area of region $C$ and $n \propto T$ is the factor increase in counts, which scales linearly with exposure time $T$. A fallacy here is that the likelihood function must be zero exterior of $C^\dagger(X = 1)$, which is a region within prior support, lest not all credible regions shrink by this factor—the latter is realistic, however, and can be viewed as consistent with numerical operation in finite-sample context.

In this work we simulated increased exposure time by assuming a scaling of $A^\dagger/A \approx 1/2$. Numerically, given a set $\{(s, w)_i\}_{i=1...N}$ of importance samples with weights $\{w_i\}$ from the density $p(M, R|d)$, we: (i) define a fiducial vector $l$ as that of the sample reporting the highest nuisance-marginalized likelihood value; (ii) calculate $\forall i, s_i^\dagger = (s_i - l)/\sqrt{2} + l$; and (iii) define $\{(s^\dagger, w)_i\}_{i=1...N}$ as a set of importance samples from density $p^\dagger(M, R|d)$.

## ORCID iDs

G. Raaijmakers https://orcid.org/0000-0002-9397-786X

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

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

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

T. Hinderer https://orcid.org/0000-0002-3394-6105

S. Nissanke https://orcid.org/0000-0001-6573-7773

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

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

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

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

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

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

## 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. 2019, PhRvX, 9, 011001

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

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018b, LRR, 21, 3

Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, PhRvC, 58, 1804

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

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

Baldo, M., Bombaci, I., & Burgio, G. F. 1997, A&A, 328, 274

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

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

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

Bogdanov, S. 2016, EPJA, 52, 37

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

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

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

Cromartie, H., Fonseca, E., Ransom, S. M., et al. 2019, NatAs, in press (doi:10.1038/s41550-019-0880-2)

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

Drischler, C., Hebeler, K., & Schwenk, A. 2019, PhRvL, 122, 042501

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

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

Flanagan, E. E., & Hinderer, T. 2008, PhRvD, 77, 021502

Forum, M. P. 1994, MPI: A Message-Passing Interface Standard, Tech. rep., Knoxville, TN, USA

Fraga, E. S., Kurkela, A., & Vuorinen, A. 2014, ApJL, 781, L25

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

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

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

Harding, A. K., & Muslimov, A. G. 2002, ApJ, 568, 862

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

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

Hebeler, K., Holt, J. D., Menéndez, J., & Schwenk, A. 2015, ARNPS, 65, 457

Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11

Hebeler, K., & Schwenk, A. 2010, PhRvC, 82, 014314

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

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/

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

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

Kullback, S., & Leibler, R. A. 1951, Ann. Math. Statist., 22, 79
Lattimer, J. M., & Prakash, M. 2016, PhR, 621, 127

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

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

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

Lynn, J. E., Tews, I., Carlson, J., et al. 2016, PhRvL, 116, 062501

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

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

Miller, M. C., Chirenti, C., & Lamb, F. K. 2019a, ApJ, in press (arXiv:1904.08907)

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

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

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

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

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

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

Pandharipande, V. R., & Smith, R. A. 1975, NuPhA, 237, 507

Paschalidis, V., & Stergioulas, N. 2017, LRR, 20, 7

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

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

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

Raithel, C. A., Özel, F., & Psaltis, D. 2018, ApJL, 857, L23

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

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

Shoemaker, D. 2019, arXiv:1904.03187

Skilling, J. 2006, BayAn, 1, 833

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

Stergioulas, N., & Friedman, J. L. 1995, ApJ, 444, 306

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

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

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

Watts, A. L. 2019, arXiv:1904.07012

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