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

Our understanding of the dense matter equation of state (EOS) of neutron stars has made significant progress over the last few years due to the arrival of new avenues to measure observables like mass, radius, and tidal deformability, that connect to the behavior of matter at supranuclear densities. Recently NASA's X-ray timing telescope, the Neutron Star Interior Composition Explorer (NICER), has delivered the first joint measurement of mass and radius through pulse profile modeling of the millisecond pulsar PSR J0030+0451 (Miller et al. 2019; Riley et al. 2019). The impact of this measurement on the dense matter EOS has been extensively studied in various EOS frameworks (see, e.g., Miller et al. 2019; Raaijmakers et al. 2019, 2020; Dietrich et al. 2020; Jiang et al. 2020; Landry et al. 2020; Essick et al. 2020b; Al-Mamun et al. 2021), including EOS with phase transitions to quark matter (see, e.g., Alvarez-Castillo et al. 2020; Blaschke et al. 2020; Li et al. 2020; Tang et al. 2021; Xie & Li 2021) and models that explore the possibility of there being two stable neutron star branches (Christian & Schaffner-Bielich 2020).

Concurrently, the second and third observing runs of LIGO/Virgo have so far resulted in the confirmed gravitational wave detections of two (most-likely) binary neutron star mergers: GW170817 (Abbott et al. 2017c, 2019a) and GW190425 (Abbott et al. 2020a). By accurately measuring the gravitational wave phase, limits can be put on the EOS-dependent tidal deformability of the neutron stars (Flanagan & Hinderer 2008; Hinderer et al. 2010). While for GW170817 the tidal deformability could be measured within a 90% highest posterior density interval when adopting low-spin priors (see, e.g., Abbott et al. 2018, 2019b), the low signal-to-noise ratio of GW190425 resulted in only weak upper limits on the tidal deformability even when assuming low spins (Abbott et al. 2020a). We consider the ∼2.6 M_⊙ secondary object in GW190814 (Abbott et al. 2020c) to be a black hole (Nathanail et al. 2021), and will therefore not use this third event in our analysis.

At nuclear densities, the EOS is well constrained by nuclear theory and experiments (see, e.g., Tsang et al. 2012; Lattimer & Lim 2013; Huth et al. 2021). In particular, many-body calculations based on chiral effective field theory (EFT) interactions have enabled systematic predictions for the neutron matter EOS up to nuclear saturation density including theoretical uncertainties (see, e.g., Hebeler et al. 2013; Tews et al. 2013; Lynn et al. 2016; Drischler et al. 2019, 2020). Up to saturation density, the resulting symmetry energy and pressure of neutron matter are also consistent with extractions from nuclear experiments (Lattimer & Lim 2013), including from measurements of the dipole polarizability of neutron-rich nuclei (Roca-Maza et al. 2015; Birkhan et al. 2017; Kaufmann et al. 2020). Taking these results at nuclear densities, combined with standard crust EOS, different extrapolations to high densities have been found to lead to neutron star radii consistent with all multimessenger observations (see, e.g., Annala et al. 2020; Dietrich et al. 2020; Raaijmakers et al. 2020; Essick et al. 2020b; Biswas et al. 2021). Recently, the results of PREX-II have pointed to higher pressures (Adhikari et al. 2021; Reed et al. 2021), but with very large uncertainties, so that in a combined analysis with astrophysical and chiral EFT constraints, the overall consistency still persists (Essick et al. 2021).

NICER data have now enabled a joint estimate of the mass and radius of the high-mass, rotation-powered millisecond pulsar PSR J0740+6620. Since PSR J0740+6620 (unlike PSR J0030+0451) is in a binary with an inclination that allows measurement of the Shapiro delay, its mass can be measured independently via radio timing. Cromartie et al. (2020) reported a mass of ${2.14}_{-0.09}^{+0.10}$ M_⊙, and a joint campaign by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) and the Canadian Hydrogen Intensity Mapping Experiment (CHIME)/Pulsar collaborations has now resulted in an updated mass of 2.08 ± 0.07 M_⊙ (Fonseca et al. 2021).

Riley et al. (2021) have used this mass measurement as an informative prior for pulse-profile modeling analysis that is joint over the phase-resolved spectroscopic data from NICER and phase-averaged data from the XMM-Newton European Photon Imaging Camera. The inclusion of the smaller XMM-Newton (hereafter XMM) data set allows for better constraints on the proportion of the X-ray emission that is attributable to background rather than PSR J0740+6620, ultimately acting to cut out solutions with high compactness. This results in an inferred radius of ${12.39}_{-0.98}^{+1.30}\,\mathrm{km}$, and a mass of ${2.072}_{-0.066}^{+0.067}$ M_⊙ that is little changed from the radio prior. For a full description of the methodology employed in the mass–radius inference we refer the reader to Riley et al. (2021).

In this Letter, we use the mass and radius from Riley et al. (2021) for PSR J0740+6620 as input for inferring the dense matter EOS, combining it with other constraints from nuclear theory and multimessenger observations. It should be considered as a follow-up to our previous work that built on NICER's results for PSR J0030+0451 (Raaijmakers et al. 2019, 2020), where in this work we explore also a broader range of multimessenger constraints. As the high-density constraints from astrophysical observations get more precise, with the new NICER results and future LIGO/Virgo measurements, it will be intriguing to see them play out with the present nuclear constraints. In this Letter, we also explore this for the new NICER results and how they constrain the EOS above nuclear densities starting from different chiral EFT calculations. ¹⁰

