FROM NEUTRON STAR OBSERVABLES TO THE EQUATION OF STATE. I. AN OPTIMAL PARAMETRIZATION

A key goal of observations of neutron stars is to determine the equation of state (EOS) of cold, ultradense matter. The EOS at high densities results from the interactions of nucleons, quarks, and possibly other constituents and fully characterizes the microphysics of the neutron star interior. Over the years, the sample of proposed EOS has incorporated vastly different physics and, accordingly, has predicted a variety of different global properties for neutron stars (see Özel & Freire 2016 for a recent review). Some early studies employed a purely nucleonic framework (e.g., Baym et al. 1971; Friedman & Pandharipande 1981; Akmal et al. 1998; Douchin & Haensel 2001), while others explored the role of hyperons (e.g., Balberg & Gal 1997), pion condensates (e.g., Pandharipande & Smith 1975), and kaon condensates (e.g., Kaplan & Nelson 1986). It is even possible, albeit contrived, to construct EOS that generate parallel stable branches of neutron stars with the same masses, but different radii by introducing a first-order phase transition into the EOS of hybrid neutron stars (Glendenning & Kettner 2000; Blaschke et al. 2013). More recently, there have been studies that incorporate the expected quark degrees of freedom at high densities, either from phenomenological models or derived from early results of lattice QCD (e.g., Alford et al. 2005, 2013; Kojo et al. 2015).

Nuclear physics experiments are able to constrain these various EOS only up to densities near the nuclear saturation density, ${\rho }_{\mathrm{sat}}$ ∼ 2.7 × 10¹⁴ g cm⁻³ (see Lattimer 2012). In this regime, there have been attempts to convert experimental results more directly into physical constraints. Even though, at these densities, the interactions between particles can be formally described in terms of static, few-body potentials (see, e.g., Akmal et al. 1998; Morales et al. 2002; Gandolfi et al. 2012), there are large uncertainties arising from extrapolations to β-equilibrium and low-temperature matter. Furthermore, at densities well above ${\rho }_{\mathrm{sat}}$, the interactions between particles can no longer be expanded in these terms.

At densities beyond ${\rho }_{\mathrm{sat}}$, experimentally constraining the dense-matter EOS can be accomplished using observations of neutron stars and exploiting the direct mapping between the EOS and the mass–radius relation (Lindblom 1992). This can be achieved by comparing measured masses (Demorest et al. 2010; Antoniadis et al. 2013) and radii (Guillot et al. 2013; Guillot & Rutledge 2014; Heinke et al. 2014; Nättilä et al. 2015; Bogdanov et al. 2016; Özel et al. 2016; see Özel & Freire 2016 for a recent review) to the predictions of particular EOS. However, this approach is limited in scope, as it serves to constrain the parameters of already formulated EOS frameworks and may not necessarily span the full range of possibilities. In other words, it may not be capable of recreating the EOS from observations in a model-independent way.

A separate approach is to empirically infer the EOS from the combination of these radii and mass measurements by utilizing a model-independent parametrization of the EOS. To date, several parametrizations have been proposed. The EOS can be written as a spectral expansion in terms of the enthalpy (Lindblom & Indik 2012, 2014). Alternatively, the EOS can be represented as a discrete number of segments that are piecewise polytropic or linear (Özel & Psaltis 2009; Read et al. 2009a; Hebeler et al. 2010; Steiner et al. 2016).

Despite these various proposals, only limited systematic optimizations have been undertaken to determine the ideal number of segments and functional forms for a parametric EOS. Moreover, these earlier studies used an existing sample of theoretical EOS as benchmarks and sought to represent only this subset. Finally, many of these early studies only sought to reproduce the predicted radii to rather large uncertainties by today's standards, reflecting the available data at the time. However, the current sample of mass and radius measurements, as well as the anticipated moment of inertia measurement from the binary system PSR J0737–3039, which is expected to be measured with 10% accuracy in the next five years (Lyne et al. 2004; Kramer & Wex 2009), require reconstructions of the EOS to a much higher accuracy than ever before. Any useful parametrization must be able to recreate masses and radii to high precision with a minimum set of parameters.

