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

###### 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. The method allows any parametrization of the equation of state to be used. We use a spectral parametrization to illustrate the implications of current measurements and show how future measurements in many domains could improve our understanding of cold catalyzed matter. In particular, we find that measurements of the masses of massive neutron stars such as PSR J0740+6620 (Cromartie et al., 2019) and multiple, high-precision measurements of tidal deformabilities will provide significant information about the equation of state of matter above nuclear saturation density.

## 1 Introduction

There are several recent or upcoming astronomical observations with important implications for the properties of the cold, catalyzed matter in the cores of neutron stars. Chief among those observations are the constraints on the binary tidal deformability from the gravitational wave event GW170817 (Abbott et al., 2017, 2019) and the expected measurements of neutron star radii and masses using data from the Neutron star Interior Composition ExploreR (NICER; Gendreau et al. 2016). This information, combined with nuclear data and 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 formalism to combine EOS information from disparate sources and measurements, as well as the results obtained by using different combinations of existing and potential future measurements. We note that our statistical approach could also incorporate other types of data, for example information about the cooling of neutron stars (see Potekhin et al. 2015 and Wijnands et al. 2017 for recent reviews), but here we focus only on constraints on the equation of state. In Section 2 we discuss our general statistical methodology. In Section 3 we discuss the incorporation 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 give our results, and our conclusions are in Section 6.

## 2 Statistical approach

Suppose that different types of measurements have been made of multiple neutron stars. If we have a parameterized model of the EOS of cold catalyzed matter, how should we analyze the data to provide an estimate of the probability density of those EOS parameters? In this section we focus on our general methodology. We therefore postpone discussion of particular types of constraints to Section 3, comparison with previous work to Section 4, and specific use of existing and future data to Section 5. We give the details of our illustrative EOS model and priors in Section 5.

Suppose that some set of neutron stars has been observed. 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:

(1) |

By assumption, the true value of is the same for all neutron stars, and the true value of is fixed for a given star (and thus does not vary with the measurement ), 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 can vary from one measurement to the next of a single star (such as the surface emission pattern during a thermonuclear burst) and can certainly vary from one star to another. The measurements could be of entirely distinct types.

We are interested in the posterior probability density . We obtain this by marginalizing the full posterior probability density over the nuisance parameters (i.e., the parameters that do not depend directly on the equation of state) , , and :

(2) |

The proportionality in this expression is to remind us that we will need, as a final step, to normalize so that . The likelihood is the product of all of the individual likelihoods, so

(3) |

We make the following two simplifying assumptions:

Assumption 1: the prior in expression (2) can be represented as the product of the following factors:

(4) |

Here we write the prior on as because it is possible that the prior will depend on other parameters (for example, the maximum central density of a stable star will in general depend on ).

Assumption 2: We assume that when we break the overall likelihood into a product of the likelihoods of the individual observations, the parameters not associated with a given observation do not contribute to the likelihood of that observation. For example, we assume that the central density of one star has no influence on the likelihood of the data taken from another star. This means we can write

(5) |

Here “” means “the set of measurements of star ”.

These assumptions allow us to write the posterior probability density for as follows:

(6) |

When we compare this with the general expression we see that given our assumptions the likelihood of the full set of all data given the model and parameter values is

(7) |

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 equation of state 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 or . Similarly, when energy-dependent X-ray waveforms from NICER are analyzed, only the gravitational mass and circumferential radius depend on and . Thus in, e.g., the waveform case the marginalized likelihood associated with given and will be the same as the corresponding marginalized likelihood associated with the corresponding and , 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 by (1) determining the pressures predicted at 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 in increasing order. We then determine a given credibility quantile (e.g., the 5% quantile of the pressure at ) by summing the normalized weights of the pressures at until 5% is reached.

## 3 Incorporation of Specific Constraints

Different measurements and observations require different approaches to incorporate them into our overall constraints. Some, such as the nuclear symmetry energy, can be obtained directly from the equation of state in a broad category of nuclear models. Others, such as the binary tidal deformability measured from GW170817 or future events, require marginalization. Here we discuss each of the constraints in turn and how we use them. As we discussed in Section 2, we can obtain the likelihood from multiple independent constraints by simply multiplying the individual likelihoods.

### 3.1 Constraints not requiring marginalization