In this work we will closely follow the analysis framework developed previously in Greif et al. (2019) and Raaijmakers et al. (2019, 2020). Below, we summarize this method and highlight several updates to the framework.

We consider two EOS parameterizations: (i) a piecewise-polytropic (PP) model with three segments between varying transitions densities (Hebeler et al. 2013), and (ii) a speed-of-sound (CS) model first introduced in Greif et al. (2019). To capture the uncertainty in the EOS around nuclear saturation density (n₀ = 0.16 fm⁻³), both parameterizations are matched to a power-law fit of a range of EOS calculated from chiral EFT interactions (Hebeler & Schwenk 2010; Hebeler et al. 2013) below 1.1n₀. At densities below 0.5n₀ this power-law fit is connected to the BPS crust EOS (Baym et al. 1971).

To constrain these EOS parameterizations, governed by the EOS parameters θ , we employ Bayes' theorem and write the posterior distributions of the EOS parameters and central energy densities ε as

$$\begin{eqnarray}&&p({\boldsymbol{\theta }},{\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{d}},{\mathbb{M}})\propto p({\boldsymbol{\theta }}\,| \,{\mathbb{M}})\,p({\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{\theta }},{\mathbb{M}})\,p({\boldsymbol{d}}\,| \,{\boldsymbol{\theta }},{\mathbb{M}}),\end{eqnarray} \tag{ 1 }$$

where ${\mathbb{M}}$ denotes the model including all assumed physics and d the data set used to constrain the EOS, consisting of, e.g., radio, X-ray, and gravitational wave data. When assuming each of these data sets to be independent of each other, we can separate the likelihoods and write

$$\begin{eqnarray}&&\begin{array}{l}p({\boldsymbol{\theta }},{\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{d}},{\mathbb{M}})\propto p({\boldsymbol{\theta }}\,| \,{\mathbb{M}})\,p({\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{\theta }},{\mathbb{M}})\\ \quad \times \,\prod _{i}p({{\rm{\Lambda }}}_{1,i},{{\rm{\Lambda }}}_{2,i},{M}_{1,i},{M}_{2,i}\,| \,{{\boldsymbol{d}}}_{\mathrm{GW},\ {\rm{i}}}(,{{\boldsymbol{d}}}_{\mathrm{EM},\ {\rm{i}}}))\\ \quad \times \,\prod _{j}p({M}_{j},{R}_{j}\,| \,{{\boldsymbol{d}}}_{\mathrm{NICER},{\rm{j}}})\\ \quad \times \,\prod _{k}p({M}_{k}\,| \,{{\boldsymbol{d}}}_{\mathrm{radio},{\rm{k}}}).\end{array}\end{eqnarray} \tag{ 2 }$$

Here the products run over the number of different observed stars, or mergers, in the case of the gravitational wave data. Furthermore, in Equation (2) we have equated the nuisance-marginalized likelihoods to the nuisance-marginalized posterior distributions derived in Riley et al. (2019, 2021), Abbott et al. (2019a, 2020a), and Fonseca et al. (2021). This approximation is justifiable when the priors used in estimating these nuisance-marginalized posterior distributions are uninformative, which for simplicity we will assume to be a uniform prior in this case. The posterior distributions derived by Riley et al. (2019, 2021) already use a jointly uniform prior in mass and radius. The posterior distributions derived by Abbott et al. (2019a, 2020a) use a jointly uniform prior in the tidal deformabilities of the two components Λ_i within the range Λ_i ⊂ [0, 5000] (for GW190425 the upper bound of Λ₂ was set to 10⁴.). The prior on the detector frame masses, which are redshifted with respect to the source frame masses (${M}_{\det }={M}_{i}(1+z)$), is uniform within the range ${M}_{\det }\subset [0.5,7.7]$ and ${M}_{\det }\subset [1,5.31]$ for GW170817 and GW190425 respectively. However, the posterior distribution on component masses from gravitational waves is highly degenerate because of the accurately measured chirp mass ${{ \mathcal M }}_{c}={\left({M}_{1}{M}_{2}\right)}^{3/5}/{\left({M}_{1}+{M}_{2}\right)}^{1/5}$. To speed up the convergence of our parameter estimation, we therefore transform the gravitational wave posterior distributions to include the two tidal deformabilities, chirp mass and mass ratio q, while reweighing such that the prior distribution on these parameters is uniform. Further, we also fix the chirp mass to its median value, since the small uncertainty in this parameter does not affect the EOS parameter estimation (see Raaijmakers et al. 2020), and thus have