We present here an optimization of a piecewise-polytropic parametric EOS. In order to test our parametrization as generally as possible, we generate mock EOS that incorporate physically extreme behavior. We then apply our parametrization to these mock EOS and determine the differences in masses, radii, and moments of inertia between the parametrized and the full EOS. We optimize the number of segments included in the parametric EOS by requiring these differences to be smaller than the expected accuracy of observations. We find that sampling the EOS with five fiducial densities that are evenly spaced in the logarithm of density (i.e., using five polytropes) recreates the radii of the assumed EOS to within ≲0.5 km for the extreme cases and to within ≲0.12 km for 95% of a sample of 42 proposed EOS with a wide range of input physics. Such a parametrization is also able to reproduce the maximum mass to within ≲0.04 ${M}_{\odot }$ and the moment of inertia to within ≲10% for 95% of the proposed EOS sample.

We parametrize the EOS in terms of n piecewise polytropes,¹ spaced between two densities, ρ_min and ρ_max. We define the dividing density and pressure between each piecewise polytrope to be ρ_i and P_i, respectively. Each polytropic segment is then given by

$$\begin{eqnarray}&&P={K}_{i}{\rho }^{{{\rm{\Gamma }}}_{i}}\quad ({\rho }_{i-1}\leqslant \rho \leqslant {\rho }_{i}),\end{eqnarray} \tag{ 1 }$$

where the constant, K_i, is determined from the pressure and density at the previous fiducial point according to

$$\begin{eqnarray}&&{K}_{i}=\displaystyle \frac{{P}_{i-1}}{{\rho }_{i-1}^{{{\rm{\Gamma }}}_{i}}}=\displaystyle \frac{{P}_{i}}{{\rho }_{i}^{{{\rm{\Gamma }}}_{i}}}\end{eqnarray} \tag{ 2 }$$

and the polytropic index for the segment, Γ_i, is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i}=\displaystyle \frac{\,{\mathrm{log}}_{10}({P}_{i}/{P}_{i-1})}{\,{\mathrm{log}}_{10}({\rho }_{i}/{\rho }_{i-1})}.\end{eqnarray} \tag{ 3 }$$

Figure 1 shows an example of piecewise polytropes over three density segments, with various values of their polytropic indices, Γ, to illustrate the behavior of Equations (1)–(3). Our primary goal in optimizing the parametrization is to reduce the errors in the prediction of observables (i.e., mass, radius, and moment of inertia) below a threshold that is comparable to the uncertainties in present or upcoming observations, while keeping the number of polytropic segments to a minimum. For the optimization process, we produce extreme, albeit physically allowed EOS between 1–8 ${\rho }_{\mathrm{sat}}$ to test how well our parametrizations reproduce observables with various numbers of polytropic segments included. Once our parametrization is optimized, we then apply it to more reasonable, physically motivated EOS to test its ability to recreate those as well.

Zoom In Zoom Out Reset image size

In addition to the number of polytropes to include in the parametrization, there are two other variables that we have to optimize: the density at which the parametrization should start and the spacing of the polytropic segments. For the question of where to start the parametrization, we explored starting at ρ₀ = 10¹⁴ g cm⁻³ as well as at ρ₀ = ${\rho }_{\mathrm{sat}}$. It is typically assumed that the EOS is known up to ${\rho }_{\mathrm{sat}}$. However, Lattimer & Prakash (2001) showed that for a sample of around 30 proposed EOS, the predicted pressures vary by a factor of 5 over the range 0.5 ${\rho }_{\mathrm{sat}}$ < ρ < ${\rho }_{\mathrm{sat}}$, even though these EOS are all meant to be consistent with nuclear physics experiments in this density regime. This is because the extrapolation of pressures from symmetric nuclear matter to neutron-rich matter is poorly constrained. Meanwhile, densities below ∼0.5 ${\rho }_{\mathrm{sat}}$ do not significantly affect the global properties of the star. Therefore, we allowed our parametrization to start at both 0.5 ${\rho }_{\mathrm{sat}}$ and ${\rho }_{\mathrm{sat}}$, in order to explore these two limits. As for the question of how to space the polytropic segments, the preferred option is to space the segments evenly in the logarithm of the density. A logarithmic spacing more finely samples the low-density region of the EOS, which is the region that most affects the resulting neutron star observables (Lattimer & Prakash 2001; Özel & Psaltis 2009; Read et al. 2009a). For completeness, we also explored a second possibility: spacing the fiducial densities between each polytropic segment linearly.

