# Neutron Star Radii, Universal Relations, and the Role of Prior Distributions

###### Abstract

We investigate constraints on neutron star structure arising from the assumptions that neutron stars have crusts, that recent calculations of pure neutron matter limit the equation of state of neutron star matter near the nuclear saturation density, that the high-density equation of state is limited by causality and the largest high-accuracy neutron star mass measurement, and that general relativity is the correct theory of gravity. We explore the role of prior assumptions by considering two classes of equation of state models. In a first, the intermediate- and high-density behavior of the equation of state is parameterized by piecewise polytropes. In the second class, the high-density behavior of the equation of state is parameterized by piecewise continuous line segments. The smallest density at which high-density matter appears is varied in order to allow for strong phase transitions above the nuclear saturation density. We critically examine correlations among the pressure of matter, radii, maximum masses, the binding energy, the moment of inertia, and the tidal deformability, paying special attention to the sensitivity of these correlations to prior assumptions about the equation of state. It is possible to constrain the radii of neutron stars to a be larger than 10 km, even without consideration of additional astrophysical observations, for example, those from photospheric radius expansion bursts or quiescent low-mass X-ray binaries. We are able to improve the accuracy of known correlations between the moment of inertia and compactness as well as the binding energy and compactness. We also demonstrate the existence of a correlation between the neutron star binding energy and the moment of inertia.

###### pacs:

97.60.JdNeutron stars and 26.60.-cnuclear matter aspects of neutron stars and 97.80.JpX-ray binaries and 21.65.Cdnuclear matter## 1 A Brief History of Neutron Star Radii

Theoretical estimates of radii of spherically-symmetric non-rotating neutron stars came directly from analytic and numerical solutions of the Tolman-Oppenheimer-Volkov equations Tolman39ss (); Oppenheimer39om (). At this time, the critical role played by the maximum mass was already known. Twenty years later, improved estimates of radii came directly from improvements in nuclear physics: Ref. Cameron59ns () used a recently developed zero-range nucleon-nucleon interaction (the Skyrme interaction Skyrme59te ()) to compute neutron stars with masses up to a maximum of and radii as small as 8 km and even larger than 27 km. It was not realized at this time that the EOS was acausal at high densities, thus leading to extremely small radii. Also, the extremely large radii were found for stars with masses below what is commonly considered to be the minimum formation mass of neutron stars, about . Later calculations of the EOS of neutron star matter using both Skyrme and relativistic field-theoretical interactions constrained by nuclear experiment have found radii between about 9 and 15 km Steiner05ia (). Astronomical observations of radio pulsations in neutron stars followed 8 years later in 1967; observations of pulsars in binary systems have lead to numerous mass measurements (e.g., see Ref. Phinney94 ()). Importantly, the largest accurately measured masses Demorest10 (); Antoniadis13 () set a lower limit to the neutron star maximum mass, which, combined with general relativity and causality, constrain minimum values for neutron star radii Lattimer12 () of all masses. For a star, the current minimum maximum mass establishes that its radius km. Larger mass measurements will increase this minimum radius. An important theoretical development was the realization that neutron star radii for stars in the mass range are primarily determined by the EOS in the density range Lattimer01 (), where is the nuclear saturation density. This provides a direct link between both nuclear experiment and theory and neutron star radii (though the link is weaker if there is a strong phase transition in this density range, as demonstrated below). One of the objectives of this paper is to use recent developments in nuclear experiment and theory, coupled with the assumption that neutron stars have crusts, to establish neutron star radius constraints. In addition, we investigate how and whether astrophysical observations of neutron stars in X-rays can realistically further improve these limits.

Ref. vanParadijs79 () proposed using photospheric radius expansion (PRE) X-ray bursts to obtain simultaneous mass and radius measurements in 1979, but the method did not lead to interesting constraints until 2006 Ozel06 (). Masses and radii from PRE X-ray bursts Guver10a (); Guver10b (); Ozel09b () gave radii between 8 and 12 km and masses between 1.2 and 2 under the assumption that the photosphere at “touchdown” has receded to the neutron star surface. Under the assumption that this model describing PRE X-ray bursts was correct, large radii and EOS models with high pressures between and seemed to be ruled out.

Observations of thermal emission from accreting neutron stars in binary systems, known as quiescent low-mass X-ray binaries (qLMXBs), also provide radius measurements. The earliest such observations led to radius estimates of less than 1 km, much smaller than that predicted from theoretical models vanParadijs87sx (). For accreting neutron stars, the photosphere is expected to consist of hydrogen bildsten92 (), and Ref. Rutledge99tt () showed that models of thermal emission from pure hydrogen photospheres rajagopal96 (); zavlin96 () gave inferred radii on the order of 10 km. Neutron star radius information from quiescent low-mass X-ray binaries began to provide significant constraints on radii in 2006 (at the same time as PRE X-ray bursts), due to improved observations and a consistent treatment of the surface gravity in the neutron star atmosphere model used to fit the spectrum Heinke06 (); Webb07 (). This led to refined estimates of radii between 8 and 16 km, which, however, by themselves did not rule out a significant number of contemporary theoretical EOS models Lattimer01 ().

There were two principal difficulties which prevented these measurements from tightly constraining the EOS. The first difficulty was the fact that taking neutron star mass and radius measurements and constructing an EOS formally requires an inversion of the TOV equations. A method to invert the TOV equations was developed in 1992 Lindblom92 (), but this approach is subject to numerical difficulties (see recent work in Ref. Lindblom14 ()). The second difficulty was that the traditional approach of performing a chi-squared fit is prohibited by the large uncertainties and the potential for the mass-radius curve to be non-monotonic in both mass and radius (see a more recent summary of the possible mass-radius curve shapes in Ref. Alford13 ()). Ref. Ozel09 () showed that, if one chooses a sufficiently generic EOS parameterization, a formal inversion is unnecessary. Ref. Read09 () showed that a three-parameter model based on piecewise polytropes (commonly used, e.g., in numerical solutions of rotating relativistic stars Eriguchi85 ()) is general enough to reproduce modern theoretical EOSs to within a few percent. Thus Ref. Ozel09 () used this parameterization to show that, given simultaneous mass and radius measurements of three neutron stars, one can marginalize over the unknown neutron star masses to obtain probability distributions for the three unknown EOS parameters. This was applied in Ref. Ozel10 () to three simultaneous mass and radius constraints for three neutron stars. For the EOS near and below the nuclear saturation density (assumed to have zero uncertainty), Ref. Ozel10 () used the Skyrme interaction SLy4 that was fit to nuclei and theoretical calculations of low-density neutron matter Chabanat95ns (). It was claimed that the mass and radius measurements presented a challenge for theories of neutron star matter, but it is now understood that systematic uncertainties in the X-ray burst model are the more likely culprit for the tension that was observed.