Nuclear symmetry energy.—For the discussion here we assume that the nuclear symmetry energy is the difference in energy per nucleon between pure neutron matter and symmetric nuclear matter, at nuclear saturation density . Our EOS is for catalyzed matter at each density, which means that the energy density we have is for matter that is not pure neutrons, but as pointed out by Lattimer & Prakash (2016) the proton fraction at is only %, which is a small enough correction to be neglected. Thus under these assumptions , where MeV (Tsang et al., 2012). If the measured value of is and the predicted value for EOS parameters is , then the likelihood factor associated with the symmetry energy is simply

(8) |

Maximum gravitational mass.—The maximum gravitational mass is a function only of the equation of state in the limit of slow rotation. It is the gravitational mass evaluated at the largest central total mass-energy density such that .

If the posterior probability distribution for the mass of star is , then the likelihood factor for EOS parameters for that star is

(9) |

We could write a similar integral to incorporate observations that disfavor large maximum masses, and in Section 5 we show the results for one such hypothetical constraint.

Moment of inertia.—For a given equation of state, the expected moment of inertia can be computed given a central density or a mass. If we assume we know the mass well (as we do for both components of the double pulsar PSR J07373039, 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

(10) |

where is the probability of observing a moment of inertia if the expected value at is for EOS parameters .

Gravitational binding energy.—Suppose that a star with a precisely measured gravitational mass is believed to have a baryonic mass (and thus a binding energy ) with some likelihood (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); see Section 5 for details and caveats). Then

(11) |

where is the baryonic mass for a gravitational mass that is predicted using the equation of state with parameters .

### 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 from gravitational wave observations of GW170817 (Abbott et al., 2019). The dimensionless form of the tidal deformability, for a star of gravitational mass and circumferential radius , is

(12) |

Here is the tidal Love number. Hinderer (2008) has a good discussion of how to compute given an EOS and central density (see also the erratum at Hinderer 2009). Gravitational wave measurements give a tighter constraint on the binary tidal deformability than on the individual tidal deformabilities of the two stars: the most easily measurable quantity for stars of masses and with tidal deformabilities and is (Wade et al., 2014)

(13) |

In such events the masses are not measured well individually, but the chirp mass is known precisely; for example, for GW170817 (Abbott et al., 2019). We note that, for fixed , is relatively insensitive to the mass ratio . For instance, using the scaling suggested by De et al. (2018), for is only larger than for .

Because only 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 space. For given EOS parameters, 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 parameters, we can compute . Thus in general we would compute this likelihood factor by integrating over both and :

(14) |

where and are the priors for and and is the three-dimensional likelihood obtained from analysis of gravitational wave data, given EOS parameters .

However, is known with such high precision and accuracy that, given a value for , is known to high accuracy. Therefore we can recast the likelihood factor as

(15) |

where is the prior probability density for at the value of implied by and , and the integral is over the probability distribution for obtained from the gravitational wave analysis.

#### 3.2.2 Radius and mass

Suppose that for a given star the likelihood of a mass and radius is . For a given stellar mass, the radius is determined precisely for given EOS parameters. Thus the likelihood factor associated with a radius measurement is

(16) |

where is the circumferential radius for a gravitational mass given EOS parameters , and is the prior on . Note that the integration is equivalent to integrating the full likelihood over the full curve predicted using a given equation of state.

#### 3.2.3 Combination of constraints

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

(17) |

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-parametric approach of Landry & Essick (2018) 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 discards important information. However, to our knowledge all previous analyses except Alvarez-Castillo et al. (2016) have used a hard lower bound on the maximum mass, in the sense that a given equation of state or parameter combination is allowed if it has a maximum mass above some specified value (often because the quoted of PSR J0348+0432 (Antoniadis et al., 2013) was the highest measured mass until the of PSR J0740+6620 (Cromartie et al., 2019)), 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 of the estimates of the mass of PSR 16142230. The first measurement, by Demorest et al. (2010), was . The second measurement, by Fonseca et al. (2016), was . The most recent measurement, by Arzoumanian et al. (2018), is . Thus in both updates the best estimate of the mass has been slightly more than one standard deviation below the previous best estimate. Thus a strict lower bound at the mass from the first measurement would be too restrictive given our current knowledge of the mass of PSR 16142230. Instead, one should use the full posterior distribution of the mass.