We found that the combination of starting the parametrization at ${\rho }_{\mathrm{sat}}$ and spacing the fiducial densities evenly in the logarithm resulted in the smallest errors in mass and radius. We, therefore, start the first polytrope at ρ₀ = ${\rho }_{\mathrm{sat}}$ and space the remaining fiducial densities evenly in the logarithm between ρ₀ and 7.4 ${\rho }_{\mathrm{sat}}$. We set the last point, ρ_n = 7.4 ${\rho }_{\mathrm{sat}}$, following the results of Read et al. (2009a) and Özel & Psaltis (2009) who found that the pressure at this density determines the neutron star maximum mass and that pressures at higher densities do not significantly affect the overall shape of the resulting mass–radius curve. We determined the pressure corresponding to each fiducial density by sampling whichever EOS we were parametrizing, i.e., P_i = P_EOS(ρ_i). For ρ ≤ ρ₀, we connected our parametrization to a low-density EOS.

We varied the total number of fiducial densities above ${\rho }_{\mathrm{sat}}$ from 3–12. Clearly, as the number of polytropes used to represent the EOS increases, the errors in the observables are expected to reduce. Our goal in the remainder of this paper is to determine the minimum number of fiducial densities required to reproduce the mass, radius, and moment of inertia of a neutron star to within the desired observational uncertainties.

In order to determine how well our parametrizations were able to reproduce the observations predicted by a given EOS, we used the Tolman–Oppenheimer–Volkoff (TOV) equations and solved them to find the mass, radius, and moment of inertia.

The TOV equations give the pressure, P, and the enclosed mass, M, of the star as a function of radius, according to

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}=-\displaystyle \frac{G}{{c}^{2}}\displaystyle \frac{(\epsilon +P)(M+4\pi {r}^{3}P/{c}^{2})}{{r}^{2}-2{GMr}/{c}^{2}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dr}}=\displaystyle \frac{4\pi {r}^{2}\epsilon }{{c}^{2}},\end{eqnarray} \tag{ 5 }$$

where the energy density, , is given by

$$\begin{eqnarray}&&d\displaystyle \frac{\epsilon }{\rho }=-{Pd}\displaystyle \frac{1}{\rho }.\end{eqnarray} \tag{ 6 }$$

To get the full relation between energy density and mass density, we can integrate Equation (6) for ${\rm{\Gamma }}\ne 1$ to

$$\begin{eqnarray}&&\epsilon (\rho )=(1+a)\rho {c}^{2}+\displaystyle \frac{K}{{\rm{\Gamma }}-1}{\rho }^{{\rm{\Gamma }}},\end{eqnarray} \tag{ 7 }$$

where a is an integration constant. Along any density section of the EOS, requiring continuity at either endpoint determines a such that Equation (7) becomes

$$\begin{eqnarray}\epsilon (\rho ) & = & \left[\displaystyle \frac{\epsilon ({\rho }_{i-1})}{{\rho }_{i-1}}-\displaystyle \frac{{P}_{i-1}}{{\rho }_{i-1}({{\rm{\Gamma }}}_{i}-1)}\right]\rho +\displaystyle \frac{{K}_{i}}{{{\rm{\Gamma }}}_{i}-1}{\rho }^{{{\rm{\Gamma }}}_{i}},\\ & & ({\rho }_{i-1}\leqslant \rho \leqslant {\rho }_{i})\end{eqnarray} \tag{ 8 }$$

where K_i and Γ_i are determined as in Equations (2) and (3).

Similarly, for the case of Γ = 1, Equation (6) becomes

$$\begin{eqnarray}\epsilon (\rho ) & = & \displaystyle \frac{\epsilon ({\rho }_{i-1})}{{\rho }_{i-1}}\rho +{K}_{i}\mathrm{ln}\left(\displaystyle \frac{1}{{\rho }_{i-1}}\right)\rho -{K}_{i}\mathrm{ln}\left(\displaystyle \frac{1}{\rho }\right)\rho \\ & & {\rho }_{i-1}\leqslant \rho \leqslant {\rho }_{i}.\end{eqnarray} \tag{ 9 }$$

We used Equations (8) or (9) to relate an EOS to the energy density, and then used that energy density to integrate the TOV equations outwards from the center of the star. The radius at which the pressure becomes negligible gives the total mass and radius of the star.

In order to calculate the moment of inertia, we simultaneously solved Equations (4) and (5) with two coupled differential equations for the relativistic moment of inertia,