$$\begin{eqnarray}&&\begin{array}{l}p({\boldsymbol{\theta }},{\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{d}},{\mathbb{M}})\propto p({\boldsymbol{\theta }}\,| \,{\mathbb{M}})\,p({\boldsymbol{\varepsilon }}\,| \,{\boldsymbol{\theta }},{\mathbb{M}})\\ \quad \times \,\prod _{i}p({{\rm{\Lambda }}}_{1,i},{{\rm{\Lambda }}}_{2,i},{q}_{i}\,| \,{{ \mathcal M }}_{c},{{\boldsymbol{d}}}_{\mathrm{GW},\ {\rm{i}}}(,{{\boldsymbol{d}}}_{\mathrm{EM},{\rm{i}}}))\\ \quad \times \,\prod _{j}p({M}_{j},{R}_{j}\,| \,{{\boldsymbol{d}}}_{\mathrm{NICER},{\rm{j}}})\\ \quad \times \,\prod _{k}p({M}_{k}\,| \,{{\boldsymbol{d}}}_{\mathrm{radio},{\rm{k}}}).\end{array}\end{eqnarray} \tag{ 3 }$$

Fixing the chirp mass means that the vector ε only contains one central density per merger, where the tidal deformability of the second component is now set by Λ₂ = Λ₂( θ ; q). If a gravitational wave event has an associated electromagnetic (EM) counterpart, the likelihood for that event becomes a product of the nuisance-marginalized posterior distribution from the gravitational wave data and the nuisance-marginalized posterior distribution from the EM analysis, such that

$$\begin{eqnarray}&&\begin{array}{l}p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},q\,| \,{{ \mathcal M }}_{c},{{\boldsymbol{d}}}_{\mathrm{GW}},{{\boldsymbol{d}}}_{\mathrm{EM}})\propto p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},q\,| \,{{ \mathcal M }}_{c},{{\boldsymbol{d}}}_{\mathrm{GW}})\\ \quad \times \,p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},q\,| \,{{ \mathcal M }}_{c},{{\boldsymbol{d}}}_{\mathrm{EM}}).\end{array}\end{eqnarray} \tag{ 4 }$$

Obtaining the posterior distribution $p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},q\,| \,{{ \mathcal M }}_{c},{{\boldsymbol{d}}}_{\mathrm{EM}})$ is discussed in Section 3.2.1 for the specific case of AT2017gfo, the kilonova associated with GW170817 (see, e.g., Arcavi et al. 2017; Chornock et al. 2017; Cowperthwaite et al. 2017; Coulter et al. 2017; Kasliwal et al. 2017; Nicholl et al. 2017; Tanvir et al. 2017; Abbott et al. 2017a, 2017b).

We then sample from the posterior distribution $p({\boldsymbol{\theta }},{\boldsymbol{\varepsilon }}\,| {\boldsymbol{d}},{\mathbb{M}})$, compute the corresponding M, R, and Λ, and then evaluate the likelihood by applying a kernel density estimation to the posterior distributions from Riley et al. (2019, 2021) and Abbott et al. (2019a, 2020a) using the nested sampling software MultiNest. The same prior distribution $p({\boldsymbol{\theta }}\,| \,{\mathbb{M}})$ is used as in previous work; we refer the reader to Section 2.3 of Raaijmakers et al. (2020) and references therein for a more detailed description.

In this Section we investigate the impact of the Riley et al. (2021) mass–radius measurement for PSR J0740+6620 on the dense matter EOS, both separately and when combined with previous constraints.

3.1. Radio Mass Measurement of PSR J0740+6620

First, we constrain the EOS using the updated mass measurement of 2.08 ± 0.07 M_⊙ for PSR J0740+6620 derived using radio timing (Fonseca et al. 2021), and compare this to the constraints from the previously published mass of ${2.14}_{-0.09}^{+0.1}$ M_⊙ (Cromartie et al. 2020). In Figure 1 we show the posterior distribution on EOS parameters θ when transformed to the mass–radius parameter space. We note that, as expected, the slightly lower updated mass measurement shifts the posterior distributions to lower maximum neutron star masses and lower radii, although the effect is almost negligible. Since the radio timing mass measurement is already incorporated in the joint mass–radius estimate from NICER we will not use this measurement in the remainder of this work.

Zoom In Zoom Out Reset image size

3.2. GW170817 and GW190425

The gravitational wave events GW170817 and GW190425 have so far been the only confirmed neutron star binary mergers during the recent observing runs of the LIGO/Virgo collaboration (Abbott et al. 2021). Although both events have a non-negligible chance of being neutron star–black hole mergers (see, e.g., Yang et al. 2018, Ascenzi et al. 2019, Coughlin & Dietrich 2019, and Hinderer et al. 2019 for GW170817 and, e.g., Kyutoku et al. 2020 and Han et al. 2020 for GW190425), in the following we will assume both objects to be neutron stars. We use the low-spin ¹¹ posterior distributions on tidal deformability and mass ratio with the IMRPhenomPv2_NRTidal ¹² waveform model (Hannam et al. 2014; Khan et al. 2016; Dietrich et al. 2019) for GW170817 and GW190425. Furthermore we use the median chirp mass values of ${{ \mathcal M }}_{c}=1.186$ M_⊙ for GW170817 ¹³ and ${{ \mathcal M }}_{c}=1.44$ M_⊙ for GW190425. ¹⁴