The second reason that this approach is suboptimal is that there is, after all, uncertainty in the mass measurements. If we accept as the mass estimate for PSR J0740+6620, then using the hard bound approach an EOS with a maximum mass of is just as viable as an EOS with a maximum mass of . But if we assume that the measurement has only Gaussian statistical uncertainties, there is an 84% probability that the mass of PSR J0740+6620 is greater than . Thus in reality the EOS with is considerably more consistent with the data than the EOS with . 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 equation of state are from PSR J0740+6620, from PSR J0348+0432, and from PSR 16142230. An EOS with is disfavored at the level from PSR J0740+6620 alone, but at more than when constraints from 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 is excluded much more strongly based on 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 , then using the hard bound method it would not contribute at all to EOS constraints, whereas in reality it would make low- 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 in constraints, rather than integrating in the full posterior space. Similarly, numerous papers use only the maximum-likelihood or minimum- point along the curve implied by a given EOS, whereas instead integration should 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. However, this 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.

However, even setting aside the fundamentally statistical nature of the problem, in realistic situations attempts to invert observed quantities to determine equation of state parameters fails because the inversion is singular. For example, if two 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). They point out that some neutron stars might not have a central density large enough to reach the highest density in the equation of state 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 a 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 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 the corresponding circumferential radiii 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 fm), where all the pressures agree by construction because up to that density we use the SLy (Douchin & Haensel, 2001) equation of state. We then construct the cumulative probability distribution for the pressure at progressively higher densities. We also plot the curves that follow from the equations of state implied by the 5%, 50%, and 95% curves of pressure versus density. Currently, constraints on the EOS are relatively weak, which means that a large fraction of our EOS parameter combinations end up with high likelihood, and thus we do not need to perform sophisticated searches through parameter space.

Our method can be used with any parametrization of the equation of state. We assume that the pressure is a function only of the density, i.e., that the EOS is barotropic. We therefore do not have an explicit dependence of pressure on the proton fraction, because we assume that the matter is in beta equilibrium. We follow Abbott et al. (2018a) in using the spectral parametrization introduced by Lindblom (Lindblom, 2010, 2018), in which the free parameters are spectral indices that represent the adiabatic index (where is the pressure and is the full mass-energy density) using the expansion

(18) |

where and is the pressure at half of nuclear saturation density. We also follow previous work (e.g., Carney et al. 2018; Abbott et al. 2018a) by using an expansion up to with the following uniform priors on the coefficients : , , , and . We do not additionally require, as some papers have, that at all densities; the parametrization itself guarantees that , which is needed to enforce thermodynamic stability. We also require that the adiabatic speed of sound be less than the speed of light.

Once the equation of state is chosen, 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 parametrization). 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 our constraints 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 on the symmetry energy). We then present results assuming only measurements of (1) and (2) (with a standard precision on the symmetry energy and different potential measurements on the mass), and so on. This makes it possible to see how additional measurements would progressively improve the precision of our understanding of the equation of state 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 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 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 parameters (we assume that the central density can with equal probability be anywhere between the density that would produce a 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 current constraints, although given our assumptions our 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 in the context of 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 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:

(19) |

where is the symmetry energy predicted using specified values of the equation of state parameters . For our standard constraint we choose MeV and 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.

### 5.2 Maximum mass

A viable EOS must be able to support a gravitational mass 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 for PSR J16142230 (original mass measurement in Demorest et al. 2010 and current mass measurement in Arzoumanian et al. 2018), for PSR J0348+0432 (Antoniadis et al., 2013), and 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 (van Kerkwijk et al., 2011), and another black widow system, PSR 13113430, has an estimated mass of (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 rapidly to a black hole, then , 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 that have been obtained via comparison of simulations with the kilonova that followed GW170817 (e.g., Shibata et al. 2017; Rezzolla et al. 2018; Coughlin et al. 2018).

We adopt as our standard maximum mass constraint the combination of three mass measurements: for PSR J0740+6620, for PSR J0348+0432 and for PSR J16142230. We also explore the constraints we would obtain if there is a future mass measurement of , or a future mass measurement of , or a confirmed upper limit of . 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 from an event such as GW170817 depend on the waveform model, with a spread of % among models used thus far (Abbott et al., 2019). Bearing this caveat in mind, the middle 90% of the posterior credible range for has been reported as (Abbott et al., 2019). 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. 2018b), 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 samples at this site as input for a kernel density estimate of the marginalized posterior in 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 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 (e.g., Radice et al. 2018 find a lower limit for GW170817), but we have not included such limits in our analysis.

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 (2015), and see additional caveats related to our uncertainty about the equation of state 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 from 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 than assumed in the fitting, 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 Raithel et al. 2018), 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 has been obtained for a given star. The posterior need not be a product of independent posteriors in and , or independent posteriors in and ; 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:

(20) |

where we explore the consequences of selecting and 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, fractional precision of % for a single star is necessary to add significantly to our information. As a check of our method, we find, as we must, in Figure 8 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 J07373039 (Burgay et al., 2003), it was pointed out that in principle spin-orbit coupling could be measured 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 J07373039A (which has a mass of : Kramer et al. 2006), within a few years (Lattimer & Schutz, 2005; Kramer & Wex, 2009). 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 %. For our illustrative constraint we select g cm, with Gaussian uncertainties, because this is consistent with the other real and hypothetical constraints we are applying.

We add our second hypothetical constraint in Figure 9: moment of inertia measurements for the neutron star PSR J07373039A. Measurements to the hoped-for precision of % 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 a star would have little influence on estimates of the radii of stars with masses around , but would significantly improve the estimates of the radii of 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 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 (Nomoto, 1984; Podsiadlowski et al., 2004, 2005). 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 the 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 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 J07373039B, originated from an electron-capture supernova, and that its gravitational mass of should therefore be identified with .

Thus when we incorporate this hypothetical factor into our analysis we do so by assuming that the baryonic mass corresponding to a gravitational mass is with a Gaussian likelihood.

## 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 measurements of large neutron star masses, the measurements of tidal deformability, and the future measurements of neutron star radii and masses. We find that different types of measurements will play significantly different roles in constraining the equation of state 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 dramatically enhanced understanding of the nature of dense matter.

## References

- Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Physical Review Letters, 119, 161101
- Abbott et al. (2018a) —. 2018a, Physical Review Letters, 121, 161101
- Abbott et al. (2018b) —. 2018b, arXiv e-prints, arXiv:1811.12907
- Abbott et al. (2019) —. 2019, Physical Review X, 9, 011001
- Agathos et al. (2015) Agathos, M., Meidam, J., Del Pozzo, W., et al. 2015, Phys. Rev. D, 92, 023012
- Alvarez-Castillo et al. (2016) Alvarez-Castillo, D., Ayriyan, A., Benic, S., et al. 2016, European Physical Journal A, 52, 69
- Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Physical Review Letters, 120, 172703
- Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448
- Arzoumanian et al. (2018) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47
- Bauswein et al. (2013) Bauswein, A., Baumgarte, T. W., & Janka, H.-T. 2013, Physical Review Letters, 111, 131101
- Burgay et al. (2003) Burgay, M., D’Amico, N., Possenti, A., et al. 2003, Nature, 426, 531
- Carney et al. (2018) Carney, M. F., Wade, L. E., & Irwin, B. S. 2018, Phys. Rev. D, 98, 063004
- Coughlin et al. (2018) Coughlin, M. W., Dietrich, T., Margalit, B., & Metzger, B. D. 2018, arXiv e-prints, arXiv:1812.04803
- Cromartie et al. (2019) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2019, arXiv e-prints, arXiv:1904.06759
- De et al. (2018) De, S., Finstad, D., Lattimer, J. M., et al. 2018, Physical Review Letters, 121, 091102
- DeDeo & Psaltis (2003) DeDeo, S., & Psaltis, D. 2003, Physical Review Letters, 90, 141101
- Demorest et al. (2010) Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081
- Douchin & Haensel (2001) Douchin, F., & Haensel, P. 2001, AAP, 380, 151
- Fonseca et al. (2016) Fonseca, E., Pennucci, T. T., Ellis, J. A., et al. 2016, ApJ, 832, 167
- Freire (2009) Freire, P. C. C. 2009, arXiv e-prints, arXiv:0907.3219
- Fryer et al. (2015) Fryer, C. L., Belczynski, K., Ramirez-Ruiz, E., et al. 2015, ApJ, 812, 24
- Gamba et al. (2019) Gamba, R., Read, J. S., & Wade, L. E. 2019, arXiv e-prints, arXiv:1902.04616
- Gendreau et al. (2016) Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, in Proc. SPIE, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, 99051H
- Hartle (1967) Hartle, J. B. 1967, ApJ, 150, 1005
- Hinderer (2008) Hinderer, T. 2008, ApJ, 677, 1216
- Hinderer (2009) —. 2009, ApJ, 697, 964
- Kramer & Wex (2009) Kramer, M., & Wex, N. 2009, Classical and Quantum Gravity, 26, 073001
- Kramer et al. (2006) Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2006, Science, 314, 97
- Lackey & Wade (2015) Lackey, B. D., & Wade, L. 2015, Phys. Rev. D, 91, 043002
- Landry & Essick (2018) Landry, P., & Essick, R. 2018, arXiv e-prints, arXiv:1811.12529
- Lattimer & Prakash (2016) Lattimer, J. M., & Prakash, M. 2016, Phys. Rep., 621, 127
- Lattimer & Schutz (2005) Lattimer, J. M., & Schutz, B. F. 2005, ApJ, 629, 979
- Lawrence et al. (2015) Lawrence, S., Tervala, J. G., Bedaque, P. F., & Miller, M. C. 2015, ApJ, 808, 186
- Lindblom (2010) Lindblom, L. 2010, Phys. Rev. D, 82, 103011
- Lindblom (2018) —. 2018, Phys. Rev. D, 97, 123019
- Lo et al. (2013) Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19
- Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, ApJL, 850, L19
- Martinez et al. (2015) Martinez, J. G., Stovall, K., Freire, P. C. C., et al. 2015, ApJ, 812, 143
- Miller (2013) Miller, M. C. 2013, arXiv e-prints, arXiv:1312.0029
- Miller & Lamb (2015) Miller, M. C., & Lamb, F. K. 2015, ApJ, 808, 31
- Nomoto (1984) Nomoto, K. 1984, ApJ, 277, 791
- Oppenheimer & Volkoff (1939) Oppenheimer, J. R., & Volkoff, G. M. 1939, Physical Review, 55, 374
- Orellana et al. (2013) Orellana, M., García, F., Teppa Pannia, F. A., & Romero, G. E. 2013, General Relativity and Gravitation, 45, 771
- Özel et al. (2016) Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28
- Podsiadlowski et al. (2005) Podsiadlowski, P., Dewi, J. D. M., Lesaffre, P., et al. 2005, MNRAS, 361, 1243
- Podsiadlowski et al. (2004) Podsiadlowski, P., Langer, N., Poelarends, A. J. T., et al. 2004, ApJ, 612, 1044
- Potekhin et al. (2015) Potekhin, A. Y., Pons, J. A., & Page, D. 2015, Space Sci. Rev., 191, 239
- Raaijmakers et al. (2018) Raaijmakers, G., Riley, T. E., & Watts, A. L. 2018, MNRAS, 478, 2177
- Radice et al. (2018) Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018, ApJL, 852, L29
- Raithel et al. (2018) Raithel, C. A., Özel, F., & Psaltis, D. 2018, ApJL, 857, L23
- Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Phys. Rev. D, 79, 124032
- Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, ApJL, 852, L25
- Riley et al. (2018) Riley, T. E., Raaijmakers, G., & Watts, A. L. 2018, MNRAS, 478, 1093
- Romani et al. (2012) Romani, R. W., Filippenko, A. V., Silverman, J. M., et al. 2012, ApJL, 760, L36
- Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012
- Steiner et al. (2010) Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33
- Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170
- Tolman (1939) Tolman, R. C. 1939, Physical Review, 55, 364
- Tooper (1965) Tooper, R. F. 1965, ApJ, 142, 1541
- Tsang et al. (2012) Tsang, M. B., Stone, J. R., Camera, F., et al. 2012, Phys. Rev. C, 86, 015803
- van Kerkwijk et al. (2011) van Kerkwijk, M. H., Breton, R. P., & Kulkarni, S. R. 2011, ApJ, 728, 95
- Wade et al. (2014) Wade, L., Creighton, J. D. E., Ochsner, E., et al. 2014, Phys. Rev. D, 89, 103012
- Wijnands et al. (2017) Wijnands, R., Degenaar, N., & Page, D. 2017, Journal of Astrophysics and Astronomy, 38, 49
- Yagi & Yunes (2013) Yagi, K., & Yunes, N. 2013, Phys. Rev. D, 88, 023009
- Yunes et al. (2016) Yunes, N., Yagi, K., & Pretorius, F. 2016, Phys. Rev. D, 94, 084002