$$\begin{eqnarray}&&\displaystyle \frac{{dI}}{{dr}}=\displaystyle \frac{8\pi }{3}\displaystyle \frac{(\epsilon +P)}{{c}^{2}}\displaystyle \frac{{{fjr}}^{4}}{1-2\tfrac{{GM}}{{{rc}}^{2}}},\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left({r}^{4}j\displaystyle \frac{{df}}{{dr}}\right)+4{r}^{3}\displaystyle \frac{{dj}}{{dr}}f=0,\end{eqnarray} \tag{ 11 }$$

where $f(r)\equiv 1-\tfrac{\omega (r)}{{\rm{\Omega }}}$, $j\equiv {e}^{-\nu /2}{(1-2{GM}/{{rc}}^{2})}^{1/2}$, ω(r) is the rotational frequency of the local inertial frame at radius r, and Ω is the spin frequency of the star. The boundary conditions for the second-order partial differential Equation (11) are

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{df}}{{dr}}\right]}_{r=0}=0\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&f(r={R}_{\mathrm{NS}})=1-2\displaystyle \frac{G}{{c}^{2}}\displaystyle \frac{I}{{R}_{\mathrm{NS}}^{3}}.\end{eqnarray} \tag{ 13 }$$

To solve these coupled equations, we integrated Equations (10) and (11) outwards from the center of the star, using Equation (12) as one boundary condition, and iterated it to find the value of f₀ for which Equation (13) is valid.

In this way, we determined the mass, radius, and moment of inertia for a given EOS.

In order to be as general as possible and go beyond the current sample of proposed EOS, we tested our parametrization on a sample of mock EOS that incorporated extreme but physically allowed behavior, with the hypothesis that if our parametrization could accurately capture these extreme cases, then it would be able to reproduce more reasonable EOS as well.

For our mock EOS, we assumed the true EOS of neutron stars to be well known up to the nuclear saturation density. We then divided the high-density regime, from 1–8 ρ_sat, into 15 segments that were evenly spaced in the logarithm of the density, so that each segment would be small relative to the overall range of densities. The logarithmic sampling also sampled the lower-density region better and thus allowed more variability in the region that most affects the global properties of the star (Lattimer & Prakash 2001; Özel & Psaltis 2009; Read et al. 2009a). In each of the 15 segments, we allowed our mock EOS to be one of two extremes, which we describe below (see Section 5 for the addition of phase transitions to these mock EOS).

We set both the upper and lower extremes of each segment to be polytropes. To determine the polytropic index of each, we looked at a sample of 49 proposed EOS. This sample of EOS was compiled in order to include a wide variety of physics and calculation methods, as in Read et al. (2009a). Using the tabular data for each of these EOS, we calculated the polytropic indices for all adjacent sets of pressure and density using Equation (3) for two density ranges. The cumulative distribution of the resulting polytropic indices is shown in Figure 2. We found that the majority of polytropic indices fell between Γ = 1 and 5. Therefore, we set our nominal upper and lower extremes to have these indices and allowed any sequence of polytropic indices across our density range.

Zoom In Zoom Out Reset image size

There is, however, an absolute upper bound on the allowed polytropic index that is set by the condition of causality, which requires that the pressure gradient inside a neutron star obey the relation

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{d\epsilon }\equiv {\left(\displaystyle \frac{{c}_{s}}{c}\right)}^{2}\leqslant 1,\end{eqnarray} \tag{ 14 }$$

where c_s is the local sound speed.

We can use this relation to derive the corresponding bound on Γ by noting that we can write Equation (14) as

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{d\epsilon }=\displaystyle \frac{{dP}}{d\rho }{\left(\displaystyle \frac{d\epsilon }{d\rho }\right)}^{-1}={\left(\displaystyle \frac{{c}_{s}}{c}\right)}^{2}.\end{eqnarray} \tag{ 15 }$$

For a polytrope, dP/dρ = ΓP/ρ, and the mass density-energy density relation of Equation (6) can be expanded to

$$\begin{eqnarray}&&\displaystyle \frac{d\epsilon }{d\rho }=\displaystyle \frac{1}{\rho }(\epsilon +P).\end{eqnarray} \tag{ 16 }$$

The polytropic index can, therefore, be written as

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{c}_{s}^{2}}{{c}^{2}}\displaystyle \frac{\epsilon +P}{P}\leqslant \displaystyle \frac{\epsilon +P}{P}\equiv {{\rm{\Gamma }}}_{\mathrm{luminal}}.\end{eqnarray} \tag{ 17 }$$

In order to ensure that our upper extreme did not violate causality, we set our upper polytropic index to be the minimum of Γ = 5 and Γ_luminal.

