© ESO, 2016

1. Introduction

The equation of state (EoS) of the cold dense matter inside neutron stars (NS) has remained a mystery for decades. Experiments on Earth and theoretical many-body calculations have constrained the pressure-density relation of matter near the nuclear saturation densities. Recently, progress has also been made in measuring the NS radii (for a review, see Miller 2013; Özel 2013; Suleimanov et al. 2015) that allow us to constrain the behavior of the EoS at higher densities by inverting the Tolman-Volkoff-Oppenheimer (TOV; Tolman 1939; Oppenheimer & Volkoff 1939) structure equations. Furthermore, these measurements probe the phase diagram of dense quantum chromodynamics at lower temperatures and higher baryon densities than the measurements of, for example, ultrarelativistic heavy-ion collisions inside earthly laboratories (see, e.g., Lattimer 2012). One of the most promising candidates for obtaining accurate astrophysical mass-radius (M − R) measurements has been the thermal emission originating from NS surface layers. One possibility is to use the cooling of NS surface during type-I X-ray bursts from low-mass X-ray binary (LMXB) systems, where the cooling tail is shown to follow theoretical model predictions surprisingly well (Poutanen et al. 2014). In these systems, the NS is accompanied by a lighter star, usually a main sequence or evolved late-type star, that fills its Roche lobe and transfers material through an accretion disk onto the NS. After accumulating enough material, the fuel is rapidly burned in a thermonuclear explosion occurring below the surface in the NS ocean. Some of these bursts can be so energetic that the Eddington limit is reached, causing the entire NS photosphere to expand. These photospheric radius expansion (PRE) bursts can then be used to obtain mass and radius measurements by comparing the cooling tail of the burst to accurate theoretical predictions (for early work, see, e.g., Damen et al. 1990; van Paradijs et al. 1990; Lewin et al. 1993).

Recent studies (Suleimanov et al. 2011a; Poutanen et al. 2014; Kajava et al. 2014) have demonstrated that the X-ray burst cooling properties heavily depend on the accretion rate and spectral state of the source. The key finding was that care must be taken to select only those bursts that show “passive cooling”, meaning the ones that occur at the hard spectral state and at small accretion rate, where the extra heating from the in-falling material appears to be negligible. Kajava et al. (2014) also showed that the evolution of the blackbody normalization can be used as a trace to pin down the passively cooling bursts used in the M − R measurements. The soft state bursts, on the other hand, show only weak or completely non-existent evolution of the normalization that is in contradiction with the theoretical atmosphere model predictions.

In addition to only using the passively cooling bursts, it is possible to improve the analysis by using the information from the whole cooling tail by applying the so-called cooling tail method. In this recently developed method, the observed cooling track in the blackbody normalization K vs. the flux F (or rather K^{− 1 / 4} vs. F) plane is compared to the theoretical model evolution of the color-correction factor versus the luminosity (in units of the Eddington), f_c − L/L_Edd, that is the so-called color-correction curve. By comparing the whole cooling track to the models (with variable f_c) – in contrast, for example, to the “touchdown method” where f_c is assumed to be constant – we can infer more robust constraints from the data. This also allows us to circumvent the problematic issue of deciding where the Eddington limit is reached, and when the photospheric radius coincides with the NS radius (see, e.g., Steiner et al. 2010). In fact, because it is the curvature of the evolving color-correction factor that is used to constrain the Eddington flux, the method is valid even for bursts that do not reach the Eddington limit (see, e.g., Zamfir et al. 2012). So far, no work exists where the cooling tail method has been applied with this kind of strict hard-state burst selection criteria for many sources simultaneously. We therefore now pay special attention to the burst selection and choose only the most well-behaved hard-state PRE bursts for our analysis. We then use these bursts to constrain the mass and radii of three different NSs by applying the cooling tail method to them.

Using these constraints we can then go one step further and address the issue of unknown EoS of the cold dense matter inside neutron stars. To do this, we use Bayesian inference to derive empirical pressure-density and mass-radius relations based on our burst results. Using a Markov chain Monte Carlo (MCMC) algorithm, we fit the three M − R constraints from each source simultaineously and in parallel allowing us to put astrophysical constraints on some of the nuclear physics parameters, such as the symmetry energy and the pressure of neutron-rich matter at the saturation density ℒ (Lattimer & Steiner 2014). In addition, the combination of cooling tail observations and parameterized EoS allows us to make more accurate mass and radius measurements for each of the sources, indicating a new way of probing individual NS characteristics.

The paper is structured as follows: in Sect. 2, we present the methods used for the data reduction of the bursts. In Sect. 3, we use this data to obtain separate mass, radius, and distance constraints for the three sources in our sample. In the second part of the paper, in Sect. 4, we use Bayesian analysis to obtain the parameterized EoS. Finally, in Sect. 5, we discuss the constraints and compare our results to the measurements made previously.

2. Data

In our analysis we used the data from the Proportional Counter Array (PCA; Jahoda et al. 2006) instrument on board of the Rossi X-ray Timing Explorer (RXTE) satellite. Our sample consists of three neutron stars: 4U 1702−429, 4U 1724−307, and SAX J1810.8−2609. These sources were selected because they have been known to exhibit PRE bursts in the hard state (i.e., at low accretion rate where the NS is thought to cool passively without any external heating; Kajava et al. 2014). They also show the most robust evolution of the normalization down to very low luminosities enabling us to use them as clear examples of bursts to which the cooling tail method can be applied. These bursts are visually selected, based on the evolution of their normalization, to be in good agreement with the model used. We omit one long burst from 4U 1724−307 which has already been analyzed (Suleimanov et al. 2011a,b), since there is evidence that this particular burst might have a high metallicity content in its atmosphere (see Appendix A).

RXTE/PCA data were analyzed with the heasoft package (version 6.16) and response matrices were generated using pcarsp (11.7) task of this package. All data were fitted using the xspec 12.8.1 package (Arnaud 1996) where the recommended systematic error of 0.5% was added to the spectra (Jahoda et al. 2006). We identified X-ray bursts using a similar method to that in Galloway et al. (2008). The time-resolved spectra for the bursts were extracted using an initial integration time of 0.25 s. In order to maintain approximately the same signal-to-noise ratio the integration time was doubled every time the count rate decreased by a factor of $\sqrt{\mathrm{2}}$. The exposures were dead-time corrected following the approach recommended by the instrument team¹. The correction resulted in a roughly 10–15% increase in the peak flux, with the difference decreasing quickly as the observed count rate declined. A standard method of removing a 16 s spectrum taken from prior to the burst was used to account for the possible background emission (Kuulkers et al. 2002, and references therein). This standard method assumes the background (i.e., mainly the persistent emission) to be constant during the burst, even though this might not be the case. The changes due to the background emission are, however, not significant in the cooling phase (see Fig. 6 in Worpel et al. 2013). The differences in burst characteristics with and without this background subtraction was also checked and found to be negligible at least at high burst fluxes (i.e., at fluxes larger than 20% of the peak flux). These deadtime-corrected spectra were then fitted with a blackbody model multiplied by an interstellar absorption model with constant hydrogen column density N_H (value obtained from the literature, see Table 1). The best-fit parameters are the color temperature T_c and the normalization constant K ≡ (R_bb [ km ] /D₁₀)², where D₁₀ = D/ 10 kpc. From the corresponding χ² distributions (see Fig. 1) of each source, we also conclude that the data is sufficiently well described by the blackbody model. It should also be noted that the theoretical atmosphere model spectra cannot be perfectly fit by a (diluted) blackbody model either (Suleimanov et al. 2011b, 2012), so in reality we do not even expect the observed χ² distribution to be close to ideal. The bolometric flux was estimated using the cflux-model in the range 0.01–200 keV. All error limits were obtained using error -task in xspec.

Fig. 1 Reduced χ2 distributions for the blackbody spectral fits consisting of the points used in the cooling tail analysis. The dashed curve shows the theoretical expected χ2 distribution. Both obtained and theoretical distribution are normalized so that the encapsulated area of the curves is unity.

