Questions Related to the Equation of State of High-Density Matter

: Astronomical data about neutron stars can be combined with laboratory nuclear data to give us a strong base from which to infer the equation of state of cold catalyzed matter beyond nuclear density. However, the nuclear and astrophysical communities are largely distinct; each has their own methods, which means that there is often imperfect communication between the communities regarding caveats about claimed measurements and constraints. Here we present a brief summary from one astronomer’s perspective of relevant observations of neutron stars, with warnings as appropriate, followed by a set of questions that are intended to help enhance the dialog between nuclear physicists and astrophysicists.


Introduction
It has been noted in innumerable papers that, in the quest to understand the properties of cold catalyzed matter, laboratory measurements and astronomical observations are complementary. Fundamentally this is because (1) laboratory experiments yield precise information on matter that is nearly symmetric (i.e., with comparable numbers of neutrons and protons) and that has a density less than or comparable to nuclear saturation density, whereas (2) astronomical observations of neutron stars give us much more indirect data about significantly denser, highly neutron-rich matter. Nuclear theories that are well tuned to agree with laboratory data diverge dramatically from each other at a few times nuclear saturation density. Not only does the equation of state P( ) (where P is the pressure and is the total mass-energy density) differ between predictions: the expected composition has substantial diversity, and includes the possibilities of mainly nucleons, quark matter of various types, hyperons, and condensates.
It is, therefore, understandable that the nuclear community eagerly anticipates new astronomical information against which they can test their theories. However, astronomical data are highly indirect compared with laboratory data; although it is probably good for humanity as a whole that we are currently unable to experiment directly on neutron stars, it is inconvenient for nuclear physicists! Thus for each particular claim we need to ask: what assumptions went into deriving that constraint, and how reliable are those assumptions?
Similarly, from the standpoint of the astronomer, most nuclear experiments and theories are black boxes. It is no doubt frustrating to nuclear theorists that their own carefully constructed, deeply physically motivated, work on dense matter is placed on the same level as their rival's shoddily assembled mishmash of phenomenology. However, astronomers do not know what questions to ask. How can we judge whether a proposed equation of state, or composition, is reasonable? Even for nuclear experiments, we do not typically understand what specific assumptions are needed to go from raw data to implications about dense matter.
Here we write from one astronomer's perspective. We therefore begin by providing a brief summary of claimed astronomical constraints on dense matter. More detailed critical appraisals of astronomical measurements of neutron stars are in, e.g., references [1][2][3][4][5][6]. We then pose some questions to nuclear physicists, to help make the greatest progress in constructing precise and reliable models of dense matter.

Astronomical Measurements and Caveats
In this section we discuss astronomical measurements of neutron stars that help place constraints on the properties of cold, dense, catalyzed matter, with caveats as appropriate. Our focus will be on structural aspects of neutron stars (mass, radius, etc.), which therefore help in the inference of the equation of state; see [7,8] for recent reviews of what can be learned about thermal transport properties from observations of cooling neutron stars. See [9] for a recent discussion of how to combine disparate laboratory and astronomical observations to constrain the equation of state.

