From Neutron Star Observables to the Equation of State. II. Bayesian Inference of Equation of State Pressures

TL;DR

This paper develops a Bayesian method to infer neutron star interior pressures at fixed densities from mass-radius data, highlighting the importance of regularization and full posterior analysis for accurate EoS characterization.

Contribution

It introduces a Bayesian framework for inferring neutron star EoS pressures at fixed densities, accounting for regularization and avoiding biases from marginalization.

Findings

Pressure at twice nuclear saturation density can be inferred within 0.3 dex.

EoS with phase transitions can be inferred within ~30%.

Full posterior analysis prevents radius bias of nearly 1 km.

Abstract

One of the key goals of observing neutron stars is to infer the equation of state (EoS) of the cold, ultradense matter in their interiors. We present here a Bayesian statistical method of inferring the pressures at five fixed densities, from a sample of mock neutron star masses and radii. We show that while five polytropic segments are needed for maximum flexibility in the absence of any prior knowledge of the EoS, regularizers are also necessary to ensure that simple underlying EoS are not over-parametrized. For ideal data with small measurement uncertainties, we show that the pressure at roughly twice the nuclear saturation density, rho_sat, can be inferred to within 0.3 dex for many realizations of potential sources of uncertainties. The pressures of more complicated EoS with significant phase transitions can also be inferred to within ~30%. We also find that marginalizing the multi-dimensional parameter space of pressure to infer a mass-radius relation can lead to biases of nearly 1 km in radius, towards larger radii. Using the full, five-dimensional posterior likelihoods avoids this bias.

Equations

$$P(\rho )=K_i{\rho }^{{\Gamma }_i},{\rho }_{i-1}<\rho <{\rho }_i,$$

$$K_i=\frac{P_{i-1}}{{\rho }_{i-1}^{{\Gamma }_i}}=\frac{P_i}{{\rho }_i^{{\Gamma }_i}}$$

$${\Gamma }_i=\frac{{\log }_{10}(P_i/P_{i-1})}{{\log }_{10}({\rho }_i/{\rho }_{i-1})}.$$

$$\frac{\text{d}P}{\text{d}r}=-\frac{G}{c^2}\frac{(ϵ+P)(M+4\pi r^{3}P/c^2)}{r^2-2GMr/c^2}$$

$$\frac{\text{d}M}{\text{d}r}=\frac{4\pi r^2ϵ}{c^2},$$

$$d\frac{ϵ}{\rho }=-Pd\frac{1}{\rho }.$$

$$P(EoS|data)=P(P_1,P_2,P_3,P_4,P_5|\mathrm{data}).$$

$$\begin{gathered} P(P_1,P_2,P_3,P_4,P_5|\mathrm{data})= \\ CP(\mathrm{data}|P_1,...,P_5)\times P_p(P_1,...,P_5), \end{gathered}$$

$$P(\mathrm{data}|P_1,...,P_5)=\prod _{i=1}^N{P}_i(M_i,R_i|P_1,...,P_5)$$

$$P_i(M_i,R_i|P_1,...,P_5)=P_{\max }(M_i,R_i|P_1,...,P_5,{\rho }_c),$$

$$P_{i}(M_{i},R_{i}|P_{1},...,P_{5},\rho_{c})=\frac{1}{2\pi\sigma_{R_{i}}\sigma_{M_{i}}}\\ \exp\Bigl{\{}-\frac{[M_{i}-M_{\mathrm{EoS}}(\rho_{c})]^{2}}{2\sigma_{M_{i}}^{2}}-\frac{[R_{i}-R_{\mathrm{EoS}}(\rho_{c})]^{2}}{2\sigma_{R_{i}}^{2}}\Bigr{\}},$$

$$P_i\le P_{i+1}.$$

$$\frac{dP}{dϵ}=\frac{c_s^2}{c^2}\le 1,$$

$$\xi =\exp \left[-\frac{(d^2(\ln P)/d(\ln \rho {)}^2{)}^2}{2{\lambda }^2}\right],$$

$$\begin{gathered} \int \int \int \int \int P(P_1,...,P_5|\mathrm{data})d{P}_{1}d{P}_{2}d{P}_{3}d{P}_{4}d{P}_5 \\ =0.68, \end{gathered}$$

$$R=R({\rho }_c,P_1,...,P_5)$$

