THE EQUATION OF STATE FROM OBSERVED MASSES AND RADII OF NEUTRON STARS

Observational information for three Type-I X-ray bursters with PRE bursts (Özel et al. 2009; Güver et al. 2010a, 2010b) is presented in Table 1, along with associated uncertainties. We have reexamined the uncertainties for each object as described in the Appendix; in some cases our values for the uncertainties differ from those used previously.

We first fix f_c = 1.4 and X = 0 as was done in previous work for EXO 1745–248 (Özel et al. 2009). We observe in this case from Table 1 that no real-valued solution of Equations (7)–(9) is possible for any of the sources because α>1/8 for the central values of the observables. Moreover, it is not consistent with X-ray burst models to fix the color-correction factor to f_c = 1.4 and the hydrogen abundance to X = 0, and indeed these restrictions were relaxed in subsequent analyses of 4U 1608–522 and 4U 1820–30 (Güver et al. 2010a, 2010b). Atmosphere models (Madej et al. 2004) suggest a range 1.33 < f_c < 1.81; these are consistent with earlier calculations (London et al. 1986; Ebisuzaki & Nakamura 1988). The largest values of f_c are approached as the flux approaches the Eddington limit. In the tail of the burst, the flux is lower and f_c does not vary strongly. We choose to select the value of f_c in our MC analysis from a boxcar distribution centered at f_c = 1.4 with an uncertainty of 0.07. This selection is comparable to the boxcar distribution with f_c = 1.35 ± 0.05 used previously (Güver et al. 2010a, 2010b). (In the following, we denote a boxcar uncertainty in X as ΔX and a Gaussian uncertainty with σ_X.) In addition, as described in Appendices A.1 and A.2, we find insufficient information from the observations to exclude any possible values of X for EXO 1745–248 and 4U 1608–522, so we choose a uniform distribution with 0 < X < 0.7, which is also assumed by Güver et al. (2010a, 2010b). The source 4U 1820–30 is an ultracompact and accretes He-rich fuel (see Appendix A.3).

If values of f_c and X or the observables F_TD,∞, A, or D are selected at random from their respective probability distributions, some combinations will satisfy the condition α < 1/8. In this way, using only these acceptable combinations, the most probable values $\hat{\alpha }$ and $\hat{\gamma }$ can be established using MC sampling. (Here and below, we define $\hat{X}$ to be the most probable value of X obtained from an MC simulation, and the uncertainties on this quantity are determined by selecting the region surrounding the most probable value that includes 68% of the total MC weight.) Values for M and R, and their uncertainties, can then be determined from $\hat{\alpha }$ and $\hat{\gamma }$ using Equations (7)–(9). But the fraction of accepted realizations is then very small as shown in Table 2. The entries in the first group of this table give the most probable values $\hat{\alpha }$, $\hat{\gamma }$, and $\hat{R}_{\infty }$ using r_ph = R and values for F_TD,∞, A, and D from their probability distributions summarized in Table 1 and, for f_c and X, from previous discussion. For each quantity, the uncertainty is determined by creating a histogram for the indicated MC realizations, sorting the bins by decreasing weight, and selecting the range which encloses 68% of the sum of all bins. We identify $R_{\infty } = [F_{\infty }/(4\pi \sigma T_{\mathrm{eff,}\infty }^{4})]^{1/2} = R/\sqrt{1-2\beta } = \alpha \gamma$ and impose no a priori restriction for the value of α. The most probable value of α, $\hat{\alpha }$, is observed to be approximately a factor $1+\bar{X}$ larger than the value of α in Table 1, where $\bar{X}$ is the average value of X in its assumed range, i.e., 1.35 for EXO 1745–248 and 4U 1608–522, and 1.0 for 4U 1820–30. For all three sources, $\hat{\alpha }> 1/8$ by several standard deviations.

Table 2. Comparison of Monte Carlo Realizations for PRE Bursts Employing Data from Table 1

Quantity	EXO 1745–248	4U 1608–522	4U 1820–30
FTD,∞ = FEdd[1 − 2β(R)]1/2; α unrestricted
$\hat{\alpha }$	0.165+0.045−0.028	0.229+0.054−0.038	0.167+0.018−0.017
$\hat{\gamma }(\rm {km})$	70.4+19.2−14.7	80.2+18.8−15.5	80.9+21.2−3.56
${\hat{R}_{\infty }}(\rm {km})$	13.2+1.7−1.6	20.5+5.7−5.5	15.3+1.7−1.5
FTD,∞ = FEdd[1 − 2β(R)]1/2; α < 1/8 restriction
$\hat{\alpha }$	0.122+0.003−0.007	0.123+0.003−0.005	0.123+0.003−0.004
$\hat{\gamma }(\rm {km})$	107.9+14.2−14.2	127.1+8.0−7.9	109.3+5.3−8.1
${\hat{R}_{\infty }}(\rm {km})$	12.7+1.6−1.6	15.2+0.85−0.81	13.2+0.78−1.1
Points accepted	4.4%	0.24%	0.44%
FTD,∞ = FEdd[1 − 2β(R)]1/2; α < 1/8; original inputs from Ozel et al.
$\hat{\alpha }$	0.122+0.003−0.003	0.122+0.003−0.003	0.124+0.001−0.001
$\hat{\gamma }(\rm {km})$	109.0+4.2−4.2	113.1+5.0−5.0	104.0+1.8−1.8
${\hat{R}_{\infty }}(\rm {km})$	13.3+0.4−0.4	13.8+0.5−0.5	12.9+0.2−0.2
Points accepted	13%	0.18%	1.5 × 10−6%
FTD,∞ = FEdd; α < 3−3/2 = 0.192 restriction
$\hat{\alpha }$	0.165+0.026−0.018	0.180+0.012−0.023	0.167+0.016−0.016
$\hat{\gamma }(\rm {km})$	79.0+17.2−12.4	97.0+15.4−14.4	91.0+13.8−10.6
${\hat{R}_{\infty }}(\rm {km})$	13.1+1.8−1.6	15.3+2.8−1.3	15.0+1.8−1.4
Points accepted	65%	15%	91%
FEdd[1 − 2β(R)]1/2 < FTD,∞ < FEdd; q3 + p2 < 0 restriction
$\hat{\alpha }$	0.143+0.022−0.018	0.155+0.018−0.019	0.158+0.014−0.015
$\hat{\gamma }(\rm {km})$	87.1+17.5−13.8	103.0+16.6−13.6	104.8+5.6−17.0
${\hat{R}_{\infty }}(\rm {km})$	13.0+1.8−1.6	15.3+2.3−1.4	15.0+1.6−1.6
Points accepted	32%	5.8%	37%