The upper panels of Figure 3 show the posterior distributions on the EOS for both events in the mass–radius space. We note that the constraints on tidal deformability from GW170817 give more support to softer EOS, although the 95% credible region spans a relatively large range of radii. GW190425 only led to weak upper limits on the tidal deformability due to its low signal-to-noise ratio and single-detector detection. The EOS is however constrained as a result of the high mass of the primary component (with 95% credible range 1.60–1.87 M_⊙), excluding EOS that do not support these masses.

3.2.1. AT2017gfo

Following the detection of GW170817 an EM counterpart was observed across the frequency spectrum; see, e.g., Abbott et al. (2017a, 2017b) and references therein, Coulter et al. (2017), Chornock et al. (2017), Drout et al. (2017), Hallinan et al. (2017), Kasliwal et al. (2017, 2019), Margutti et al. (2017), Pian et al. (2017), Smartt et al. (2017), and Troja et al. (2017). Of particular interest here is the thermal infrared–optical–ultraviolet transient powered by radioactive decay of r-process nucleosynthesis in the neutron-rich material ejected during merger: the so-called kilonova or macronova (e.g., Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010). The kilonova properties depend on the mass, velocity, and composition of the ejected material, which in turn depend on the binary progenitor parameters such as the tidal deformability of the neutron stars. Using this connection it is possible to constrain the EOS from the kilonova light curve (see, e.g., Coughlin et al. 2018; Hinderer et al. 2019; Radice & Dai 2019; Capano et al. 2020; Dietrich et al. 2020).

Here we analyze the bolometric luminosity of GW170817 (as compiled in Kasliwal et al. 2017) via the new Bayesian framework outlined in Raaijmakers et al. (2021b). We consider a two-component kilonova model, where the first component, the dynamical ejecta, is associated with material ejected through tidal forces and the shock interface between the two neutron stars (see, e.g., Radice et al. 2018 and references therein). The second component is associated with neutrino-driven winds or material ejected through viscous forces. We connect the outflow properties of these components to the binary progenitor properties using the formulae presented in Krüger & Foucart (2020) for dynamical ejecta and disk mass, which are fitted to numerical simulations of compact mergers. The velocity of the dynamical ejecta is calculated using the formula in Coughlin et al. (2019), while the velocity of the disk wind ejecta is left as a free parameter. The dynamical ejecta includes both material ejected through tidal forces and material ejected through shocks on the contact interface between the stars (see, e.g., Sekiguchi et al. 2016; Dietrich & Ujevic 2017; Tanaka et al. 2020; Nedora et al. 2021). To distinguish these we consider two different opacities in the dynamical ejecta, corresponding to tidal tail and shock ejecta, where the latter is less neutron-rich compared to the tidal tail and thus has a lower opacity (see Table 1). For simplicity we take a single opacity for the disk wind ejecta. The outflow properties are then connected to a bolometric luminosity through the semi-analytic light-curve model of Hotokezaka & Nakar (2020). The priors on all parameters are shown in Table 1.

Table 1. Parameters Used in the Model Described in Section 3.2.1 and Their Prior Support in the Analysis of AT2017gfo

Parameters	Prior Density and Support
Binary properties
${{ \mathcal M }}_{c}$ (M⊙)	∼U(1.18, 1.2)
q	∼U(0.2, 1)
Λ1	∼U(0, 2500)
Λ2	∼U(0, 2500)
Ejecta and light-curve properties
Mdyn (M⊙)	Equation (2), Raaijmakers et al. (2021b)
vdyn (c)	Equation (D5), Coughlin et al. (2019)
${v}_{\min ,\ \mathrm{dyn}}$ (c)	∼U(0.1, 1.0) vdyn
${v}_{\max ,\ \mathrm{dyn}}$ (c)	∼U(1.5, 2.5) vdyn
vκ (c)	∼U(v${}_{\min ,\ \mathrm{dyn}}$, ${{\rm{v}}}_{\max ,\ \mathrm{dyn}}$)
κlow (cm2 g−1)	∼U(0.1, 5)
κhigh (cm2 g−1])	∼U(5, 30)
Mwind (M⊙)	Equation (10), Raaijmakers et al. (2021b)
vwind (c)	∼U(0.03, 0.15)
${v}_{\min ,\ \mathrm{wind}}$ (c)	∼U(0.1, 1.0) vwind
${v}_{\max ,\ \mathrm{wind}}$ (c)	∼U(1.5, 2.0) vwind
κwind [cm2 g−1]	∼U(0.1, 5)

Note. The notation U(a, b) here means uniformly drawn between boundaries a and b.

Download table as:  ASCIITypeset image