Soon after, Ref. Steiner10te () used Bayesian inference to combine mass and radius constraints from both qLMXBs and PRE X-ray bursts to obtain constraints on the EOS. This work used an alternate nuclear physics-based parameterization for matter near the nuclear saturation density and piecewise polytropes at higher densities (in a slightly different form than that presented in Ref. Read09 ()). In the case that the prior distributions for the polytrope parameters are uniform, the Bayesian inference-based method to obtain the EOS parameters is similar to the method developed in Ref. Ozel09 (), except that Ref. Steiner10te () additionally obtained probability distributions for the mass-radius curve by marginalizing over the posteriors for the radius as a function of mass. Also, the direct use of Bayesian inference in Ref. Steiner10te () allowed the use of more parameters which more fully explored the uncertainties in the EOS near saturation densities and led to constraints on the density dependence of the nuclear symmetry energy. Finally, Ref. Steiner10te () showed that the assumption that the photosphere at touchdown is coincident with the surface (i.e., not extended) is not consistent with the data (within the context of the model being used for PRE X-ray bursts in Refs. Steiner10te () and Ozel10 ()). The final result was a radius range of 10.7 to 12.5 km for a 1.4 neutron star. Radii smaller than 10.7 km were ruled out, in contrast to earlier results based on PRE X-ray bursts. These improved radius constraints came principally as a result of the marginalization from Bayesian inference. However, it was reported in Ref. Steiner10te () that there were several systematic uncertainties which potentially invalidated this result.

One important systematic uncertainty was the potential for accretion to muddle the interpretation of PRE X-ray bursts. In 2011, Ref. Suleimanov11 () used longer bursts and a more complete model of the neutron star atmosphere to obtain a radius greater than 14 km for the single source studied. For the qLMXB sources, the composition of the atmosphere plays an important role. Ref. Servillat12 () showed that a helium (rather than hydrogen) atmosphere changes the inferred radius range for the neutron star in M28 from 6 to 11.5 km to 7 to 17 km. Both of these systematic uncertainties continue to play a role in interpreting mass and radius measurements. At around the same time, work progressed on constraining the nature of the nucleon-nucleon interaction from neutron star mass and radius observations Hebeler10b (); Steiner12cn (), showing that current observations suggested a relatively weak density dependence in the nuclear symmetry energy, predicting radii less than 14 km.

In 2013, Ref. Steiner13tn () re-analyzed the qLMXB and PRE data including some systematic uncertainties ignored in previous works, including variations in the color correction factor used to infer masses and radii from the PRE X-ray bursts. This systematic uncertainty extended the radius range from 10.4 to 12.9 km. Another critical systematic included in Ref. Steiner13tn () is the ambiguity in the choice of the prior distribution. Marginalizations over model parameters include a choice of prior distribution, whether implicit as in Refs. Ozel09 (); Ozel10 () or explicit as in Ref. Steiner10te (). The ambiguity from the choice of prior is not important if the number of parameters is sufficiently small. Also, prior assumptions are not important if uncertainty in the parameter that represents the independent variable of the model function (e.g., the mass of the neutron star) is small. If both of these conditions held, one could perform the classical chi-squared fit and directly obtain the mass-radius curve and the associated EOS. However, these conditions did not then (and still do not) hold for neutron star mass and radius observations.

One way the choice of prior distribution is important is the extent to which it favors or disfavors the presence of phase transitions where the pressure is nearly flat with increasing density. Polytropic EOSs, when used with uniform prior distributions for the polytropic exponent, naturally disfavor phase transitions. The reason for this is that all polytropes go through the origin, thus the polytropic exponent must be very small in order to produce a nearly flat pressure as a function of energy density (which would result from a strong phase transition). Ref. Steiner13tn () found that this prior ambiguity was also important and also that the 14 km radius obtained in Ref. Suleimanov11 () was difficult to reconcile with the smaller radii obtained from qLMXB sources.

Ref. Guillot13 () combined new qLMXB observations and a constant radius model (motivated by the vertical shape of the mass-radius curves obtained in Ref. Steiner13tn ()) to obtain a very small preferred neutron star radius with rather small error, km. However, when individually analyzed, the sources produced predicted radii in a wide range (7 to 20 km). Ref. Lattimer14co () concluded that two systematic uncertainties were responsible for this result: (i) some of the small radius neutron stars may have helium rather than hydrogen atmospheres (as discussed above) and (ii) uncertainties in the hydrogen column density may give incorrect radii. Ref. Heinke14 () confirmed that the hydrogen column density inferred for one neutron star in the data set was indeed systematically shifted. Ref. Heinke14 () also showed that the choice of galactic abundance model was an important systematic.

Ref. Lattimer14co () also introduced the use of the Bayes factor
^{1}^{1}1The Bayes factor is sometimes referred to as the “odds
factor” or “likelihood ratio”. The latter term, however, is often
used to describe the ratio of the maximum likelihoods, rather than
the ratio of the integrals. to decide between various models. In
the case of a chi-squared fit, the model with the smaller value of
chi-squared gives a better fit. In the frequentist approach, an
equivalent method is the maximum likelihood test: the
preferred model is the one with the largest maximum likelihood.
A test using the Bayes factor is the Bayesian analog of the maximum
likelihood test. The Bayes factor for model A versus model B is a
ratio of the “evidence”. The evidence, in turn, is defined as the
integral over the posterior (see Ref. Lattimer14co () for a brief
review and Ref. Steiner15mb () for additional discussion).