Notes. In the third group, we use the values and distributions of F_TD, A, D, f_c, and X from Özel et al. (2009), Güver et al. (2010a), and Güver et al. (2010b). For all other entries, we use 0 < X < 0.7 for the hydrogen mass fraction, except for 4U 1820–30, for which X = 0 and we use 1.33 < f_c < 1.47 for the color-correction factor.

Download table as:  ASCIITypeset image

In the second group of Table 2, we repeat the analysis, but this time only accept realizations for which α ⩽ 1/8. Note that at most 4% of the realizations are accepted. This implies, at a minimum, that the assumptions in this model are incomplete. The most probable value for α decreases and now lies approximately 1σ below the value 1/8, and the size of confidence interval for $\hat{\alpha }$ is now significantly reduced compared to the unrestricted case, by factors of 4–18. The net effect is to pre-select the value 1/8 for $\hat{\alpha }$ irrespective of the observed values of F_TD,∞, A, and D. Another consequence of the selective rejection of most of the MC realizations is to increase the most probable value of γ and decreases the most probable value of R_∞, and to greatly decrease their uncertainties; the latter result was already noted by Özel et al. (2010). In addition, the resulting values of $\bar{\gamma }$ and $\bar{R}_\infty$ for the different sources are also "herded" into similar values (Özel et al. 2009; Güver et al. 2010a, 2010b). The third group summarizes the MC results obtained using the exact same observational inputs found in those references. Note that for 4U 1820–30, only about 1 realization out of 100 million is now accepted! The larger acceptance rate for EXO 1745–248 is due to the assumption that X = 0 which gives the smallest possible α. Analyses of X-ray bursts with this model result in certain values for M and R that are nearly independent of α and γ, and are therefore almost completely independent of the observables F_TD,∞, A, and D, which seems unrealistic. The fact that only a small fraction of realizations are accepted may indicate that the underlying model contains systematic errors or is unphysical, a point we shall explore in the next section.

Figure 1 displays the MC probability distributions for the masses and radii of EXO 1745–248, 4U 1608–522, and 4U 1820–30 as determined from α and γ following Equations (7)–(9), and with the restriction α ⩽ 1/8. There are two peaks in the distributions because of the quadratic in Equation (7). We also impose that β < 1/2.94 ≈ 0.34 (Glendenning 1992), which is based on requiring a subluminal sound speed throughout the star (see also Lattimer & Prakash 2007). This causes the rejection of a small fraction of realizations from the higher-redshift region. For all three sources, our confidence intervals are larger than computed in previous works (Özel et al. 2009; Güver et al. 2010a, 2010b). The origin of this distinction is because we use larger variations f_c and, for EXO 1745–248, in X, and we employ a different uncertainty in D for 4U 1820–30. A detailed discussion of the treatment of the uncertainties in the observables for each object is given in the Appendix.

Zoom In Zoom Out Reset image size

When using a Monte Carlo scheme to determine uncertainties, we find that the two error ranges in Figure 1 are equally populated, while Özel et al. (2009), Güver et al. (2010a, 2010b), and Özel et al. (2010), transforming probabilities using a Jacobian, indicate a much higher probability for the higher redshift solution. This difference is caused by a numerical error in their Jacobian (F. Özel 2010, private communication); when the correct Jacobian is used, both integration methods agree.

The lack of real-valued solutions for M and R for the most probable values of the observables, and indeed, the rejection of the vast majority of MC realizations, motivates us to consider another possibility; namely, that at "touchdown" the photosphere is still extended. To explore this situation, we first consider the extreme possibility that r_ph ≫ R at the point identified as touchdown from the maximum in T_bb,∞. In this case, F_TD,∞ = F_Edd and is independent of the stellar radius R. With this assumption, α (Equation (3)) and γ (Equation (4)) are now related to β and R via

$$\begin{eqnarray} \alpha &=& \beta \sqrt{1-2\beta }, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \gamma &=& \frac{R}{\beta (1-2\beta)}. \end{eqnarray} \tag{ 11 }$$

Defining a quantity θ = cos⁻¹(1 − 54α²), the expressions for the compactness, radius, and mass are then

$$\begin{eqnarray} \beta _1 &=& {1\over 6}[1+\sqrt{3}\sin (\theta /3)- \cos (\theta /3)], \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \beta _2 &=& {1\over 6}\left[1+2\cos (\theta /3)\right], \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} R &=& \alpha \gamma \sqrt{1-2\beta _{1,2}}, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} M &=& \frac{c^2}{G}\alpha ^2\gamma. \end{eqnarray} \tag{ 15 }$$

Obviously, M is real for all α, γ>0. For α < 3^−3/2 ≃ 0.192, θ is real and there are three real roots for β and R. Only two of these, 0 ⩽ β₁ ⩽ 1/3 and 1/3 ⩽ β₂ ⩽ 1/2, are physically meaningful: the other root is negative. When α>3^−3/2 there is one real root for β and two imaginary ones, but the real root is negative.