The fit to the bolometric luminosity of AT2017gfo using the data compiled in Kasliwal et al. (2017) is shown in Figure 2, showing all data points to be contained within the 95% credible region of the posterior distribution. In the lower panels of Figure 3 we show the updated prior distribution for the EOS with GW170817 and with the inclusion of AT2017gfo. The EM data give more posterior support to stiffer over softer EOS, due to the estimated ejected mass requiring a neutron star with larger tidal deformability. The estimated radius of a 1.4 M_⊙ neutron star for the PP and CS models is ${12.12}_{-1.44}^{+1.10}\,\mathrm{km}$ and ${11.53}_{-1.15}^{+1.16}\,\mathrm{km}$, respectively, which is broadly consistent with multimessenger constraints obtained by other works (see, e.g., Coughlin et al. 2019; Capano et al. 2020; Dietrich et al. 2020; Breschi et al. 2021; Nicholl et al. 2021). Important to note is that the EM modeling of the kilonova here is simplified and relies on a few assumptions that are known to affect results, such as spherical ejecta geometry (see, e.g., Heinzel et al. 2021; Korobkin et al. 2021), fixed nuclear heating rate (see, e.g., Barnes et al. 2020), and an incomplete mapping between properties of the binary system and the ejecta outflows. It is also dependent on the choice of light-curve modeling, where the distinction can be made between semi-analytic modeling (such as in this work and, e.g., Breschi et al. 2021; Nicholl et al. 2021) and interpolating between radiative transfer simulations (e.g., Coughlin et al. 2018; Dietrich et al. 2020). We use a semi-analytical model from Hotokezaka & Nakar (2020), which for the current statistical uncertainty in gravitational wave parameter estimation and uncertainty in light-curve observations produces consistent results to full radiative transport models (see, e.g., Coughlin et al. 2020a, 2020b), although this will change in the future with improved gravitational wave detectors and optical telescopes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. NICER Mass–Radius and Multimessenger Constraints

Next we study the constraints on the EOS from the new mass–radius estimate of PSR J0740+6620 using data from NICER and XMM, presented in Riley et al. (2021). They find a radius of ${12.39}_{-0.98}^{+1.30}\,\mathrm{km}$ and a mass of ${2.072}_{-0.066}^{+0.067}$ M_⊙, where the upper and lower limits bound the 68% credible regions. The EOS results are shown in Figure 4, both in energy density–pressure and mass–radius space. From the Kullback–Leibler divergence (Kullback & Leibler 1951) plotted as a function of energy density in the upper insets, we find that, especially at higher energy densities, there is a significant information gain from prior-to-posterior. Note that similar but, especially for the CS model, broader constraints are found for the posterior distribution when only using the radio mass measurement of PSR J0740+6620, as indicated by the orange dashed-dotted lines. This is a result of the mass–radius estimate of PSR J0740+6620 being very consistent with our prior ranges informed by low-density chiral EFT calculations. The chiral EFT calculations do exclude, however, stiffer EOS with radii >14 km, where the mass–radius posterior of PSR J0740+6620 has non-negligible posterior support. For the CS model this effect is stronger as additional constraints on the speed of sound at 1.5n₀ in the CS model lead to overall smaller radii than in the PP model (see Section 2.3 of Raaijmakers et al. 2020).

Zoom In Zoom Out Reset image size

In the mass–radius space we also plot the EOS constraints given the joint NICER mass–radius estimate excluding the XMM data. For this analysis Riley et al. (2021) report a value of ${11.29}_{-0.81}^{+1.20}\,\mathrm{km}$ for the radius and ${2.078}_{-0.063}^{+0.066}\,{M}_{\odot }$ for the mass. As this joint mass–radius estimate has slightly more posterior support for lower radii, the corresponding EOS constraints suggest a softening of the EOS at high densities. These results should be interpreted with caution, however, because the NICER-only analysis leads to an under-prediction of the background (the contribution from instrumental or astrophysical background to the unpulsed component of the pulse profile). This results in more of the unpulsed component being attributed to the hot regions via high-compactness solutions. The XMM data show that a larger component of the unpulsed emission must come from true background, eliminating these high-compactness solutions and increasing the inferred radius in the joint NICER–XMM analysis (see also Section 4.2 in Riley et al. 2021).

Finally, in Figure 5 we show the constraints on the EOS from PSR J0740+6620, PSR J0030+0451 (first derived in Raaijmakers et al. 2019, but here no information on high-mass pulsars is included) and the combination of the two pulsars. Note that for the combined constraints, most of the information comes from PSR J0740+6620, since the 68% credible region of the mass–radius posterior of PSR J0030+0451 covers a broad range in radii that are consistent with the EOS constraints from PSR J0740+6620.

Zoom In Zoom Out Reset image size

In the lower panels of Figure 5 we show the combined constraints on the EOS including mass–radius estimates from PSR J0740+6620 (Riley et al. 2021), PSR J0030+0451 (Riley et al. 2019) and mass–tidal deformability estimates from GW170817 (Abbott et al. 2019a) and GW190425 (Abbott et al. 2020a), and the kilonova AT2017gfo (Kasliwal et al. 2017). We find that, in particular, the pulsar mass–radius estimates by NICER favor stiffer EOS, as well as GW170817 when the associated kilonova AT2017gfo (Kasliwal et al. 2017) is included. The weak constraints from GW190425 on the tidal deformability are also broadly consistent with the constraints coming from the other sources. As a comparison we show the posterior distribution when combining all analyses excluding the mass–radius estimate of PSR J0740+6620, but with the radio mass measurement of Fonseca et al. (2021). We note that the additional radius information on PSR J0740+6620 constrains the softer EOS, especially for the CS model.