Bayesian inference continues to be an important tool for analyzing neutron star mass and radius observations, in part because the uncertainties are still large and the problem of determining the mass-radius curve (or the EOS) is underconstrained. For this reason, the Bayes factor remains one of the best tools to compare the evidence for (or against) the various models and assumptions which are employed. Models with tightly constrained posteriors are not preferred because they correspond to small values of the evidence. One way to think about this result is to note that models with tightly constrained posteriors are too finely-tuned because there are only a very few subset of model parameters for which the model is not ruled out.

More recently, Ref. Lattimer14co () showed that Bayes factors imply that an extended photosphere at touchdown is preferred in PRE X-ray bursts (the posteriors for an extended photosphere are broad and thus that the evidence for the model is large) and helium, rather than hydrogen, atmospheres are favored in some qLMXBs. However, this work did not rule out the potential for other interpretations and Ref. Ozel15a () has since suggested that rotation and temperature corrections to the Eddington limit may also play a role.

Nevertheless, Ref. Guillot14 (), implementing new data which resolved some difficulties in the previous study Guillot13 (), but retaining their assumptions of a common radius and hydrogen atmospheres for all sources, continued to obtain a small radius ( km).

## 2 The Role of Prior Distributions

The choice of prior distribution continues to play an important role in interpreting observations and theory. This statement trivially holds true, as in the context of Bayesian inference all model assumptions can be viewed as a particular choice of prior distribution. Nevertheless, it is important when one can quantify the effect that these assumptions have. A recent analysis of both the qLMXB and PRE data in Ref. Steiner15un () gives two different ranges for the radius of a 1.4 solar mass neutron star under different prior assumptions. Model A (based on polytropes) gave a radius range of 10.8 to 12.4 km while Model C (which allows for stronger phase transitions) gave a radius range of 10.2 to 11.9 km (this distinction was first observed in Ref. Steiner13tn ()). Since the Bayes factor for Model A versus Model C is not significantly different from 1, one must presume that either model could be correct, and the full range of allowed radii for a 1.4 neutron star is between 10.2 and 12.4 km. The prior ambiguity increases the uncertainty by at least 30%. The constraints on the pressure at twice the saturation density, as inferred from neutron star radius measurements, are also sensitive to the prior distribution: Model A gives a range between 9.1 and 23.0 while Model C gives a range between 2.3 and 17.0 Steiner15un (). Had one considered only a polytrope-based parameterization of the EOS one would have overestimated the lower limit on the pressure by a factor of four. Recently, Ref. Nattila15 () found a similar variation in the radius when analyzing X-ray bursts in the hard state. The final result was 11.3 to 12.8 km for a 1.4 neutron star for Model A and 10.5 to 12.5 km for Model C.

Ref. Ozel15b () compared two different methods for obtaining masses and radii from PRE X-ray bursts. The first method is basically the method given in Ref. Ozel10 () and Ref. Steiner10te (). In this first method, the marginalization over model parameters includes a Jacobian factor which transforms between a two-dimensional space over touchdown flux and normalization to a two-dimensional space over mass and radius. In the second method of Ref. Ozel15b (), this Jacobian is removed. Ref. Ozel15b () confusingly refers to their first method as the “frequentist” approach and the second method as a “Bayesian” approach, even though the presence of the Jacobian is simply a transformation of variables (and is widely used by Bayesian practitioners Gelman14 ()). In this article, we employ a language that is more typical in the statistics literature. In Ref. Ozel15b (), these two methods are compared by a qualitative examination of “bias” in the method. Bias is measured by comparing the posterior in the mass-radius space distribution to the prior distribution in mass-radius space.

A quantitative approach which is more similar to the language of the statistics literature would compute the evidence for these two methods and form the ratio to compute the Bayes factor. The Bayes factor is not the same as the “bias” used in Ref. Ozel15b (): it measures the integral under the posterior rather than its location in parameter space relative to the prior. In the case that the relative Bayes factor between these two models is near unity both methods would need to be considered, just as was done with Models A and C in Ref. Steiner15un () and also below.

Prior distributions are also relevant for determining the nature of
“correlations” or “universal relations”. There have been several
correlations observed within nuclear structure and between
nuclear structure and neutron stars: the correlation between the
EOS of neutron matter and the neutron skin thickness of
lead Typel01 (), the correlation between the pressure of
neutron-rich matter and neutron star radii Lattimer01 (), and the
correlation between the neutron skin thickness of lead and neutron
star radii Horowitz01 (). These correlations are principally
explored by studying a range of several model parameterizations. They
reflect the nature of the uncertainties in two quantities vis à
vis the current knowledge of the nucleon-nucleon interaction. In the
context of Bayesian inference, these correlations describe the nature
of the posterior distribution of two quantities given a particular
prior distribution. The shape of the correlation is thus dependent on
the prior distribution^{2}^{2}2Alternatively, in the context of a
hierarchical Bayesian analysis, the shape of the correlation is
dependent on the hyper-prior distribution of hyper-parameters..

In the past several years, other kind of correlations within neutron stars, which have been referred to as “universal relations”, e.g., the correlation between the moment of inertia and the Love number of a neutron star Yagi13 (). The uncertainties in these universal relations are also dependent on the prior assumptions, though the magnitude of this prior dependence can be very different from one correlation (or one universal relation) to another.

## 3 The Equation of State

### 3.1 Neutron Star Crusts and the Low-Density EOS

We assume that neutron stars have crusts, i.e., they have a surface region with densities less than approximately g cm composed largely of nuclei, neutrons and electrons in beta equilibrium Baym71 (). The pressure in this region is largely due to relativistic, degenerate electrons with at most a 5% contribution from nuclei and neutrons. Since the nuclei are in pressure equilibrium with the neutrons, they individually contribute almost no pressure since their baryon density is close to the nuclear saturation density fm or g cm where uniformly dense symmetric matter has zero pressure. The major contribution of baryons to the pressure is from the collective Coulomb pressure due to the nuclear lattice, and is therefore largely independent of uncertainties in the equation of state of nuclear matter. Various calculations indicate that the transition from crustal material to uniform nuclear matter, , occurs in the range .