From the top group in Table 2, we can see that both 4U 1820–30 and EXO 1745–248 have most probable values of α that satisfy the constraint for positive real values of β and R, and 4U 1608–522 lies less than 1σ above this limit. Therefore, we now find that a much larger fraction of the MC realizations are accepted when selecting values for F_TD,∞, A, D, f_c, and X from their probability distributions. This is shown in the fourth group of Table 2, for which we use the model given in Equations (10) and (11) and impose the restriction α < 0.192. In contrast to the case in which r_ph = R, the uncertainties in $\hat{\alpha }$ and $\hat{\gamma }$ are not as strongly diminished. Moreover, the values for $\hat{\alpha }$ are no longer nearly the same for the three sources. While in principle each accepted realization results in two M, R values, one corresponding to β₁ and the other to β₂, nearly all the β₂ realizations are rejected on the basis of causality when β₂>1/2.94. Figure 2 displays the probability distributions for the fourth group of runs listed in Table 2, as well as their 68% and 95% contours. The error contours are larger than those shown in Figure 1 and considerably larger than determined previously (Özel 2006; Özel et al. 2009; Güver et al. 2010a, 2010b).

Zoom In Zoom Out Reset image size

It is perhaps unphysical to make the extreme assumption r_ph ≫ R when T_eff reaches a maximum. But because of the relative consistency of the solutions and the much larger fraction of MC points accepted in this case, it seems reasonable to consider the possibility that at touchdown, that is, when T_bb,∞ reaches a maximum and the normalization reaches a minimum, the photospheric radius may in fact be larger than the quiescent stellar radius. We can write $\alpha = \beta \sqrt{1-2\beta }\sqrt{1-\beta _{\mathrm{ph}}}$, where β_ph = GM/(r_phc²) < β. Defining the quantity h = 2R/r_ph, the quantities α and γ now become

$$\begin{eqnarray} \alpha &=& \frac{F_{\mathrm{Edd}}}{\sqrt{A}} \frac{\kappa D}{c^{3}f_{c}^{2}} = \beta \sqrt{1-2\beta }\sqrt{1-h\beta }, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \gamma &=& \frac{Ac^{3}f_{c}^{4}}{F_{\mathrm{Edd}}\,\kappa } = \frac{R}{\beta (1-2\beta)\sqrt{1-h\beta }}. \end{eqnarray} \tag{ 17 }$$

Note that R_∞ = αγ remains unchanged. Equation (16) is a quartic equation for β. Two positive real solutions for β exist when q² + p³ < 0, where

$$\begin{eqnarray} p &=&{1\over 6h}\left(\alpha ^2-{1\over 24h}\right)\,, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} q &=&-{1\over 12h}\left({1\over 144h^2}+3a\alpha ^2\right)\,, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} a &=& {1\over 6h}\left(1-{3(2+h)^2\over 16h}\right)\,; \end{eqnarray} \tag{ 20 }$$

otherwise the only real solution is negative. Defining the quantities

$$\begin{eqnarray} \theta &=& \cos ^{-1}\left[\left({q\over -p}\right)^{3/2}\right]\,, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} v &=& 2\sqrt{-p}\cos ({\theta /3})\,, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} w &=& \sqrt{2v-2a}\,, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} z_\pm &=& -\left(2v+4a\pm {2b\over w}\right)\,, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} b &=& {(2+h)^2\over 8h^2}\left[1-{2+h\over 8h}\right]\,, \end{eqnarray} \tag{ 25 }$$

the two solutions for β with values in the range 0 ⩽ β ⩽ 0.5 are

$$\begin{equation} \beta _{\pm } ={2+h\over 8h}\pm {w-\sqrt{z_\pm }\over 2}\,, \end{equation} \tag{ 26 }$$

where the sign has the same sense in the two occurrences.

Figure 3 displays the resulting probability distributions assuming that h is uniformly distributed in the range 0 < h < 2, i.e., R < r_ph < ∞. The error contours are much larger than those determined in the h = 2 (r_ph = R) case and slightly larger than those in the h = 0 (r_ph ≫ R) cases. Most of the solutions for β₊ are rejected because they lead to acausal combinations of M and R. The fifth group in Table 2 presents the associated values of α, γ, and R_∞. These results are nearly identical to those of Figure 2 because small values of the photospheric radius are strongly disfavored by the requirement that the masses and radii are real valued. We will therefore use the r_ph ≫ R probability distributions in M and R for the Bayesian estimation of the EOS in Section 4, and note that these results will be essentially identical to what one would obtain without assuming a particular value for the radius of the photosphere, as long as it is not a priori assumed to be equal to R.

Zoom In Zoom Out Reset image size

Using the number of accepted MC realizations as a guide, we can also obtain an estimate for the lower limit of the location of the photosphere at touchdown. As is common in Bayesian analysis, we use the ratio 0.1 as a guide; models for which the fraction of accepted realizations is less than 0.1 are rejected. This implies that r_ph>5.0R for 4U 1608–522, r_ph>1.1R for EXO 1745–248, and r_ph>1.4R for 4U 1820–30. This analysis appears to strongly disfavor an interpretation of the touchdown radius r_ph = R and we find it likely that the photospheric radius is extended at touchdown. At the same time, there is little difference between the results assuming h = 0 or a uniform range 0 < h < 2. It would then seem justified that we, in our remaining analysis, reject the interpretation r_ph = R, and make use of the r_ph ≫ R, h = 0 interpretation. Nonetheless, we will also include results for the r_ph = R interpretation in the work below for comparison to previous studies. We will show that some aspects of the EOS constraints are dependent on whether one assumes that r_ph = R or r_ph ≫ R.

Bayes theorem can be formulated as (see, e.g., Grinstead & Snell 1997)

$$\begin{equation} P({\cal M}|D) = \frac{P(D|{\cal M}) P({\cal M})}{P(D)} \,, \end{equation} \tag{ 27 }$$

where $P({\cal M})$ is the prior probability of the model ${\cal M}$ without any information from the data D, P(D) is the prior probability of the data D, $P(D|{\cal M})$ is the conditional probability of the data D given the model ${\cal M}$, and $P({\cal M}|D)$ is the conditional probability of the model ${\cal M}$ given the data D. This latter quantity, $P({\cal M}|D)$ is what we want to obtain, namely, the probability that a given model is correct given the data.