To investigate the impact of the EOS constraints from nuclear physics we compare our analysis of PSR J0740+6620 using four different chiral EFT uncertainty bands. All bands are based on microscopic calculations for pure neutron matter, which are then extended to neutron star matter in beta-equilibrium using the formalism discussed in Hebeler et al. (2013). In order to improve the description of all employed EOS, we generalized the density dependence of the energy-density functional (see Equation (2) in Hebeler et al. 2013) by enlarging the range of the exponent γ to $\gamma \in \left[1.2,2.5\right]$.

The results from Hebeler et al. (2013) formed the basis of our previous studies (Raaijmakers et al. 2019, 2020). The calculations for pure neutron matter were initially performed in Hebeler & Schwenk (2010) using many-body perturbation theory, while the uncertainty band results mainly from variations of the couplings involved in three-nucleon interactions. Second, in Tews et al. (2013) the calculations for neutron matter were improved by including for the first time all two-, three-, and four-neutron interactions to next-to-next-to-next-to-leading order (N³LO), which are predicted in a parameter-free way for neutron matter (see, e.g., Hebeler et al. 2015; Hebeler 2021 for reviews). Third, in Drischler et al. (2019) the calculations were further optimized by improving the treatment of three-nucleon interactions and extending the many-body expansion to higher orders. In addition, the EOS uncertainty bands also include effects from variations of regulator scales in state-of-the-art nucleon–nucleon and three-nucleon interactions. In this work, we use the combined 450 MeV and 500 MeV N³LO uncertainty bands from Drischler et al. (2019). Finally, we include results of Lynn et al. (2016). These were obtained by non-perturbative quantum Monte Carlo simulations of neutron matter at next-to-next-to-leading order (N²LO). This represents a completely different many-body method than those used for the other three bands, and the results of Lynn et al. (2016) are also based on a different set of local two- and three-nucleon interactions derived from chiral EFT.

Similar to Raaijmakers et al. (2020) we approximate the EOS within these bands with a single polytrope P = Nn^Γ. However, to obtain a better fit to the additional bands considered here, we vary the polytropic index Γ as a function of the normalization N,

$$\begin{eqnarray}&&{\rm{\Gamma }}(N)=\displaystyle \frac{(N-{N}_{\min })}{({N}_{\max }-{N}_{\min })}({{\rm{\Gamma }}}_{\max }-{{\rm{\Gamma }}}_{\min })+{{\rm{\Gamma }}}_{\min },\end{eqnarray} \tag{ 5 }$$

where ${N}_{\min /\max }$ and ${{\rm{\Gamma }}}_{\min /\max }$ are determined by fitting a polytrope to the lower and upper bound of the band. In Figure 6 we show the four different bands for the pressure of neutron star matter with an example of the fit through each band. This shows the consistency of these different chiral EFT calculations, with different methods, interactions, and approximations. The first point of the band where n/n₀ > 0.5 is matched to the BPS crust EOS at 0.5n₀ via a linear interpolation.

Zoom In Zoom Out Reset image size

We study the dependence of the EOS constraints on the different chiral EFT bands by inferring the EOS from the mass–radius estimate of PSR J0740+6620 using each band and both high-density parameterizations. The results are shown in Figure 7. We also show the 95% credible region of the updated prior distribution when directly joining the PP or CS high-density parameterization to the crust EOS at 0.5ρ_ns. As expected the chiral EFT calculations mostly exclude stiffer EOS. While the different chiral EFT bands yield very good agreement on the upper bound of the radius estimates, the lower bound on the radius does slightly depend on the chiral EFT band used, especially at lower neutron star masses, depending on how soft the chiral EFT band is (see Figure 6).

Zoom In Zoom Out Reset image size

In the lower panels of Figure 7 we also show the posterior distributions on the pressure at densities n = 1.5n₀ and n = 2n₀ above the chiral EFT bands. These results demonstrate that the PSR J0740+6620 mass–radius measurement systematically prefers higher pressures at these densities compared to the corresponding prior distributions of each chiral EFT band. Furthermore, the posteriors at n = 2n₀ agree very well for all chiral EFT bands and are peaked around P ∼ 10^34.5dyn cm⁻² ∼ 20 MeV fm⁻³.

In this Letter, we have investigated the constraints on the EOS posed by the new joint mass–radius estimate from NICER–XMM data (Riley et al. 2021), and compared and combined with multimessenger EOS constraints from radio timing, gravitational wave mergers, and their counterparts, and the previous PSR J0030+0451 mass–radius estimate by NICER. In Table 2 we summarize the results obtained in Sections 3 and 4 for the constraints on the radius of a 1.4, 1.6, and 1.8 M_⊙ neutron star, as well as ΔR = R₂ − R_1.4, and the maximum mass of a non-rotating neutron star M_TOV, as well as the constraints on the central energy density and pressure for PSR J0740+6620.

Table 2. Key Quantities from the Posterior Distributions Obtained in Sections 3 and 4: Radius of a 1.4 M_⊙, 1.6 M_⊙, and 1.8 M_⊙ Neutron Star, as Well as ΔR = R₂ − R_1.4, and the Maximum Mass of a Non-rotating Neutron Star M_TOV

	PSR J0740+6620, NICER x XMM				PSR J0030	GW170817	GW170817	Combined	Combined
	Heb 13	Tews 13	Lynn 16	Dri 19	+0451	+ GW190425	+ AT2017gfo	without	with