$$M=M({\rho }_c,P_1,...,P_5).$$

$$R=R(P_1,...P_5;M),$$

$$\begin{gathered} P(R;M)=\int P\left[R(P_2,...,P_5;M),P_2,...,P_5\right] \\ \times \hat{J}\left(\frac{R}{P_1}\right)d{P}_2...d{P}_5, \end{gathered}$$

Full text

1. Introduction

Observations of neutron stars provide a unique way to study the equation of state (EoS) of cold, ultradense matter at high densities. Many theoretical EoS have been proposed, incorporating a wide variety of physics and calculation methods. The sample of proposed EoS ranges from purely nucleonic models (e.g., Baym et al. 1971; Friedman & Pandharipande 1981; Akmal et al. 1998; Douchin & Haensel 2001), to models including hyperons (e.g., Balberg & Gal 1997), pion condensates (e.g., Pandharipande & Smith 1975), or kaon condensates (e.g., Kaplan & Nelson 1986). A number of newer studies incorporate the effects of the additional quark degrees of freedom expected at high densities, either from phenomenological models or from the early results of lattice QCD (e.g., Alford et al. 2005, 2013; Kojo et al. 2015).

While nuclear physics experiments can be used to inform and constrain these models at low densities, experimental data become sparse even in the vicinity of the nuclear saturation density, $\rho_{\rm sat}~{}\sim~{}2.7~{}\times~{}10^{14}$ g cm*-3*, for neutron-rich matter (see Lattimer 2012). Neutron stars, with central densities up to $8~{}\rho_{\rm sat}$, offer a way to constrain the EoS at higher densities. Any particular EoS predicts a unique combination of allowed neutron star masses and radii. Observations of these quantities can, therefore, be compared to the predictions of various EoS models to provide constraints.

The masses of numerous neutron stars have already been measured, of which at least a dozen have simultaneous radii measurements (Guillot et al. 2013; Guillot & Rutledge 2014; Heinke et al. 2014; Nättilä et al. 2015; Bogdanov et al. 2016; Özel et al. 2016; for a recent review, see Özel & Freire 2016). Some of these measurements already place stringent constraints on the EoS: for example, the observations of two, roughly 2 $M_{\odot}$ neutron stars (Demorest et al. 2010, see Fonseca et al. 2016 for the revised mass; Antoniadis et al. 2013) imply that the EoS must be relatively stiff in order to support the mass against gravitational collapse. However, this approach of using observations to constrain individual EoS models is inherently limited in scope: it is possible that the current sample of proposed EoS does not probe the full range of physical possibilities.

One alternate approach to constraining individual EoS models is to infer, instead, their functional form directly from the observations, by exploiting the one-to-one mapping of the EoS to the mass-radius curve. It has been shown that it is formally possible to invert that one-to-one mapping to recover the EoS, given sufficient mass and radius observations (Lindblom, 1992). Such an inversion requires observations that completely populate the mass-radius curve. However, because there is no known mechanism for producing neutron stars below $\sim 1M_{\odot}$ and none have been observed, performing a full inversion is observationally unrealistic.

A second approach is to infer a simpler functional form, using a parametric EoS. This requires fewer observations, allows for a more limited data range, and, at the same time, allows for different types of data, such as the moment of inertia, to be combined. Many parametrizations have been proposed. It has been shown, for example, that the EoS can be written as a spectral expansion in terms of the enthalpy (Lindblom & Indik, 2012, 2014). It has also been shown that the EoS can be represented as a discrete number of segments that are piecewise polytropic or linear (Read et al., 2009; Özel & Psaltis, 2009; Hebeler et al., 2010; Steiner et al., 2016). In Raithel et al. (2016), we performed an optimization of the segmented, parametric EoS and found that using five piecewise polytropic segments above the nuclear saturation density was necessary to reproduce the mass, radius, and moment of inertia for a variety of diverse EoS to within the expected uncertainties of next-generation experiments, i.e., to within 0.5 km, 0.1 $M_{\odot}$, and $\sim 10\%$, respectively.

Given a parametric EoS, the parameters can be determined from a Bayesian inference method. Such methods have been previously developed in Steiner et al. (2010, 2016) and Özel et al. (2016). While a Bayesian method has the potential of being robust, important statistical considerations have yet to be explored. For example, it is not obvious what types of prior distributions should be used; the effects of measurement uncertainties have not yet been quantified; and biases in currently-used methods of marginalization persist. These uncertainties are particularly important when there is an absence of data at low masses to constrain the solution.