With this set of steps, we obtained the pressure, mass density, and energy density for each segment of our mock EOS. The mock EOS starts at ${\rho }_{\mathrm{sat}}$, so once the pressure at this density is determined, the above relations will uniquely determine the rest of the behavior of each mock EOS. We introduced further freedom in our mock EOS by allowing two significantly different pressures at the starting point, ${\rho }_{\mathrm{sat}}$. This is motivated by the fact that such a bifurcation is also seen among the set of proposed EOS. We used the EOS SLy (Douchin & Haensel 2001) and PS (Pandharipande & Smith 1975) to determine the lower and higher starting pressures, respectively.

With 15 density segments, two options for the polytropic behavior along each segment, and two options for the starting pressure, our algorithm produced 2 × 2¹⁵ = 65,536 different mock EOS against which we could test our parametrization. However, we excluded any mock EOS that were too soft to produce a 1.9 ${M}_{\odot }$ neutron star (Demorest et al. 2010; Antoniadis et al. 2013), reducing our final sample size to 53,343.

4.1. The Mock EOS

The 53,343 extreme mock EOS that reached 1.9 ${M}_{\odot }$ are shown in Figure 3. All possible trajectories through this grid are indeed included, with the constraint that pressure must always increase with density (i.e., only monotonic behavior in this grid is allowed). The broadening of the mock EOS at high pressures is a result of our requirement that the upper polytropic limit must be the minimum of Γ = 5 and Γ_luminal.

Zoom In Zoom Out Reset image size

The corresponding mass–radius curves, calculated according to the method described in Section 3, are shown in Figure 4. As Figure 4 demonstrates, starting the mock EOS at the two different pressures corresponding to the two families of proposed EOS allowed us to fully span the range of reasonable radii. With this choice, we achieved a dense sampling of mass–radius curves that span radii from ∼9–17 km for M ≲ 1 ${M}_{\odot }$. Furthermore, the mock mass–radius curves include curves that shallowly slope upwards, that are nearly vertical, and that bend backwards, indicating that we have sampled a wide range of possible underlying behavior.

Zoom In Zoom Out Reset image size

The trend of increasing maximum mass with radius is a result of our causality constraint. Lindblom (1984) derived the maximum gravitational redshift of a neutron star as a function of mass by assuming that the EOS is trusted up to 3 × 10¹⁴ g cm⁻³ and configuring the resulting mass of the star to maximize the redshift. As in our analysis, that study required that the relationship between pressure and density inside the star not violate causality. After converting their relationship between the mass and maximum gravitational redshift to a mass–radius relationship, we find that the corresponding relationship of ${M}_{\max }\sim 0.24R(\,{M}_{\odot }/\mathrm{km})$, shown as the dashed line in Figure 4, is very close to what is seen in our mock EOS. The small differences between this relationship and that in our mock EOS can be attributed to the different EOS assumed up to the first fiducial density: Lindblom (1984) assumed the EOS BPS, while we assumed either SLy or PS. Furthermore, they assumed their EOS up to 3 × 10¹⁴ g cm⁻³, while we assumed SLy or PS only up to ${\rho }_{\mathrm{sat}}$ ∼ 2.7 × 10¹⁴ g cm⁻³.

4.2. Determining the Goodness of the Parametric Representation

4.3. Results of Parametrization of Mock EOS

Figure 5 shows the cumulative distribution of residuals in radius and mass for when various numbers of polytropic segments were included in the parametrization for our full sample of 53,343 extreme mock EOS. We found that our goal of ΔR < 0.5 km was achieved by sampling the EOS with five fiducial densities above ${\rho }_{\mathrm{sat}}$ (i.e., including 5 polytropes). Specifically, using five fiducial densities, we found that for 95% of the mock EOS, ΔR = 0.50, 0.44, and 0.48 km at M = 1.4, 1.6 and 1.8 ${M}_{\odot }$, respectively. In addition, a parametrization with five fiducial densities reproduced M_max to within 0.12 ${M}_{\odot }$ in 95% of the cases.

Zoom In Zoom Out Reset image size