PP model
R1.4	${12.56}_{-0.91}^{+0.80}$	${12.85}_{-0.95}^{+0.77}$	${12.35}_{-0.98}^{+0.96}$	${12.87}_{-0.98}^{+0.85}$	${12.35}_{-1.99}^{+0.99}$	${11.51}_{-1.47}^{+1.51}$	${12.12}_{-1.44}^{+1.10}$	${12.30}_{-0.76}^{+0.72}$	${12.33}_{-0.81}^{+0.76}$
R1.6	${12.60}_{-1.00}^{+0.87}$	${12.87}_{-1.05}^{+0.87}$	${12.40}_{-1.04}^{+1.03}$	${12.90}_{-1.08}^{+0.94}$	${12.50}_{-2.08}^{+0.96}$	${11.43}_{-1.53}^{+1.68}$	${12.10}_{-1.69}^{+1.23}$	${12.32}_{-0.99}^{+0.92}$	${12.35}_{-0.90}^{+0.83}$
R1.8	${12.62}_{-1.19}^{+0.98}$	${12.86}_{-1.26}^{+1.00}$	${12.42}_{-1.19}^{+1.13}$	${12.89}_{-1.29}^{+1.08}$	${12.68}_{-1.99}^{+0.94}$	${11.65}_{-1.80}^{+1.64}$	${12.22}_{-1.91}^{+1.26}$	${12.29}_{-1.19}^{+1.06}$	${12.33}_{-1.06}^{+0.97}$
ΔR	$-{0.24}_{-1.04}^{+0.65}$	$-{0.17}_{-1.21}^{+1.26}$	$-{0.22}_{-1.05}^{+0.60}$	$-{0.45}_{-1.07}^{+1.14}$	$-{0.13}_{-1.02}^{+0.76}$	$-{0.35}_{-1.09}^{+0.80}$	$-{0.26}_{-1.14}^{+0.77}$	$-{0.30}_{-1.06}^{+0.64}$	$-{0.29}_{-0.98}^{+0.61}$
MTOV	${2.26}_{-0.23}^{+0.15}$	${2.33}_{-0.30}^{+0.14}$	${2.22}_{-0.21}^{+0.19}$	${2.33}_{-0.31}^{+0.18}$	${1.74}_{-0.57}^{+0.66}$	${1.84}_{-0.17}^{+0.51}$	${1.96}_{-0.44}^{+0.42}$	${2.23}_{-0.23}^{+0.15}$	${2.23}_{-0.23}^{+0.14}$
${\mathrm{log}}_{10}({\varepsilon }_{c})$	${14.99}_{-0.16}^{+0.27}$	${14.99}_{-0.18}^{+0.28}$	${15.00}_{-0.16}^{+0.26}$	${14.99}_{-0.19}^{+0.28}$	${14.86}_{-0.13}^{+0.28}$	⋯	⋯	⋯
${\mathrm{log}}_{10}({{\rm{P}}}_{c})$	${35.39}_{-0.24}^{+0.39}$	${35.37}_{-0.26}^{+0.41}$	${35.41}_{-0.25}^{+0.37}$	${35.37}_{-0.28}^{+0.43}$	${34.92}_{-0.21}^{+0.30}$	⋯	⋯	⋯
K	1.00	0.89	1.00	0.85	⋯	⋯	⋯	⋯
CS model
R1.4	${12.27}_{-0.90}^{+0.54}$	${12.49}_{-0.87}^{+0.49}$	${12.16}_{-0.97}^{+0.63}$	${12.56}_{-0.92}^{+0.51}$	${11.51}_{-1.90}^{+1.12}$	${11.18}_{-1.51}^{+1.33}$	${11.53}_{-1.15}^{+1.16}$	${11.98}_{-0.71}^{+0.63}$	${12.18}_{-0.79}^{+0.56}$
R1.6	${12.25}_{-0.94}^{+0.59}$	${12.43}_{-0.92}^{+0.55}$	${12.16}_{-0.99}^{+0.66}$	${12.50}_{-0.96}^{+0.55}$	${11.48}_{-1.92}^{+1.20}$	${10.92}_{-1.56}^{+1.58}$	${11.33}_{-1.41}^{+1.38}$	${11.91}_{-0.94}^{+0.78}$	${12.14}_{-0.84}^{+0.61}$
R1.8	${12.14}_{-1.05}^{+0.69}$	${12.27}_{-1.04}^{+0.66}$	${12.08}_{-1.07}^{+0.74}$	${12.33}_{-1.06}^{+0.65}$	${11.52}_{-1.77}^{+1.23}$	${10.85}_{-1.47}^{+1.73}$	${11.34}_{-1.65}^{+1.44}$	${11.72}_{-1.06}^{+0.96}$	${12.00}_{-0.96}^{+0.74}$
ΔR	$-{0.69}_{-1.02}^{+1.10}$	$-{0.72}_{-1.08}^{+1.12}$	$-{0.58}_{-1.08}^{+1.03}$	$-{1.06}_{-0.83}^{+1.46}$	$-{0.93}_{-0.86}^{+1.31}$	$-{0.93}_{-0.83}^{+1.36}$	$-{0.81}_{-0.92}^{+1.22}$	$-{0.91}_{-0.85}^{+1.15}$	$-{0.74}_{-0.95}^{+1.09}$
MTOV	${2.13}_{-0.16}^{+0.33}$	${2.13}_{-0.18}^{+0.29}$	${2.14}_{-0.17}^{+0.34}$	${2.12}_{-0.16}^{+0.31}$	${1.46}_{-0.42}^{+0.82}$	${1.81}_{-0.15}^{+0.45}$	${1.85}_{-0.30}^{+0.56}$	${2.09}_{-0.15}^{+0.26}$	${2.11}_{-0.16}^{+0.29}$
${\mathrm{log}}_{10}({\varepsilon }_{c})$	${15.19}_{-0.20}^{+0.21}$	${15.19}_{-0.20}^{+0.20}$	${15.18}_{-0.21}^{+0.21}$	${15.20}_{-0.20}^{+0.19}$	${15.03}_{-0.21}^{+0.33}$	⋯	⋯	⋯
${\mathrm{log}}_{10}({{\rm{P}}}_{c})$	${35.61}_{-0.27}^{+0.30}$	${35.62}_{-0.28}^{+0.30}$	${35.60}_{-0.28}^{+0.30}$	${35.63}_{-0.29}^{+0.30}$	${35.05}_{-0.24}^{+0.30}$	⋯	⋯	⋯
K	1.00	1.04	0.92	1.05	⋯	⋯	⋯	⋯