In this paper, we perform the natural next step to Raithel et al. (2016), i.e., using Bayesian statistical techniques to infer the pressures of our optimized five-polytrope parameterization given samples of simulated astrophysical measurements. We review our optimal parametrization in $\S$2 and introduce the Bayesian inference method in $\S$3. In $\S$4, we test the Bayesian inference on a number of equations of state, ranging from simple polytropes to complex EoS with significant phase transitions.

We find that, if one wants to be fully agnostic about the physics of the EoS at high densities, the full five-polytrope parametrization must be used. This is the only way to ensure that complex structure and/or significant phase transitions are allowed in the inferred EoS. However, if the true EoS is relatively simple or even a single polytrope, the five-polytrope model will lead to over-parametrization. We find that a regularizer can help reduce the effects of the over-parametrization, while still allowing the complexity that is the benefit of the five-polytrope model.

Finally, in either case, it is important to use the full five- (or three-) dimensional posterior likelihoods rather than a marginalization, which, as we demonstrate in $\S$5, can introduce biases as large as 1 km.

2. The Optimal Parametric Equation of State

In order to infer the functional form of the dense-matter equation of state, we use a parametric EoS composed of piecewise polytropes. As we showed in Raithel et al. (2016), a parametric EoS with five polytropic segments is optimal for reproducing the observable properties of neutron stars for a diverse set of EoS, to within expected observational uncertainties.

We begin the parametric EoS at the nuclear saturation density, $\rho_{\rm sat}$, below which we assume a low-density EoS. We space the five polytropic segments evenly in the logarithm of density, so that the fixed densities separating the polytropes are at 1, 1.4, 2.2, 3.3, 4.9, and 7.4 $\times\rho_{\rm sat}$. The EoS across each segment is given by

[TABLE]

where the constant $K_{i}$ and the polytropic index, $\Gamma_{i}$ are determined by the pressure and density of the previous fiducial point, i.e.,

[TABLE]

and

[TABLE]

The above EoS can be mapped to a mass-radius curve using the Tolman-Oppenheimer-Volkoff (TOV) equations. The TOV equations relate the pressure, $P$, and enclosed mass, $M$, of the star as a function of radius, according to

[TABLE]

and

[TABLE]

where the energy density, $\epsilon$, is given by

[TABLE]

The TOV equations provide a one-to-one mapping of the EoS to mass-radius space. The goal of this paper is to demonstrate the reverse process: the inference of our parametric EoS, given simultaneous observations of masses and radii.

3. Bayesian Inference of EoS Pressures

In order to infer the pressures of our parametric EoS, we follow the Bayesian approach of Özel et al. (2016), which we recreate below. (For a similar analysis, see Steiner et al. 2010). The pressures of interest are those at our five fiducial densities, which, as described in $\S$ 2, completely determine our piecewise-polytropic EoS.

The posterior that a particular realization of our parametric EoS correctly describes a set of data can be written as

[TABLE]

By Bayes’ theorem, we can rewrite this as

[TABLE]

where $C$ is a normalization constant, $P_{p}(P_{1},...,P_{5})$ is the prior on the set of five pressures, and

[TABLE]

is the likelihood of a particular realization of $N$ total mass-radius observations, given a set of EoS parameters.

To calculate the likelihood of observing a particular value of $(M,R)$ given an EoS, we compute the probability that the observation is consistent with each point along the predicted M-R curve, and then take the maximum likelihood. That is, we calculate

[TABLE]

where we have used the central density, $\rho_{c}$, to parametrize the mass-radius curve. In the limit of small errors in either $M$ or $R$, this method is equivalent to taking the “closest approach” of the curve to the data point.

Finally, the likelihood that an observation of $(M_{i},R_{i})$ is consistent with a point on the mass-radius curve is given by

[TABLE]

where $\sigma_{R_{i}}$ and $\sigma_{M_{i}}$ are the measurement uncertainties associated with the radius and mass, respectively. Here, $R_{\rm EoS}(\rho_{c})$ and $M_{\rm EoS}(\rho_{c})$ are the radius and mass predicted by the set of pressures, $(P_{1},...,P_{5})$, that comprise our parametrized EoS for a particular central density, $\rho_{c}$. In order to populate the five-dimensional posterior of eq. (8), we use Markov-Chain Monte Carlo simulations following the Metropolis-Hastings algorithm.