Mass Measurements
The most reliable, and most precise, structural measurements of neutron stars are of the masses of these stars when they are in binaries. If a neutron star in a binary appears to us as a pulsar, then precise pulsar timing can yield both classical parameters such as the orbital period and line-of-sight velocity, and "post-Keplerian" parameters including the Shapiro delay, advance of pericenter, the Einstein delay, and the detection of orbital decay due to gravitational radiation (see [10] for a particularly clear discussion, with equations, of all of these observables). The combination of these measurements leads to a full solution of the orbital system including both masses. The Shapiro delay, which is the magnitude and orbital phase dependence of the relative time of arrival of pulses delayed by propagation through the gravitational field of the companion, is particularly useful because the delay depends only on the mass of the companion rather than on its nature (in contrast, for example, tides on a non-pointlike companion such as a main sequence star can contribute to pericenter precession).
Indeed, the current record-holder for the highest mass established for a neutron star is M = 2.17 +0.11 −0.10 M for PSR J0740+6620, which was determined using Shapiro delay [11]. Shapiro delay was also used to measure the mass of first neutron star, PSR J1614−2230, with a gravitational mass M definitely in the ∼2 M range (original reference [12]; the best current estimate of the mass is M = 1.908 ± 0.016 M [13]). The third of the well-established > ∼ 2 M neutron stars is PSR J0348+0432. Its mass of M = 2.01 ± 0.04 M [14] was determined using a different method involving spectroscopy of the companion white dwarf.
All three of these measurements are definitive. The observations are clean, the underlying theory is well-understood and well-verified, and there are no significant residuals to the fits. We can therefore use these mass measurements, with confidence, to constrain equations of state.
Stronger but less reliable constraints follow from more astrophysically complicated analyses of additional neutron stars. It has, for example, been noted for several years that the best current estimates of the masses of some neutron stars in "black widow" systems are quite high, e.g., 2.40 ± 0.12 M for PSR B1957+20 [15], and 2.68 ± 0.14 M for PSR 1311−3430 [16]. Clearly, if the masses of these stars are as high as the best estimates, then there are important implications for the high-density equation of state. However, there are significant residuals and systematics in these fits that suggest caution in adopting those masses [15][16][17]. More recently, arguments for a relatively low maximum mass have been made based on the properties of short gamma-ray bursts [18][19][20], and when those arguments are applied to the double neutron star merger event GW170817 they suggest a maximum mass of M ≈ 2.17 M for a nonrotating star [21]. These arguments, however, assume that in such events a black hole forms rapidly. There are plausible arguments supporting that hypothesis (see [21] for such arguments), but there is no direct confirmation of the formation of a black hole and thus this limit is also subject to question. Arguments for a low upper limit to the maximum mass have also been made based on the outflow properties of the matter combined with numerical simulations; for example, [22] suggest an upper limit of 2.15-2.25 M , [23] find an upper limit of 2.16 +0. 17 −0.15 M , and [24] argue for an upper limit of 2.16-2.28 M . However, the simulations and comparisons with observations are complex and thus these limits also cannot be taken as definitive.
In sum, it is reasonable to take into account the well-established M > ∼ 2 M neutron stars when constraining equations of state. There are, for example, several papers that argue that these masses demonstrate that the conformal limit on the speed of sound, c s < c/ √ 3, is violated inside neutron stars (e.g., [25][26][27][28]). There are also long-standing nuclear arguments that M max < 2.9-3.2 M [29,30]. However, more stringent constraints are not currently robust.

Tidal Deformability Measurements
The most intriguing recent addition to neutron star measurements is the wealth of data and inference that accompanied the GW170817 double neutron star coalescence. In addition to the astrophysically uncertain mass limits discussed above, a key output of the analysis is a combination of the tidal deformabilities of the two stars involved in the merger.
The basic idea is straightforward: given that neutron stars are not point masses, they exert mutual tidal effects as they spiral close to each other. Tidal distortions require energy, which ultimately comes from the orbit, which means that two neutron stars will coalesce in a slightly shorter time than would two black holes of the same mass as the neutron stars. The more easily the stars are deformed, the greater will be the change in the coalescence time and thus in the gravitational waveform. Current analyses proceed on the assumption that tidal effects can be encapsulated in one parameter, the (dimensionless) tidal deformability where k 2 is the tidal Love number and R is the circumferential radius (see [31] for how to compute Λ). It turns out that for binary coalescences, the most tightly constrained tidal parameter is a combination of the individual tidal deformabilities [32]. This binary tidal deformability is where m 1 and m 2 ≤ m 1 are the gravitational masses, and Λ 1 and Λ 2 are the tidal deformabilities, for stars 1 and 2. The quoted limits onΛ from GW170817 depend on both assumptions and priors. Notable among these are: • Because tidal effects are strongly dependent on the separation between the stars, waveform differences increase strongly with increasing orbital frequency. This is precisely where waveform models, and specifically those that include tidal deformability as a parameter, are least certain. As a result, there is a systematic error inherent in uncertainties regarding the correct waveform model. The analysis papers therefore display tidal deformability posteriors from multiple different waveform families, which differ by ∼10% at a few hundred Hertz [33].