Some of these bursts show typical characteristics of a PRE: A peak in the normalization after a few seconds of the ignition. The evolution of the observed temperature should also show the characteristic double-peaked structure, arising from the cooling of the photosphere when it expands and the subsequent heating when it collapses back toward the surface due to the changing radiation pressure. The aforementioned signs of the expansion also indicate that the flux has reached (or exceeded) the Eddington limit during the burst. The moment when the atmosphere collapses back to the NS surface – that is the normalization reaches its minimum value K_td and the temperature its second peak – is defined as the touchdown. This also marks the beginning of the cooling phase where the subsequent evolution is dependent on the spectral state of the source, meaning that if the burst occurred in the hard state or in the soft state. In the hard state the normalization rises to a nearly constant level while the flux and the temperature continue to decrease for the rest of the burst. This increase of normalization is due to the changing color-correction factor f_c as we approximate the emerging spectrum with a diluted blackbody model as ${\mathit{F}}_{\mathit{E}}\mathrm{\approx }\frac{\mathrm{1}}{{\mathit{f}}_{\mathrm{c}}^{\mathrm{4}}}{\mathit{B}}_{\mathit{E}}\mathrm{\left(}{\mathit{f}}_{\mathrm{c}}{\mathit{T}}_{\mathrm{eff}}\mathrm{\right)}\mathit{,}$(1)where B_E is the blackbody function and T_eff is the effective temperature that is connected to the observed blackbody color temperature as T_c = f_cT_eff(1 + z)^-1 , where 1 + z = (1−2GM/Rc²)^−1/2 is the redshift. We stress here that the decrease in the color-correction factor during the cooling is a feature predicted by numerous atmosphere model computations (London et al. 1986; Lapidus et al. 1986; Suleimanov et al. 2011b, 2012). Consequently, the decrease of the color correction then leads into an increase in the (observed) normalization because ${\mathit{K}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{=}{\mathit{f}}_{\mathrm{c}}\mathit{A,}\mathit{A}\mathrm{=}\mathrm{\left(}{\mathit{R}}_{\mathrm{\infty }}\mathrm{\right[}\mathrm{km}\mathrm{]}\mathit{/}{\mathit{D}}_{\mathrm{10}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(2)where R_∞ = R(1 + z) is the apparent NS radius (Penninx et al. 1989; van Paradijs et al. 1990). On the other hand, for the soft-state bursts the normalization is nearly constant, contrary to the theory². It is also crucial to notice that the normalization value in the tail of the soft-state bursts is different to that observed for the hard-state bursts. In addition, the touchdown flux can also vary³ (see Fig. 1 in Kajava et al. 2014). Because of these differences, the burst selection becomes extremely important as our model assumptions are only valid if the NS surface alone is emitting. This seems to be valid only in the hard state (see Poutanen et al. 2014; Kajava et al. 2014,for more information about the soft vs. hard-state burst selection). The increasing emission area during the PRE phase (i.e., increase in the blackbody normalization K) before the touchdown is mostly related to the increase of the photospheric radius. In the cooling tail, one believes that the evolution of K happens just because of varying f_c with the constant actual photospheric radius equal to the NS radius. Therefore, the ratio of the normalizations in the expansion phase and the cooling tail is $\frac{{\mathit{K}}_{\mathrm{e}}}{{\mathit{K}}_{\mathrm{t}}}\mathrm{=}{\left(\frac{{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{t}}}{{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{e}}}\right)}^{\mathrm{4}}{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathrm{e}}\mathrm{)}{\mathit{R}}_{\mathrm{e}}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathrm{t}}\mathrm{)}{\mathit{R}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\mathit{,}$(3)where indices e and t refer to the expansion and tail, respectively. By taking f_c,t ≈ 1.4 in the tail and f_c,e ≳ 2 during the expansion (Pavlov et al. 1991; Suleimanov et al. 2012) and demanding that R_e>R_t (which also means (1 + z_e)R_e> (1 + z_t)R_t for $\mathit{R}\mathit{>}\frac{\mathrm{3}}{\mathrm{2}}\frac{\mathrm{2}\mathit{GM}}{{\mathit{c}}^{\mathrm{2}}}$), we end up with a simple PRE condition $\frac{{\mathit{K}}_{\mathrm{e}}}{{\mathit{K}}_{\mathrm{t}}}\gtrsim \frac{\mathrm{1}}{\mathrm{4}}\mathrm{·}$(4)When the normalizations are equal K_e = K_t, we get R_e ≳ 2 R_t (note also that z_t>z_e). What is remarkable with this condition is that the observed “expansion”, can be less than unity (compared to the tail) in order for the burst to have a PRE episode. The PRE condition can be transformed into a requirement that the observed peak normalization at the expansion phase K_e must be larger than the normalization at the touchdown K_td as both of them should have similar values of the color-correction factors. But this is equivalent to the standard criterion that K should have a local minimum (at the touchdown).

Fig. 2 Bolometric flux, temperature and blackbody normalization evolution during the hard-state PRE bursts. The black line shows the estimated bolometric flux (left-hand y-axis) in units of 10-7 erg cm-2 s-1. The blue ribbon shows the 1σ limits of the blackbody normalization (outer right-hand y-axis) in (km / 10 kpc)2. Similarly, the dashed orange line shows the color-corrected angular size (R∞/D10)2 using the same axis. The red ribbon show the 1σ errors for blackbody color temperature (inner right-hand y-axis) in keV. Highlighted gray area marks the region of the cooling tail used in the fitting procedures.

Time-resolved spectral parameters for the bursts in our sample are presented in Fig. 2. Additionally, we show the color-corrected angular size with the assumption f_c = 2 for the evolution before the touchdown. For the f_c values after the touchdown, we use the cooling tail model fits from Sect. 3 to correct for the varying color-correction factor. Because of the new PRE criteria we choose also to keep the bursts that show only modest photospheric expansion in our sample (bursts #1 and #3 from 4U 1702−429 and burst #1 from SAX J1810.8−2609). Even though the expansion phase in these bursts is not very long, it is clear from Fig. 2 that the subsequent cooling phase is still similar (compare, for example, the bursts #1 and #2 from 4U 1702−429).

3. The Bayesian cooling tail method

For our mass and radius analysis we use the cooling tail method (see Suleimanov et al. 2011a,b and Appendix A of Poutanen et al. 2014). With this method the information from the whole cooling track after the peak of the burst is used and the observed cooling is compared to the theoretical evolution predicted by passively cooling NS atmosphere models. To relate the observed data and the individual masses and radii of the NSs, we use Bayesian analysis (see also Özel & Psaltis 2015). Bayes’ theorem can be presented simply as (see, e.g., Grinstead & Snell 1997) $\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left(}\mathrm{ℳ}\mathrm{\right|}\mathrm{𝒟}\mathrm{)}\mathrm{=}\frac{\mathrm{Pr}\mathrm{\left(}\mathrm{𝒟}\mathrm{\right|}\mathrm{ℳ}\mathrm{\left)}\mathrm{Pr}\mathrm{\right(}\mathrm{ℳ}\mathrm{)}}{\mathrm{Pr}\mathrm{\left(}\mathrm{𝒟}\mathrm{\right)}}\mathit{,}$(5)where Pr(ℳ) is the prior probability of the model ℳ without any additional information from the data , is the prior probability of the data , is the conditional probability of the data given the model ℳ, and is the conditional probability of the model ℳ given the data . Here the last quantity is what we want and it gives us the probability that a given model is correct, given the data. We can extend Bayes theorem further by having many non-overlapping models ℳ_i, which exhaust the total model space ℳ. Then the relation can be written as $\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left(}{\mathrm{ℳ}}_{\mathit{i}}\mathrm{\right|}\mathrm{𝒟}\mathrm{)}\mathrm{=}\frac{\mathrm{Pr}\mathrm{\left(}\mathrm{𝒟}\mathrm{\right|}{\mathrm{ℳ}}_{\mathit{i}}\mathrm{\left)}\mathrm{Pr}\mathrm{\right(}{\mathrm{ℳ}}_{\mathit{i}}\mathrm{)}}{{\sum }_{\mathit{j}}\mathrm{Pr}\mathrm{\left(}\mathrm{𝒟}\mathrm{\right|}{\mathrm{ℳ}}_{\mathit{j}}\mathrm{\left)}\mathrm{Pr}\mathrm{\right(}{\mathrm{ℳ}}_{\mathit{j}}\mathrm{)}}\mathrm{·}$(6)In the cooling tail method the model space consists of four parameters: mass M and radius R of the NS, hydrogen mass fraction X in the atmosphere, and the distance D to the source. For the distance D a uniform flat prior distribution is assumed without any restrictions. Similarly, a uniform two-dimensional prior distribution is assumed for (M,R) space. We also take into account the causality requirement, R> 2.824GM/c² (Lattimer 2012) and limit the mass so that it lies between 0.8 M_⊙<M< 2.5 M_⊙. For the hydrogen fraction X, a Gaussian prior distribution is used and is discussed in more detail later on.