Recent calculations of the properties of pure neutron matter have produced estimates of the pressure-energy density relation up to about . Matter in neutron stars at densities between and is extremely neutron-rich because it is in beta equilibrium. This condition is equivalent to minimization of the total energy per baryon with respect to the charge fraction where and are the neutron and proton baryon densities, respectively, and . The difference between the energy of pure neutron matter and symmetric matter (with equal numbers of neutrons and protons) is called the nuclear symmetry energy , and the energies of intermediate proton fractions can be approximated with a quadratic interpolation between these extremes:

(1) |

where is the energy per baryon of symmetric matter. A crude estimate for the symmetry energy near is

(2) |

nuclear experimental information and neutron matter calculations indicate that and . The pressure corresponding to Eqs. (1) and (2) is

(3) | |||||

where is the pressure of symmetric matter. Note that, by definition, ; to leading order, near , the symmetric matter pressure increases linearly with density, , where MeV is the nuclear incompressibility parameter.

As stated above, matter in neutron stars is in beta equilibrium:

(4) |

where is the electron energy per baryon assuming relativistic degenerate electrons, and the ’s are chemical potentials. This is equivalent to

(5) |

This equation can be solved as a cubic equation for at a specific density with the approximate solution

(6) |

which has the value when . The pressure of pure neutron matter with ansatz Eq. (2), at , is In beta equilibrium, the approximate neutron star pressure at is

(7) |

The correction term in Eq. (7) is of order 1.4%, and can be ignored to good approximation. At higher densities, the proton fraction and the correction term generally increase due to the increasing symmetry energy. There is also a contribution from . However, for densities up to the neutron star matter pressure is essentially equivalent to pure neutron matter pressure.

Therefore, for densities between the core-crust transition density and about , experimental information about the symmetry energy and calculations of pure neutron matter pressures offer viable constraints on the equation of state. Experimental information concerning the symmetry energy is usually encoded in the parameters and , defined as

(8) | ||||

(9) |

For the symmetry energy of Eq. (2), one finds so that .

Model | ||||||||
---|---|---|---|---|---|---|---|---|

MeV | MeV | MeV | MeV | MeV fm | ||||

GCR 0 | 12.7 | 0.49 | 1.78 | 2.26 | 30.5 | 31.3 | 7.272 | 2.113 |

GCR 1 | 12.7 | 0.48 | 3.45 | 2.12 | 32.1 | 30.8 | 10.402 | 2.335 |

GCR 2 | 12.8 | 0.488 | 3.19 | 2.20 | 32.0 | 40.6 | 10.537 | 2.343 |

GCR 3 | 13.0 | 0.475 | 3.21 | 2.47 | 32.0 | 44.0 | 13.274 | 2.487 |

GCR 4 | 12.6 | 0.475 | 5.16 | 2.12 | 33.7 | 51.5 | 14.304 | 2.533 |

GCR 5 | 13.0 | 0.50 | 4.71 | 2.49 | 33.8 | 56.2 | 18.678 | 2.700 |

GCR 6 | 13.4 | 0.514 | 5.62 | 2.436 | 35.1 | 63.6 | 20.933 | 2.770 |

DSS 0 | 10.94 | 0.459 | 4.106 | 1.977 | 31.1 | 39.4 | 8.125 | 2.182 |

DSS 1 | 11.00 | 0.460 | 4.425 | 1.947 | 31.4 | 41.0 | 8.453 | 2.206 |

DSS 2 | 11.95 | 0.495 | 3.493 | 2.632 | 31.4 | 45.3 | 13.760 | 2.509 |

DSS 3 | 11.02 | 0.460 | 4.683 | 1.935 | 31.7 | 42.4 | 8.768 | 2.229 |

DSS 4 | 10.95 | 0.454 | 5.158 | 1.972 | 32.1 | 45.4 | 9.676 | 2.290 |

DSS 5 | 10.34 | 0.429 | 4.954 | 2.024 | 31.3 | 43.4 | 9.180 | 2.258 |

DSS 6 | 10.29 | 0.433 | 7.227 | 1.842 | 33.5 | 53.3 | 11.241 | 2.384 |

### 3.2 Piecewise Polytropes For the High-Density EOS

In our first model, we assume the crust-core transition density is , and use the crust EOS from Ref. Baym71 (), for which the pressure is MeV fm, the energy density is MeV fm, and the internal energy per baryon is MeV. For densities above the core-crust transition density, a scheme similar to a three piecewise polytrope scheme explored by Ref. Read09 () is used. There are 7 parameters. Four parameters correspond to boundary density points (), and three correspond to the polytropic exponents () of the region between the densities and . Read et al. Read09 () discovered that a variety of equations of state were successfully approximated by piecewise polytropes with a common set of boundary densities and . Assuming that the core-crust transition pressure and energy density are fixed by the crustal equation of state, and the core-crust transition density is chosen to be , fixes two more parameters. In this case, the three remaining parameters are the polytropic exponents . We choose, equivalently, the pressures and at the boundaries as parameters, following Ref. Ozel09 (). The polytropic exponents are given by

(10) |

for . Assuming continuity of both energies and pressures at the boundary points , the density and energy density within , for , and for , for , are

(11) | |||||

(12) |

where and are the energy density and energy per baryon at the point , respectively.

The boundary is sufficiently close to that neutron matter calculations, which are claimed to be reliable to such densities, offer a viable method of estimating and , and therefore . We shall adopt the approach that is bracketed by the range of neutron matter calculations performed by Refs. Gandolfi12 () and Drischler14 (). As noted in Ref. Lattimer01 (), the neutron star matter pressure between and essentially determines the radii of neutron stars in the mass range to (as long as one assumes no strong phase transitions in this region, see e.g. Fig. 7). We therefore expect that the value of will play the same role in this EOS. Refs. Gandolfi12 () showed that the neutron matter energy for densities less than about was adequately approximated by the double power law

(13) |

