Constraining the equation of state of high-density cold matter using nuclear and astronomical measurements

TL;DR

This paper introduces a Bayesian framework that integrates diverse nuclear and astronomical data to constrain the equation of state of cold dense matter, enhancing our understanding of neutron star properties across various densities.

Contribution

It presents a practical Bayesian method that incorporates full posterior distributions and can be applied with different EOS parameterizations, improving constraints from current and future measurements.

Findings

Different measurements constrain different density ranges of the EOS.

Better symmetry energy data impacts sub-nuclear densities.

Precise radius and tidal deformability measurements improve overall EOS understanding.

Abstract

The increasing richness of data related to cold dense matter, from laboratory experiments to neutron-star observations, requires a framework for constraining the properties of such matter that makes use of all relevant information. Here, we present a rigorous but practical Bayesian approach that can include diverse evidence, such as nuclear data and the inferred masses, radii, tidal deformabilities, moments of inertia, and gravitational binding energies of neutron stars. We emphasize that the full posterior probability distributions of measurements should be used rather than, as is common, imposing a cut on the maximum mass or other quantities. Our method can be used with any parameterization of the equation of state (EOS). We use both a spectral parameterization and a piecewise polytropic parameterization with variable transition densities to illustrate the implications of current measurements and show how future measurements in many domains could improve our understanding of cold catalyzed matter. We find that different types of measurements will play distinct roles in constraining the EOS in different density ranges. For example, better symmetry energy measurements will have a major influence on our understanding of matter somewhat below nuclear saturation density but little influence above that density. In contrast, precise radius measurements or multiple tidal deformability measurements of the quality of those from GW170817 or better will improve our knowledge of the EOS over a broader density range.

Equations

$$\begin{array}[]{rl}{\vec{\alpha}}&\qquad{\rm Equation~{}of~{}state~{}parameters}\\ {\rho_{c,i}}&\qquad{\rm Central~{}density~{}of~{}star}~{}i\\ {\vec{\beta}}_{i}&\qquad{\rm Other~{}parameters~{}fixed~{}for~{}star}~{}i\\ {\vec{\gamma}}_{i,j}&\qquad{\rm Parameters~{}that~{}could~{}vary~{}between~{}measurements}~{}j~{}{\rm of~{}star}~{}i\\ q(\ldots)&\qquad{\rm Prior~{}probability~{}density}\\ P(\ldots)&\qquad{\rm Posterior~{}probability~{}density}\\ {\cal L}(\ldots)&\qquad{\rm Likelihood~{}of~{}the~{}data~{}given~{}a~{}model~{}with~{}associated~{}parameter~{}values.}\\ \end{array}$$

$$P(\vec{\alpha })\propto \int q(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})\mathcal{L}(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})d{\rho }_{c,i}d{\vec{\beta }}_{i}d{\vec{\gamma }}_{i,j}.$$

$$\mathcal{L}(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})=\prod _{i,j}{\mathcal{L}}_{i,j}(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j}),$$

$$q(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})=q(\vec{\alpha })\left[\prod _{i}q({\rho }_{c,i}|\vec{\alpha })\right]\left[\prod _{i}q({\vec{\beta }}_i)\right]\left[\prod _{i,j}q({\vec{\gamma }}_{i,j})\right].$$

$$\begin{array}[]{rl}{\cal L}({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})&=\prod_{i,j}{\cal L}_{i,j}({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})\\ &=\prod_{i}\left[\prod_{j|i}{\cal L}_{i,j}({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})\right]\;.\\ \end{array}$$

$$P(\vec{\alpha })\propto q(\vec{\alpha })\prod _i\left[\int q({\rho }_{c,i}|\vec{\alpha })q({\vec{\beta }}_i)\prod _{j|i}\left(\int q({\vec{\gamma }}_{i,j}){\mathcal{L}}_{i,j}(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})d{\vec{\gamma }}_{i,j}\right)d{\rho }_{c,i}d{\vec{\beta }}_i\right].$$

$$\mathcal{L}(\vec{\alpha })=\prod _i\left[\int q({\rho }_{c,i}|\vec{\alpha })q({\vec{\beta }}_i)\prod _{j|i}\left(\int q({\vec{\gamma }}_{i,j}){\mathcal{L}}_{i,j}(\vec{\alpha },{\rho }_{c,i},{\vec{\beta }}_i,{\vec{\gamma }}_{i,j})d{\vec{\gamma }}_{i,j}\right)d{\rho }_{c,i}d{\vec{\beta }}_i\right].$$

$${\mathcal{L}}_S(\vec{\alpha })=\mathcal{L}(S_0|S(\vec{\alpha })).$$

$${\mathcal{L}}_{M_j}(\vec{\alpha })={\int }_0^{M_{\max }(\vec{\alpha })}P(M_j)dM.$$

$${\mathcal{L}}_I(\vec{\alpha })=\mathcal{L}(I_{\mathrm{obs}}|I(\vec{\alpha },M=M_0)),$$

$${\mathcal{L}}_{E_{\mathrm{bind}}}(\vec{\alpha })=\mathcal{L}(M_{\mathrm{bary},0}|M_{\mathrm{bary}}(\vec{\alpha },M_0)),$$

$$\Lambda =\frac{2}{3}{k}_2{\left(\frac{R{c}^2}{GM}\right)}^5.$$

$$\tilde{\Lambda }=\frac{16}{13}\frac{(M_1+12M_2)M_1^4{\Lambda }_1+(M_2+12M_1)M_2^4{\Lambda }_2}{(M_1+M_2{)}^5}.$$

$${\mathcal{L}}_{\Lambda }(\vec{\alpha })=\int d{M}_1\int q(M_1)q(M_2)\mathcal{L}(M_1,M_2,\tilde{\Lambda }|\vec{\alpha })d{M}_2,$$

$${\mathcal{L}}_{\Lambda }(\vec{\alpha })=\int d{M}_{1}q(M_1)\int q(M_2|M_{\mathrm{ch}},M_1)\mathcal{L}(M_1,M_2,\tilde{\Lambda }|\vec{\alpha })d{M}_{\mathrm{ch}},$$

$${\mathcal{L}}_R(\vec{\alpha })=\int dMq(M){\mathcal{L}}_l(M,R(M,\vec{\alpha })),$$

$$\mathcal{L}(\vec{\alpha })=\left[\prod _i{\mathcal{L}}_{S,i}(\vec{\alpha })\right]\left[\prod _j{\mathcal{L}}_{M_j}(\vec{\alpha })\right]\left[\prod _k{\mathcal{L}}_{\Lambda ,k}(\vec{\alpha })\right]\left[\prod _l{\mathcal{L}}_{R,l}(\vec{\alpha })\right]\left[\prod _m{\mathcal{L}}_{I,m}(\vec{\alpha })\right]\left[\prod _n{\mathcal{L}}_{E_{\mathrm{bind}},n}(\vec{\alpha })\right].$$