We also calculated the difference in the moment of inertia for a neutron star with the same mass as Pulsar A in the double pulsar system, i.e., M_A = 1.338 ${M}_{\odot }$. The cumulative distribution of these residuals as a function of the number of polytropes included in the parametrizations is shown in Figure 6. By sampling the EOS at five fiducial densities, we found that the residuals in the moment of inertia are less than 0.17 × 10⁴⁵ g cm² in 95% of the cases. The moment of inertia for Pulsar A is expected to be on the order of 10⁴⁵ g cm² and is expected to be measured with 10% accuracy (Kramer & Wex 2009). As a result, depending on the exact value of the forthcoming moment of inertia measurement, 5–6 fiducial densities may be required to recreate the moment of inertia to the 10% accuracy level.

Zoom In Zoom Out Reset image size

Even a parametrization with just three fiducial densities reproduced the radii of ∼80% of our extreme mock EOS to within 0.5 km. However, requiring ΔM_max < 0.1 ${M}_{\odot }$ and Δ I/I ≲ 10% requires more points. We, therefore, conclude that, given the most recent observational uncertainties and the continuing prospects for even smaller errors in the near future, a parametrization that samples the EOS at five fiducial densities is optimal. Specifically, we recommend spacing the five fiducial densities evenly in the logarithm of the density, such that (ρ₀, ρ₁, ρ₂, ρ₃, ρ₄, ρ₅) = (1.0, 1.4, 2.2, 3.3, 4.9, 7.4) ${\rho }_{\mathrm{sat}}$.

We also considered more diverse EOS by allowing there to be a first-order phase transition in the mock EOS described in Section 4. We allowed only one phase transition per mock EOS, but allowed the phase transition to start in any of our 15 density segments and to last anywhere between 1 and 15 segments. For the remaining segments, the mock EOS was polytropic with an index of Γ = 1 or the minimum of Γ = 5 and Γ_luminal, as above. With the addition of these phase transitions, there are now ${2}^{N}\times N(N+1)/2$ possibilities for N segments of the EOS, for each possible starting pressure. For 15 segments and our two starting pressures (corresponding to the EOS SLy and PS at ${\rho }_{\mathrm{sat}}$), this corresponds to 7,864,320 possibilities for the mock EOS.

We randomly sampled ∼175,000 of these mock EOS, with roughly half starting at each of the initial pressures. We then applied our parametrization to each mock EOS. As we did for the original set of mock EOS, we varied the number of polytropic segments in the parametrization between 3 and 12 and calculated the resulting residuals in radius, mass, and moment of inertia. The mass and radius residuals are shown in Figure 7. With five fiducial densities above ${\rho }_{\mathrm{sat}}$, i.e., the optimal number of polytropic segments found in Section 4.3, we found that 95% of the radius residuals were less than 0.48, 0.43, and 0.46 km at 1.4, 1.6, and 1.8 ${M}_{\odot }$, respectively. These errors are comparable to those from the mock EOS without phase transitions. We also found that 95% of the differences in maximum mass were less than 0.22 ${M}_{\odot }$, which is a larger error than in our previous, less extreme sample of mock EOS. However, 70% of the errors in maximum mass were still less than 0.1 ${M}_{\odot }$ for this sample, indicating that this parametrization reasonably recreated the maximum mass for many of our most extreme sample of mock EOS. The error distribution for the moment of inertia was almost identical to the distribution for the sample without phase transitions. We found that a parametrization with five fiducial densities reproduced the moment of inertia to within 0.17 × 10⁴⁵ g cm² for 95% of the mock EOS with phase transitions.

Zoom In Zoom Out Reset image size

The tail end of all four distributions in Figure 7 extended to higher errors than did the tail in Figure 5. This is because the arbitrarily placed phase transitions in this sample of mock EOS can make the resulting mass–radius curve have very sharp turn-overs. If a sharp turn-over occurs near one of the fiducial masses at which we measure radius residuals, we will infer an artificially high error. However, the distributions in Figure 7 show that the probability of the true EOS falling in this tail is small (≲5% with five fiducial densities). We, therefore, find that, for the vast majority of cases, a parametrization that samples the EOS at five fiducial densities is sufficient to reproduce the radius and maximum mass observables to within 0.5 km and 0.1–0.2 ${M}_{\odot }$, even if the EOS contains a first-order phase transition at an arbitrary density and over an arbitrary range.