where and are parameters. For this energy formula, and assuming the validity of a quadratic symmetry energy interpolation (Eq. 1), we immediately find that and , where MeV is the bulk binding energy of symmetric matter. Table 1 displays parameter values found by Ref. Gandolfi12 () for quantum Monte Carlo neutron matter calculations (see also Ref. Steiner12cn ()). We have also displayed parameter values that fit the neutron matter results of Ref. Drischler14 () up to densities , and we will assume these calculations can be extrapolated to the slightly higher density . For each set, the corresponding values of , , and the pressure at have been tabulated. In the quadratic approximation for the isospin dependence of the nucleon energies, , so and are highly correlated (Fig. 1). The piecewise polytrope approximation for the EOS, assuming a value for , explicitly predicts a value for ,

(14) |

with given by Eq. (10). These values are close to, but generally larger than, the corresponding neutron matter predictions for a given value of , indicating that the piecewise polytropic EOS is a reasonable approximate in this density range.

Since the parameter set GCR 0 corresponds to an EOS with no three-body interactions, the range

(15) |

is predicted from realistic neutron matter studies. Assuming no strong phase transitions at lower densities, we expect the minimum radius limit, to neutron stars, which primarily depends on , from neutron matter constraints to be significantly larger than the absolute minimum limit, km Lattimer12 () established from the “maximally compact” EOS of Ref. Koranda97 () and an observed neutron star. will also play an indirect role in determining , the maximum neutron star mass, in that stars with larger radii at intermediate densities typically support larger masses. On the other hand, we expect that and will play little role in determining neutron star radii but will play significant roles in determining .

Hydrodynamic stability, causality, and observed neutron star masses impose important constraints on the choices of parameters. In this scheme the sound speed increases monotonically with density within each polytrope, so that causality within regions 1 to 3 requires

(16) |

which gives implicit equations for since

(17) |

Therefore, depends upon and , and depends upon and , but depends on . In general, is so much larger than the realistic ranges of established from neutron matter studies that the causality condition is not important to our studies ( MeV fm from causality considerations (Fig. 2)). depends upon and , and therefore also upon . However, the dependence upon is weak, as shown in Fig. 3. In addition, the central density of the star is, in many cases, larger than , in which case the limiting value is smaller than the limit given by Eq. (16). The actual limit must be found numerically from TOV integrations of the star’s structure requiring that the maximum sound speed at the center of the maximum mass star be smaller than .

Minimum values for , and/or also exist in order to satisfy hydrodynamic stability, which requires that . This is obviously satisfied by any realistic value of . Parameter ranges allowed by hydrodynamical stability and causality are portrayed in Figs. 2 and 3 as the white regions. For and , more restrictive minima can result from the constraint that the maximum mass exceeds the largest well-measured neutron star mass, which we take to be , the lower limit to the measured mass of PSR J0548+0432 Antoniadis13 (). (Note that there is no minimum value for based on this condition, due to the presence of the polytropic regions 2 and 3.) These lower limits also must be found numerically from TOV integrations, which indicates that the effective lower limit to is approximately 100 MeV fm for virtually all realistic choices of (Fig. 2).

The result of each TOV integration with a different EOS (i.e., different combinations of and ) is indicated by a symbol in Figs. 2 and 3 (many parameter combinations yield nearly identical configurations and cannot be distinguished). Solid circles show parameters that support causal configurations with, respectively, (black), (orange), (blue) and (red). Parameters that yield acausal configurations, those in which the sound speed at the center exceeds , or configurations incapable of supporting at least , are shown as teal crosses. It is clear that causal configurations capable of supporting must have MeV fm, and if , MeV fm. On the other hand, the specific value of plays relatively little role as long as lead to causal configurations of the required mass.

We note that the parameter boundaries due to causality we find are in disagreement with the results of Ref. Ozel15a () in spite of the fact that the parameterized EOSs in the two studies are identical. For example, when it is assumed that , Ref. Ozel15a () indicate that , while we find ; when it is assumed that , Ref. Ozel15a () indicates that , while we find . We also note that the solution preferred by the Bayesian analysis of Ref. Ozel15a () has optimum parameters and that lie very near the acausal boundary, and that their parameter region with finite likelihoods extends beyond this boundary.

### 3.3 Alternative Models for High-Density Matter

Choosing different models can be thought of as choosing different prior distributions, as discussed in Ref. Lattimer14co (). Thus, in order to understand how our results vary with different prior distributions, we also employ two models from Ref. Steiner15un (). The crust EOS is described in the same way using the results of Ref. Baym71 () but is supplemented with the results of Ref. Negele73 () at the highest densities in the crust. From to , the GCR EOS is used (see above). At higher densities, either Model A or Model C from Ref. Steiner13tn () is used. Model A is built on piecewise polytropes and gives very similar results to our first model described above. Model C is constructed using line segments in the - plane and prefers stronger phase transitions (this model is called “qmc_fixp” in Ref. bamr ()). In particular, Model C allows strong phase transitions just above . This can make a significant difference in the results, as shown below. It is clear that exotic matter cannot appear at since laboratory nuclei are known to consist only of neutrons and protons. The minimum possible density at which a phase transition may appear, however, is not well-determined from either theory or experiment.

We can also investigate the observed differences between Models A and C within the piecewise polytropic method described in Section 3.2 by explicitly creating a first-order phase transition between the fiducial densities and and allowing those densities to vary. One requires that and , where is the chemical potential. The maximum effect of the phase transition results if the EOS above is given by the causal limit, . In this causality-limited case, one can show that the phase transition strengths and are related by:

(18) |

## 4 Results

### 4.1 Masses, Radii, and the EOS

For the polytropic EOS from section 3.2, TOV integrations were computed with a grid of values for , and . Generally, we assume to be limited by the range described above for realistic neutron matter studies. However, we slightly extended both the lower and upper limits of to be conservative. Following Ref. Ozel15a (), we choose the value MeV fm, which proves crucial to the computation of the minimum realistic neutron star radius. We increased the upper limit to be about 1.6 times larger than , which proves crucial in establishing a maximum realistic neutron star radius. Values of were chosen to be smoothly distributed in between MeV fm (to avoid unnecessary computations of configurations with ) and (determined from causality considerations). Values for were taken to be smoothly distributed in between (hydrodynamic stability) and (causality up to ).