The model parameters can be combined into two new parameters, related more closely to the color-correction curve fitting. The first one is the Eddington flux ${\mathit{F}}_{\mathrm{Edd}}\mathrm{=}\frac{\mathit{GMc}}{{\mathit{D}}^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{e}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{,}$(7)where κ_e = 0.2(1 + X) cm² g^-1. The second parameter is related to the apparent (non-color-corrected) angular size (2). These parameters then relate our observed flux to the (relative) luminosity as F/F_Edd ∝ L/L_Edd (where L_Edd is the Eddington luminosity) and observed blackbody normalization to the color-correction as K^{− 1 / 4} = f_cA. We also note here that it is possible to assume uniform priors for F_Edd and A (in contrast to assuming uniform flat distribution for (M,R) space; see Appendix B) as in Suleimanov et al. (2011b, 2012), Poutanen et al. (2014).

As our actual model, we use the recently computed hot neutron star atmosphere models (Suleimanov et al. 2012) that account for the Klein-Nishina reduction of the electron scattering opacity using an exact relativistic Compton-scattering kernel (Poutanen & Svensson 1996). These models give us the color-correction as a function of relative luminosity, f_c(ℓ ≡ L/L_Edd). Model uncertainty is taken into account by considering a boxcar distribution of a width of (1 ± ϵ) × f_c (where ϵ = 0.03) centered around the “real” value (see Suleimanov et al. 2011b, 2012 for a discussion of model uncertainties). Compositions considered are a pure helium (He) and a solar composition of H and He with sub-solar metal abundance of Z = 0.01 Z_⊙ (SolA001). It seems that Z< 0.1 Z_⊙ in the surface layers of the NS, because in the opposite case the atmosphere model predicts a drop of around 20% in the f_c (and correspondingly in K^{− 1 / 4}) at F ~ 0.3F_Edd (Suleimanov et al. 2011b, 2012), which is not observed. Alternatively, the absence of the drop in the low luminosities might be due to an extra heating because of the accretion that starts again after the burst. In any case, owing to these uncertainties in the low luminosity regime, we neglect this area from the fit and consider only the regime where F> 0.2F_td, where F_td is the flux at the touchdown point.

We also relax the assumption of having a fixed hydrogen mass fraction X (X = 0.738 for SolA001 and X = 0 for He models) and use Gaussian priors around the exact values with one sigma tails of 0.05. Both compositions are tested for each source and the physically more justified value is selected. Note also that this selection is simple as a wrong composition gives R ≲ 6 km or R ≳ 18 km. Strictly speaking this should be taken into account by using atmosphere models that are computed with the corresponding hydrogen fractions but the models (i.e., color-correction factors f_c) depend so weakly on this value (as our one sigma limits were $\mathit{X}\mathrm{=}{\mathrm{0}}_{-0.0}^{\mathrm{+}\mathrm{0.05}}$ or X = 0.738 ± 0.05) that it is possible to neglect the effect that this has on the F_Edd and A (see, e.g., Suleimanov et al. 2012,where the difference is relatively small even for X = 0 compared to X = 1). For the M, R, and D this, however, has some non-negligible effects that introduces a small scatter of about 5 per cent to the posterior distributions around the “exact” value. The uncertainty in the hydrogen fraction is also backed up by physical arguments because for hydrogen-poor companions (in the case of X = 0, i.e., He models) the evolutionary models do not exclude the possibility of the envelope containing some fraction of H (X ≲ 0.1; Podsiadlowski et al. 2002). The value of solar ratio of H/He is relatively accurately measured but here the uncertainties are possible and due to the possible stratification on top of the NS and/or because of the light ashes from the previous bursts that may stratify and accumulate to the surface. In the end, however, one should remember that the value of the hydrogen mass fraction for each model is still just a model assumption. By selecting a Gaussian prior around the presumed value we do weaken the effect that this selection has, but we are unable to remove it completely. If no assumption for the hydrogen mass fraction were to be made, we could not infer the radius at all. Reassuringly, however, the end results do seem to gather around similar radii, which means that our assumed X values were close to the real values.

In our method the data is constructed as a set of N points $\mathrm{\left(}{\mathit{F}}_{\mathit{i}}\mathit{,}{\mathit{K}}_{\mathit{i}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{\right)}$ obtained from the blackbody fits, starting from the touchdown (i = 1) and continuing down to 0.2F_td (i = N). The lower limit here is selected so that we can maximize the available data (in contrast to 1 / e F_td used in previous work) as the theoretical models nicely follow the data. Below the 0.2F_td limit the background emission can start to play too important a role so we choose to leave it out even though some of the bursts might follow the model even beyond this. These data points are then transformed into two-dimensional probability density distributions by assuming a Gaussian measurement error model. Next we implicitly assume that all of the data distributions are independent of each other and also independent of the model assumptions and prior distributions. We then assume that the conditional probability of the data given the model, , is proportional to the product of every individual probability $\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left(}\mathrm{𝒟}\mathrm{\right|}\mathrm{ℳ}\mathrm{)}\mathrm{\propto }\underset{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{S}}}{\mathrm{􏽙}}\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{(}{\mathrm{𝒟}}_{\mathit{i}}\mathrm{\left|}\mathrm{ℳ}\mathrm{\right)}\mathit{.}$(8)Each separate probability is assumed to be proportional to the path-integral evaluated through the two-dimensional “data space”, (F,K^{− 1 / 4}), along the color-correction curve as $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left[}{\mathrm{𝒟}}_{\mathit{i}}\mathrm{\right|}\mathrm{ℳ}\mathrm{\left(}\mathit{M,R,D,X}\mathrm{\right)}\mathrm{]}\mathrm{\propto }\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array} \frac{\mathrm{1}}{\mathrm{𝒩}}{\mathrm{\int }}_{{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{lo}}}^{{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{hi}}}\mathrm{d}{\mathit{f}}_{\mathrm{c}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}{\mathrm{\int }}_{{\mathrm{ℱ}}_{\mathrm{c}}}{\mathrm{𝒟}}_{\mathit{i}}\mathrm{\left(}\mathit{F,}{\mathit{K}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{\right)}\hspace{0.17em}\mathit{J}\left(\frac{\mathit{\ell ,}{\mathit{f}}_{\mathrm{c}}}{\mathit{F,}{\mathit{K}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}}\right)\mathrm{d}\mathit{s}\end{array}& & \end{array}$(9)where f_c,lo,hi = (1 ± ϵ) × f_c(ℓ) are the lower and upper limits of the prior boxcar distribution of the color-correction f_c evaluated at relative luminosity ℓ and where ϵ = 0.03, ℱ_c = ℱ_c(F_Edd,A) is the color-correction curve in (ℓ,f_c) space, $\mathit{J}{}^{\mathrm{\right(}}\frac{\mathit{\ell ,}{\mathit{f}}_{\mathrm{c}}}{\mathit{F,}{\mathit{K}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}}^{\mathrm{\left)}}$ is the Jacobian transforming the path from the model space to the observed data space (using Eqs. (2) and (7)) and ds is the line element of ℱ_c. The path-integral is also area-normalized (or length-normalized if ϵ = 0) with the factor that is defined as the aforementioned integral (9) without the data . This normalization takes into account the variable maximum ℓ that evolves as a function of log g. We also note that the presented Bayesian path-integral formalism is related to the two-dimensional frequentist formulation of the normalized minimum distance (see Eq. (2) in Poutanen et al. 2014).

In Bayesian interference, model parameters are determined using a marginal estimation where the posterior probability for a model parameter p_j is given by $\begin{array}{ccc}& & \begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left[}{\mathit{p}}_{\mathit{j}}\mathrm{\right|}\mathrm{𝒟}\mathrm{\left]}\mathrm{\right(}{\mathit{p}}_{\mathit{j}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{V}}\mathrm{\int }\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{[}\mathrm{𝒟}\mathrm{\left|}\mathrm{ℳ}\mathrm{\right]}\\ & & \mathrm{\times }\mathrm{d}{\mathit{p}}_{\mathrm{1}}\mathrm{d}{\mathit{p}}_{\mathrm{2}}\mathit{...}\mathrm{d}{\mathit{p}}_{\mathit{j}\mathrm{-}\mathrm{1}}\mathrm{d}{\mathit{p}}_{\mathit{j}\mathrm{+}\mathrm{1}}\mathit{...}\mathrm{d}{\mathit{p}}_{{\mathit{N}}_{\mathit{P}}}\mathit{,}\end{array}$(10)where N_P = 4 (corresponding to M, R, D, and X) and $\mathit{V}\mathrm{=}\mathrm{\int }\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left[}\mathrm{𝒟}\mathrm{\right|}\mathrm{ℳ}\mathrm{\left]}\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\right[}\mathrm{ℳ}\mathrm{]}{\mathrm{d}}^{\mathit{N}}\mathrm{ℳ}\mathit{,}$(11)without the model priors that determine the integration limits. Then the one-dimensional function represents the probability that the jth parameter will take a particular value given the observational data .

Fig. 3 Left panel: combined cooling tail in the F ∝ L/LEdd vs. K− 1 / 4 ∝ fc plane with 1σ error limits presented by crosses. Gray crosses show the burst evolution before the touchdown. Best-fit theoretical atmosphere models are shown by blue (SolA001) or red (He) curves. Right panel: temperature evolution of the bursts. Blackbody color temperature Tc is shown for each cooling tail with black crosses. Red (He) or blue (SolA001) crosses show the color-corrected temperatures Teff(1 + z)-1. Dashed lines show a powerlaw with an index of 4 corresponding to the F ∝ T4 relation. Highlighted gray area marks the region of the cooling tail used in the fitting procedures.

The best-fit atmosphere models are presented in Fig. 3 (left panel). In addition, the right panel of the figure depicts the observed color temperatures T_c and the corrected effective temperatures T_eff(1 + z)^-1 for a distant observer. Here the temperature is seen to follow the $\mathit{L}\mathrm{\propto }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}$ law, that is, it shows passive cooling.