For many non-overlapping models ${\cal M}_i$ which exhaust the total model space ${\cal M}$, this relation can be rewritten

$$\begin{equation} P({\cal M}_i|D) = \frac{P(D|{\cal M}_i) P({\cal M}_i)}{\sum _j \! P(D|{\cal M}_j) P({\cal M}_j)} \,. \end{equation} \tag{ 28 }$$

For our problem, the model space consists of all of the parameters for the EOS, $p_{i=1,\ldots,N_p}$, plus values for all of the masses of the neutron stars for which we have data, $M_{i=1,\ldots,N_M}$. From the parameters, p_i, we can construct the EOS and solve the TOV equations to get a radius R_i for each of the neutron star masses M_i. To be more concise in the following, we refer to our model ${\cal M}(p_1,p_2,\ldots,p_{N_p},M_1,M_2,\ldots,M_{N_M})$ as ${\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})$. Applied to this specific problem, Equation (28) is

$$\begin{eqnarray} &&\lefteqn{P[{\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})|D] = P[D|{\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})] }\nonumber \\ &&\quad\times \! P[{\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})] \left\lbrace \int P[D|{\cal M}] P[{\cal M}] d^{N} {\cal M} \right\rbrace ^{-1},\nonumber\\ \end{eqnarray} \tag{ 29 }$$

where N = N_p + N_M is the dimensionality of our model space. The total number of EOS parameters is N_p = 8 and the total number of neutron stars in our data set is N_M = 6.

We construct our data D as a set of N_M probability distributions, ${\cal D}_i(M,R)$ in the (M, R) plane, which are all normalized to unity, i.e.,

$$\begin{equation} \int _{M_{\mathrm{low}}}^{M_{\mathrm{high}}} {\it dM} \int _{R_{\mathrm{low}}}^{R_{\mathrm{high}}} {\it dR} {\cal D}_i(M,R) = 1 \forall i \,. \end{equation} \tag{ 30 }$$

This normalization ensures that the data set for each neutron star is treated on an equal footing. We choose M_low = 0.8 M_☉ because current core-collapse supernova simulations fail to generate lower neutron star masses. The remaining limits, M_high = 2.5 M_☉, R_low = 5 km, and R_high = 18 km are extreme enough to ensure that they have no impact on our final results. None of the probability distributions ${\cal D}$ inferred from the data described above have a significant probability for neutron stars outside of these ranges. Note also that the results for the neutron stars in M13 and ω Cen are cut off in Webb & Barret (2007) for radii below 8 km so our distributions also exhibit the same cutoff.

In order to apply Equation (29) to our problem, we assume that $P(D|{\cal M})$, the conditional probability of the data given the model, is proportional to the product over the probability distributions ${\cal D}_i$ evaluated at the masses which are chosen in the model and evaluated at the radii which are determined from the model, i.e.,

$$\begin{equation} P[D|{\cal M}(p_{1{\ldots }N_p}, \! M_{1{\ldots }N_M})] \propto \prod _{i=1,\ldots,N_M} \left. \!{\cal D}_i(M,R) \right|_{M=M_i,R=R(M_i)}. \end{equation} \tag{ 31 }$$

This implicitly assumes that all of the data distributions ${\cal D}_i$ are independent of each other and also independent of the model assumptions and prior distributions. Another required input for Equation (29) is the prior distribution. We assume that the prior distribution is uniform in all of the n_p + n_M model parameters, except for a few physical constraints on the parameter space described below. Taking a uniform distribution just means that the $P({\cal M})$ terms cancel from Equation (29) and the integration limits become the corresponding prior parameter limits.

In the maximum likelihood formalism, the function $P[D|{\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})]$ is equivalent to the likelihood function, and maximizing the likelihood function gives the best fit to the data. In the case where the probability distributions represent Gaussian uncertainties, then the best fit is equivalent to a least-squares fit (Bevington & Robinson 2002). In Bayesian inference, model parameters are determined using marginal estimation, where the posterior probability for a model parameter p_j is given by

$$\begin{eqnarray} P[p_j|D](p_j) &=& \frac{1}{V} \! \int \! P[D|{\cal M}] {\it d p}_1 {\it d p}_2 \ldots {\it d p}_{j-1} {\it d p}_{j+1} \ldots {\it d p}_{N_p}\nonumber\\ &&\times {\it d M}_1 {\it d M}_2 \ldots {\it d M}_{N_M}, \end{eqnarray} \tag{ 32 }$$

where V is the denominator in Equation (29), without the model priors that determine the integration limits. The one-dimensional function P[p_j|D](p_j) represents the probability that the jth parameter will take a particular value given the observational data. Our problem thus boils down to computing integrals of the form in Equation (32). We use the Metropolis–Hastings algorithm to construct a Markov chain to simulate the distribution $P[D|{\cal M}(p_{1{\ldots }N_p},M_{1{\ldots }N_M})]$. For each point, we generate the EOS parameters p_i and the neutron star masses M_i from a uniform distribution within limits that are described below. The TOV equations are solved and this generates an R(M) curve and the radii for each neutron star, R_i. From these six masses and radii, the weight function $P\left[D|{\cal M}\right]$ is computed from Equation (31) and the point is rejected or accepted using the Metropolis algorithm. In order to compute posterior probabilities P[p_j|D](p_j), we construct several histograms of the integrand $P[D|{\cal M}]$ from all of the points that were accepted. To construct the 1σ regions, we sort the histogram bins by decreasing probability, and select the first N bins which exhaust 68% of the total weight. A similar procedure is used for the 2σ regions. In order to constrain the full EOS of neutron star matter, we histogram the value of the pressure predicted by each EOS in the Markov chain on a fixed grid of energy density. To create the predicted curve, R(M), we construct a histogram of the predicted radii for each EOS in the Markov chain on a fixed mass grid.