•
The inferred constraints on tidal deformability also depend on assumptions regarding the dimensionless spin parameter χ ≡ cJ/(GM 2 ) (where J is the angular momentum of a neutron star). A "high-spin" analysis assumes that the prior probability for χ is uniform between 0 and 0.89, which is the highest possible for a maximally-spinning star with a plausible equation of state [33]. A "low-spin" analysis assumes that χ < 0.05, which is consistent with observed pulsars in double neutron star binaries [33]. Fortunately, the tidal deformability does not depend strongly on the spin parameter prior (see [33] and the discussion below).

•
The large number of gravitational wave cycles in GW170817 means that there is an extremely precise measurement of the "chirp mass", which is M ch = η 3/5 M tot , where M tot = m 1 + m 2 and η = m 1 m 2 /M 2 tot . For GW170817, M ch = 1.186 ± 0.001 M at 90% credibility [33]. However, the mass ratio is not well constrained; at 90% credibility, using a low-spin prior, m 2 /m 1 = 0.73-1.0 [33]. This means that the primary mass m 1 could range from 1.36 M to 1.60 M , and the secondary mass m 2 could range from 1.16 M to 1.36 M [33]. These ranges are wide enough that their values depend on the prior probability distribution assumed for the masses of neutron stars. For example, should one assume that the mass has equal prior probability from some minimum mass (say, 1.0 M ) to the maximum allowed mass? Or should one place a mass prior based on the observed (and highly incomplete!) set of measured neutron star masses, or based on a different criterion?
In addition, many analyses assume that both stars are in the same family of compact objects (e.g., that both are hadronic stars), but the possibility of different families must be considered (e.g., [34][35][36]); among other possibilities, transitions between families could help explode M > 50 M stars [37] and might eventually be detected in the peak frequency of the postemerger gravitational wave signal [38]. With all of these caveats, the quoted 90% credibility range forΛ is (0, 630) for the high-spin prior and (70, 720) for the low-spin prior [33]. When combined with a M max > 1.97 M restriction (although note that, properly, the full mass posterior should be used; see [9]), the tidal deformability measurement implies that the radii of both stars are R = 11.9 ± 1.4 km at 90% credibility [39]. This range rules out especially hard equations of state, at least for the family of stars involved in the merger (if there are multiple families then other families could be hard; see, e.g., references [40,41]).
The prospects for improvement based on future observations are good, with some reservations. GW170817 produced a tremendously strong signal; its signal to noise ratio of ∼33 was the highest of any gravitational wave event yet observed [42]. This stemmed from its extreme closeness (just ∼40 Mpc), which in turn meant that the host galaxy was readily identified. Most future detections of double neutron star coalescences will not be so fortunate. However, the improved sensitivity from additional runs, including the ongoing LIGO/Virgo O3 run, will certainly help, as will joint analysis of multiple events [43,44]. In particular, if the high frequency sensitivity is improved as expected by the use of squeezed light the effects of tidal deformation will be more significant, although as a consequence the waveform models will be expected to diverge more and thus systematic errors are likely to become more prominent. Overall, as gravitational wave astronomy progresses we can anticipate major improvements in our understanding of the equation of state of dense matter.