$$\Gamma (p)=\exp \left(\sum _k{\gamma }_k{x}^k\right),$$

$${\mathcal{L}}_S(\vec{\alpha })=\frac{1}{(2\pi {\sigma }_S^2{)}^{1/2}}{e}^{-(S(\vec{\alpha })-S_0{)}^2/2{\sigma }_S^2},$$

$$\mathcal{L}(M,R)\propto e^{-(M-1.4M_{\odot }{)}^2/2{\Delta }_M^2}{e}^{-(R-12\mathrm{km}{)}^2/2{\Delta }_R^2},$$

Full text

dense matter, equation of state, neutron star

††software: SuperMongo (https://www.astro.princeton.edu/$\sim$rhl/sm/sm.html)

1 Introduction

Several recent or upcoming astronomical measurements have or will have important implications for our understanding of the properties of the cold, catalyzed matter in the cores of neutron stars. Chief among them are the measurements of the binary tidal deformability from the gravitational-wave event GW170817 (Abbott et al., 2017, 2019b) and the expected measurements of neutron-star radii and masses using the Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016). This information, combined with nuclear data and other astronomical constraints, such as the high measured masses of a few neutron stars (Demorest et al., 2010; Antoniadis et al., 2013; Cromartie et al., 2019), opens up new opportunities to constrain the equation of state (EOS) of cold high-density matter.

Here, we present a rigorous and practical Bayesian procedure that can be used to combine information from different types of nuclear measurements and observations of neutron star systems to constrain the EOS of high-density cold matter. Our procedure can also be used to constrain other properties of neutron stars. For example, data on the cooling of neutron stars could be used to constrain the composition of the interior of neutron stars (see Potekhin et al. 2015 and Wijnands et al. 2017 for recent reviews), providing information that would be complementary to constraints on the EOS. Here, however, we focus only on constraints on the EOS. In Section 2 we discuss our general statistical methodology. In Section 3 we discuss the use of particular types of data, such as the highest measured masses of neutron stars and tidal deformabilities from individual events. In Section 4 we compare our methodology to previous work on constraining the high-density EOS. In Section 5 we discuss our assumed priors and present the results for our illustrative EOS models. We summarize our conclusions in Section 6.

2 Statistical approach

In this section we describe our methodology. We suppose that different types of observations have been made of a set of neutron stars and that we are considering one or more parameterized models of the EOS of neutron-star matter. How should we analyze these observations to correctly use all the available information to obtain estimates of the posterior probability densities of the parameters in these EOS models?

We assume that $n$ neutron stars have been observed and denote a property of star $i$ by a subscript $i$ on the symbol for that property. The observations can be of very different types, e.g., separate measurements of different stars could inform us about their masses, or masses and radii, or moments of inertia, or tidal deformabilities. Our notation is:

[TABLE]

Note that given ${\vec{\alpha}}$ and a rotation rate, $\rho_{c,i}$ determines the mass $M_{i}$ of star $i$. By assumption, the true value of ${\vec{\alpha}}$ is the same for all neutron stars, and the true value of $\rho_{c,i}$ is fixed for a given star (and thus does not vary with the measurement $j$), but can vary from one star to another. Examples of other parameters that are fixed for a given star are the observer inclination and distance to the star; those parameters can, of course, vary from one star to another. The other parameters ${\vec{\gamma}}_{i,j}$ (such as the surface emission pattern during a thermonuclear burst) can vary from one measurement to the next of a single star and can vary from one star to another. The measurements could be of entirely distinct types.

We are interested in the posterior probability density $P({\vec{\alpha}})$. We obtain this by marginalizing the full posterior probability density $P({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})\propto q({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j}){\cal L}({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})$ over the nuisance parameters (i.e., the parameters that do not depend directly on the EOS) $\rho_{c,i}$, ${\vec{\beta}}_{i}$, and ${\vec{\gamma}}_{i,j}$:

[TABLE]

The proportionality in this expression is to remind us that we will need, as a final step, to normalize $P({\vec{\alpha}})$ so that $\int P({\vec{\alpha}})d{\vec{\alpha}}=1$. The likelihood ${\cal L}$ is the product of all of the individual likelihoods, so

[TABLE]

where ${\cal L}_{i,j}({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})$ is the likelihood of measuring data set $j$ from star $i$ given the model under consideration with parameter values ${\vec{\alpha}}$, $\rho_{c,i}$, ${\vec{\beta}}_{i}$, and ${\vec{\gamma}}_{i,j}$.

We make the following two simplifying assumptions:

Assumption 1: the prior $q({\vec{\alpha}},\rho_{c,i},{\vec{\beta}}_{i},{\vec{\gamma}}_{i,j})$ in expression (2) can be represented as the product of the following factors:

[TABLE]

Thus, we assume that the priors are independent of each other, with the exception of the prior on the central density. We write the prior on $\rho_{c,i}$ as $q(\rho_{c,i}|{\vec{\alpha}})$ because it is possible that the prior will depend on other parameters (for example, for a nonrotating star, the maximum central density of a stable star will often depend only on ${\vec{\alpha}}$, but in general, the maximum stable density will also depend on the rotation rate). In principle, other parameters could also be codependent, e.g., if one of our parameters is the rotation frequency, then its maximum value depends on both $\rho_{c,i}$ and ${\vec{\alpha}}$. However, for the cases we consider here, the rotation frequency is small enough that it is unimportant.

Assumption 2: We assume that when we break the overall likelihood into a product of the likelihoods of the individual data sets given the model and parameter values, the parameters not associated with a given observed quantity do not affect the likelihood of the measured value of that quantity. For example, for a given distribution of central densities, we assume that the central density of one star has no influence on the likelihoods of the measured values of the parameters that describe another star. This means we can write

[TABLE]

Here, “$j|i$” means “the set $j$ of measurements of star $i$”.

These assumptions allow us to write the posterior probability density for ${\vec{\alpha}}$ as follows:

[TABLE]

When we compare this with the general expression $P({\vec{\alpha}})\propto q({\vec{\alpha}}){\cal L}({\vec{\alpha}})$ we see that, given our assumptions, the likelihood of the full set of all data given the model and parameter values ${\vec{\alpha}}$ is

[TABLE]

Loosely speaking, this approach assigns the likelihood of each set of values of the measured quantities, given the data, to all combinations of the model parameter values that yield these values of the measured quantities. To see why this is appropriate, note that when we analyze particular neutron-star data, we find that the central density and EOS parameters only influence a subset of the parameters that are used to describe the data. For example, the distance, direction, and orientation of a merging binary do not depend on either $\rho_{c,i}$ or ${\vec{\alpha}}$. Similarly, when energy-dependent X-ray waveforms from NICER are analyzed, only the gravitational mass $M$ and circumferential radius $R$ depend on $\rho_{c,i}$ and ${\vec{\alpha}}$. Thus, in, e.g., the waveform case, the marginalized likelihood associated with given $\rho_{c,i}$ and ${\vec{\alpha}}$ will be the same as the corresponding marginalized likelihood associated with the corresponding $M$ and $R$, where the marginalization is performed over all of the other parameters that describe the particular data set.