3.1. Priors on the pressures

For the priors on $(P_{1},...,P_{5})$ in eq. (8), we employ constraints from physical principles, laboratory nuclear physics experiments, and astrophysical observations. Specifically, we require that:

(1.) The EoS be microscopically stable, i.e.,

[TABLE]

(2.) The EoS remain causal between the fiducial densities, i.e.,

[TABLE]

where $c_{s}$ is the local sound speed.

(3.) Each EoS produce a neutron star with a mass of at least 1.97 $M_{\odot}$, in order to be within $1~{}\sigma$ of the mass measurement of the most massive neutron stars (Antoniadis et al., 2013; Fonseca et al., 2016).

(4.) Pressures $P_{1}\geq 3.60$ MeV/fm3 and $P_{2}\geq 11.70$ MeV/fm3, in order to be consistent with nuclear physics experiments.

Constraint (4) provides lower limits on the pressures at the first two fiducial densities. As noted in Özel et al. (2016), important constraints on the EoS in the density regime around $\rho_{\rm sat}$ are obtained through nucleon-nucleon scattering experiments at energies below 350 MeV and from the properties of light nuclei. Results from such experiments can be extended by assuming two- and three-body potentials at densities near $\rho_{\rm sat}$ (Akmal et al., 1998). However, the interactions at higher densities cannot be written in terms of static few-body potentials. We, therefore, impose this nuclear physics constraint only on the pressures at our two lowest fiducial densities: $\rho_{1}=1.4~{}\rho_{\rm sat}$ and $\rho_{2}=2.2~{}\rho_{\rm sat}$. Because the three-nucleon interaction is always repulsive, the most model-independent lower limit uses only the two-nucleon interaction. Using the Argonne AV8 two-nucleon pressure as calculated in Gandolfi et al. (2014), we find $P(\rho_{1})$=3.60 MeV/fm3 and $P(\rho_{2})$=11.70 MeV/fm3. The AV8 potential is a simplified version of the Argonne AV18 potential (Wiringa et al., 1995); however, as noted in Özel et al. (2016), the two-nucleon interaction pressures are approximately the same for either version of the potential.

3.2. Regularizers

We also include a regularizer in our prior distributions, in order to reduce the tendency of our model to over-parameterize simple EoS. The regularizer, $\xi$, is Gaussian over the second logarithmic derivative of the EoS, i.e.,

[TABLE]

where $\lambda$ is the characteristic scale. We determine a suitable value for $\lambda$ by calculating $d^{2}(\ln{P})/d(\ln{\rho})^{2}$ at our fiducial densities for a sample of 49 proposed EoS. This sample of EoS was compiled from the literature in order to incorporate a wide variety of physics and calculation methods, as in Read et al. (2009). The cumulative distribution of second logarithmic derivatives for this sample are shown in Fig. 1. From this cumulative distribution, we find that 95% of the derivatives are $\lesssim~{}2$. We, therefore, use a characteristic scale of 4 $\times$ this value, resulting in $\lambda=8$ as our Gaussian regularizer. Such a regularizer will apply, at most, a penalty of $\sim$ 3% to the likelihood that we would calculate for an EoS that has second derivatives that occur in our sample of physically-motivated EoS.

4. Testing the Bayesian Inference with Mock Data

In this section, we test the Bayesian inference method described in $\S$3 using different sets of mock data. For most of the simulations described below, we assume an underlying EoS and create a realization of a sample of mass-radius data that are equidistant in mass between 1.2 and 2.0 $M_{\odot}$. We assume Gaussian measurement uncertainties with the same dispersions among the data points, which we denote by $\sigma_{\rm M}$ and $\sigma_{\rm R}$. We use a Monte Carlo method to draw a particular realization of simulated measurements from these distributions and apply our Bayesian inference method to this data set. We explore below how well the Bayesian inference method works for different types of underlying EoS, as well as for different number and quality of data points.

In Fig. 2, we show the result for one realization of mock data drawn from the EoS SLy, with uncertainties of $\sigma_{\rm R}=0.5$ km and $\sigma_{\rm M}=0.1~{}M_{\odot}$. EoS SLy is particularly challenging for a parametrization like ours that is optimized for potentially more complex EoS because it is practically a single polytrope in the density range of interest. We, therefore, use this example to explore the strengths and limits of the inference as well as of the regularizer.