We also use Model C EOS as described in section 3.3. The calculation proceeds as for the polytropic EOS described above, except that the results for Monte Carlo realizations are histogrammed and binned (using Ref. o2scl ()). In effect, because Model C prefers stronger phase transitions, comparing with the polytropic EOS compares two different prior assumptions. The correlation between the radius of 1.4 and the pressure, is plotted for the polytropic EOS in Fig. 4, and for Model A and Model C in Fig. 5. As expected, a large degree of correlation exists between the radii of configurations, , and in either case. Eliminating acausal configurations, or configurations incapable of supporting , which are shown in Fig. 4 as green circles, the correlation becomes more robust. The correlation between and is best described with a quadratic contribution, as can be seen in Figs. 4 and 5. The correlation between and is

(19) |

where is in units of MeV fm.

Ref. Lattimer01 () had shown the existence of a phenomenological relation between and , which was later modified Lattimer13 () to consider only EOSs capable of supporting :

(20) |

The polytrope results are consistent with this relation (Fig. 6), where we approximated from the piecewise polytropic EOSs, as can readily be expected given the high degree of correlation between and or (Fig. 1). The tendency for our predicted radii to be slightly smaller than those predicted by Eq. (20) is due to the fact that both the lower and upper limits to are smaller than for the range of the EOS samples considered by Refs. Lattimer10 () and Lattimer13 (). However, this correlation is strongly sensitive to the prior assumptions, i.e. the possible presence of a phase transition just above the saturation density. This is demonstrated by the stark difference between the polytrope-based models and Model C (compare Figs. 6 and 7.)

The stark difference in permitted values of is readily understood if one permits strong phase transitions and a high-density causality-limited EOS in the piecewise polytrope model by permitting and to vary, as described in Sec. 3.3. There it was shown, Eq. (18), that the strength of the phase transition is described by , or, equivalently, by . The range of permitted values of is extremized when . The largest value, km, results when . The smallest value, km, results when , and is nearly as small as allowed by general relativity and causality, 8.15 km Lattimer12 (), in the case when . The lower limit to increases for larger values of (accompanied by a decrease in ). Furthermore, a phase transition is only allowed, assuming , if ; further decreases in permitted values of result for larger values of . These conclusions are similar to those of Ref. Alford15 ().

There is a lack of correlation between and . Instead, a high degree of correlation exists between and (Figs. 8 and 9). The presence of this correlation appears to be relatively prior-independent. Model C, which allows for stronger phase transitions, also tends to allow EOSs which have higher pressures at high densities, and thus gives relatively larger maximum masses. That is, the posterior distribution of is strongly prior dependent.

Ref. Lattimer12 () determined the minimum radii of stars as a function of the minimum value of using the maximally compact EOS () from Ref. Koranda97 (). The relation between and for the piecewise polytropic EOS is shown in Fig 10. Configurations that are causal and capable of supporting at least have km (10.85 km for the minimum value of from neutron matter calculations). As the maximum mass limit is raised, the lower limit for in causal configurations slowly increases. For , it is necessary that km and MeV fm. The relation between and for Models A and C are plotted in Fig. 11, demonstrating the dependence on the prior distribution.

Note that causality and the constraint that strongly limits configurations with radii less then 10 km, as shown in Figs. 10 and 11. In the polytrope model, it is important to emphasize the parameters leading to these acausal configurations satisfy causality at the boundary points and . Acausality occurs at densities larger than but below the central density of maximum mass configurations. Thus, the acausality boundary shown for as a function of indicated in Fig. 3 is effectively extended to lower values of at lower values of , and also implies a lower limit to when combined with the constraint .

The minimal radius limit from the maximally compact EOS () is shown as the solid curve in Fig. 10 which excludes the green shaded region. This boundary approximately represents an extreme limit because most EOSs have larger radii. It was recently demonstrated that existence of neutron stars implies that the sound speed must be larger than Alford13 (); Bedaque15sv () (the dotted line in Fig. 10), and most of our parameterized EOSs also exceed this limit.

In the polytrope model, the realistic upper limit to from neutron matter theory sets an interesting upper limit to of about 13 km. With the present value of , neutron matter constraints on restrict neutron star radii for stars to lie in the narrow range 11 km km. It is interesting to examine the frequency distribution of among the models that satisfy (Fig. 12). Although the minimum radius is about 10.6 km, assuming that the parameters and have equal likelihoods within their ranges implies that it is highly unlikely that km. On the other hand, Model C which allows strong phase transitions also allows the pressure to increase quickly above the saturation density. In turn, this allows radii as large as 14.5 km. The corresponding cumulative probability distribution is given in Fig. 13. It should be noted however, it is unclear what physical mechanism would give rise to a strong increase in the pressure just above the nuclear saturation density. A strongly repulsive four-neutron force, for example, seems unlikely.

The restrictions of causality and large maximum masses severely restrict the allowed EOSs. Fig. 14 shows boundaries in the pressure-energy density plane with different assumptions for permitted by causality and the assumed low-density EOS for the crust and for neutron matter. For , the maximum uncertainty in pressure for a given energy density is no larger than a factor of 3 (which occurs near ), and is slightly larger than a factor of 2 near the central densities of maximum mass stars. The corresponding regions in the mass-radius plane that can be populated by EOSs satisfying the , causality and the low-density EOS constraints are shown in Fig. 15. These figures show the importance of neutron star mass measurements: the larger the minimum value of , the more restricted are the ranges of and and the more accurately the EOS can be predicted. The corresponding regions in the mass-radius phase for Model C are plotted in Fig. 16, where larger radii are allowed as discussed above.