Even though we optimized our parametrization using EOS that span a much wider range of possibilities than the currently proposed ones, we nevertheless explored how well this parametrization reproduced the physically motivated EOS found in the literature. To this end, we applied our optimized parametrization (five fiducial densities above ${\rho }_{\mathrm{sat}}$) to a sample of 42 proposed EOS,² which incorporate a variety of different physical possibilities and calculation methods. (The tabular data for these EOS are compiled in Cook et al. 1994; Lattimer & Prakash 2001; Read et al. 2009a; and Özel & Freire 2016). Our sample included purely nucleonic EOS, such as: relativistic (BPAL12 and ENG) and nonrelativistic (BBB2) Brueckner–Hartree–Fock EOS; variational-method EOS (e.g., FPS and WFF3); and a potential-method EOS (SLy). We also include models which incorporate more exotic particles, including, for example, a neutron-only EOS with pion condensates (PS), a relativistic mean-field theory EOS with hyperons and quarks (PCL2), and an effective-potential EOS with hyperons (BGN1H1).

In applying our parametrization to these proposed EOS, we no longer connected to SLy or PS for ρ < ${\rho }_{\mathrm{sat}}$. Instead, we assumed that each EOS is known up to ${\rho }_{\mathrm{sat}}$, and we, therefore, used the full EOS that we were parametrizing for the low-density regime. For ρ > ${\rho }_{\mathrm{sat}}$, we applied our parametrization as above. After applying our parametrization, we calculated the resulting residuals in radii at the three fiducial masses,³ the maximum mass, and the moment of inertia. The residuals for our optimized, 5-polytrope parametrization are shown in Figure 8.

Zoom In Zoom Out Reset image size

The errors in applying the parametrization to the proposed, physically motivated EOS were much lower than for the more extreme, mock EOS. Our optimized parametrization reproduced the radii of ∼95% of the proposed EOS to within 0.10, 0.12, and 0.09 km (at 1.4, 1.6, and 1.8 ${M}_{\odot }$, respectively) and the maximum masses of 95% of the EOS to within 0.04 ${M}_{\odot }$. As examples, we show in Figure 9 the full mass–radius relation as well as the one calculated from our parametrized EOS for several proposed EOS: AP4 (nucleonic), ALF4 (quark hybrid), SLy (nucleonic), and BGN1H1 (includes hyperons). As can be seen here, the differences in mass–radius space are extremely small.

Zoom In Zoom Out Reset image size

In this paper, we investigated an optimal parametrization of the neutron star EOS that can be used to interpret neutron star observations. We found that a parametric EOS with five polytropic segments evenly spaced in logarithm between 1 and 7.4 ${\rho }_{\mathrm{sat}}$ was sufficient to reproduce the radii of the proposed EOS to within 0.12 km and the maximum mass to within 0.04 ${M}_{\odot }$ in 95% of cases. This parametrization was also able to reproduce the radii of our more extreme, mock EOS to within 0.5 km, suggesting that even if a more extreme EOS is proposed or realized in nature, our parametrization will be robust enough to reproduce it well.

The radii of approximately 15 neutron stars have already been measured, for most of which the masses are also known (Guillot et al. 2013; Guillot & Rutledge 2014; Heinke et al. 2014; Nättilä et al. 2015; Bogdanov et al. 2016; Özel et al. 2016). Even though the uncertainties in the individual radius measurements are in the order of ∼2 km, combining all the measurements leads to a rather narrow range of predicted neutron star radii. Moreover, these spectroscopic measurements will be independently checked in the near future by, e.g., waveform modeling with NICER (Gendreau et al. 2012). The forthcoming measurement of the moment of inertia, which is a higher-order moment of the mass distribution of the star (Lyne et al. 2004; Kramer & Wex 2009), and the prospect of eventual measurements of other higher-order moments using gravitational wave observations from coalescing neutron stars (e.g., Read et al. 2009b; Del Pozzo et al. 2013; Agathos et al. 2015) will provide additional constraints. Finally, the recent mass measurements of approximately two-solar-mass neutron stars (Demorest et al. 2010; Antoniadis et al. 2013) already place strong priors on empirically inferred EOS.

Previous empirical inferences of the EOS based on existing observations typically relied on functional forms with 3–4 parameters (e.g., Özel et al. 2010, 2016; Steiner et al. 2010, 2016). It is difficult to compare the errors produced in these parametrizations directly to our own, as they are not always quantified in mass–radius space. Read et al. (2009a) report errors in their parametrization in pressure versus density. They find that a four-parameter EOS produced average errors in $\mathrm{log}P$ of 0.013, which corresponds to a fractional error in pressure of ∼3%. Similarly, Lindblom & Indik (2012) report errors in their spectral parametrization in energy density versus enthalpy space. Their errors ranged from 1%–15% when they used 4–5 spectral parameters. Finally, Steiner et al. (2013) performed a preliminary comparison between five different types of models (e.g., a parametrization of two polytropes and one with four line segments) and found the largest difference in the predicted radii between any two of their models to be 0.8 km. However, these models were not necessarily optimized in terms of the number of parameters.