Future Radius Measurements
A precise and reliable radius measurement for a neutron star of known mass would have major implications for our understanding of neutron star structure. Our opinion, which we have given in much greater detail in [3,4], is that at present systematic errors in the interpretation of (mostly X-ray) data are large enough that none of the published estimates yield reliable constraints on models for the equation of state. Here we will simply give one indication of how biases can enter current analyses, and will then talk about the brighter prospects for the future.
Neutron star radius estimates currently concentrate on X-ray observations, because there is a reasonable argument that X-rays from neutron star surfaces can be understood in the context of thermal mechanisms. For neutron stars with surface magnetic fields B < 10 9 G, and a known surface composition, spectra and beaming patterns can then be computed with high precision (e.g., references [45][46][47] and used in the analysis of data to measure the mass and sometimes also the radius. The estimates of radius have focused on the analysis of neutron stars cooling after thermonuclear X-ray bursts (starting with the pioneering work of Jan van Paradijs; see [48] for a good review of some of the early ideas) and the so-called "quiescent low-mass X-ray binaries" (qLMXBs; see [49,50] for some early work), which are thought to be non-accreting neutron stars in binaries. A standard assumption in such analyses is that the entire surface of the star emits uniformly. This cannot be true to arbitrary precision, so we can explore what happens to our fits when the temperature varies over the surface.
In Figure 1 we display the result of such a fit, for a nonrotating neutron star. We assume that the spectrum everywhere is a blackbody, with an effective temperature that ranges smoothly from kT eff = 2 keV at the equator to kT eff = 0.2 keV at the poles. Assuming a perfect detector response and no interstellar absorption (both assumptions make it easier to detect deviations from a blackbody), we fit the resulting spectrum with a single-temperature blackbody. The radius used to construct the synthetic data was 12 km, and we generated roughly N = 3 × 10 4 counts, which is comparable to the number seen in a typical thermonuclear X-ray burst. The fitted radius was R fit = 10.45 ± 0.14 km, and χ 2 /dof = 180.8/173 for the single-temperature fit. Figure 1. Demonstration that if there is a temperature gradient over a surface, a single-temperature fit can appear to be statistically good but nonetheless produce strongly biased radius measurements. Here the synthetic data have ∼3 × 10 4 counts, which is typical of individual thermonuclear X-ray bursts, and were constructed using a large temperature gradient: the effective temperature is kT eff = 2 keV at the equator but is kT eff = 0.2 keV at the poles. The inferred radius of R fit = 10.45 ± 0.14 km strongly excludes the radius R = 12 km used to construct the data, and yet the formally good fit χ 2 /dof = 180.8/173 (solid line) gives no hint that the model and inference are incorrect.
The takeaways here are that (a) the quality of the fit is excellent, which means that we have no indication that the model is incorrect, and (b) the fitted radius is highly biased; the true value is excluded at high confidence. This is just one indication of the potential issues with current fits of X-ray data. For a comprehensive discussion of the effects that distance, composition, and other uncertainties have on radius estimates for qLMXBs, see [51]. For burst sources, many thermonuclear bursts from neutron stars show periodic modulation of their X-ray intensity, which implies that their surfaces do not emit uniformly, and even those that show no modulation often do not behave in accordance with the models that are commonly used to interpret them (this was originally pointed out by [52], and explored in much greater detail with more bursts by [53]).
The crux of the problem is that for these analyses it is possible to get formally acceptable fits, which yield apparently precise measurements of the radii (and sometimes the masses) of the neutron stars, but which are actually biased significantly because the models are incorrect. Thus we could be wrong, but not know it.
There is reason to hope that the estimates that emerge from observations using the Neutron star Interior Composition ExploreR (NICER, see [54]) will not suffer from that deficiency. NICER has accumulated more than a million seconds each on a small number of weakly magnetic, rapidly-rotating neutron stars, and the data that are being analyzed are the energy-dependent X-ray waveforms. Hotter regions on the surface, which rotate with the star, modulate the observed spectrum and flux. Models involving one or more hotter regions are then fit to the data and yield, among other parameters, the mass and radius of the star.
The models are necessarily simplified and usually assume, for example that the hotter regions have uniform temperature (although overlap of different spots could simulate a temperature gradient). However, when [55,56] used these simplified models to analyze synthetic data that were constructed using different assumptions (e.g., oval spots, temperature gradients in the spots, different beaming patterns, and different spectra) than were used in the parameter inference, they found that if the fit is formally statistically good (as assessed, e.g., using a chi squared test) then the estimated parameters are not significantly biased. This leads to cautious optimism that analysis of NICER data will yield reliable as well as interestingly precise constraints on the radius and mass of a few neutron stars.
We conclude this section by noting that, given a specified equation of state family, tidal deformability measurements by themselves can be used to constrain the radius (e.g., [57][58][59]). Given additional assumptions about merger observations (in particular the question of whether, and how rapidly, a black hole forms), somewhat more model-dependent radius estimates are also possible [60][61][62][63][64] (see some important caveats at [65]).

Future Moment of Inertia Measurements
The rotation of neutron stars in double neutron star binaries can produce potentially observable effects depending on the orientation of the system. One effect, sometimes known as geodetic precession, is that in a general case of misalignment between the rotation axes and the orbital axis, all axes will precess around the direction of the total angular momentum. If the rate of orbital precession is measured precisely enough, and if the orbital angular momentum as well as the rotation rate and orientation of the pulsar axes are known well, then the rotational angular momentum and thus the moment of inertia of the neutron stars could be established. In combination with well-known masses, this would provide a complementary way to get at the equation of state.
A second method involves the advance of the pericenter of the system, due to spin-orbit coupling. This effect is expected to be measurable in the double pulsar system PSR J0737−3039, which has a primary that has a rotation period of 22.7 ms [66]. Some estimates are that the moment of inertia will eventually be determined to a precision of ∼10% [67,68], although the analysis has proven to be more challenging than originally envisioned. If such precision can be realized in practice then the moment of inertia plus the excellent knowledge of the masses (better than 0.1% precision; see [66] and many later papers) will yield a strongly competitive constraint on the equation of state.

Future Constraints Based on Gravitational Binding Energy
An astrophysically uncertain constraint on relatively low-mass neutron stars comes from the suggestion that a particular type of core-collapse supernova, called an "electron-capture supernova", can occur only if the baryonic rest mass of the core of a massive star, just prior to collapse, is in the mass range M bary = 1.366-1.375 M [69][70][71]. If the subsequent collapse contains all of those baryons and no more (so that matter is neither driven away nor added by later fallback), and if the neutron stars resulting from this mechanism can be identified and their gravitational mass established, then for this particular mass we know the gravitational binding energy of the star. The difficulties of applying this argument are multiple: for example, it could be that the mechanism is not as clear as it appears to be, or that matter is expelled or accreted. Moreover, the most precisely determined masses (from binary pulsar observations) do not pile up at a particular value that might make us suspect that several of those neutron stars formed by the electron capture mechanism (see https://www.stellarcollapse.org/nsmasses for an updated list of neutron star masses). Indeed, if the M = 1.174 ± 0.004 M companion to PSR J0453+1159 [72] is a neutron star then it is too light to have formed via an electron-capture supernova and other mechanisms such as ultra-stripped supernovae [73] would need to be invoked. However, future theoretical developments and the march of observational results could bring this constraint into play.

Questions for Nuclear Physicists
In the previous section, we presented an overview of existing astronomical constraints on neutron star matter, along with caveats as appropriate. The fundamental recommendation is: when a claimed astronomical constraint is proposed, nuclear physicists should examine the assumptions carefully and look for the possibility of systematic error. The most insidious systematic errors are those that can bias the results without leading to an obviously bad fit, because then it is possible to unknowingly draw incorrect inferences. It is useful to talk with multiple specialists in the relevant astrophysical field, to get a sense for whether the results are reliable and believed throughout the field.
In this section, we go the other way, by discussing some questions that astrophysicists could ask of nuclear physicists. We again focus on the equation of state, but similar questions could be asked regarding transport properties.
What model-dependence is there in quotations of results?-Various quantities, such as the symmetry energy and its density derivative, or the neutron skin thickness of 208 Pb or 48 Ca, have been advertised as having important implications for the equation of state of nuclear matter. However, in detail there are some complexities. For example, the quoted range of the symmetry energy depends on an experimental quantity (the binding energy of symmetric nuclear matter at nuclear saturation density) and a theoretical quantity (the energy per nucleon of pure neutron matter at the same density). How should we judge the uncertainties? As another example, the neutron skin thickness of lead promises to tell us about very neutron-rich matter, albeit at densities much less than nuclear saturation (e.g., [74]). Will that have immediate implications for the radii of neutron stars, or does that association require the assumption of a particular model framework?
What astronomical information is most important to make progress?-Suppose that you, as a nuclear physicist, could get one outstandingly precise piece of information about a single neutron star, but you could know nothing else about the star. What would tell you the most? The radius? The tidal deformability? The mass? The moment of inertia? In some cases, it would surely matter what the answer is. For example, a single neutron star with a precisely known mass of 2.4 M would have important implications. What if instead of a single precisely known value, you had fuzzier information about multiple quantities, or about multiple stars? What combination, and to what precision, would be the most informative for you?
What parameterized equation of state model should we use?-When we are interested in the structural properties of neutron stars, we could use specifically proposed equations of state, but given that none of them is likely to be exactly right, it seems reasonable to use a parameterization. However, several forms exist in the literature, including piecewise polytropic equations of state (e.g., [75,76]), spectral models of the adiabatic polytropic index (e.g., [77]), and polynomial forms based on quantum Monte Carlo techniques [78]. Are any of these preferred over any other, or should we simply try a few possibilities to determine how much the results depend on our choice?
What priors should astronomers apply?-Ideally, data analysis will be largely independent of any prior probability distributions that are applied because the data will be highly informative. Realistically, however, other than some neutron star mass measurements, astronomical observations of neutron stars have not reached that ideal. Thus priors matter to some degree. What, then, should be applied? Should we have uniform priors on the mass? Mass priors corresponding to the neutron star masses we have measured (which, as a reminder, are a tiny fraction of all neutron stars)? Priors on the radius or tidal deformability? If we apply priors to our equation of state parameters, how do we include all reasonable possibilities while excluding the ones that are unreasonable?
How can we judge the physical realism of specific equation of state models?-From the typical astronomer's perspective, one equation of state model is as good as another; all are thrown into the mix to determine whether new observations can provide discriminatory power. However, many of those models are decades old and do not satisfy the latest laboratory data, or contain assumptions that might range from implausible to impossible. What should we look for? Is there, for example, a standard set of experimental and sanity tests that all proposals must pass? If so, we can be alert to whether a particular proposal satisfies those criteria.
Can your model be ruled out, realistically?-One of the frustrations of astronomy is that theorists can be endlessly slippery. For example, suppose that a dark matter model is proposed in which the matter has an arbitrarily adjustable interaction strength. Then a detection confirms the model and leads to a trip to Stockholm, but a nondetection simply means that the interaction strength has been constrained to be below some value. What is the case for your model of cold, catalyzed matter at up to several times nuclear saturation density? If any realistic observation would merely constrain your parameters rather than potentially ruling out the entire model, which would be disappointing to astronomers because it would mean that your model is not falsifiable and thus, ultimately, does not have much explanatory power. By "realistic observation" we take into account that, sure, a neutron star with M = 1.4 M and R = 5 km might kill your model, but we are thinking about possibilities within the bounds of current observations. If in contrast there is some plausible observation that might completely kill your model, then as a result your model will attract enhanced attention and will receive greater credibility if it survives.

Conclusions
The increasingly rapid flow of information about dense matter from experiments and observations makes this a heady time for the development of this branch of nuclear physics. In order to make optimal progress, nuclear physicists and astrophysicists will need to work together closely, with critical assessment of each inference, to assure that the reasoning is as sound as possible. It will be exciting to follow and participate in that partnership! Funding: This research received no external funding