This analysis is easily extensible to a different number of EOS parameters or a different number of neutron star data sets. The only issue is that of computer time: the TOV equations must be solved for each point in the model space, and enough points must be selected to cover the model space fully. Another advantage is the explicit presence of the prior distributions, $P({\cal M})$. Although we have set these distributions to unity for this work, future work will utilize these terms to examine the impact of constraints on the EOS from terrestrial experiments.

We divide the EOS into four energy density regimes. The region below the transition energy density ε_trans ≈ ε₀/2 is the crust, for which we use the EOS of Baym et al. (1971) and Negele & Vautherin (1983). Here, ε₀ is the nuclear saturation energy density; it is convenient to remember that the nuclear saturation baryon density 0.16 fm^-3 corresponds to an energy density of ≈160 MeV fm^-3 and a mass density of ≈2.7 × 10¹⁴ g cm^-3. For ε_trans < ε < ε₁, we use a schematic expression representing charge-neutral uniform baryonic matter in beta equilibrium that is compatible with laboratory data. Finally, two polytropic pressure–density relations are used in the regions ε₁ < ε < ε₂ and ε>ε₂. The densities ε₁ and ε₂ are themselves parameters of the model. The schematic EOS for ε_trans < ε < ε₁ is taken to be

$$\begin{eqnarray} \varepsilon &=& n_B\left\lbrace m_B+B+\frac{K}{18}(u-1)^2+ \frac{K^{\prime }}{162}(u-1)^3\right.\nonumber\\ &&\left.+ \ (1- 2 x)^2 [S_k u^{2/3}+ S_p u^{\gamma }]+{3\over 4}\hbar c x(3\pi ^2n_bx)^{1/3}\right\rbrace,\qquad \end{eqnarray} \tag{ 33 }$$

where n_B is the baryon number density, m_B is the baryon mass, u = n_B/n₀, and x is the proton (electron) fraction. The saturation number density, n₀, is fixed at 0.16 fm⁻³, the binding energy of saturated nuclear matter, B, is fixed at −16 MeV, and the kinetic part of the symmetry energy, S_k, is fixed at 17 MeV. The compressibility K, the skewness K', the bulk symmetry energy parameter S_v ≡ S_p + S_k (where S_p is the potential part of the symmetry energy), and the density dependence of the symmetry energy γ, are parameters which we constrain to lie within the ranges specified in Table 3. These limits operate as constraints in our otherwise trivial prior distributions, $P({\cal M})$. To avoid bias in our results and to ensure that our model space is not overconstrained, we have intentionally made these ranges larger than normally expected from modern models of the EOS of uniform matter which are fit to laboratory nuclei. While in principle the crust EOS for each set of EOS parameters could be different, as described in Steiner (2008), in practice the masses and radii are not strongly affected by changes in the crust at this level. The transition between the crust EOS and the low-density EOS is typically around half of the nuclear saturation density and is determined for each parameter set by ensuring that the energy is minimized as a function of the number density. We have opted, at this stage, not to include correlations between parameters that have been shown to exist from nuclear systematics or neutron matter calculations. For example, the values of S_v and γ (or, equivalently, S_v and S_s, the surface symmetry parameter) are highly correlated (Steiner et al. 2005) in liquid drop mass formula fits to nuclear masses. Such correlations will be considered in a future publication.

The last term in Equation (33) is due to electrons. The proton fraction x is determined as a function of density by the condition of beta equilibrium

$$\begin{equation} {\partial \varepsilon \over \partial x}=\hbar c(3\pi ^2n_Bx)^{1/3}-4 [S_ku^{2/3}+S_pu^\gamma ]\,(1-2x)=0\,, \end{equation} \tag{ 34 }$$

which has the solution

$$\begin{equation} x={1\over 4}[(\sqrt{d+1}+1)^{1/3}- (\sqrt{d+1}-1)^{1/3}]^3\,, \end{equation} \tag{ 35 }$$

where

$$\begin{equation} d={\pi ^2n_B\over 288}\left({\hbar c\over [S_ku^{2/3}+S_pu^\gamma ]}\right)^3\,. \end{equation} \tag{ 36 }$$

We also include muons in our EOS, but they are not included in the above expressions for clarity.

For lack of a clear theoretical understanding of the nature of matter above saturation densities, we parameterize the high-density EOS by a pair of polytropes. Read et al. (2009) have shown that such a parameterization can effectively model a wide variety of theoretical model predictions in this density range. Specifically, for energy densities above a parameterized value ε₁, the EOS is

$$\begin{equation} P = K_1 \varepsilon ^{1+1/n_1}, \end{equation} \tag{ 37 }$$

where P is the pressure, n₁ is the polytropic index, and K₁ is a coefficient. Note that we are parameterizing the high-density EOS as polytropes in the total energy density ε rather than the number density, and that this EOS is assumed to be in beta equilibrium so it automatically includes leptonic contributions. Above a parameterized energy density, ε₂, we use a second polytrope

$$\begin{equation} P = K_2 \varepsilon ^{1+1/n_2} \,. \end{equation} \tag{ 38 }$$

We choose the transition densities ε₁ and ε₂ and the polytropic indices n₁ and n₂ as parameters. The coefficients K₁ and K₂ are determined by pressure and energy density continuity at the transition densities, with values limited such that 1600 MeV fm⁻³>ε₂>ε₁>150 MeV fm^-3. Note that 1600 MeV fm^-3 is either larger than or nearly as large as the central energy density of most configurations. In addition, we limit ε₁ < 600 MeV fm^-3, so that the parameters of the schematic EOS maintain a close connection to their usual definitions in terms of properties near the saturation density. Finally, we limit the polytropic indices with 0.2 < n₁ < 1.5 and 0.2 < n₂ < 2.0; our results are insensitive to this choice of limits.

While phase transitions are included in our models because they will appear like successive polytropes with different indices, some parameter sets do not imply any phase transition, as they simply model equations of state which have effective polytrope indices which vary with density. Models with more than two strong phase transitions will not be well reproduced by our parameterization, but it is not clear that these models are particularly realistic.