The black lines in Fig. 2 represent the EoS SLy, while the magenta lines show the most-likely inferred EoS found with our Bayesian method. The gray bands represent the 68% credibility regions. For five-dimensional likelihoods, the 68% credibility region is defined as the region where

[TABLE]

exactly analogous to the lower-dimensional case. It should be noted that these credibility regions show the spread of possible solutions only, and should not be over-interpreted. That is, there are many curves that may be drawn through these regions that either have very small probabilities or violate one of our priors and are therefore unphysical. For example, the curve that could be drawn along the very edge of any of the bands is much less likely than those truly contained in the 68% credibility region. To emphasize this point further, we show in Fig. 3 several of the individual mass-radius curves that contribute to the shape of the 68% credibility region in the top panel of Fig. 2. The jagged edge of the contour is created by multiple different mass-radius curves. Indeed, the curve that follows either edge in its entirety is not included in the 68% interval. The same can be said for all of the grayed 68% credibility regions shown throughout this paper, and thus interpretations of those regions should be made with caution.

In the top panel of Fig. 2, we show the results of the inversion using uniform priors in pressure in the absence of any regularizers. The inferred EoS contains several sharp transitions between different polytropic indices. The middle panel shows an inversion for identical data, but in which the Gaussian regularizer has been added to the uniform priors. The stepped behavior that was shown in the top panel is effectively eliminated by the regularizer. With the addition of the regularizer, our inferred, most-likely EoS closely follows SLy. The errors in pressure for our most-likely EoS are all less than 30% for this realization of mock data, while the errors for $P_{2}-P_{5}$ are 3$-$7%.

Because there is no physical motivation to assume that the prior distribution of the pressures is uniform, we also tested the inversion with priors that are uniform in the logarithm of pressure. The results of this test, with identical data as above and the Gaussian regularizer included, are shown in the bottom panel of Fig. 2.

The top panel of Fig. 2 shows that assuming a uniform prior introduces a preference towards higher pressures. On the other hand, assuming that prior is uniform in the logarithm of pressures biases the results toward lower pressures. The inversion is particularly sensitive to this bias in the low-mass/low-density region of the EoS, where we lack data. Our Gaussian regularizer helps reduce this bias, as shown in the middle and bottom panels of Fig. 2, in which the results are similar for the two types of priors when the regularizer is also included.

In the absence of a regularizer, the freedom introduced by using five polytropes in the parametrization combined with this sensitivity to prior distributions could lead to significantly skewed results, or even the false inference of a phase transition with a perceived high statistical confidence. It is, therefore, important to use a Gaussian regularizer on the second derivative, with characteristic scale $\lambda=8$, to avoid this sensitivity to over-fitting.

In Fig. 4, we explore potential correlations between the inferred parameters for this test. Specifically, we show the correlations for the inversion using uniform priors and the Gaussian regularizer (i.e., the middle panel of Fig. 2). There are slight anti-correlations between adjacent inferred pressures, due to fact that, even with the addition of the regularizer, we are still over-parametrizing our model while trying to fit the effectively single-polytrope EoS SLy. Overall, however, the pressures that are not adjacent are uncorrelated with one another. Özel & Psaltis (2009) showed that, for a three-polytrope parametrization, parametrizing with the pressures as free parameters instead of the polytropic indices reduced the correlations between the inferred values. Figure 4 shows that the low levels of correlations are maintained here.

In order to ensure that the Bayesian inference works well for other underlying EoS and that the regularizer does not adversely limit our ability to detect potential phase transitions, we tested the method on a number of different EoS as well. As an example,111See discussion around Fig. 8 for other examples of EoS with a wide range of predicted radii. we show in Fig. 5 the results of our method obtained for an underlying EoS with significantly more structure than SLy; specifically, we generated a mock EoS with an extreme change in the polytropic index (from $\Gamma=1$ to $\Gamma=5$) that occurs in between two of our fiducial densities. This EoS is shown in the black solid line of Fig. 5. This was again designed to challenge the inversion procedure, but in the opposite extreme from EoS SLy. Even in this case, the most likely solution still recovers all of the pressures to within $\sim$30%, and recovers $P_{2}$ to within 11%.