Once we have the posterior density at each of a large number of EOS parameter combinations, we compute the posterior density in pressure at a specific density $\rho_{0}$ by (1) determining the pressures predicted at $\rho_{0}$ using each parameter combination, (2) assigning a statistical weight to each pressure that is the same as the posterior density for the parameter combination, and then (3) sorting the predicted pressures at $\rho_{0}$ in increasing order. We then determine a given credibility quantile (e.g., the 5% quantile of the pressure at $\rho_{0}$) by summing the normalized weights of the pressures at $\rho_{0}$ until 5% is reached.

3 Using Different Types of Measurements

Different measurements and observations require different approaches to use them in our statistical procedure for constraining the EOS. Some, such as the nuclear symmetry energy, can be computed directly from the EOS for a broad category of nuclear models. Others, such as the binary tidal deformability measured for GW170817 or future events, require marginalization. We now discuss illustrative measurements and observations and how they can be used in our statistical procedure. As we discussed in Section 2, we can obtain the likelihood from multiple independent measurements by simply multiplying their individual likelihoods. We also note that additional measurements beyond what we consider here, such as measurements of neutron-star quasinormal modes (e.g., Kokkotas & Schmidt 1999), may be available in the future.

3.1 Constraints Not Requiring Marginalization

3.1.1 Nuclear Symmetry Energy

In our discussion here, we assume that the nuclear symmetry energy $S$ is the difference in the energy per nucleon between pure neutron matter (which at density $n$ we denote by $\epsilon(n)/n$) and symmetric nuclear matter, at the nuclear saturation density $n_{s}$. We are interested in the EOS of cold catalyzed matter, which is not purely neutrons. However, as pointed out by Lattimer & Prakash (2016) the proton fraction at $n_{s}$ is only $\sim 1$%, which is small enough to be neglected. With this approximation $S=\epsilon(n_{s})/n_{s}-E_{\rm sym}$, where $E_{\rm sym}=-16.0$ MeV is the energy per nucleon of symmetric matter at nuclear saturation density (Tsang et al., 2012). If the measured value of $S$ is $S_{0}$ and the predicted value for EOS parameters ${\vec{\alpha}}$ is $S({\vec{\alpha}})$, then the likelihood factor associated with the symmetry energy is simply

[TABLE]

3.1.2 Gravitational Mass

We can in principle obtain information about the EOS from any measurement of a neutron-star mass. For example, if an EOS has a maximum mass of $2.5~{}M_{\odot}$ but no neutron stars are found to have masses larger than $2.2~{}M_{\odot}$, that EOS could be disfavored (we thank R. Essick for emphasizing this point to us). However, the mass distribution of neutron stars depends on more than the EOS. For instance, although all equations of state allow $0.5~{}M_{\odot}$ neutron stars to exist, there are no plausible suggestions for how such stars can form. Moreover, the path to forming high-mass neutron stars is not well understood; in the example above, it could be that it is simply extremely rare that a star’s birth and subsequent accretion will produce a mass above $2.2~{}M_{\odot}$ even if significantly higher masses are allowed by the EOS. A complete analysis would take all measured masses into account using a joint, parameterized model of birth and accretion as well as the EOS, but this is not currently feasible. We therefore focus on the highest measured masses.

In the limit of slow rotation, the maximum gravitational mass is a function only of the EOS. It is the gravitational mass $M_{\rm max}$ at the largest central total mass$-$energy density $\rho_{c}$ such that $dM/d\rho_{c}\geq 0$.

If the posterior probability distribution for the mass of star $j$ is $P(M_{j})$, then the likelihood factor for the EOS parameter values ${\vec{\alpha}}$ for that star is

[TABLE]

A similar integral can take into account observations that disfavor large maximum masses, and in Section 5 we show the results for one such hypothetical constraint.

3.1.3 Moment of Inertia

For a given EOS, the expected moment of inertia can be computed given either a central density or a mass (Hartle, 1967). If we assume that we know the mass $M_{0}$ very precisely (as is the case for both components of the double pulsar PSR J0737$-$3039, which is the system of greatest promise for moment-of-inertia measurements), then when a measurement is made of the moment of inertia of the pulsar, the likelihood factor will be

[TABLE]

where ${\cal L}(I_{\rm obs}|I({\vec{\alpha}},M=M_{0}))$ is the likelihood of observing a moment of inertia $I_{\rm obs}$ if the expected value at $M=M_{0}$ is $I({\vec{\alpha}})$ for EOS parameter values ${\vec{\alpha}}$.

3.1.4 Gravitational Binding Energy