We thus have eight total EOS parameters, K, K', S_v, γ, n₁, n₂, ε₁, and ε₂, with limiting values summarized in Table 3. In addition to these parameter limits, some combination of parameters must be rejected because they are unphysical. Unphysical combinations include those in which

• 1.  
  the maximum mass is smaller than 1.66 M_☉, which is 2σ below the mass of PSR J1903+0327, 1.74 ± 0.04 M_☉ (Champion et al. 2008);
• 2.  
  the EOS becomes acausal below the central density of the maximum mass star;
• 3.  
  the EOS is anywhere hydrodynamically unstable, i.e., has a pressure that decreases with increasing density; and
• 4.  
  the maximum mass star has a maximum stable rotation rate less than 716 Hz, the spin frequency of the fastest known pulsar, Ter 5AD (Hessels et al. 2006). The spin frequency at which equatorial mass shedding commences is given to within a few percent by Haensel et al. (2009):

  $$\begin{equation} f_K \simeq 1.08 \left(\frac{M}{ M_{\odot }}\right)^{1/2} \left({10{\,\rm km}\over R}\right)^{3/2}{\rm kHz}\,. \end{equation} \tag{ 39 }$$

  There has been claimed possible evidence for a higher spin frequency in XTE J1239−285 (Kaaret et al. 2007), but this observation does not have strong statistical significance and has not been confirmed by subsequent observations.

In addition to these criteria, during the MC generation of neutron star mass sets for the Bayesian analysis, if one of the seven masses is larger than the maximum mass for the selected EOS, that realization is discarded and a new one is selected.

To show that this model accurately represents significantly different EOSs, we fit it to the Skyrme model SLy4 which gives relatively small neutron star radii (of order 10 km for M = 1.4 M_☉), and the field-theoretical model NL3 which gives rather larger radii, of order 15 km, for the same mass. These fits are shown in Figure 5, illustrating that the associated EOSs are reproduced to within a few percent. Our parameterization includes many extreme models, including those like NL3 with rather large neutron star radii. We will find below that such models are ruled out, however, by the observational data.

Zoom In Zoom Out Reset image size

We first analyze the EOS implied by observations for the case in which X-ray bursts are modeled assuming r_ph ≫ R. The histograms for the EOS parameters are given in Figures 6–7. The double-hatched (red) and single-hatched (green) regions outline the 1σ and 2σ confidence regions which are suggested by the simulation. The 1σ and 2σ parameter limits are summarized in Table 4 along with representative values obtained from terrestrial experiments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is remarkable that the inferred ranges for the schematic EOS parameters K, K', S_v, and γ are consistent with the values derived from nuclear experiments. Especially significant is the inferred low value for γ, a parameter which controls the pressure of the EOS in the range of 1–3 ε₀ and therefore the radii of neutron stars in the mass range 1 < M/M_☉ < 1.5 (Lattimer & Prakash 2001). The value of γ also controls the energy and pressure of pure neutron matter. Hebeler & Schwenk (2010) have shown that several recent studies indicate a convergence in predictions for neutron matter. For the saturation density used there, n₀ = 0.16 fm⁻³, the mean values of the neutron matter energy and pressure are 16.3 MeV and 2.5 MeV fm⁻³, which imply S_v = 32 MeV and γ = 0.28 with the parameterization of Equation (33), well within our predicted ranges.

One must take care in comparing experimental constraints directly to our results, however. The parameters determined experimentally are often "local" quantities, in the sense that they are properties of the EOS only at densities close to saturation. Also, the results from Tsang et al. (2009) mostly constrain the EOS in region near half the saturation density. We have used the schematic EOS parameters over a somewhat larger range of densities, from 1/2 to, typically, up to 2 or 3 times the saturation density.

The predicted skewness for the schematic EOS has a relatively small magnitude, centered at K' = −280 MeV. The most probable value for the transition between the schematic EOS and the first polytrope is around twice the saturation density, ε₀. The index of the first polytrope is sharply peaked around 0.5, which corresponds to a polytropic exponent γ₁ = 1 + 1/n₁ ≃ 3 which is rather large. Coupled with the small magnitude of the skewness, this implies a quite stiff EOS at supernuclear densities. Finally, note that, typically, n₂>n₁, indicating that the EOS softens at high densities. A summary of the results of the EOS parameters for both the schematic EOS and the polytropes is given in Table 4.

We next consider the results of parameter fitting when the X-ray burst data are modeled assuming r_ph = R. Results are summarized and compared to the previous case in Table 4. Although we do not show histograms for the parameters in this case, their behavior is similar to those of Figures 6 and 7 modulo the different means and variations of the two cases. It is significant that the small value for γ found for the case r_ph ≫ R is duplicated for r_ph = R, supporting a weak density dependence for the symmetry energy and a consequent small estimate for neutron star radii in the mass range 1–1.5 M_☉. We note that Özel et al. (2010) also concluded, on the basis of observations of X-ray bursts, that the implied neutron star radii were relatively small. Özel et al. (2010) obtained radii approximately 1 km smaller than we do, however, for r_ph = R (which are already 0.5–1 km smaller than those we obtain for r_ph ≫ R) because their analysis favors the high-redshift solution, and because they do not impose the same causality constraint on their MC sampling.

Although the low-density EOS is not strongly affected by the assumption that r_ph = R, at higher densities this choice leads to a softer EOS: the magnitude of the skewness K' and first polytropic index n₁ are both larger for the case r_ph = R (Table 4). The differences between these predicted EOS can be easily seen by referring to Figure 8. Each panel displays an ensemble of one-dimensional histograms over a fixed grid in one of the axes (note that this is not quite the same as a two-dimensional histogram). The bottom panels present the ensemble of histograms of the pressure for each energy density. This was computed in the following way: for each energy density, we determined the histogram bins of pressure which enclose 68% and 95% of the total MC weight. The locations of those regions for each one-dimensional histogram are outlined by dotted and solid lines, respectively, and these form the contour lines. These 1σ and 2σ contour lines give constraints on the pressure of the EOS as a function of the energy density as implied by the six neutron star data sets and are presented in Tables 5 and 6 along with the most probable EOS for the cases r_ph ≫ R and r_ph = R, respectively. The softer nature of the EOS in the r_ph = R case is most apparent at the highest energy density, for which the pressure is 10% less than in the r_ph ≫ R case.