Some previous studies (e.g., Steiner et al. 2016) have suggested that parametrizing with polytropic segments disfavors phase transitions because polytropes naturally go through the origin. However, continuity between segments, as required in any reasonable parametrization, implies that the power law segments are never required to go through the origin. Small values of the exponent, and hence phase transitions, are thus fully allowed. Moreover, in Fig. 5, we show that our parametrization is able to recover a phase transition, even when the phase transition is offset from the fiducial densities.

Finally, in order to quantify the range of uncertainties in the inferred pressures for different statistical realizations of the mock measurements, we generate a large number of mock data sets drawn from the EoS SLy (as in Figs. 2$-$4) and applied our method to each set. We summarize these results in Fig. 6, which shows the cumulative distribution of errors in the most likely inferred pressures at each of our five fiducial densities. We find that the pressure at $\rho_{2}=2.2~{}\rho_{\rm sat}$ is the best constrained, with errors less than 15% in 95% of the realizations. The other four pressures have errors less than $\sim$20% in approximately half of the realizations, and are correct to within less than 0.3 dex in every realization.

5. Biases due to marginalization

The Bayesian inference scheme described in $\S$3 provides five-dimensional posteriors, $P(P_{1},...,P_{5}|\mathrm{data})$. While one might explore the marginalized distributions over each pressure, $P_{i}$, this can easily lead to misinterpretations. It is possible, for example, that just considering the most likely pressures from each marginalized distribution of $P_{i}$ will produce an EoS that violates one of our priors. The pressures are coupled to one another and marginalizing over any one dimension removes that dependence. For this reason, we exclusively use the five-dimensional likelihoods to interpret our results.

Another approach that is taken in some earlier studies (e.g., Steiner et al. 2010) is to marginalize the output of the Bayesian inference in mass-radius space, rather than over the pressures. This is done by creating a one-dimensional histogram of the radii over a fixed grid of masses, for all possible mass-radius curves. Analytically, we can derive such a marginalization by first writing the radius and mass as functions of the five pressures and a central density, i.e.,

[TABLE]

For a fixed mass, we can write the radius equation instead as

[TABLE]

which we can invert to recover $P_{1}$ in terms of $R$ and $(P_{2},...,P_{5})$ for a fixed mass. One can then express the radius as a function of only the mass by marginalizing across the other four pressures, i.e.,

[TABLE]

where $\hat{J}(R/P_{1})$ is the Jacobian that transforms from $P_{1}$ to $R$. This is the analytical equivalent of taking the one-dimensional histogram of radii over a grid of masses.

Equation (18) highlights the issues that marginalizing introduces. If the full posterior, $P\left[R(P_{2},...,P_{5};M),P_{2},...,P_{5}\right]$, is relatively flat, then the highly non-linear Jacobian will dominate the resulting marginalization. The marginalization is particularly sensitive to this bias when the data are sparse and have large errors. However, the marginalization can be skewed for any data set, if the posterior distribution is not sharply peaked enough to overcome the influence of the Jacobian.

Indeed, in Fig. 7, we show this effect for two different sets of data: one that is clustered around $R~{}\sim~{}10$ km and one that is clustered around $R~{}\sim~{}12$ km. Both inversions use 5 simulated $(M,R)$ data points, with masses that are spaced evenly between 1.2 and 1.8 $M_{\odot}$ and with measurement uncertainties of $\sigma_{R}$=1 km and $\sigma_{M}=0.1~{}M_{\odot}$. For the smaller-radii dataset, the marginalized solution is offset by $\gtrsim~{}0.5$ km from the data at all masses (upper left panel). The most likely solution, in contrast, goes right through the data (lower left panel). For the larger-radii dataset, the effect is less extreme but still present: the marginalized solution is shifted approximately 0.4 km to the right of the data. The most likely solution, again, goes through every data point. While the marginalized solution is indeed within the error bars of the data, the most likely solution recreates the data almost perfectly for a variety of data. Current radius data have even larger and often overlapping error contours that will make this effect hard to identify by eye. It is, therefore, extremely important that only the full five-dimensional likelihoods be used to identify EoS constraints.