Corresponding best-fit values of the model parameters F_Edd and A are listed in Table 2 along with the 1σ and 2σ confidence limits of the posterior distributions. After marginalizing over the radius R, mass M, and hydrogen mass fraction X we get the posterior distribution for the distance D too, that is also listed in Table 2. Figure 4 shows the two-dimensional mass and radius probability posterior distributions of our analysis. The obtained contours are arched and elongated along the curves of constant Eddington temperature⁴${\mathit{T}}_{\mathrm{Edd}\mathit{,}\mathrm{\infty }}\mathrm{=}{\left(\frac{\mathit{gc}}{{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{\kappa }}_{\mathrm{e}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{z}}\mathrm{=}\mathrm{1.14}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\hspace{0.17em}\mathit{A}\hspace{0.17em}{\mathit{F}}_{-7}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{K}\mathit{,}$(12)where g is the surface gravity, σ_SB is the Stefan-Boltzmann constant and F_-7 = F_Edd/ 10^-7 erg cm^-2 s^-1. The width of the contours are defined by the errors in T_Edd,∞ that originate from the uncertainty in F_Edd, A, and X. Our location on this curve is defined by the distance D. Because the distance is free to vary in our analysis, we end up with the aforementioned arched posteriors.

Contours from 4U 1724−307 and SAX J1810.8−2609 are seen to be almost identical with the radius constrained between about 11−13 km (for a more strict lower mass limit of M> 1.1 M_⊙) where the largest scatter is being produced by the unknown distance. Both of these sources are also best explained by a solar-like composition with an almost zero metallicity (SolA001 model). For 4U 1702−429 the radius is seen to be in the same range if mass ≲1.8 M_⊙ is assumed. The model in this case consists of a pure helium composition.

Fig. 4 Mass-radius constraints for the sources from the hard state PRE bursts. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upper-left region is excluded by constraints from the causality and general relativistic requirements (Haensel et al. 2007; Lattimer & Prakash 2007).

The choice of pure helium composition causes the Bayesian model to favor larger mass. This happens because F_Edd ∝ M/ (1 + X), so when the hydrogen fraction increases or decreases (because of the possible uncertainty in the X parameter) the change can be balanced by also increasing or decreasing M. With the solar composition of H and He our X priors are symmetrical and this effect is canceled out. In the case of He models the lower-limit of X = 0 makes the hydrogen fraction prior distribution asymmetrical and hence the bias for larger mass is present. We note, however, that in addition to this bias there is a slight preference for the 4U 1702−429 to favor larger log g values and hence larger masses. A similar effect is also present in the M vs. R posteriors because of our choice of flat distance prior. As F_Edd ∝ M/D² we end up oversampling the small distance values of our flat prior that can then create a bias that favors smaller mass. This effect is visible as an increased probability density around M ≲ 1.5 for 4U 1724−307 and SAX J1810.8−2609. Because of these model biases one should be careful not to give too much emphasis to the specific values of the M vs. R distributions presented in Fig. 4, but to focus more on the overall structure of the posteriors given by the confidence contours.

It is also possible to set some constraints to the unknown distance by using the cooling tail method. From the values of F_Edd and A, inferred from the data, we can derive the maximum distance where we still have M and R solutions (see Appendix A of Poutanen et al. 2014, for the full set of equations) as ${\mathit{D}}_{\mathrm{10}}\mathrm{\le }{\mathit{D}}_{\mathrm{10}\mathit{,}\max }\mathrm{=}\frac{\mathrm{1.77}\mathrm{\times }{\mathrm{10}}^{-2}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{X}\mathrm{)}{\mathit{A}}^{\mathrm{2}}{\mathit{F}}_{-7}}\mathrm{·}$(13)On the other hand, our lower-limit of the mass prior distribution (M_min = 0.8 M_⊙) also sets a lower-limit on the distance when combined with F_Edd and A. From these two constraints we are then also able to put some limits on the distance to the NS.

4. EoS constraints

The final goal of the mass and radius measurements is to constrain the pressure-density relationship of the cold dense matter. Here we use these new mass and radius constraints from the three NSs to probe the EoS by applying a Bayesian analysis to the data.

Here the model space consists of EoS parameters p_{i = 1,...,Np} in addition to the values of neutron star masses M_{i = 1,...,NM} with a total dimensionality of our model space as N = N_p + N_M. The total number of neutron stars in our sample is N_M = 3 and the number of EoS parameters N_p depends on our initial choice of the model (see Sect. 4.1).

The data is now constructed as a set of N_M probability distributions, in the (M,R) plane obtained from the cooling tail posteriors presented in Sect. 3. All of these distributions are normalized to unity by computing the integral ${\mathrm{\int }}_{{\mathit{M}}_{\mathrm{low}}}^{{\mathit{M}}_{\mathrm{high}}}\mathrm{d}\mathit{M}{\mathrm{\int }}_{{\mathit{R}}_{\mathrm{low}}}^{{\mathit{R}}_{\mathrm{high}}}\mathrm{d}\mathit{R}\hspace{0.17em}{\mathrm{𝒟}}_{\mathit{i}}\mathrm{\left(}\mathit{M,R}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{\forall }\mathit{i}\mathit{.}$(14)This normalization ensures that each source is treated equally in the analysis. As integration limits we use the same constraints as in the cooling tail method analysis where M_low = 0.8 M_⊙, M_high = 2.5 M_⊙, R_low = 5 km and R_high = 18 km.

Next we implicitly assume again that all of the new data distributions are independent of each other and also independent of the model assumptions and prior distributions. We then assume that the conditional probability of the data given the model, , is proportional to the product of the probability distributions evaluated at some mass M_i and radius (determined from the model) R(M_i) so that $\begin{array}{ccc}\begin{array}{c}\mathrm{Pr}\end{array}\mathrm{\left[}\mathrm{𝒟}\mathrm{\right|}\mathrm{ℳ}& & \mathrm{\left(}{\mathit{p}}_{\mathrm{1}\mathit{,}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{P}}}\mathit{,}{\mathit{M}}_{\mathrm{1}\mathit{,}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{M}}}\mathrm{\right)}\mathrm{]}\mathrm{\propto }\\ & & \underset{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{M}}}{\mathrm{􏽙}}{\mathrm{𝒟}}_{\mathit{i}}\mathrm{\left(}\mathit{M,R}\mathrm{\right)}{\mathrm{|}}_{\mathit{M}\mathrm{=}{\mathit{M}}_{\mathit{i}}\mathit{,R}\mathrm{=}\mathit{R}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathrm{\right)}}\mathit{.}\end{array}$(15)For the model parameters n_p + n_M we assume a uniform distribution (i.e., weakly informative physical priors) except for a few physical constraints described below.

In order to obtain all of the posterior probability distributions for the parameters we use the publicly available bamr code (Steiner 2014a). The TOV solver and data analysis routines were obtained from the O₂scl library (Steiner 2014b). To solve the integral (10)the code uses the Metropolis-Hastings algorithm to construct a Markov chain to simulate the distribution . For each point, an EoS parameter p_i and the neutron star mass M_i is generated. A corresponding radius curve R(M) (and radius R_i) is then obtained by solving the TOV equations. From these three masses and radii, the weight function is computed and the point is either accepted or rejected according to the Metropolis algorithm.

4.1. EoS parameterization

When building our EoS we follow the work by Steiner et al. (2015) and separate our density into three effective regimes: the crust, and regions below and above the nuclear saturation density n₀. We assume nuclear binding energy of E_nuc(n₀) = −16 MeV and saturation density of 0.16 fm^-3, typical values obtained from Kortelainen et al. (2010). In mass densities these values correspond to about 2.7 × 10¹⁴ g cm^-3. The nuclear symmetry energy (the difference between energy per baryon of neutron matter and that of the nuclear matter)⁵ is denoted as , where n_B is the baryon number density and . The pressure of neutron-rich matter at the saturation density n₀ is denoted by . In addition, the nuclear masses and theoretical models imply a correlation between ℒ and (Lattimer & Steiner 2014) and thus we additionally restrict these parameters as enclosing the aforementioned constraints⁶. The transition density from the crust to core is fixed to be at nuclear baryon density of 0.08 fm^-3. We note that the crust model has almost no effect to the resulting radii (nor does the fixed transition density) as the results are much more dependent on the high density behavior of our EoS.