Ultimately, it is the inverse process of what we have shown that will be of interest: the inference of a parametric EOS from astrophysical measurements. We will address this aspect of the problem in a forthcoming paper. However, our work here suggests that the uncertainties in the inference of the EOS will be dominated by the quality of the data and the uncertainties in the inversion process itself, rather than those introduced by the parametrization of the EOS.

We thank Gordon Baym for useful comments on the manuscript. We gratefully acknowledge support from NASA Grant NNX16AC56G.

For completeness, we also explored a parametrization that uses linear segments between a number of density points to represent the EOS. As in the case of our polytropic parametrization, we started this parameterization at ${\rho }_{\mathrm{sat}}$ and spaced the segments evenly in the logarithm of the density. The EOS along each linear segment is given by

$$\begin{eqnarray}&&P={m}_{i}\rho +{b}_{i}\quad \quad ({\rho }_{i-1}\leqslant \rho \leqslant {\rho }_{i}),\end{eqnarray} \tag{ 18 }$$

where continuity at the endpoints implies

$$\begin{eqnarray}&&{m}_{i}=\displaystyle \frac{{P}_{i}-{P}_{i-1}}{{\rho }_{i}-{\rho }_{i-1}}\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&{b}_{i}={P}_{i-1}-\left(\displaystyle \frac{{P}_{i}-{P}_{i-1}}{{\rho }_{i}-{\rho }_{i-1}}\right){\rho }_{i-1}.\end{eqnarray} \tag{ 20 }$$

Using this linear relationship for pressure to integrate the differential Equation (6), we find that the energy density in this case is given by

$$\begin{eqnarray}&&\epsilon (\rho )=(1+a)\rho {c}^{2}+m\rho \mathrm{log}\rho -b,\end{eqnarray} \tag{ 21 }$$

where a is an integration constant. By requiring continuity in the energy density at the endpoints of each segment, we can solve for the integration constant such that

$$\begin{eqnarray}\epsilon (\rho ) & = & \left(\displaystyle \frac{{\epsilon }_{i-1}+{b}_{i}}{{\rho }_{i-1}}-{m}_{i}\mathrm{log}{\rho }_{i-1}\right)\rho +{m}_{i}\rho \mathrm{log}\rho -{b}_{i},\\ & & ({\rho }_{i-1}\leqslant \rho \leqslant {\rho }_{i}).\end{eqnarray} \tag{ 22 }$$

We used Equation (22) to relate the linear EOS to the energy density and then used that energy density to integrate the TOV equations and solved for the total mass, radius, and moment of inertia of the neutron star.

We applied a five-segment parametrization to the ∼53,000 mock EOS shown in Figure 4 and calculated the errors in radius at 1.4, 1.6, and 1.8 ${M}_{\odot }$ as well as the errors in the maximum mass. In this way, we can directly compare these results with those from the five-polytrope optimal parametrization in Figure 10. We find that the polytropic parametrization performs modestly better. However, the linear parametrization is still able to recreate the radii of the full EOS to within ≲0.5 km for ∼80% of the extreme, mock EOS.

Zoom In Zoom Out Reset image size

Even though our mock EOS are composed of polytropic segments, there are multiple mock EOS segments per parametrization segment. Moreover, the mock EOS segments are offset from the parametrization segments. We, therefore, do not expect that the mock EOS should significantly bias the performance of a polytropic parametrization over a linear parametrization.

In order to see how well the linear parametrization performs for more physically motivated EOS, we also applied it to the sample of 42 proposed EOS from Section 6. We compare these results to those of the polytropic parametrization in Figure 11. We find again that the polytropic parametrization performs better than the linear parametrization, although the differences between the two parametrizations are most significant at radius and maximum mass errors that are well below observational uncertainties.

Zoom In Zoom Out Reset image size

It is likely that the same level of errors could be achieved by the linear parametrization if more than five segments were included. However, for five segments, the polytropic parametrization performs modestly better in most cases. A polytropic parametrization is also the more natural choice for the neutron star EOS. We, therefore, recommend a five-segment polytropic parametrization over a linear one.