Zoom In Zoom Out Reset image size

The upper panels of Figure 8 display the ensemble of histograms of the ratio of the pressure of the EOS to the pressure of SLy4 over a grid of energy density. The choice of the SLy4 EOS here is essentially arbitrary, and just assists in plotting the results since the pressure varies over a couple orders of a magnitude. The inferred pressure ratio from several Skyrme models is also shown. These Skyrme models were selected in Steiner & Watts (2009) because they are a representative set of the recommended models from Stone et al. (2003) which have symmetry energies which do not become negative at high density. In both the r_ph = R and r_ph ≫ R cases, there appears to be a softening of the EOS at low densities which is incompatible with some of the Skyrme models. Some models are apparently ruled out independently of assumptions about the radius of the photosphere in X-ray bursts: field-theoretical models like NL3 (Lalazissis et al. 1997), which have stiff symmetry energies; FSUGold (Todd-Rutel & Piekarewicz 2005), which has a softer symmetry energy but becomes too soft at high density; and the Akmal, Pandharipande, and Ravenhall EOS (APR; Akmal et al. 1998), which becomes too stiff at high densities. The nearly vertical nature of APR for energy densities just over 200 MeV fm^-3 is due to the phase transition in APR to matter which includes a π⁰ condensate. Even after applying a Gibbs phase construction, the pressure increases only very slowly with density in the small mixed phase region.

The upper panels in Figure 9 present our results for the predicted mass–radius relation according to our two assumptions regarding the photospheric radius for X-ray bursts. They give the ensemble of histograms of the radius over a fixed grid in neutron star mass. The width of the contours at masses below 1 M_☉ tends to be large because the available neutron star mass and radius data generally imply larger masses. In general, the implied M–R curve suggests relatively small radii, which is consistent with our conclusions regarding the softness of the nuclear symmetry energy in the vicinity of the nuclear saturation density. Tables 7 and 8 summarize the 1σ and 2σ contour lines from these panels, and give as well the most probable M–R curve. The assumption that r_ph = R implies smaller radii for neutron star masses greater than about 1.3 M_☉ and perhaps very small radii for masses in excess of 1.5 M_☉. This is suggestive of the onset of a phase transition above the nuclear saturation density or perhaps simply the approach to the neutron star maximum mass limit in this case. The mass and radius contour lines, determined in the same way as for Figure 8, are also given in the figure, and summarized in Tables 7 and 8. The radii in the r_ph = R case average about 1 km less than in the r_ph ≫ R case, except around 1.4 M_☉, where the difference is about 0.4 km.

Zoom In Zoom Out Reset image size

The lower panels in Figure 9 summarize the output probability distributions for M and R for the six neutron stars which were used in this analysis. Note that the scales have been modified, and these output probability distributions are much smaller than the input distributions given in Figures 1, 2, and 4. The combination of several neutron star mass and radius measurements with the assumption that all neutron stars must lie on the same mass–radius curve puts a significant constraint on the mass and radius of each object. Note that these output probability distributions closely match the implied M versus R curves presented in the upper panels of Figure 9. The tendency for smaller radii when M>1.3–1.4 M_☉ is apparent.

The M and R constraints for each object are given in Table 9, with their corresponding 1σ uncertainties. The three PRE burst sources suggest masses near 1.5 M_☉ and the radii near 11 km in the case where r_ph = R. We have already observed that this agreement is due to the extreme restrictions placed on the acceptability of points during the MC sampling. In the case r_ph ≫ R, the PRE burst sources are predicted to have a wider range of masses, from 1.3 to 1.6 solar masses. The quiescent LMXB masses tend to be smaller, but are strongly dependent on assumptions about the radius of the photosphere in the PRE burst sources. Particularly uncertain is the mass for X7. The observations imply a rather large value of R_∞, which is compatible with a small radius if the mass is large (r_ph ≫ R), or a large radius if the mass is small (r_ph = R). The stars in M13 and ω Cen show the opposite trend. They have small values of R_∞ which is compatible with small radii in the r_ph ≫ R case if the mass is relatively small. In general, the predicted radii of all stars range from about 10 km to 12.5 km, with the exception of X7 which may have a large radius. Note that a 13 km radius for an 0.8 M_☉ star is beyond the limit implied by rotation at 716 Hz and thus if the neutron star in X7 is observed to spin rapidly enough, a much larger mass and smaller radius would be implied instead.

The largest difference between the predicted equations of state between the r_ph = R and r_ph ≫ R cases is the high-density behavior. This leads to large differences in the predicted maximum masses. The probability distributions for the maximum neutron star mass are given in Figure 10, along with the associated 1σ confidence regions. The two probability distributions are arbitrarily normalized so that their peak is unity. These results are strongly dependent on assumptions of the photospheric radius at touchdown. The two-peaked behavior in the case r_ph = R suggests a possible phase transition could match the data and implies a maximum mass very close to the observed limit of 1.66 M_☉. This result is similar to that claimed by Özel et al. (2010), but there it is stated that the results are incompatible with a nucleonic EOS. Our results do not support this extreme interpretation. Although the neutron star radii implied by this analysis are small, they are not small enough to require a phase transition; for neutron stars of mass 1.4 M_☉, radii smaller than 10 km can be generated by purely nucleonic equations of state (Steiner et al. 2005).

Zoom In Zoom Out Reset image size