Suppose that a star with a precisely measured gravitational mass $M_{0}$ is thought to have a baryonic mass $M_{\rm bary,0}$ (and thus a binding energy $M_{\rm bary,0}c^{2}-M_{0}c^{2}$) with some likelihood ${\cal L}(M_{\rm bary,0}|M_{0})$ (one such possible scenario is if there is evidence that the neutron star was formed in an electron-capture supernova (Nomoto 1984; Podsiadlowski et al. 2004, 2005; Zha et al. 2019; see Section 5 for details and caveats). Then

[TABLE]

where $M_{\rm bary}({\vec{\alpha}},M_{0})$ is the baryonic mass for a gravitational mass $M_{0}$ that is predicted using the EOS with parameter values ${\vec{\alpha}}$.

3.2 Constraints Requiring Marginalization

3.2.1 Binary Tidal Deformability in Neutron-star Mergers

The newest category of EOS-relevant neutron-star observations is the constraint on the tidal deformability of neutron stars that has been obtained using gravitational-wave observations of GW170817 (Abbott et al., 2019b). The dimensionless form of the tidal deformability, for a star of gravitational mass $M$ and circumferential radius $R$, is

[TABLE]

Here, $k_{2}$ is the tidal Love number. Hinderer (2008) has a good discussion of how to compute $\Lambda$ given an EOS and the central density (see also the erratum at Hinderer 2009). Gravitational-wave measurements give a tighter constraint on the binary tidal deformability than on the tidal deformabilities of the two stars individually: indeed, at least for the Taylor family of post-Newtonian waveforms, the most easily measurable quantity for stars of masses $M_{1}$ and $M_{2}\leq M_{1}$ with tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ is (Wade et al., 2014)

[TABLE]

In such events, the masses are not measured well individually, but the chirp mass $M_{\rm ch}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$ is known precisely; for example, for GW170817, $M_{\rm ch}=1.186\pm 0.001~{}M_{\odot}$ (Abbott et al., 2019b). We note that, for fixed $M_{\rm ch}$, ${\tilde{\Lambda}}$ is relatively insensitive to the mass ratio $M_{2}/M_{1}$. For instance, using the scaling $\Lambda\propto M^{-6}$ suggested by De et al. (2018), ${\tilde{\Lambda}}$ for $M_{2}/M_{1}=0.6$ is only $\sim 5\%$ larger than ${\tilde{\Lambda}}$ for $M_{2}/M_{1}=1$.

Because only $M_{\rm ch}$ is measured precisely, we need to marginalize over the masses. We approach this marginalization problem by assuming that gravitational-wave data analysis has given us a full posterior in $(M_{1},M_{2},{\tilde{\Lambda}})$ space. For given EOS parameter values, the prior probability distribution for the masses is set by the prior probability distribution for the central densities (or by equivalent criteria). For fixed EOS parameter values, we can compute ${\tilde{\Lambda}}={\tilde{\Lambda}}(M_{1},M_{2},{\vec{\alpha}})$. Thus, in general, we would compute this likelihood factor by integrating over both $M_{1}$ and $M_{2}$:

[TABLE]

where $q(M_{1})$ and $q(M_{2})$ are the priors for $M_{1}$ and $M_{2}$ and ${\cal L}(M_{1},M_{2},{\tilde{\Lambda}}|{\vec{\alpha}})$ is the three-dimensional likelihood obtained from the analysis of the gravitational-wave data, given EOS parameter values ${\vec{\alpha}}$.

However, $M_{\rm ch}$ is known with such high precision and accuracy that, given a value for $M_{1}$, $M_{2}$ is known to high accuracy. Therefore, we can recast the likelihood factor as

[TABLE]

where $q(M_{2}|M_{\rm ch},M_{1})$ is the prior probability density for $M_{2}$ at the value of $M_{2}$ implied by $M_{\rm ch}$ and $M_{1}$, and the integral is over the probability distribution for $M_{\rm ch}$ obtained from the gravitational-wave analysis. Note that even if $M_{\rm ch}$ is known with high precision, we cannot write the second integral as a delta function. This is because the total probability in the narrow range of $M_{2}$ allowed for a given $M_{\rm ch}$ and $M_{1}$ depends on the EOS and the prior for the masses (or central densities). As a result, this factor must be calculated directly for each EOS.

3.2.2 Radius and Mass

Suppose that for a given star the likelihood of a mass $M$ and radius $R$ is ${\cal L}(M,R)$. For a given stellar mass, the radius $R$ is determined precisely for given EOS parameter values. Thus, the likelihood factor associated with a radius measurement is

[TABLE]

where $R(M,{\vec{\alpha}})$ is the circumferential radius for a gravitational mass $M$ given EOS parameter values ${\vec{\alpha}}$, and $q(M)$ is the prior on $M$. Note that the integration is equivalent to integrating the full $(M,R)$ likelihood over the full $(M,R)$ curve predicted using a given EOS.

3.3 Combination of Constraints

Under the assumption of independent measurements that we described earlier, we can determine the final likelihood ${\cal L}({\vec{\alpha}})$ at a given set of values of the EOS parameters ${\vec{\alpha}}$ by simply setting it equal to the product of the individual likelihoods. Thus, if there is some set $i$ of independent symmetry energy measurements, some set $j$ of neutron-star mass measurements high enough to be constraining (noting that here the use of $i$ and $j$ is different than it was in Section 2), some set $k$ of binary tidal deformability measurements, some set $l$ of mass$-$radius pairs, some set $m$ of moments of inertia, and some set $n$ of gravitational binding energies, then the final likelihood is

[TABLE]

We stress that this expression implicitly assumes that systematic errors can be neglected. If they cannot, then — as always — there is the prospect for significant bias.

4 Comparison with Previous Approaches

In this section, we compare our statistical method with EOS constraint methods in the literature. In Section 5 we will discuss specific inferences of masses, radii, etc. that are then used to constrain the EOS. Here, we focus on the statistical approaches themselves. Our method is generally consistent with other methods that are fully Bayesian, e.g., among recent papers Lackey & Wade (2015), Agathos et al. (2015), Alvarez-Castillo et al. (2016), and Riley et al. (2018). The non-parameteric approach of Landry & Essick (2019) is also worth consideration.

4.1 Use of Bounds in Mass or Other Quantities

As we have emphasized, in a fully consistent Bayesian analysis, a given observation needs to be incorporated using a likelihood-based procedure. Imposing a strict bound of any kind, other than bounds stemming from fundamental physical laws, may discard important information. However, to our knowledge, all previous analyses except that of Alvarez-Castillo et al. (2016) have used a hard lower bound on the maximum mass, in the sense that a given EOS or parameter combination is allowed if it has a maximum mass above some specified value (often $1.97~{}M_{\odot}$, because the $M=2.01\pm 0.04~{}M_{\odot}$ mass reported by Antoniadis et al. 2013 for PSR J0347+0432 was the highest reported mass until the $M=2.14^{+0.10}_{-0.09}~{}M_{\odot}$ mass reported by Cromartie et al. 2019 for PSR J0740+6620), and disallowed if the maximum mass is below the bound. A similar approach is taken commonly, but not as universally, with the tidal deformability measurement from GW170817.

The first reason that this is incorrect is illustrated nicely by the progression in time of the estimates of the mass of PSR 1614$-$2230. The first measurement, by Demorest et al. (2010), was $M=1.97\pm 0.04~{}M_{\odot}$. The second measurement, by Fonseca et al. (2016), was $M=1.928\pm 0.017~{}M_{\odot}$. The most recent measurement, by Arzoumanian et al. (2018), is $M=1.908\pm 0.016~{}M_{\odot}$. Thus, the best estimate of the mass in both updates is slightly more than one standard deviation lower than the previous best estimate. Thus, a strict lower bound at the $-1\sigma$ mass $M=1.93~{}M_{\odot}$ from the first measurement would be too restrictive given our current knowledge of the mass of PSR 1614$-$2230. Instead, one should use the full posterior distribution of the mass.

The second reason why this approach is suboptimal is that there is, after all, uncertainty in the mass measurements. If we accept $M=2.14^{+0.10}_{-0.09}~{}M_{\odot}$ as the mass estimate for PSR J0740+6620, then using the hard-bound approach, an EOS with a maximum mass of $2.05~{}M_{\odot}$ is just as viable as an EOS with a maximum mass of $2.14~{}M_{\odot}$. But if we assume that the measurement has only Gaussian statistical uncertainties, there is an $\sim$84% probability that the mass of PSR J0740+6620 is greater than $2.05~{}M_{\odot}$. Thus, in reality, the EOS with $M_{\rm max}=2.14~{}M_{\odot}$ is considerably more consistent with the data than the EOS with $M_{\rm max}=2.05~{}M_{\odot}$. Applying a lower bound is not a statistically appropriate approach.

The third reason that strict bounds should not be used is that this approach does not allow the incorporation of information from multiple stars. For example, at the moment, the only published masses that pose significant constraints to the EOS are $M=2.14^{+0.10}_{-0.09}~{}M_{\odot}$ for PSR J0740+6620, $M=2.01\pm 0.04~{}M_{\odot}$ for PSR J0348+0432, and $M=1.908\pm 0.016~{}M_{\odot}$ for PSR 1614$-$2230. An EOS with $M_{\rm max}=1.8~{}M_{\odot}$ is disfavored at the $3.8\sigma$ level for PSR J0740+6620 alone, but at more than $9\sigma$ when measurements of the masses of all three pulsars are included (using the simple assumption that the uncertainties are exactly Gaussian, which is unlikely to be true at several standard deviations). Thus, $M_{\rm max}=1.8~{}M_{\odot}$ is excluded much more strongly based on the data from all three stars than it would be using just the most massive of the three. If a future star is discovered with, say, a mass $M=2.01\pm 0.05~{}M_{\odot}$, then using the hard-bound method, it would not contribute at all to EOS constraints, whereas in reality, it would make low-$M_{\rm max}$ EOSs significantly less probable.

4.2 Lack of Marginalization

It is common, although not universal, for post-GW170817 EOS constraint papers to use an estimate of the tidal deformability parameter at $1.4~{}M_{\odot}$ in constraints, rather than integrating over the full posterior space. Similarly, numerous papers use only the maximum-likelihood or minimum-$\chi^{2}$ point along the $R(M)$ curve implied by a given EOS, whereas the integration should instead be performed over the whole curve (see, for example, Steiner et al. 2010; Özel et al. 2016).

4.3 Attempts to Invert Measurements to Obtain the EOS

Early papers on inference of the EOS from neutron-star measurements often presented EOS determination as an inversion of neutron-star measurements, sometimes using a Jacobian formalism to map neutron-star observables into EOS parameters. Such an approach misses the fact that this is intrinsically a measurement problem, not a problem of inverting a mathematical relation, and thus must be approached statistically. Not approaching the analysis as a measurement problem can lead to fundamental difficulties.

Even setting aside for the moment the fundamentally statistical nature of the problem, in realistic situations, attempts to invert observed quantities to determine EOS parameter values fail because the inversion is singular. For example, if two $M(R)$ curves obtained from different equations of state cross, then the inversion is clearly singular at the crossing point. Another difficulty with this approach has been emphasized by Riley et al. (2018) and Raaijmakers et al. (2018), in the context of EOS models that have separately parameterized segments at different densities, such as models that use a sequence of polytropes. They point out that some neutron stars might not have a central density large enough to reach the highest density in the EOS model. In that case, the parameters describing higher densities have no influence on the mass and radius of that star, and thus nothing can be inferred about those parameters (Raaijmakers et al., 2018). A further difficulty with approaching EOS parameter estimation as a mathematical inversion problem rather than as statistical inference is that a one-to-one mapping requires that the number of EOS parameters be equal to the number of observables. Of course, the hope is that there are many more observations than model parameters!

For these reasons, most papers in the last decade have approached this problem correctly, as a statistical inference problem, rather than as a problem of inverting a map between observations and model parameters.

5 Results

In this section, we present the 5%, 50%, and 95% credibility quantiles for the pressure at a set of densities and for the circumferential radii at a set of gravitational masses, obtained using progressively more restrictive data. The densities start at half of nuclear saturation density (i.e., at 0.08 baryons per fm3), where the pressures of all of our EOS models agree by construction, because up to that density we use the SLy (Douchin & Haensel, 2001) EOS. We then construct the cumulative probability distribution for the pressure at progressively higher densities. We also plot the $M-R$ curves that bound the region that makes up 90% of the total probability. Currently, constraints on the EOS are relatively weak, which means that most EOS parameter combinations have high likelihoods, and thus, we do not need to perform sophisticated searches through parameter space.

Our method can be used with any parameterization of the EOS. We assume that the pressure is a function only of the density, i.e., that the EOS is barotropic. The pressure does not depend explicitly on the temperature or the proton fraction, because we assume that the matter is in beta equilibrium. For our primary parameterization we follow Abbott et al. (2018) in using the spectral parameterization introduced by Lindblom (Lindblom, 2010, 2018), in which the free parameters are spectral indices $\gamma_{k}$ that represent the adiabatic index $\Gamma(p)=[(\rho+p)/p](dp/d\rho)$ (where $p$ is the pressure and $\rho$ is the total mass$-$energy density) using the expansion

[TABLE]

where $x\equiv\log(p/p_{0})$ and $p_{0}$ is the pressure at half of nuclear saturation density. We also follow previous work (e.g., Abbott et al. 2018; Carney et al. 2018) by using an expansion up to $x^{3}$ with the following uniform priors on the coefficients $\gamma_{k}$: $\gamma_{0}\in[0.2,2]$, $\gamma_{1}\in[-1.6,1.7]$, $\gamma_{2}\in[-0.6,0.6]$, and $\gamma_{3}\in[-0.02,0.02]$. We do not additionally require, as some papers have, that $\Gamma(p)\in[0.6,4.5]$ at all densities. The parameterization itself guarantees that $\Gamma(p)>0$, which is needed to enforce thermodynamical stability. We also require that the adiabatic speed of sound be less than the speed of light. In Section 5.7 we display results using an alternative parameterization which has potentially different polytropic indices at variable transition densities.

Once the EOS is chosen, then in the slow rotation limit, the mass and radius as functions of the central density, the maximum stable mass, and the gravitational binding energy for a given gravitational mass follow from the Tolman-Oppenheimer-Volkoff (TOV) equation (Tolman, 1939; Oppenheimer & Volkoff, 1939), and from the relation between baryonic mass density and total mass$-$energy density discussed in Tooper (1965). We compute the moment of inertia and spin quadrupole moment following the approach in Hartle (1967), and the tidal Love number using the development in Hinderer (2008) (see also the erratum at Hinderer 2009). We verified the accuracy of our code by comparing our outputs with those listed in Table III of Read et al. (2009) (using their equation of state rather than the spectral parameterization). We also checked that our moments of inertia, quadrupole moments, and tidal deformabilities follow closely the I-Love-Q relations (Yagi & Yunes 2013 and subsequent papers).

The order in which we add measurements is (1) symmetry energy (from laboratory measurements), (2) mass measurements, (3) tidal deformability measurements, (4) hypothetical future measurements of both radius and mass, (5) hypothetical future measurements of moments of inertia, and finally (6) hypothetical future measurements of the binding energy of stars with precisely measured gravitational masses. That is, in the first section we present results assuming only measurements of (1) (with different illustrative levels of precision for the symmetry energy). We then present results assuming only measurements of (1) and (2) (with a standard precision for the symmetry energy and different potential measurements for the mass), and so on. This makes it possible to see how additional measurements progressively improve the precision of our understanding of the EOS and, as a consequence, the neutron-star mass$-$radius relation. Note that when a new measurement is incorporated, the new EOS constraints can shift beyond the previous 5% or 95% quantile. For example, if a neutron star is measured to have a high mass then soft equations of state are disfavored, which then shifts the quantiles to higher pressure at a given density.

Whereas in Section 2 we presented our general statistical method, and in Section 3 we discussed how to apply our method to particular types of measurements, here we use both existing and potential future measurements to find credibility regions in $P(\rho)$ space. Thus we need to make choices about which measurements to use. For example, after GW170817, there have been many detailed simulations and comparisons with electromagnetic information (especially the details of the resulting kilonova) that have endeavored to constrain the maximum mass of neutron stars, or to place lower limits on the tidal deformability of neutron stars of particular masses. We also need to specify the prior on the mass or the central density for a given combination of EOS parameter values. In the results we present here, we assume that the central density can with equal probability be anywhere between the density that would produce an $M=1.0~{}M_{\odot}$ neutron star with that EOS, and the density that produces the maximum mass possible for that EOS. We stress that although we make particular choices, these are only illustrative. Our focus is not to produce our own version of the constraints, although given our assumptions, our current constraints are in the top left panels of Figure 5 and Figure 6. Instead, we make these choices to demonstrate how our method works in practice; other choices of measurements and even of EOS families and the priors on their parameters could be used straightforwardly with our method.

Our final note prior to presenting our results is a reminder that all measurements and observations have to be interpreted within a model framework, and this means that we rely on that framework to obtain quantities of interest. For example, virtually all neutron-star observations are interpreted under the assumption that general relativity properly describes extreme gravity. Many papers prior to the direct detection of gravitational waves pointed out that the mass$-$radius relation (and thus all other structural aspects of stars) could be modified considerably in different theories of gravity (see DeDeo & Psaltis 2003 and Orellana et al. 2013 for just two examples). Careful analysis of gravitational-wave data has limited the prospects for deviations from general relativity in stellar-mass objects (see Yunes et al. 2016 for an excellent summary after the first two events), but it is useful to keep an open mind.

5.1 Nuclear Symmetry Energy

Tsang et al. (2012) give the status of a number of different laboratory measurements that could constrain the nuclear symmetry energy. We treat the likelihood factor from the symmetry energy as a Gaussian:

[TABLE]

where $S({\vec{\alpha}})$ is the symmetry energy predicted using specified values of the EOS parameters ${\vec{\alpha}}$. For our standard constraint, we choose $S_{0}=32$ MeV and $\sigma_{S}=2$ MeV from a rough averaging of the various results presented in Tsang et al. (2012). Non-Gaussian likelihoods are also straightforward to include in our framework.

Figure 1 shows that more precise measurements of $S$ would strongly constrain the EOS below nuclear saturation density but would have little effect above $n_{s}$. Figure 2 shows that knowledge of $S$ has little impact on our knowledge of the radius of stars with $M>1.0~{}M_{\odot}$.

5.2 Maximum Mass

A viable EOS must be able to support a maximum gravitational mass $M$ that is at least as great as the largest reliably measured neutron-star mass. For masses, the gold standard is neutron stars in relativistic binaries, for which it is possible to measure post-Keplerian parameters such as the Shapiro delay, pericenter precession, and orbital decay due to the emission of gravitational radiation (see Freire 2009 for a good discussion of how these parameters are measured and the governing equations). The precision with which these masses can be measured, plus the reliability of the underlying theory, makes inferred masses the bedrock of astronomical constraints on the EOS of cold high-density matter. Particularly notable are the mass measurements $M=1.908\pm 0.016~{}M_{\odot}$ for PSR J1614$-$2230 (original mass measurement in Demorest et al. 2010 and current mass measurement in Arzoumanian et al. 2018), $M=2.01\pm 0.04~{}M_{\odot}$ for PSR J0348+0432 (Antoniadis et al., 2013), and $M=2.14^{+0.10}_{-0.09}~{}M_{\odot}$ for PSR J0740+6620 (Cromartie et al., 2019).

There are intriguing suggestions of even higher-mass neutron stars. For example, the “black widow” system PSR B1757+20 has an estimated mass of $2.40\pm 0.14~{}M_{\odot}$ (van Kerkwijk et al., 2011), and another black widow system, PSR 1311$-$3430, has an estimated mass of $2.68\pm 0.14~{}M_{\odot}$ (Romani et al., 2012). However, these measurements are less reliable than the relativistic binary masses because of potential systematic errors and the residuals in the fits (van Kerkwijk et al., 2011; Romani et al., 2012).

There are also arguments based on short gamma-ray bursts (Bauswein et al., 2013; Fryer et al., 2015; Lawrence et al., 2015) that were later applied to the double-neutron-star coalescence event GW170817 (Margalit & Metzger, 2017), which suggest a relatively low maximum mass. For example, Margalit & Metzger (2017) suggest that if the two neutron stars in GW170817 formed a hypermassive neutron star that collapsed within tens or hundreds of milliseconds to a black hole, then M_{\rm max}\mathrel{\raise 1.29167pt\hbox{<}\mkern-14.0mu\lower 2.58334pt\hbox{\sim}}2.17~{}M_{\odot}, which is precisely consistent with the predictions of Fryer et al. (2015) and Lawrence et al. (2015). However, there is no direct evidence that there was a collapse to a black hole. Similarly, there are various model-dependent upper limits on $M_{\rm max}$ that have been obtained via comparison of simulations with the kilonova that followed GW170817 (e.g., Shibata et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Coughlin et al. 2019).

We adopt, as our standard maximum mass constraint, the combination of three mass measurements: $M=2.14^{+0.10}_{-0.09}~{}M_{\odot}$ for PSR J0740+6620, $M=2.01\pm 0.04~{}M_{\odot}$ for PSR J0348+0432 and $M=1.908\pm 0.016~{}M_{\odot}$ for PSR J1614$-$2230. We also explore the constraints we would obtain if there is a future mass measurement of $M=2.3\pm 0.1~{}M_{\odot}$, or a future mass measurement of $M=2.4\pm 0.1~{}M_{\odot}$, or a confirmed upper limit of $M_{\rm max}=2.2\pm 0.05~{}M_{\odot}$. In all cases, we assume that the masses or mass limits have Gaussian likelihoods.

Figure 3 shows the second level of constraints, in which we consider that the probability distribution of the symmetry energy is a Gaussian with a mean of 32 MeV and a standard deviation of 2 MeV, and add the mass constraints described above. It is clear from the figure that such measurements place important constraints on the high-density EOS. Likewise, Figure 4 shows the resulting mass$-$radius constraints. We see that, as expected, lower limits on the maximum mass push radii to higher values, whereas upper limits push them to lower values.

5.3 Tidal Deformability

The limits on ${\tilde{\Lambda}}$ from an event such as GW170817 depend on the waveform model, with a spread of $\sim 10$% among models used thus far (Abbott et al., 2019b). Bearing this caveat in mind, the middle 90% of the posterior credible range for ${\tilde{\Lambda}}$ has been reported as (70,720; Abbott et al. 2019b). Future improvements in gravitational-wave sensitivity, plus the simple accumulation of observing time, are expected to yield a rapidly growing number of detected double-neutron-star coalescences, and potentially a few mergers between neutron stars and black holes. These additional observations will improve the constraints on the tidal deformability, especially given the anticipated improvements in high-frequency sensitivity due to the use of squeezed light. It is, however, worth tempering expectations for two reasons: (1) although tidal effects will be more pronounced at higher frequencies and thus constraints could in principle be improved substantially, waveform families also diverge more at higher frequencies and thus the role of systematic errors will be more prominent, and (2) GW170817 was an exceptionally strong event (its signal to noise was the largest of any event in the first two LIGO runs; see Abbott et al. 2019a), which means that future events are likely to be measured less precisely.

For GW170817, the full posterior over all model parameters is available at

https://dcc.ligo.org/LIGO-P1800061/public. We use the $\sim 4000$ samples at this site as input for a kernel density estimate (see Rosenblatt 1956; Parzen 1962; Silverman 1986 for details) of the marginalized posterior for the primary mass and binary tidal deformability, which we use in our estimates of the constraints that we display in Figure 5 and Figure 6. Here, we add to our standard $S+M_{\rm max}$ constraints information from tidal deformability measurements. We begin with the single event GW170817, and then suppose that we have a succession of identical events. From these figures, it is clear that precise tidal deformability measurements will contribute substantially to our understanding of the dense matter EOS, and to our knowledge of the radius at a wide range of masses. We also note that various groups have modeled the electromagnetic emission and have proposed other limits on ${\tilde{\Lambda}}$ (e.g., Radice et al. 2018 find a lower limit ${\tilde{\Lambda}}>400$ for GW170817), but we have not included such limits in our analysis (see Kiuchi et al. 2019 for cautionary remarks about lower limits to ${\tilde{\Lambda}}$ obtained in this manner).

Thus far, we have used existing measurements, plus plausible extrapolations. We will now explore the effect of adding additional types of constraints that could be obtained in the future.

5.4 Radius Measurements

Reliable and precise radius measurements would be extremely useful in constraining the properties of high-density matter, and much effort has been devoted to the analysis of, in particular, X-ray data from isolated and bursting neutron stars. However, there are potentially large systematic errors in current reports of neutron-star radii; for detailed discussions, see Miller (2013) and Miller & Lamb (2016), and see additional caveats related to our uncertainty about the EOS of the crust in Gamba et al. (2019).

There is optimism that systematic errors might not be significant for the results that will be obtained using Neutron Star Interior Composition Explorer (NICER) measurements of the X-ray pulse waveforms of a few non-accreting neutron stars that are pulsars. This optimism is based on studies that have been performed of the method, which involves fitting the energy-resolved X-ray waveforms to models with thermally emitting spots that rotate with the neutron star. Lo et al. (2013) and Miller & Lamb (2015) generated synthetic waveforms using various geometries and assumptions, and fit them with standard models that had uniformly emitting circular spots. Although in many cases, the generated spots were oval, or had temperature gradients, or had spectra or beaming patterns different from those assumed in the fitted model, in no case was there a statistically good fit that was significantly biased in mass or radius. This stands in strong contrast to alternative methods, for which an apparently excellent fit with large bias is possible or even likely, meaning that the fit quality alone does not give a hint that there are potential problems.

Thus, our opinion is that although current radius measurements may have significant systematic errors, future NICER measurements are promising. In addition, as was pointed out by Annala et al. (2018) (see also De et al. 2018, Raithel et al. 2018, and other papers), gravitational-wave measurements from double-neutron-star mergers can place limits on neutron-star radii, but because these are not independent from tidal deformability estimates, we have not included them separately in our constraints.

In this section, we suppose that a posterior in $(M,R)$ has been obtained for a given star. The posterior need not be a product of independent posteriors in $M$ and $R$, or independent posteriors in $M/R$ and $M$; the correlations, if any, depend on the details of the system (see Lo et al. 2013; Miller & Lamb 2015). For the purposes of illustration only, we suppose here that the posterior in mass and radius is a product of independent Gaussians:

[TABLE]

where we explore the consequences of selecting $\Delta_{M}$ and $\Delta_{R}$ equal to 20%, 10%, 5%, and 2% of the best values of the mass and radius, respectively.

In Figure 7 we show the effect of adding radius plus mass measurements as described in Equation (20). As can be seen in Figure 7, a fractional precision of \mathrel{\raise 1.29167pt\hbox{<}\mkern-14.0mu\lower 2.58334pt\hbox{\sim}}5% for a single star is necessary to add significantly to our information (although as pointed out by Miller 2016 and Weih et al. 2019, in some mass ranges, such as $M>2~{}M_{\odot}$, less-precise measurements could still be important). As a check of our method, we find in Figure 8, as we must, that improved measurement precision of mass and radius will dramatically tighten the mass$-$radius relation.

5.5 Moment of Inertia for Neutron Star of Known Mass

Shortly after the discovery of the double pulsar PSR J0737$-$3039 (Burgay et al., 2003), it was pointed out (Lattimer & Schutz, 2005; Kramer & Wex, 2009) that in principle, spin-orbit coupling could be measured within a few years from the resulting extra pericenter precession, and that this might yield an interestingly precise moment of inertia for the more rapidly rotating of the two pulsars, PSR J0737$-$3039A (which has a mass of $M=1.3381\pm 0.0007~{}M_{\odot}$: Kramer et al. 2006). The measurement has been far more challenging than originally envisioned, but there is still hope that within about a decade, the moment of inertia can be measured to within $\sim 10$%. For our illustrative constraint, we select $I_{1.338}=1.37\times 10^{45}$ g cm2, with Gaussian uncertainties, because this is consistent with the other real and hypothetical constraints we are applying and because it is consistent with the moment-of-inertia range found by Landry & Kumar (2018).

We add our second hypothetical constraint in Figure 9: moment-of-inertia measurements for the $M=1.338~{}M_{\odot}$ neutron star PSR J0737$-$3039A. Measurements to the hoped-for precision of $\sim 10$% for this single star would add significantly to the constraints, but less-precise measurements would have little effect. In Figure 10 we see that better moment-of-inertia measurements for an $M=1.338~{}M_{\odot}$ star would have little influence on estimates of the radii of stars with masses around $M=1.0~{}M_{\odot}$, but would significantly improve the estimates of the radii of $M=1.8-2~{}M_{\odot}$ stars.

5.6 Binding Energy of Neutron Stars Formed in Electron-capture Supernovae

If it were possible to know the baryonic rest mass (that is, the sum of the masses of all of the constituent particles if separated to large distance at zero speed) as well as the gravitational mass, for individual neutron stars to reasonable precision, then the resulting knowledge of the binding energy for those stars would provide another constraint on the EOS. It is not possible to make a direct measurement of the baryonic rest mass of a star, but there are suggestions that a particular type of core-collapse supernova known as an electron-capture supernova might occur when the core baryonic rest mass is in the narrow range $M_{\rm bary}\sim 1.36-1.37~{}M_{\odot}$ (Nomoto, 1984; Podsiadlowski et al., 2004, 2005; Zha et al., 2019). If there is then neither expulsion of mass nor additional fallback, and if neutron stars formed via this mechanism can be identified and their gravitational masses measured, then a constraint could be applied. There are clearly several ways in which this identification, or the estimate of the baryonic rest mass, could fail. Moreover, some objects likely to be neutron stars are too light to have formed from an electron-capture supernova, e.g., the $M=1.174\pm 0.004~{}M_{\odot}$ companion to PSR J0453+1159 Martinez et al. (2015), so other mechanisms to produce low-mass neutron stars (such as ultra-stripped supernovae; see Tauris et al. 2017) could be in play in this mass range. Notwithstanding those caveats, Podsiadlowski et al. (2005) made the interesting suggestion that the second pulsar in the double pulsar system, PSR J0737$-$3039B, originated from an electron-capture supernova, and that its gravitational mass of $M=1.2489\pm 0.0007~{}M_{\odot}$ should therefore be identified with $M_{\rm bary}=1.366-1.375~{}M_{\odot}$.

Thus, when we incorporate this hypothetical factor into our analysis, we do so by assuming that the baryonic mass corresponding to a gravitational mass $M=1.2489~{}M_{\odot}$ is $1.37~{}M_{\odot}$ with a Gaussian likelihood.

In Figures 11 and 12 we show the effect of adding this constraint. A fractional uncertainty \mathrel{\raise 1.29167pt\hbox{<}\mkern-14.0mu\lower 2.58334pt\hbox{\sim}}0.5% would improve our knowledge of the EOS below $\sim{\rm few}\times n_{s}$, and would also tighten the range of radii of low-mass neutron stars.

5.7 Selected Results for a Piecewise Polytropic Parameterization of the Equation of State

One drawback of the spectral parameterization we use, with the priors we adopt, is that this does not allow the pressure to be nearly constant over a range of densities. That is, this parameterization is poor at reproducing phase transitions. Although we again stress that the main point of our paper is our Bayesian inference framework rather than specific results, we present for comparison results for a subset of the measurements presented above, for a different equation of state.

In this parameterization, we again enforce causality ($dP/d\rho<c^{2}$) and stability ($dP/d\rho>0$) and use the SLy EOS (Douchin & Haensel, 2001) up to less than half of the nuclear saturation density $\rho_{s}$. Above $\rho_{0}=\rho_{s}/2$ we represent the EOS by a sequence of polytropes with indices that can change at transition densities that are also parameters: our priors are $\rho_{1}\in[3/4,5/4]\rho_{s}$, $\rho_{2}\in[3/2,5/2]\rho_{s}$, $\rho_{3}\in[3,5]\rho_{s}$, and $\rho_{4}\in[6,10]\rho_{s}$. Our priors on the polytropic indices are $\Gamma_{1}\in[2,3]$ from $\rho_{0}$ to $\rho_{1}$, $\Gamma_{2}\in[0.1,5]$ from $\rho_{1}$ to $\rho_{2}$, $\Gamma_{3}\in[0.1,5]$ from $\rho_{2}$ to $\rho_{3}$, $\Gamma_{4}\in[0.1,5]$ from $\rho_{3}$ to $\rho_{4}$, and $\Gamma_{5}\in[0.1,5]$ for densities higher than $\rho_{4}$. All priors are flat in the bracketed range. The limited range $[2,3]$ for $\Gamma_{1}$ is informed by the study of Hebeler et al. (2013).

The results of using this parameterization with a subset of our measurements are shown in Figure 13 and Figure 14. Here, the progressive measurements are $S=32\pm 2$ MeV; the measured masses of PSR J0740+6620, PSR J0348+0432, and PSR 1614$-$2230; the tidal deformability from GW170817; a hypothetical measurement of $(M,R)=(1.4~{}M_{\odot},12~{}{\rm km})$ with 5% precision; a hypothetical measurement of the moment of inertia of an $M=1.338~{}M_{\odot}$ star with $\Delta I_{45}=0.1$; and hypothetical knowledge to within $0.005~{}M_{\odot}$ of the baryonic rest mass of a star with a gravitational mass of $M=1.2489~{}M_{\odot}$. We see that although details of the resulting constraints are somewhat different than for the spectral parameterization, the trends are similar.

6 Conclusions

We have shown that diverse sources of both laboratory and astronomical information about cold, dense, catalyzed matter can be incorporated flexibly within a straightforward, rigorous, and practical Bayesian framework. We treat carefully the constraints that stem from the existing measurements of the symmetry energy, large neutron-star masses, and tidal deformability, the expected future measurements of neutron-star radii and masses, and the possible future measurements of the moments of inertia and gravitational binding energies of neutron stars. We find that different types of measurements will play significantly different roles in constraining the EOS in different density ranges. For example, better symmetry energy measurements will have a major influence on our understanding of matter somewhat below nuclear saturation density but little influence above that density. In contrast, precise radius measurements or multiple tidal deformability measurements of the quality of those from GW170817 or better will improve our knowledge of the equation of state over a broader density range. Of course, any of these analyses would have to be revisited if systematic errors dominate; but overall, the prospects are good in the next few years for a dramatically enhanced understanding of the nature of dense matter.

We thank Chris Pethick for advice on the nuclear physics and Paulo Bedaque, David Blaschke, Reed Essick, Will Farr, Koutarou Kyutoku, Phillipe Landry, and Victor Santos Guedes for other illuminating discussions. We also thank the anonymous referee for an unusually constructive report. M.C.M. and C.C. thank the Kavli Institute for Theoretical Physics for its hospitality during the completion of this paper. M.C.M. and C.C. are grateful for the hospitality of Perimeter Institute where part of this work was carried out. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. CC was also supported in part by the Simons Foundation through the Simons Foundation Emmy Noether Fellows Program at Perimeter Institute. This work was partially supported by the Brazilian National Council for Scientific and Technological Development (CMPq). This research was also supported in part by the National Science Foundation under grant No. NSF PHY-1748958. The authors acknowledge the use of NASA’s Astrophysics Data System (ADS) Bibliographic Services and the arXiv.