It is interesting to compare our results with those of a recent study in Ref. Chen15 (). In that study, various - curves were assumed. Beginning with a relativistic mean field EOS (FSU-Garnet) that predicts km and , new curves were generated by arbitrarily translating the original curve for (corresponding to ) to smaller radii by discrete amounts. The EOS corresponding to each newly-generated curve was deduced by the inversion technique from Ref. Lindblom92 (). It was shown that was possible only if km, a result very similar to that of the present study if phase transitions are not allowed. However, it can be argued that the prescription of Ref. Chen15 () is model-dependent, being sensitive to the choices of and the fiducial EOS (FSU-Garnet). In addition, the prescription seems not consistently applied, as an abrupt reduction of at implies a discontinuity in the - curve and a first-order phase transition in . Rather, Ref. Chen15 () alters the low-density EOS for so that the values of can be reduced from the original value ( km) to values as small as 9 km, which seemingly results in higher pressures for than the original EOS. For , the pressure can be 20–30% larger than the original EOS, and the amplification grows with decreasing density. This is incompatible with our knowledge of the crust’s EOS which allows no such deviations.

### 4.2 Universal Relations

There have been shown to exist several relatively EOS-independent relations among neutron star observables. These may be very useful to reduce degeneracies in interpretations of observations, including those from gravitational radiation Yagi13 (). Ref. Lattimer89 () found a relation between the binding energy , where we set , is the number of nucleons in the star, and the compactness parameter ; this was later improved by Ref. Lattimer01 (), who found

(21) |

for a wide variety of EOSs which could support at least . Later, Ref. Lattimer05s () found another EOS-independent relation concerning the moment of inertia for EOSs which could support approximately the same :

(22) |

These relations are compared to results from piecewise polytropic EOSs that satisfy causality and various values of in Figs. 17 and 19. The results are reasonably well approximated by Eqs. (21) and (22) except for , i.e., close to the maximum masses. We find an improved approximation for the moment of inertia, showing less uncertainty, especially for compactnesses typical of stars, is given by

(23) | ||||

(24) |

The smaller uncertainties result from assuming that . It is apparent, however, from both Figs. 17 and 19 that should further observations increase the value of uncertainties in analytic approximations for BE and can be substantially reduced.

Ref. Yagi13 () found an even more remarkable EOS-independent relation relating the moment of inertia, the tidal Love number and the quadrupole polarizability, which is now known as the I-Love-Q relation. The correlation between the dimensionless moment of inertia, , and the dimensionless tidal Love number, , is shown in Fig. 21 for the piecewise polytropic EOSs. The relation for Yagi13 () and its inverse, , are

(25) | ||||

(26) | ||||

(27) | ||||

(28) |

The deviation of piecewise polytropic EOSs from these analytic approximations is negligible, as seen in Fig. 21. Combining Eq. (28) with the approximation Eq. (22) yields an explicit approximation for , as shown in Fig. 22. The dashed lines show the combined approximations from Eqs. (28) and (22). It is apparent that uncertainties in a revised correlation would be substantially reduced if observational constraints for were taken into account.

It is interesting to examine the EOS dependence of the moment of inertia relative to the binding energy (Fig. 24). The two quantities are not as highly correlated as are and , but the correlation is still significant. The approximation shown in Fig. 24, using , is

(29) | ||||

(30) |

The improvement resulting from inclusion of maximum mass constraints, compared to the correlation inferred from Eqs. (22) and (21), is seen to be substantial, especially for compactnesses characteristic of stars.

## 5 Discussion

In the future, astronomical observations of neutron stars may be sufficiently plentiful and precise that one will be able to map out the mass-radius curve and the EOS using simpler methods like a analysis. The currently available observational data does not permit this. The problem of determining the mass-radius curve is underconstrained, in part because the data is not yet precise, and in part because there are several remaining systematic uncertainties (or several model assumptions).