Note. For the analyses of Section 4, we also show the inferred central energy density ε_c , the corresponding central pressure P_c , and the Bayes' factor K comparing with the model using the chiral EFT band from Hebeler et al. (2013). The first four column results are for the different chiral EFT bands from Hebeler et al. (2013) (Heb 13), Tews et al. (2013) (Tews 13), Lynn et al. (2016) (Lynn 16), and Drischler et al. (2019) (Dri 19), while all other results are for the baseline inference using Heb 13. The column "Combined with" refers to the NICER–XMM analysis of PSR J0740+6620, the NICER analysis of PSR J0030+0451 and multimessenger constraints combined, while in the column "Combined without" the NICER–XMM analysis of PSR J0740+6620 is replaced with just the radio mass measurement by Fonseca et al. (2021). The radii are given in km, M_TOV in M_⊙, and ε_c and P_c in g cm⁻³ and dyn cm⁻², respectively. The upper and lower values correspond to the 95% credible interval.

Download table as:  ASCIITypeset image

5.1. Implications for Nuclear Physics

5.2. Implications for Maximum Mass

An important quantity relating to the EOS is the maximum stable mass of a non-rotating neutron star, M_TOV. Accurate knowledge of M_TOV can aid in classifying compact mergers and merger remnants. In Figure 8 we show posterior distributions on M_TOV when analyzing the updated radio mass measurement of PSR J0740+6620, the joint mass–radius estimate of PSR J0740+6620, and combining GW170817, GW190425, AT2017gfo, PSR J0740+6620, and PSR J0030+0451. The last results in 95% credible ranges for ${M}_{\mathrm{TOV}}={2.23}_{-0.25}^{+0.15}$ M_⊙ and ${M}_{\mathrm{TOV}}={2.11}_{-0.16}^{+0.28}$ M_⊙ for the PP and CS model, respectively. This is in agreement with values previously found (see, e.g., Nathanail et al. 2021 and references therein) when assuming the secondary component in GW190814 was a black hole (Abbott et al. 2020c). Note that the higher end of the distribution in Figure 8 is very dependent on our choice of parameterization, as no information is included from sources with masses above 2.08 M_⊙. One could use information on the merger remnant of GW170817 to put an upper bound on M_TOV (see, e.g., Margalit & Metzger 2017; Shibata et al. 2017; Ruiz et al. 2018), but that is beyond the scope of this Letter. The lower end of the distribution, on the other hand, is strongly correlated with the radio mass measurement of PSR J0740+6620. The recently lowered mass distribution presented in Fonseca et al. (2021) results in slightly lower values for M_TOV compared to the distributions found in Raaijmakers et al. (2020).

Zoom In Zoom Out Reset image size

5.3. Systematic Uncertainties and Framework Comparisons

5.4. Summary and Future Prospects

This work was supported in part by NASA through the NICER mission and the Astrophysics Explorers Program. G.R., T.H., S.N., and T.E.R. are grateful for support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) through the VIDI and Projectruimte grants (PI Nissanke). T.H. also acknowledges support from the NWO sectorplan. T.E.R. and A.L.W. acknowledge support from ERC Consolidator grant No. 865768 AEONS (PI: Watts). S.K.G., K.H., and A.S. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 279384907—SFB 1245. W.C.G.H. acknowledges support through grant 80NSSC20K0278 from NASA. 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. We thank Sharon Morsink and Renee Ludlam for comments on a draft manuscript. We thank Mansi Kasliwal for providing us with the data for AT2017gfo.

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 (Virtanen et al. 2020), MPI (The 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), kalepy (Kelley 2021).