Normalization. The normalization factors obtained from the two PRE bursts in Özel et al. (2009) are A/(1 km² kpc^-2) = 1.04 ± 0.01 and 1.30 ± 0.01. These two measurements differ significantly, and this means that either the color-correction factors for the two bursts differ by roughly 6% or the geometry of the two bursts is different. It is not clear how to resolve this conflict. Özel et al. (2009) use a boxcar distribution with A = 1.16 km² kpc^-2 and ΔA = 0.13 km² kpc^-2, which represents the choice with the minimum possible uncertainty. We instead choose a Gaussian centered at A = 1.17 km² kpc^-2 with σ_A = 0.13 km² kpc^-2, which ensures that the two observations are included to within 1σ. This choice is consistent with the objects below, which also have normalization factors which are simulated with Gaussian distributions, but with notably smaller uncertainties.

Distance. EXO 1745–248 is located in the globular cluster Terzan 5. Ortolani et al. (2007) analyzed the distance to Terzan 5 using both the NICMOS instrumental magnitudes and the calibrations from Stephens et al. (2000) and Cohn et al. (2002). The combined distance estimation is D = 5.9 ± 0.9 kpc, where the uncertainty has been obtained from the standard deviation of the three different distance measurements, suggesting a Gaussian with the same standard deviation. On the basis that the distance measurement from the NICMOS instrumental magnitudes was independent of photometric calibrations and thus more accurate, Özel et al. (2009) used the NICMOS distance of D = 6.3 kpc with ΔD = 0.32 kpc. This choice assigns, however, a zero probability to a distance of 5.9 kpc, the central value suggested in Ortolani et al. (2007). Until this distance measurement is more clearly determined, we choose a Gaussian distribution with D = 6.3 kpc and σ_D = 0.6 kpc.

Touchdown flux. We employ the result from Özel et al. (2009), which is a Gaussian distribution with F_TD,∞ = 6.25 erg cm^-2 s⁻² and σ_F = 0.2 erg cm^-2 s⁻².

Hydrogen mass fraction. Galloway et al. (2008) noted that EXO 1745–248 has exhibited long Type-I bursts with estimated values of ∫dt F_persistent/∫dt F_burst ≈ GM/R/E_nuc ≈ 20–46. These long bursts do not always show a strong thermal evolution (D. K. Galloway 2009, private communication); but if they are indeed thermonuclear in origin, then the low ratio of persistent to burst fluence indicates H-rich fuel and larger values of X. In contrast to the determination of X from these long bursts, this object has also been identified as an ultracompact binary in Heinke et al. (2003), through a phenomenological assessment of its spectral properties, suggesting that X is small, 0 < X < 0.1. The radius expansion bursts have short durations, which would be consistent for ignition of a pure He layer (i.e., the hydrogen is consumed by the stable hot CNO cycle). We retain the full range of hydrogen composition, 0 < X < 0.7.

Normalization. The normalization factors for the four PRE bursts are A/(1 km² kpc^-2) = 3.267 ± 0.047, 3.302 ± 0.049, 3.258 ± 0.054, and 3.170 ± 0.047 (Güver et al. 2010a). We use the average from Güver et al. (2010a), 3.246 ± 0.024.

Distance. Güver et al. (2010a) give a Gaussian distribution with D = 5.8 kpc and σ_D = 2.0 kpc with a cutoff below 3.9 kpc. We also employ this result. We note that if the flux is indeed less than the Eddington value, then D < 4.36 kpc(M/1.4 M_☉)^1/2 for the central value of the touchdown flux.

Touchdown flux. Two of the four PRE bursts gave a value for the touchdown flux, F_TD,∞ = (15.58 ± 0.82) erg cm^-2 s⁻² and (15.14 ± 1.05) erg cm^-2 s⁻² (Güver et al. 2010a). The value (15.41 ± 0.65) erg cm^-2 s⁻² was obtained from the fit in Güver et al. (2010a).

Hydrogen mass fraction. The bursts in 4U 1608–522 suggest an accretion rate in of 3%–5% $\dot{M}_{\mathrm{Edd}}$, which suggests H ignition (Galloway et al. 2008); the brighter bursts from this system are of short duration, however, so it is likely that much of the hydrogen is consumed via the HCNO cycle prior to He ignition. As with EXO 1745–248, we use the full range 0 < X < 0.7.

Normalization. The three PRE bursts for which a normalization was obtained give A/(1 km² kpc^-2) = 0.8886 ± 0.0373, 0.9668 ± 0.0339, and 0.9040 ± 0.0200 (Güver et al. 2010b). We use the value quoted in Güver et al. (2010b), 0.9198 ± 0.0186.

Distance. 4U 1820–30 is in the globular cluster NGC 6624. Güver et al. (2010b) use a boxcar distribution from 6.8 kpc to 9.6 kpc, to reflect two distance measurements of (7.6 ± 0.4) kpc (Kuulkers et al. 2003) and (8.4 ± 0.6) kpc (Valenti et al. 2007). We employ a Gaussian distribution centered at 8.2 kpc with σ_D = 0.7 kpc since the 95% confidence regions for this Gaussian is the same as the range suggested by the boxcar in Güver et al. (2010b). A more recent determination (Dotter et al. 2010) gives an apparent distance of 10.2 kpc, which when corrected for extinction gives 8.1 kpc, consistent with our value described above.

Touchdown flux. Five of the bursts have a measured touchdown flux: F_TD,∞/(10⁻⁸ erg cm^-2 s⁻²) = 5.33 ± 0.27, 5.65 ± 0.20, 5.12 ± 0.15, 5.24 ± 0.19, and 5.42 ± 0.16 (Güver et al. 2010b). The fit in Güver et al. (2010b) gives (5.31 ± 0.10) × 10⁻⁸ erg cm^-2 s⁻², which we use here.

Hydrogen mass fraction. This object is likely an ultracompact binary with an H-poor donor (King & Watson 1986; Stella et al. 1987). Although evolutionary models do not exclude the possibility that the envelope could contain some H (X ≲ 0.1; Podsiadlowski et al. 2002), a comparison of burst properties with theoretical ignition models suggests H-poor fuel (Cumming 2003). We fix X = 0 for this source.