A similar difficulty exists on the side of experimental nuclear physics. Even though isospin-asymmetric matter near the saturation density is well-known, the nature of neutron-rich matter near the nuclear saturation density is still subject to uncertainties such as the nature of the three-neutron force and the highest density to which we can trust chiral effective theories (see recent progress on this front in Ref. Furnstahl15 (). Experimental observables which probe high-density matter, such as intermediate-energy heavy-ion collisions (e.g. Ref. Tsang09co ()), are also subject to strong systematics. As we have noted above, the lowest density at which a phase transition to exotic matter is possible is a critical quantity necessary for understanding neutron stars. This lowest density limit for exotic matter is not well-known.

In the context of Bayesian inference, we approached these issues by varying our prior assumptions. To the extent that we are able to express the various systematic uncertainties in terms of prior probability distributions, Bayesian inference provides us a way to quantify the uncertainty in the mass-radius curve and the EOS of dense matter. This method also allows us to critically examine the extent to which correlations between observables exist.

We found, in particular, that the possible presence of phase transitions at low-density has an important impact on the lower-limit for the radius of low-mass neutron stars. We showed that some universal relations (or correlations) are not strongly sensitive to the prior assumptions, such as the relation between the radius of a neutron star, the pressure at , and the correlation between and the pressure at . We also showed that the high degree of correlation between the moment of inertia and the tidal Love number (the I-Love correlation) is robust with respect to prior assumptions concerning the equation of state. Finally, we presented a new universal relation connecting the neutron star binding energy to its moment of inertia.

###### Acknowledgements.

AWS was supported by the U.S. DOE Office of Nuclear Physics. JML was supported by the U.S. DOE grant DE-AC02-87ER40317 and the DOE Topical Collaboration “Neutrinos and Nucleosynthesis in Hot and Dense Matter”. EFB was supported by the NSF under grant No.’s PHY-1430152 (JINA Center for the Evolution of the Elements) and AST 11-09176. This project used computational resources from the University of Tennessee and Oak Ridge National Laboratory’s Joint Institute for Computational Sciences.## References

- (1) R. C. Tolman, Phys. Rev. 55, 364 (1939).
- (2) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
- (3) A. G. W. Cameron, Astrophys. J 130, 884 (1959).
- (4) T. H. R. Skyrme, Nucl. Phys. 9, 615 (1959).
- (5) A. W. Steiner, M. Prakash, J. M. Lattimer, and P. Ellis, Phys. Rep. 411, 325 (2005).
- (6) E. S. Phinney and S. R. Kulkarni, Ann. Rev. Astron. & Astrophys. 32, 591 (1994).
- (7) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature 467, 1081 (2010).
- (8) J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, Science 340, 1233232 (2013).
- (9) J. M. Lattimer, Ann. Rev. Nucl. Part. Sci. 62, 485 (2012).
- (10) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001).
- (11) J. van Paradijs, Astrophys. J. 234, 609 (1979).
- (12) F. Ozel, Nature 441, 1115 (2006).
- (13) T. Güver, F. Özel, A. Cabrera-Lavers, and P. Wroblewski, Astrophys. J. 712, 964 (2010).
- (14) T. Güver, P. Wroblewski, L. Camarota, and F. Özel, Astrophys. J. 719, 1807 (2010).
- (15) F. Özel, T. Güver, and D. Psaltis, Astrophys. J. 693, 1775 (2009).
- (16) J. van Paradijs, F. Verbunt, R. A. Shafer, and K. A. Arnaud, Astron. & Astrophys. 182, 47 (1987).
- (17) L. Bildsten, E. E. Salpeter, and I. Wasserman, Astrophys. J. 384, 143 (1992).
- (18) R. Rutledge, L. Bildsten, E. Brown, G. Pavlov, and E. Zavlin, Astrophys. J. 514, 945 (1999).
- (19) M. Rajagopal and R. W. Romani, Astrophys. J. 461, 327 (1996).
- (20) V. E. Zavlin, G. G. Pavlov, and Y. A. Shibanov, Astron. & Astrophys. 315, 141 (1996).
- (21) C. O. Heinke, G. B. Rybicki, R. Narayan, and J. E. Grindlay, Astrophys. J. 644, 1090 (2006).
- (22) N. A. Webb and D. Barret, Astrophys. J. 671, 727 (2007).
- (23) L. Lindblom, Astrophys. J. 398, 569 (1992).
- (24) L. Lindblom and N. M. Indik, Phys. Rev. D 89, 064003 (2014).
- (25) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013).
- (26) F. Özel and D. Psaltis, Phys. Rev. D 80, 103003 (2009).
- (27) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009).
- (28) Y. Eriguchi and E. Mueller, Astron. & Astrophys. 146, 260 (1985).
- (29) F. Özel, G. Baym, and T. Güver, Phys. Rev. D 82, 101301 (2010).
- (30) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Phys. Scr. T56, 231 (1995).
- (31) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. 722, 33 (2010).
- (32) V. Suleimanov, J. Poutanen, and K. Werner, Astron. Astrophys. 527, A139 (2011).
- (33) M. Servillat, C. O. Heinke, W. C. G. Ho, J. E. Grindlay, J. Homg, M. van den Berg, and S. Bogdanov, Mon. Not. R. Astron. Soc. 423, 1556 (2012).
- (34) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Phys. Rev. Lett. 105, 161102 (2010).
- (35) A. W. Steiner and S. Gandolfi, Phys. Rev. Lett. 108, 081102 (2012).
- (36) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. Lett. 765, 5 (2013).
- (37) S. Guillot, M. Servillat, N. A. Webb, and R. E. Rutledge, Astrophys. J. 772, 7 (2013).
- (38) J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A 50, 40 (2014).
- (39) C. O. Heinke, H. N. Cohn, P. M. Lugger, N. A. Webb, W. C. G. Ho, J. Anderson, S. Campana, S. Bogdanov, D. Haggard, A. M. Cool, and J. E. Grindlay, Mon. Not. Royal Astron. Soc. 444, 443 (2014).
- (40) A. W. Steiner, J. Phys. G 42, 034004 (2015).
- (41) F. Ozel, D. Psaltis, T. Güver, G. Baym, C. Heinke, and S. Guillot, arXiv:1505.05155 (2015).
- (42) S. Guillot and R. E. Rutledge, Astrophys. J. Lett. 796, L3 (2014).
- (43) A. W. Steiner, S. Gandolfi, F. J. Fattoyev, and W. G. Newton, Phys. Rev. C 91, 015804 (2015).
- (44) J. Nättilä, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, arXiv:1509.06561 (2015).
- (45) F. Özel and D. Psaltis, Astrophys. J. 810, 135 (2015).
- (46) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, in Bayesian Data Analysis, 3rd ed. (CRC Press, Boca Raton, 2014), Chap. 1.
- (47) S. Typel and B. A. Brown, Phys. Rev. C 64, 027302 (2001).
- (48) C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001).
- (49) K. Yagi and N. Yunes, Science 341, 365 (2013).
- (50) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
- (51) S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C 85, 032801 (2012).
- (52) C. Drischler, V. Soma, and A. Schwenk, Phys. Rev. C 89, 025806 (2014).
- (53) S. Koranda, N. Stergioulas, and J. L. Friedman, Astrophys. J. 488, 799 (1997).
- (54) J. W. Negele and D. Vautherin, Nuc. Phys. A 207, 298 (1973).
- (55) A. W. Steiner, Astrophysics Source Code Library, 2014, record ascl:1408.020.
- (56) A. W. Steiner, Astrophysics Source Code Library, 2014, record ascl:1408.019.
- (57) J. M. Lattimer and Y. Lim, Astrophys. J. 771, 51 (2013).
- (58) J. M. Lattimer and M. Prakash, in From Nuclei To Stars: Festschrift in Honor of Gerald E Brown (World Scientific, Singapore, 2011), Chap. 12, arXiv:1012.3208.
- (59) M. G. Alford, G. F. Burgio, S. Han, G. Taranto, and D. Zappalà, Phys. Rev. D 92, 083002 (2015).
- (60) P. Bedaque and A. W. Steiner, Phys. Rev. Lett. 114, 031103 (2015).
- (61) W.-C. Chen and J. Piekarewicz, Phys. Rev. Lett. 115, 161101 (2015).
- (62) J. M. Lattimer and A. Yahil, Astrophys. J. 340, 426 (1989).
- (63) J. M. Lattimer and B. F. Schutz, Astrophys. J. 629, 979 (2005).
- (64) R. J. Furnstahl, D. R. Phillips, and S. Wesolowski, Jour. Phys. G. 42, 034028 (2015).
- (65) M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and A. W. Steiner, Phys. Rev. Lett. 102, 122701 (2009).