It should be noted that the choice of priors does affect the size and direction of this bias. Here, we use a prior distribution that is uniform in pressure and includes the Gaussian regularizer, while still requiring our other physical constraints (e.g., causality, a 2 $M_{\odot}$ star, etc.). Using a prior that is uniform in the logarithm of pressure pushes the bias in the other direction, i.e. toward smaller radii, and also pushes our most likely solution in that direction. The effect of the prior is stronger here than in Fig. 2 because we have more sparse data.

Figure 8 shows the bias of the marginalized solutions for inversions of still more data sets, each of which used uniform priors and the Gaussian regularizer. This figure also shows that the effect persists whether three or five polytropes are included in the parametrization. Each inversion used 5 simulated $(M,R)$ data points clustered near a radius between $\sim 10-15$ km, with masses between $1.2-1.8~{}M_{\odot}$, as in Fig. 7. We find that the bias is strongest at small radii. For data centered near 9.6 $M_{\odot}$, the bias is 1.1 km for a three-polytrope parametrization and 0.6 km for our five-polytrope parametrization. We include the results for the three-polytrope parametrization in order to emphasize that this bias is not a result of our specific choice of parametrization, but is a problem stemming from the marginalization of posteriors using sparse data. Given the large biases that can be introduced by the marginalization, it is clear that the marginalized solution should not be trusted. This is particularly true when the data have large, overlapping errors which will make this effect difficult to identify. The maximum likelihood method that we have used throughout this paper does not suffer any such biases.

6. Conclusions

In this paper, we developed a Bayesian method that can be used to infer the pressures of a parametrized EoS from a set of neutron star masses and radii. We used a parametrization containing five polytropic segments starting at $\rho_{\rm sat}$, which is the form we found to be optimal in our previous work (Raithel et al., 2016).

We investigated the influence of various priors and measurement uncertainties on the inferred pressures. We found that the freedom introduced by using five polytropic segments in the parametrization (which is necessary to recreate the EoS using next-generation data to within expected uncertainties) caused the inferred EoS to be too stepped in behavior. We, therefore, introduced a regularizer over the second derivative of pressure, to mitigate that freedom without sacrificing the ability to remain agnostic about the dense matter EoS. Combined with uniform priors over the pressures, we were able to show that measurement uncertainties expected in the near future will allow an inference of the pressure at $\rho_{2}=2.2\rho_{\rm sat}$ to within $\sim$15%, for a simple EoS. We were able to recreate the pressures at $\rho_{3}-\rho_{5}$ to within $\sim$ 20% in approximately half of the realizations, and to within 0.3 dex for all realizations. For a more complicated EoS with a significant break in the power-law indices, we were able to infer the pressures at all densities to within $\sim 30$%.

Finally, we showed that determining the posterior via marginalization in mass-radius space may lead to significant biases. We found that, for data at small radii, the marginalized mass-radius curves can be biased by nearly $+1$ km. Previously published EoS inferences from neutron star radius measurements are likely to have been affected by this bias in studies that used such a marginalization. It is better, instead, to use the maximum of the five-dimensional likelihoods computed in the Bayesian inference method. The most likely solutions do not suffer any such bias and were able to recreate the true $M-R$ curve for data at any radius.

It is difficult to compare our results to previous attempts to infer the pressures of a parametrized EoS from data because other studies do not characterize the uncertainties in the same ways. For example, Steiner et al. (2016) report that the maximum uncertainty allowed in their inversion is a factor of 3 in the pressure, purely from the requirement of causality, from their assumed crust EoS, and from the requirement that an EoS produce a star with $M\geq 1.97~{}M_{\odot}$. They report that the maximum uncertainty is closer to a factor of 2 near the central densities of the maximum mass stars. However, these maximum uncertainties are due to the priors of their model only, and would likely be smaller with the inclusion of data.

Using our approach, we find that the EoS can be inferred to high accuracy with the expected quality of next-generation data. Given that the EoS is currently poorly constrained at high densities, the possibility of constraining it to within even 0.3 dex, and possibly to within 15% at 2.2 $\rho_{\rm sat}$, will allow significant advances in our understanding of the physics at work in the ultradense regime.

Acknowledgements. CR is supported by the NSF Graduate Research Fellowship Program Grant DGE-1143953. FO gratefully acknowledges a fellowship from the John Simon Guggenheim Memorial Foundation in support of this work. DP acknowledges support from the Radcliffe Institute for Advanced Study at Harvard University.