Near the nuclear saturation density, below the core, we employ the Gandolfi-Carlson-Reddy (GCR) quantum Monte Carlo model (Gandolfi et al. 2012) that takes into account the three-body forces between the particles in the high density matter. The GCR results are accurately approximated by a two polytrope model given in terms of the energy E for the neutron matter at some nucleon number density n as $\mathit{E}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{=}\mathit{a}{\left(\frac{\mathit{n}}{{\mathit{n}}_{\mathrm{0}}}\right)}^{\mathit{\alpha }}\mathrm{+}\mathit{b}{\left(\frac{\mathit{n}}{{\mathit{n}}_{\mathrm{0}}}\right)}^{\mathit{\beta }}\mathrm{+}\mathit{m}n$(16)with coefficients (a and b) and exponents (α and β) constrained by QMC calculations and where m_n is the nucleon mass. The parameters of the first term, a and α, are mostly sensitive to the low density behavior of the EoS and are responsible for the two-nucleon part of the interaction. The limits of a and α are chosen so that we take into account all of the possible models from Gandolfi et al. (2012; see Table 3). On the other hand, the parameters of the second term, b and β, are sensitive to the high density physics and control the three-nucleon interactions. Furthermore, in our analysis we re-parameterized b and β in terms of and ℒ. Near the nuclear saturation density n₀ the symmetry energy of the neutron matter can be obtained from (16)as $\mathrm{𝒮}\mathrm{\equiv }\mathrm{𝒮}\mathrm{\left(}{\mathit{n}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathit{E}\mathrm{\left(}{\mathit{n}}_{\mathrm{0}}\mathrm{\right)}\mathrm{-}{\mathit{E}}_{\mathrm{nuc}}\mathrm{\left(}{\mathit{n}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{16}\mathrm{MeV}\mathrm{+}\mathit{a}\mathrm{+}\mathit{b,}$(17)where E_nuc(n₀) = −16 MeV is the previously mentioned nuclear binding energy at the saturation density. For the pressure at the saturation density we obtain $\mathrm{ℒ}\mathrm{\equiv }\mathrm{3}{\mathit{n}}_{\mathrm{0}}\frac{\mathrm{d𝒮}\mathrm{\left(}\mathit{n}\mathrm{\right)}}{\mathrm{d}\mathit{n}}{|}_{\mathit{n}\mathrm{=}{\mathit{n}}_{\mathrm{0}}}\mathrm{=}\mathrm{3}\mathrm{(}\mathit{a\alpha }\mathrm{+}\mathit{b\beta }\mathrm{)}\mathit{.}$(18)We also restrict the GCR model only up to nuclear saturation density as the validity of the model might not hold if a phase transition is present.

Above the saturation density n₀, a set of three piecewise polytropes are used and referred to as “model A”, similar to Steiner et al. (2013, 2015). In this way, when parameterizing the high-density EoS we introduce three continuous power laws defining the pressure as $\mathit{P}\mathrm{\propto }{\mathit{ϵ}}^{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\mathit{n}}\mathit{,}$(19)as a function of the energy density ϵ. It has been shown that it is possible to model a wide range of theoretical model predictions with these kinds of simple polytropes with a typical rms error of about 4% when compared to the actual numerical counterparts (Read et al. 2009). In practice we can mimic theoretical models with up to three phase transitions because they will appear as successive polytropes with different indices. Model A has five free parameters: the first transition energy ϵ₁ and the first polytrope index n₁, the second transition energy ϵ₂ and the second polytrope index n₂, and a third polytrope index n₃ (see Table 3 for the hard limits). Of course, we also require that ϵ₂>ϵ₁ in order to avoid double-counting of the parameter space. We have only five parameters (in contrast to six) because the transition to the first polytrope is already fixed by the EoS at the saturation density.

An alternative for the high density EoS is a piecewise linear model referred to as model C by Steiner et al. (2013, 2015). Here the low density EoS is used up to 200 MeV fm^-3, after which four line segments are considered in the ϵ vs. P plane at fixed energy densities of 400, 600, 1000, and 1400 MeV fm^-3. The linear relation between the two last regimes is extrapolated to higher densities, if needed. Model C has four free parameters: ΔP₁ = P(ϵ = 400) − P(ϵ = 200), ΔP₂ = P(ϵ = 600) − P(ϵ = 400), ΔP₃ = P(ϵ = 1000) − P(ϵ = 600), and ΔP₄ = P(ϵ = 1400) − P(ϵ = 1000). These pressure changes are always relative to the value of the pressure at the previous fixed point, effectively showing how sharply the pressure changes as a function of the energy density. This second alternative model tends to favor strong phase transitions in the core so it is interesting to study the resulting differences between it and the more smoothly evolving polytropic model.

In total our EoS models have nine (QMC + model A) or eight (QMC + model C) free parameters. In addition to the limits set on the priors (summarized in Table 3) some combinations of the parameters are rejected on a physical basis. More precisely, we ensure that all

• 1.

  mass-radius curves produce 2 M_⊙ NS in line with the recent pulsar measurements (Antoniadis et al. 2013; Demorest et al. 2010),

• 2.

  obtained EoS are causal, meaning that dP/ dϵ> 0, and

• 3.

  EoSs are hydrodynamically stable everywhere, meaning that their pressure increases with increasing energy density.

In addition, during the computations, if any of the three masses obtained is larger than the maximum mass for the selected EoS, that realization is discarded and a new one is generated.

4.2. EoS parameter results from the Bayesian analysis

Fig. 5 Histograms for the posterior distributions of the high-density model A and model C EoS parameters as implied by the data. In the first panel, the inset shows a magnified view of the histogram near n1 ≈ 0.5. Red shading corresponds to the 68% and blue to the 95% confidence regions of the parameters. Limits of the figures correspond to the hard limits set to the parameter prior distributions (see Table 3).

The most probable value and the corresponding 1σ and 2σ limits for the EoSs are summarized in Table 4 (for the computation of the confidence regions, see Lattimer & Steiner 2014). We find that the posterior distributions for a and α, corresponding to the low-density EoS behavior (which is dominated by two-body interactions), are almost flat. Thus, the neutron star observations do not constrain these parameters, as found previously (Steiner & Gandolfi 2012). We find that the derivative of the nuclear symmetry energy ℒ is only weakly constrained. However, we do find a somewhat stronger constraint on the symmetry than what has been obtained previously (Steiner et al. 2010). The origin of this constraint is the combination of the neutron star data with the correlation between and ℒ found in quantum Monte Carlo results (Gandolfi et al. 2012). It is also remarkable that with both high-density models, model A and model C, the symmetry energy is constrained around , that is in good agreement with earthly measurements (, Klüpfel et al. 2009). We note, however, that the parameters obtained here are to be considered as “local” quantities, as they are properties of the EoS only at densities close to the saturation density.

Histograms for the posterior distributions of the high-density parameters of the model A are presented in Fig. 5. The index of the first polytrope n₁ peaks sharply around 0.5 corresponding to a polytropic exponent γ₁ = 1 + 1 /n₁ ~ 3. The large value implies a rather stiff EoS at supranuclear densities. This first polytrope, corresponding to the n₁ index, is seen to continue all the way up to the first transition density at ϵ₁ ≈ 700 ± 150 MeV fm^-3, that is around four times the saturation energy density ϵ₀. In some realizations the transition occurrs as early as ϵ₁ ≈ 200 MeV fm^-3 but these cases correspond to the first sharp peak seen at n₂ ≈ 0.5, meaning that, practically, we have only two polytropes spanning our energy density range. In this case, the role of the first polytrope is superseded by the second segment corresponding to index n₂. In the opposite case, where all three polytropes span the grid the second polytrope index has values around 1.5 (polytropic exponent γ₂ ~ 1.7) with a long tail extending all the way up to around 8.0 (i.e., to the upper-limit of our prior). What this means is that the data can not constrain the high-density behavior of the EoS very well. As n₂>n₁, it, however, indicates that some softening of the EoS is present at higher densities. The third polytropic index n₃ is only weakly constrained to be ≳1 (peaking around 1.5) indicating either no phase transitions at all or additional softening as n₃>n₂>n₁.

Fig. 6 Obtained EoS constraints in the M − R (left panel) and in the P − ϵ plane (right panel). Upper panels correspond to the QMC + model A and bottom panels to the QMC + model C EoS. Red color indicates the probability density and black lines show the 68% (dotted) and 95% (solid) confidence limit contours.

Alternatively to the rather smoothly behaving polytrope model, we take the piecewise linear model C that can show strong phase transitions. Instead of extending the n₂ and n₃ parameter upper-limits of model A, that would create an apparent bias for softer EoS, we can characterize the effect of softer EoSs by applying this piecewise parameterization to the data, too. The histograms of the posterior distributions of the obtained pressure differences (at the fixed transition densities ϵ = 400, 600, 1000, and 1400 MeV fm^-3) are presented in Fig. 5. For the first segment from ϵ_trans = 200 to ϵ₁ = 400 MeV the difference is tightly constrained around ΔP₁ = 15 MeV fm^-3 with an almost symmetric Gaussian distribution with tails of 1σ ≈ 7 MeV fm^-3. This introduces a strong phase transition to the EoS as the pressure changes by only a little. After this possible transition the EoS hardens out and for ΔP₂ and ΔP₃ (at 600 and 1000 MeV fm^-3) the posterior distributions are considerably asymmetric toward larger pressure changes peaking at 180 and 350 MeV fm^-3, respectively. For the final possible line segment, the EoS appears almost unconstrained. The sharp cutoff at higher pressures for ΔP₂, ΔP₃, and ΔP₄ appears because we rule out EoSs where the speed of sound exceeds the value of the speed of light. The modest high-value tail for ΔP₄ originates from a small group of EoSs where the central value does not exceed 1000 MeV fm^-3 and hence every value of the parameter is allowed.

4.3. Predicted EoS

Fig. 7 Individual mass and radius constraints for the three neutron stars used in the analysis. The left-hand panel shows the projected mass and the middle panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence regions of these parameters. Right-hand panels show the full 2D mass and radius probability distributions. Contours of 68% (dotted, black line) and 95% (solid, black line) confidence regions are also shown.

The predicted EoS obtained from the X-ray burst data is shown in Fig. 6 in M − R and P − ϵ planes. Each panel in the figure displays an ensemble of one-dimensional (1D) histograms over a fixed grid in one of the axes (note that this is not quite the same as a 2D histogram). The right-hand 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 location of those regions for each 1D histogram are outlined by dotted and solid curves, respectively, and these form the contour lines. These 1σ and 2σ contour lines give constraints on the pressure as a function of the energy density as implied by the three NS data sets (see also Tables C.1 and C.2). Very high density behavior between the models A and C are seen to be similar as both of the models appear rather soft in this regime. At these very high energy density regions it is actually the maximum mass requirement that constraints the pressure evolution (see Steiner et al. 2013, for more extensive discussion). Most striking difference occur at lower energy densities where the sharp phase transition in the QMC + model C EoS is seen to produce a large scatter in the pressure at supranuclear densities (around ϵ ≈ 400 MeV fm^-3).

Similarly, the left-hand panels of Fig. 6 present our results for the predicted mass-radius relations. These panels present the ensemble of histograms of the radius over a fixed grid in neutron star mass with 1σ and 2σ constraints presented with dotted and solid lines, respectively, similarly to the right-hand panels. See also Tables C.3 and C.4 that summarize these contour lines and give the most probable M − R curve. The width of the contours at masses higher than about 1.8 M_⊙ tends to be large because the available NS mass and radius data in our sample generally imply smaller masses which in turn leads to weaker constraints. The obtained EoS for the QMC + model A has a predicted radius that is almost constant over the whole range of viable masses. The radius is constrained between 11.3−12.8 km for M = 1.4 M_⊙ (2σ confidence limits). Constraints this strong are obtained because the combination of weak (or non-existent) phase transitions and the NS mass and radius measurements from the cooling tail method compliment each other well: cooling tail measurements are elongated along the constant Eddington temperature curve that stretches from small mass and small radius to large mass and large radius. On the other hand, the assumption of weak phase transitions in the EoS forces the radius to be almost independent of mass. This assumption of constant radius then eliminates some of the uncertainties present in the cooling tail measurements (mostly due to the unknown distance) as each individual measurement is required to have (almost) the same radius.

With the QMC + model C, on the contrary, the first phase transition at supranuclear densities produce a slightly skewed mass-radius curve to compensate the cooling tail burst data that is elongated along the constant Eddington temperature. With this possible phase transition present in the EoS the mass-radius curve is then able to support high-mass NSs with radius of about R ≈ 11.6 km and low-mass stars with smaller radii of around R ≈ 11.3 km simultaneously. The phase transition also causes large scatter on the radius below 1 M_⊙ as the exact location of the turning point where the radius starts to increase again, cannot be constrained from the available data. Because of these uncertainties originating from the possible phase transition, model C shows a much larger scatter in the predicted radii at small masses. The radius is constrained between 10.5−12.5 km for M = 1.4 M_⊙ (2σ confidence limits).

4.4. Individual mass and radius results for the NSs

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 also puts significant constraints on the mass and radius of each object. The resulting M and R constraints for each object are given in Fig. 7 and summarized in Table 5. In general, the QMC + model A EoS tends to favor slightly larger masses and larger radii than compared to the QMC + model C. With the polytropic model A, the resulting mass is tightly constrained around M ≈ 1.5 M_⊙ for 4U 1724−307 and for SAX J1810.8−2609. Slightly larger mass of around M ≈ 1.8 M_⊙ is obtained for the remaining 4U 1702−429. Resulting radii are constrained to be around R ≈ 12.0 km for each source. With the piecewise linear model C the mass is about M ≈ 1.3 M_⊙ for the 4U 1724−307 and SAX J1810.8−2609 and, again, around M ≈ 1.8 M_⊙ for 4U 1702−429. The obtained radii are located around R ≈ 11.4 km. Because of the uncertainties from the possible phase transition occurring in model C, the resulting mass and radius constraints for each source are also much more loose.

5. Discussion

In this paper we have used the cooling tail method to constrain the mass and radius of three NS X-ray bursters: 4U 1702−429, 4U 1724−307, and SAX J1810.8−2609. Special care was taken to use only the passively cooling bursts as theoretical calculations of the color-correction factor f_c affecting the emerging spectra do not take any external heating into account. In practice this means that the blackbody normalization K is required to evolve during the burst (because K^{− 1 / 4} ∝ f_c; see, e.g., Poutanen et al. 2014; Kajava et al. 2014) and is indeed what we observed for the bursts in our sample.

First we introduced a new Bayesian cooling tail method and assumed uniform M vs. R priors in our analysis (instead of uniform F_Edd and A priors). By marginalizing over the M, R, and X prior distributions we also got distance estimates for our sources. One advantage here is that measurements like this are basically done in the X-ray band where interstellar extinction does not play such an important role, hence reducing possible model dependencies originating, for example, from selection of the interstellar extinction model. One should, however, note that in our case completely different kinds of model dependencies are present, related, for example, to the uncertain composition of the accreted material (i.e., the value of the hydrogen mass fraction) or to the X-ray burst selection. Unfortunately, distance measurements that we can compare against are available only for 4U 1724−307, as it is located in the globular cluster Terzan 2. Distance estimates to this source range from 5.3 to 7.7 kpc (Ortolani et al. 1997, using an extinction model valid for red stars or an average value from some large sample, respectively) in addition to the more recent measurement of 7.4 ± 0.5 kpc using near-IR observations of red giant branch stars (Valenti et al. 2010). We note that our distance constraints are consistent with the lower-limit end of the measurements as D_max ≈ 6 kpc. This value is, however, dependent on our selection of hydrogen mass fraction and should not be interpreted as a strict limit. If we would decrease the hydrogen mass fraction (and hence increase the D_max value) our resulting radii would also rapidly increase, creating tension between our other measurements. Interestingly, if we would impose a cut at around 5.3 kpc into our distance prior, our M vs. R results would be tightly constrained around R = 12.0 ± 0.3 km and M = 1.5 ± 0.2 M_⊙ as the new distance prior would remove the low-mass solutions. For 4U 1702−429 and SAX J1810.8−2609 we constrained the distance to be around $\mathrm{5.}{\mathrm{6}}_{-0.9}^{\mathrm{+}\mathrm{0.6}}\mathrm{kpc}$ and $\mathrm{4.}{\mathrm{3}}_{-1.0}^{\mathrm{+}\mathrm{0.6}}\mathrm{kpc}$ (2σ confidence limits), respectively.

Mass and radius constraints of 4U 1724−307 and SAX J1810.8−2609 are found to be almost identical, with the radius confined between about 11−13 km. Both of these sources are also best modeled by a solar-like composition with almost zero metallicity (SolA001 model). The best-fit model for the third, remaining source 4U 1702−429 consists of pure helium. This implies a hydrogen-poor companion like a white dwarf, that in turn, implies a compact binary system in order for the accretion to proceed via Roche lobe overflow. Best-fit values for the radius of this source (with X = 0) give R ≈ 13 km at around M = 1.5 M_⊙, a slightly larger value compared to the two other sources.

Some physical uncertainties are also still present in the X-ray burst M − R measurements. For example, no rotation is taken into account in our current work. Rotation affects the emerging spectrum because of the various special relativistic effects (see, e.g., Bauböck et al. 2015) but most importantly because the radius of the star increases at the equator and decreases at the poles as the star becomes oblate (Morsink et al. 2007; Bauböck et al. 2013) It, however, also introduces two new unknown free parameters to the fits: the spin period and the inclination. In many cases these parameters are not known a priori (especially the inclination). In addition, the flux distribution over the surface of the NS is unknown. The rotation does not, however, have a considerable effect on the M and R constraints when the spin frequency is moderate (≲400 Hz). Luckily, most of the X-ray bursters detected seem to have a frequency around this range. In our case, only 4U 1702−429 has a measured spin period of 329 Hz, meaning that it is a slow rotator and no major corrections are expected. For the two other sources, no burst oscillations nor persistent millisecond pulsations were detected and so the spin period is unknown. Other uncertainties may arise from the composition that can have an impact on the color-correction factors (Nättilä et al. 2015). Stratification may create an almost pure hydrogen layer or, in the opposite situation, strong convection might mix up the whole photosphere thoroughly, including the burning ashes consisting of heavier elements (see e.g., Weinberg et al. 2006; Malone et al. 2011, 2014). There are, however, convincing implications that these do not affect our current measurements since the observations of the normalization K do really seem to follow the theoretical atmosphere models with the simple solar-like (or pure He) compositions very closely. Additional confirmation is also obtained from the corrected temperatures that follow the $\mathit{F}\mathrm{\propto }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}$ law implying that the values of f_c used are correct. In the end, the distance is still our biggest source of uncertainty in the measurements.

Similarly to the physical uncertainties, some technical sources of systematic errors are present. For example data selection for individual bursts plays a role: fitting bursts from the touchdown down to only half of the touchdown flux (instead of the 20% used here) increases the uncertainties, as one would expect, because we use fewer data, but also increases the radius by as much as about 800 m. Also by refining the cooling tail method (where f_c − ℓ is used in the fitting procedures) into a more accurate two parameter treatment (where our model is ${\mathit{w}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{-}\mathrm{\left(}\mathit{w}{\mathit{f}}_{\mathrm{c}}^{\mathrm{4}}\mathrm{\right)}\mathit{\ell }$ and w is the dilution factor that was previously assumed to be $\mathit{w}\mathrm{\approx }{\mathit{f}}_{\mathrm{c}}^{-4}$) the radius is seen to increase by about 500 m (Suleimanov et al., in prep.). Because both of these effects act to increase our radius we can understand our current measurements as a lower limit of the real radius.

The results presented here constitute the first observational NS M − R constraints for 4U 1702−29 and SAX J1810.8−2609 using RXTE/PCA data. In addition, we constrained the compactness for the NS in 4U 1724−307 that has already been previously analyzed by Suleimanov et al. (2011a) and Özel et al. (2016). These three measurements were then used to create our parameterized EoS of cold dense matter.

In general, our EoS M − R constraints are slightly different from those found by Özel et al. (2016) as our radii tends to be somewhat larger, between about 10.6−12.4 km for model C and 11.2−12.7 km for model A, compared to their measurement of 10.1−11.1 km for M = 1.5 M_⊙. Note also that the parameterization of the EoS in Özel et al. (2016) is closer to our QMC + model A polytropic formalism. One possible cause for the difference might be data selection: the PRE burst we analyzed in this paper occurred in the hard spectral state with a low persistent flux level, whereas most of the bursts analyzed in Özel et al. (2016) occur in the soft spectral state and at higher persistent flux. As mentioned in the introduction, hard state X-ray bursts – such as the ones analyzed in this paper – tend to follow the NS atmosphere model predictions, whereas the soft state bursts never follow them (see Fig. 3 and also Suleimanov et al. 2011a; Poutanen et al. 2014; Kajava et al. 2014, for discussion). We therefore argue that in the soft state bursts there is an additional physical process (i.e. the spreading boundary layer) acting on the burst emission, that causes the assumptions on the visibility of the entire NS to break down and the value of the color-correction factor to be different from what is predicted by the passively cooling atmosphere models. One should also notice the completely different shape of the M − R contours obtained in Özel et al. (2016) where all of the M − R points are close to the region where only one mass-radius solution exists (see, e.g., Poutanen et al. 2014, for more in depth discussion about this), indicating that the lower-limit for the distance is close to the maximum distance D_max. This also creates tension for the solutions to exists on higher masses as the one-solution point lies on the R = 4GM/c² line (see, e.g., Özel & Psaltis 2015). In our case, no such a tension exists, leading, in turn, to much bigger M − R errors.

Specifically for the 4U 1724−307 we obtain R ≈ 12.0 ± 0.5 km whereas Özel et al. (2016) obtains ~12.2 km. Here the difference is already well inside error-limits. Suleimanov et al. (2011a,b) have also analyzed 4U 1724−307 using a different hard-state burst than what is in our sample and the resulting radii are considerably different (R> 14 km). There is, however, a possibility that the atmosphere during this burst might not consist of hydrogen and helium only, but is enriched with the nuclear burning ashes. These ashes would then affect the measured color-correction factor considerably (see Nättilä et al. 2015, and Appendix A for more in-depth discussion about this burst). In the recent analysis of 4U 1608−52 by Poutanen et al. (2014) they also used only the hard-state bursts from the source to constrain the mass and the radius. The resulting radii were, again, somewhat larger at around R ≳ 13 km for M = 1.4 M_⊙. We note, however, that 4U 1608−52 is a fast rotator with an observed spin frequency of 620 Hz (largest observed for an X-ray burster; Muno et al. 2002). This means that for an exact and reliable M − R analysis, a rotation-modified cooling tail method should be used (Suleimanov et al., in prep. ⁷; see also Bauböck et al. 2015).

In the second part of the paper, we combined our separate cooling tail measurements to obtain a parameterized EoS of the cold dense matter inside NSs. Using a Bayesian framework, we studied the possible constraints we can set on the EoS parameters such as the symmetry energy and its density derivative ℒ. We also relaxed our EoS priors by studying two different high-density models: model A consisting of polytropes and model C with piecewise linear segments in the ϵ − P plane. Here, model C with linear segments supports strong phase transitions in the supranuclear densities and indeed our resulting EoS is seen to have a large transition around ϵ ~ 400 MeV fm^-3 to support both large radius for large-mass and smaller radius for small-mass stars. Model A, on the other hand, gives much tighter constraints as the pressure in the core evolves rather smoothly and hence the resulting radius is almost independent of the mass. This restriction of having a constant radius then asserts tight limits to the obtained M vs. R results.

Nevertheless, all of the obtained constraints seem to favor an EoS that produces a radius of 10.5−12.8 km (95% confidence limits) for a mass of 1.4 M_⊙. One should, however, keep in mind that in reality the aforementioned systematic errors might increase these limits slightly. It is also interesting to note that we do not necessarily need a far better precision for the M − R measurements as interesting constraints can already be obtained from the existing observations. Encouragingly, the astrophysical constraints also seem to agree with nuclear physics experiments, theoretical studies, and heavy-ion experiments of neutron matter (Lattimer & Lim 2013). As well as the EoS constraints of the super-dense matter we get constraints for individual mass and radius measurements of the single NSs. By assuming that all of the sources must lie on the same M vs. R curve we can set stronger constraints than what is possible with only the cooling tail method alone. Particularly our results from the QMC + model A show how the unknown distance uncertainties (i.e., elongated M vs. R contours along the constant Eddington temperature) can be reduced by combining different sources with this kind of joint EoS fit. Finally, the mass, for example, can be constrained to an accuracy of about 0.2−0.3 M_⊙ (1σ confidence level) from the original limits spanning from M = 0.8 to 2.2 M_⊙.

Future prospects include an extended study of uninformative (for example, Jeffrey’s) priors in M − R space, both for the EoS and cooling tail fit parameters. This should be done to ensure that we do not implicitly and unknowingly transfer information to the mass and radius posteriors that already have many uncertainties present due to the poor measurements. Another obvious prospect is the addition of all possible M − R measurements including, but not restricted to, other X-ray bursting sources, quiescent LMXBs, thermally emitting isolated NSs, neutron star seismology results, pulse profiles in X-ray pulsars, moment of inertia, and crust thickness measurements.

Acknowledgments

The authors would like to thank M. C. Miller and J. Lattimer for valuable discussions and comments. We also thank the anonymous referee for an extremely thorough inspection of the manuscript and helpful comments that have significantly improved the paper. This research was supported by the Väisälä Foundation and by the University of Turku Graduate School in Physical and Chemical Sciences (JN). V.S. was supported by the German Research Foundation (DFG) grant WE 1312/48-1. J.J.E.K. acknowledges support from the ESA research fellowship programme. J.P. acknowledges funding from the Academy of Finland grant 268740. J.N. and J.J.E.K. acknowledge support from the Faculty of the European Space Astronomy Centre (ESAC). This research was undertaken on Finnish Grid Infrastructure (FGI) resources.

References

Appendix A: Two hard-state bursts from 4U 1724–307

The cooling tail method has been previously applied to another hard-state burst (obsid: 10090-01-01-021) from 4U 1724−307 (Suleimanov et al. 2011a,b). In this case, a rather large radius in excess of 13 km was obtained. Since then, there has been a second hard-state burst from the same source that does not, however, follow the same track in the K^{− 1 / 4}–F plane (this new burst is the one we use in our analysis; see Fig. A.1). We also note that the burst data presented in Suleimanov et al. (2011b, 2012) was reduced using an older RXTE data reduction pipeline without deadtime correction of the detectors. Since 2010, the RXTE team also introduced a correction to the effective area of the instrument that changed the measured flux by up to 10%. With our current (and final) RXTE data reduction methods (see Sect. 2) and with the new color-correction factors from Suleimanov et al. (2012) the older burst results in a radius of around 15 km at 1.4 M_⊙ when the cooling tail method is applied with the solar composition model (SolA001). With pure hydrogen composition the radius is brought down to the ~14 km range. We also note that in Suleimanov et al. (2011b, 2012) the first five most luminous points near the touchdown are ignored in the fit. Ignoring these points, however, leads to smaller Eddington flux and, hence, an even larger radius of around 16−17 km at 1.4 M_⊙ with our current data and new analysis methods.

The old burst also has poor χ² values at the low luminosity tail (below about <0.5F_td) originating from an additional powerlaw-like distribution of high-energy photons that cannot be described by the blackbody law. This, in turn, might imply additional surface heating by infalling material from the hot accretion flow. The extra heating does not, however, explain the discrepancy between the two bursts at higher luminosities. At these higher luminosities, it is possible that the NS atmosphere does not consist of hydrogen and helium only, but is enriched with the nuclear burning ashes (Weinberg et al. 2006). The presence of these ashes can then have a substantial impact on the color-correction factor reducing it by as much as 40% (Nättilä et al. 2015). This is also supported by the exceptionally long duration of the burst (~100 s compared to ~30 s for the other bursts in our sample) that might drive up the convection. Furthermore, the temperature evolution of this burst also shows some deviations from the L ∝ T⁴ law (see right-hand panel of Fig. A.1) that might imply inconsistencies with pure solar composition. Because to these uncertainties, we choose to leave this burst out from our sample and only use the newer hard-state burst from the source that is consistent with passive cooling.

Fig. A.1 Left-hand panel: cooling tail in the F ∝ L/LEdd vs. K− 1 / 4 ∝ fc plane with 1σ errors presented by black crosses for a hard-state PRE burst from 4U 1724−307 analyzed by Suleimanov et al. (2011b). Best-fit theoretical atmosphere model is shown with the blue curve (SolA001). A PRE burst from the same source (used in the this paper) is also presented with gray crosses and gray curve. Right-hand panel: temperature evolution of the Suleimanov et al. (2011b) burst. Blackbody temperature Tbb is shown for the cooling tail with black crosses. Blue crosses show the color-corrected temperatures Teff(1 + z)-1. The straight inclined lines show a powerlaw with an index of 4 corresponding to the F ∝ T4 relation.

Appendix B: Uniform F_Edd and A priors

Instead of assuming uniform M vs. R priors we can assume uniform F_Edd and A priors (and compute M and R) as in our previous work (see, e.g., Poutanen et al. 2014). This kind of selection of priors, however, turns out to be problematic in many ways.

First of all, we cannot truly introduce log g dependency into to the models as we know neither M nor R beforehand. In this case, we are left to iteratively select the best-fit surface gravity from the three basic values (log g = 14.0, 14.3, and 14.6) instead of interpolating between the models, as in the current work.

Secondly, there are some issues with the transformation from F_Edd and A parameters⁸ to M and R since the determinant of the Jacobian is $\begin{array}{c}\det \end{array}\left\{\mathit{J}\left(\frac{{\mathit{F}}_{\mathrm{Edd}}\mathit{,A}}{\mathit{M,R}}\right)\right\}\mathrm{=}\frac{\mathit{c}\hspace{0.17em}\mathit{G}\mathrm{|}\mathrm{1}\mathrm{-}\frac{\mathrm{4}\mathit{GM}}{\mathit{R}{\mathit{c}}^{\mathrm{2}}}\mathrm{|}}{\mathrm{2}{\mathit{D}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\kappa }}_{\mathrm{es}}{\mathit{R}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{GM}}{\mathit{R}{\mathit{c}}^{\mathrm{2}}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(B.1)From here we can see that we end up ignoring the so-called critical line where R = 4GM/c², as was shown by Özel & Psaltis (2015). Hence, some information of M and R is already present in (M,R) space owing to our choice of uniform F_Edd and A priors. The effect from this is shown visibly as a splitting of our M vs. R posteriors into two separate regions (see Fig. B.1). For clarity, no scatter in X or error in f_c is taken into account in this case. The lack of (M,R) solutions near this critical line also ends up affecting our final joint-fit results as can be seen in Fig. B.2. Because of these aforementioned issues we choose uniform priors in M and R instead, as they are final parameters that we want to study. We note, however, that in the cases analyzed before by Suleimanov et al. (2011b) and Poutanen et al. (2014) the error is not critical because the solution on M − R plane is constrained to be far from the R = 4GM/c² line because of the distance constraints.

Fig. B.1 Mass-radius constraints for the sources from the hard state PRE bursts assuming uniform FEdd and A priors. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upper-left region is excluded by constraints from the causality and general relativistic requirements (Haensel et al. 2007; Lattimer & Prakash 2007).

Fig. B.2 Individual mass and radius constraints for the three neutron stars used in the analysis with uniform FEdd and A priors. Left-hand panel shows the projected mass and the middle panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence regions of the parameters. Right-hand panels shows the full 2D mass and radius probability distributions. Contours of 68% (dotted, black line) and 95% (solid, black line) confidence regions are also shown.

Appendix C: EoS Tables and parameter correlations

The results for the two possible parameterized EoSs are tabulated here for easier access. Tables C.1 and C.2 list the pressure vs. density relation for the QMC + model A and QMC + model C EoSs, respectively. Similarly, Tables C.3 and C.4 summarize the resulting mass vs. radius relation for the two aforementioned models. Additionally, we show all possible EoS parameter correlations in Figs. C.1 and C.2.

Fig. C.1 QMC + model A EoS parameter correlations against each other. Red coloring shows the probability density of model parameters obtained from our Bayesian analysis with the cooling tail data.

Fig. C.2 QMC + model C EoS parameter correlations against each other. Symbols are the same as in Fig. C.1.

All Tables

All Figures

	Fig. 1 Reduced χ2 distributions for the blackbody spectral fits consisting of the points used in the cooling tail analysis. The dashed curve shows the theoretical expected χ2 distribution. Both obtained and theoretical distribution are normalized so that the encapsulated area of the curves is unity.
In the text

Fig. 2 Bolometric flux, temperature and blackbody normalization evolution during the hard-state PRE bursts. The black line shows the estimated bolometric flux (left-hand y-axis) in units of 10-7 erg cm-2 s-1. The blue ribbon shows the 1σ limits of the blackbody normalization (outer right-hand y-axis) in (km / 10 kpc)2. Similarly, the dashed orange line shows the color-corrected angular size (R∞/D10)2 using the same axis. The red ribbon show the 1σ errors for blackbody color temperature (inner right-hand y-axis) in keV. Highlighted gray area marks the region of the cooling tail used in the fitting procedures.

In the text

Fig. 3 Left panel: combined cooling tail in the F ∝ L/LEdd vs. K− 1 / 4 ∝ fc plane with 1σ error limits presented by crosses. Gray crosses show the burst evolution before the touchdown. Best-fit theoretical atmosphere models are shown by blue (SolA001) or red (He) curves. Right panel: temperature evolution of the bursts. Blackbody color temperature Tc is shown for each cooling tail with black crosses. Red (He) or blue (SolA001) crosses show the color-corrected temperatures Teff(1 + z)-1. Dashed lines show a powerlaw with an index of 4 corresponding to the F ∝ T4 relation. Highlighted gray area marks the region of the cooling tail used in the fitting procedures.

In the text

	Fig. 4 Mass-radius constraints for the sources from the hard state PRE bursts. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upper-left region is excluded by constraints from the causality and general relativistic requirements (Haensel et al. 2007; Lattimer & Prakash 2007).
In the text

Fig. 5 Histograms for the posterior distributions of the high-density model A and model C EoS parameters as implied by the data. In the first panel, the inset shows a magnified view of the histogram near n1 ≈ 0.5. Red shading corresponds to the 68% and blue to the 95% confidence regions of the parameters. Limits of the figures correspond to the hard limits set to the parameter prior distributions (see Table 3).

In the text

	Fig. 6 Obtained EoS constraints in the M − R (left panel) and in the P − ϵ plane (right panel). Upper panels correspond to the QMC + model A and bottom panels to the QMC + model C EoS. Red color indicates the probability density and black lines show the 68% (dotted) and 95% (solid) confidence limit contours.
In the text

Fig. 7 Individual mass and radius constraints for the three neutron stars used in the analysis. The left-hand panel shows the projected mass and the middle panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence regions of these parameters. Right-hand panels show the full 2D mass and radius probability distributions. Contours of 68% (dotted, black line) and 95% (solid, black line) confidence regions are also shown.

In the text

Fig. A.1 Left-hand panel: cooling tail in the F ∝ L/LEdd vs. K− 1 / 4 ∝ fc plane with 1σ errors presented by black crosses for a hard-state PRE burst from 4U 1724−307 analyzed by Suleimanov et al. (2011b). Best-fit theoretical atmosphere model is shown with the blue curve (SolA001). A PRE burst from the same source (used in the this paper) is also presented with gray crosses and gray curve. Right-hand panel: temperature evolution of the Suleimanov et al. (2011b) burst. Blackbody temperature Tbb is shown for the cooling tail with black crosses. Blue crosses show the color-corrected temperatures Teff(1 + z)-1. The straight inclined lines show a powerlaw with an index of 4 corresponding to the F ∝ T4 relation.

In the text

	Fig. B.1 Mass-radius constraints for the sources from the hard state PRE bursts assuming uniform FEdd and A priors. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upper-left region is excluded by constraints from the causality and general relativistic requirements (Haensel et al. 2007; Lattimer & Prakash 2007).
In the text

Fig. B.2 Individual mass and radius constraints for the three neutron stars used in the analysis with uniform FEdd and A priors. Left-hand panel shows the projected mass and the middle panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence regions of the parameters. Right-hand panels shows the full 2D mass and radius probability distributions. Contours of 68% (dotted, black line) and 95% (solid, black line) confidence regions are also shown.

In the text

	Fig. C.1 QMC + model A EoS parameter correlations against each other. Red coloring shows the probability density of model parameters obtained from our Bayesian analysis with the cooling tail data.
In the text

	Fig. C.2 QMC + model C EoS parameter correlations against each other. Symbols are the same as in Fig. C.1.
In the text