Issue 
A&A
Volume 591, July 2016



Article Number  A25  
Number of page(s)  23  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201527416  
Published online  06 June 2016 
Equation of state constraints for the cold dense matter inside neutron stars using the cooling tail method
^{1}
Tuorla ObservatoryDepartment of Physics and Astronomy, University of
Turku,
Väisäläntie 20,
21500
Piikkiö,
Finland
email:
joonas.a.nattila@utu.fi
^{2}
Department of Physics and Astronomy, University of
Tennessee, Knoxville,
TN
37996,
USA
^{3}
European Space Astronomy Centre (ESA/ESAC), Science Operations
Department, 28691 Villanueva de la
Cañada, Madrid,
Spain
^{4}
Institut für Astronomie und Astrophysik, Kepler Centre for Astro
and Particle Physics, Universität Tübingen, Sand 1, 72076
Tübingen,
Germany
^{5}
Astronomy Department, Kazan (Volga region) Federal
University, Kremlyovskaya str.
18, 420008
Kazan,
Russia
^{6}
Nordita, KTH Royal Institute of Technology and Stockholm
University, Roslagstullsbacken
23, 10691
Stockholm,
Sweden
Received: 22 September 2015
Accepted: 15 March 2016
The cooling phase of thermonuclear (typeI) Xray bursts can be used to constrain neutron star (NS) compactness by comparing the observed cooling tracks of bursts to accurate theoretical atmosphere model calculations. By applying the socalled cooling tail method, where the information from the whole cooling track is used, we constrain the mass, radius, and distance for three different NSs in lowmass Xray binaries 4U 1702−429, 4U 1724−307, and SAX J1810.8−260. Care is taken to use only the hard state bursts where it is thought that the NS surface alone is emitting. We then use a Markov chain Monte Carlo algorithm within a Bayesian framework to obtain a parameterized equation of state (EoS) of cold dense matter from our initial mass and radius constraints. This allows us to set limits on various nuclear parameters and to constrain an empirical pressuredensity relationship for the dense matter. Our predicted EoS results in NS a radius between 10.5−12.8 km (95% confidence limits) for a mass of 1.4 M_{⊙}, depending slightly on the assumed composition. Because of systematic errors and uncertainty in the composition, these results should be interpreted as lower limits for the radius.
Key words: dense matter / stars: neutron / Xrays: binaries / Xrays: bursts
© 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 manybody calculations have constrained the pressuredensity 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 TolmanVolkoffOppenheimer (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 heavyion collisions inside earthly laboratories (see, e.g., Lattimer 2012). One of the most promising candidates for obtaining accurate astrophysical massradius (M − R) measurements has been the thermal emission originating from NS surface layers. One possibility is to use the cooling of NS surface during typeI Xray bursts from lowmass Xray 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 latetype 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 Xray 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 infalling 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 nonexistent 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 socalled 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 colorcorrection factor versus the luminosity (in units of the Eddington), f_{c} − L/L_{Edd}, that is the socalled colorcorrection 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 colorcorrection 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 hardstate burst selection criteria for many sources simultaneously. We therefore now pay special attention to the burst selection and choose only the most wellbehaved hardstate 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 pressuredensity and massradius 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 neutronrich 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
Xray bursts used in the MR analysis.
In our analysis we used the data from the Proportional Counter Array (PCA; Jahoda et al. 2006) instrument on board of the Rossi Xray 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 Xray bursts using a similar method to that in Galloway et al. (2008). The timeresolved spectra for the bursts were extracted using an initial integration time of 0.25 s. In order to maintain approximately the same signaltonoise ratio the integration time was doubled every time the count rate decreased by a factor of . The exposures were deadtime corrected following the approach recommended by the instrument team^{1}. 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 deadtimecorrected 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 bestfit parameters are the color temperature T_{c} and the normalization constant K ≡ (R_{bb} [ km ] /D_{10})^{2}, where D_{10} = D/ 10 kpc. From the corresponding χ^{2} 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 χ^{2} distribution to be close to ideal. The bolometric flux was estimated using the cfluxmodel 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 doublepeaked 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 colorcorrection factor f_{c} as we approximate the emerging spectrum with a diluted blackbody model as (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_{c}T_{eff}(1 + z)^{1} , where 1 + z = (1−2GM/Rc^{2})^{−1/2} is the redshift. We stress here that the decrease in the colorcorrection 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 (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 softstate bursts the normalization is nearly constant, contrary to the theory^{2}. It is also crucial to notice that the normalization value in the tail of the softstate bursts is different to that observed for the hardstate bursts. In addition, the touchdown flux can also vary^{3} (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. hardstate 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 (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 ), we end up with a simple PRE condition (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 colorcorrection 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 hardstate PRE bursts. The black line shows the estimated bolometric flux (lefthand yaxis) in units of 10^{7} erg cm^{2} s^{1}. The blue ribbon shows the 1σ limits of the blackbody normalization (outer righthand yaxis) in (km / 10 kpc)^{2}. Similarly, the dashed orange line shows the colorcorrected angular size (R_{∞}/D_{10})^{2} using the same axis. The red ribbon show the 1σ errors for blackbody color temperature (inner righthand yaxis) in keV. Highlighted gray area marks the region of the cooling tail used in the fitting procedures. 
Timeresolved spectral parameters for the bursts in our sample are presented in Fig. 2. Additionally, we show the colorcorrected 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 colorcorrection 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) (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 nonoverlapping models ℳ_{i}, which exhaust the total model space ℳ. Then the relation can be written as (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 twodimensional prior distribution is assumed for (M,R) space. We also take into account the causality requirement, R> 2.824GM/c^{2} (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 colorcorrection curve fitting. The first one is the Eddington flux (7)where κ_{e} = 0.2(1 + X) cm^{2} g^{1}. The second parameter is related to the apparent (noncolorcorrected) 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 colorcorrection as K^{− 1 / 4} = f_{c}A. 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 KleinNishina reduction of the electron scattering opacity using an exact relativistic Comptonscattering kernel (Poutanen & Svensson 1996). These models give us the colorcorrection 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 subsolar 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., colorcorrection factors f_{c}) depend so weakly on this value (as our one sigma limits were 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 nonnegligible 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 hydrogenpoor 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 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 twodimensional 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 (8)Each separate probability is assumed to be proportional to the pathintegral evaluated through the twodimensional “data space”, (F,K^{− 1 / 4}), along the colorcorrection curve as (9)where f_{c,lo,hi} = (1 ± ϵ) × f_{c}(ℓ) are the lower and upper limits of the prior boxcar distribution of the colorcorrection f_{c} evaluated at relative luminosity ℓ and where ϵ = 0.03, ℱ_{c} = ℱ_{c}(F_{Edd},A) is the colorcorrection curve in (ℓ,f_{c}) space, 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 pathintegral is also areanormalized (or lengthnormalized 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 pathintegral formalism is related to the twodimensional 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 (10)where N_{P} = 4 (corresponding to M, R, D, and X) and (11)without the model priors that determine the integration limits. Then the onedimensional function represents the probability that the jth parameter will take a particular value given the observational data .
Results of the Bayesian cooling tail analysis.
Fig. 3 Left panel: combined cooling tail in the F ∝ L/L_{Edd} vs. K^{− 1 / 4} ∝ f_{c} plane with 1σ error limits presented by crosses. Gray crosses show the burst evolution before the touchdown. Bestfit theoretical atmosphere models are shown by blue (SolA001) or red (He) curves. Right panel: temperature evolution of the bursts. Blackbody color temperature T_{c} is shown for each cooling tail with black crosses. Red (He) or blue (SolA001) crosses show the colorcorrected temperatures T_{eff}(1 + z)^{1}. Dashed lines show a powerlaw with an index of 4 corresponding to the F ∝ T^{4} relation. Highlighted gray area marks the region of the cooling tail used in the fitting procedures. 
The bestfit 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 law, that is, it shows passive cooling.
Corresponding bestfit 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 twodimensional mass and radius probability posterior distributions of our analysis. The obtained contours are arched and elongated along the curves of constant Eddington temperature^{4}(12)where g is the surface gravity, σ_{SB} is the StefanBoltzmann 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 solarlike 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 Massradius 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 upperleft 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 lowerlimit 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^{2} 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 (13)On the other hand, our lowerlimit of the mass prior distribution (M_{min} = 0.8 M_{⊙}) also sets a lowerlimit 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 pressuredensity 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 (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 (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_{2}scl library (Steiner 2014b). To solve the integral (10)the code uses the MetropolisHastings 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_{0}. We assume nuclear binding energy of E_{nuc}(n_{0}) = −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^{14} g cm^{3}. The nuclear symmetry energy (the difference between energy per baryon of neutron matter and that of the nuclear matter)^{5} is denoted as , where n_{B} is the baryon number density and . The pressure of neutronrich matter at the saturation density n_{0} 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^{6}. 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.
Prior limits for the EoS parameters.
Most probable values and 68% and 95% confidence limits for the EoS parameters.
Near the nuclear saturation density, below the core, we employ the GandolfiCarlsonReddy (GCR) quantum Monte Carlo model (Gandolfi et al. 2012) that takes into account the threebody 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 (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 twonucleon 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 threenucleon interactions. Furthermore, in our analysis we reparameterized b and β in terms of and ℒ. Near the nuclear saturation density n_{0} the symmetry energy of the neutron matter can be obtained from (16)as (17)where E_{nuc}(n_{0}) = −16 MeV is the previously mentioned nuclear binding energy at the saturation density. For the pressure at the saturation density we obtain (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_{0}, 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 highdensity EoS we introduce three continuous power laws defining the pressure as (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 ϵ_{1} and the first polytrope index n_{1}, the second transition energy ϵ_{2} and the second polytrope index n_{2}, and a third polytrope index n_{3} (see Table 3 for the hard limits). Of course, we also require that ϵ_{2}>ϵ_{1} in order to avoid doublecounting 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_{1} = P(ϵ = 400) − P(ϵ = 200), ΔP_{2} = P(ϵ = 600) − P(ϵ = 400), ΔP_{3} = P(ϵ = 1000) − P(ϵ = 600), and ΔP_{4} = 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.
massradius 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 highdensity 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 n_{1} ≈ 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 lowdensity EoS behavior (which is dominated by twobody 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 highdensity 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 highdensity parameters of the model A are presented in Fig. 5. The index of the first polytrope n_{1} peaks sharply around 0.5 corresponding to a polytropic exponent γ_{1} = 1 + 1 /n_{1} ~ 3. The large value implies a rather stiff EoS at supranuclear densities. This first polytrope, corresponding to the n_{1} index, is seen to continue all the way up to the first transition density at ϵ_{1} ≈ 700 ± 150 MeV fm^{3}, that is around four times the saturation energy density ϵ_{0}. In some realizations the transition occurrs as early as ϵ_{1} ≈ 200 MeV fm^{3} but these cases correspond to the first sharp peak seen at n_{2} ≈ 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_{2}. In the opposite case, where all three polytropes span the grid the second polytrope index has values around 1.5 (polytropic exponent γ_{2} ~ 1.7) with a long tail extending all the way up to around 8.0 (i.e., to the upperlimit of our prior). What this means is that the data can not constrain the highdensity behavior of the EoS very well. As n_{2}>n_{1}, it, however, indicates that some softening of the EoS is present at higher densities. The third polytropic index n_{3} is only weakly constrained to be ≳1 (peaking around 1.5) indicating either no phase transitions at all or additional softening as n_{3}>n_{2}>n_{1}.
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_{2} and n_{3} parameter upperlimits 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 ϵ_{1} = 400 MeV the difference is tightly constrained around ΔP_{1} = 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_{2} and ΔP_{3} (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_{2}, ΔP_{3}, and ΔP_{4} appears because we rule out EoSs where the speed of sound exceeds the value of the speed of light. The modest highvalue tail for ΔP_{4} 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 lefthand 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. Righthand 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 Xray burst data is shown in Fig. 6 in M − R and P − ϵ planes. Each panel in the figure displays an ensemble of onedimensional (1D) histograms over a fixed grid in one of the axes (note that this is not quite the same as a 2D histogram). The righthand 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 lefthand panels of Fig. 6 present our results for the predicted massradius 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 righthand 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 nonexistent) 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 massradius 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 massradius curve is then able to support highmass NSs with radius of about R ≈ 11.6 km and lowmass 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
Most probable values for masses and radii for NSs constrained to lie on one mass versus radius curve.
The combination of several neutron star mass and radius measurements with the assumption that all neutron stars must lie on the same massradius 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 Xray 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 colorcorrection 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 Xray 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 Xray 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 nearIR observations of red giant branch stars (Valenti et al. 2010). We note that our distance constraints are consistent with the lowerlimit 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 lowmass solutions. For 4U 1702−429 and SAX J1810.8−2609 we constrained the distance to be around and (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 solarlike composition with almost zero metallicity (SolA001 model). The bestfit model for the third, remaining source 4U 1702−429 consists of pure helium. This implies a hydrogenpoor companion like a white dwarf, that in turn, implies a compact binary system in order for the accretion to proceed via Roche lobe overflow. Bestfit 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 Xray 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 Xray 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 colorcorrection 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 solarlike (or pure He) compositions very closely. Additional confirmation is also obtained from the corrected temperatures that follow the 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 and w is the dilution factor that was previously assumed to be ) 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 Xray 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 colorcorrection 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 massradius solution exists (see, e.g., Poutanen et al. 2014, for more in depth discussion about this), indicating that the lowerlimit 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 onesolution point lies on the R = 4GM/c^{2} 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 errorlimits. Suleimanov et al. (2011a,b) have also analyzed 4U 1724−307 using a different hardstate 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 colorcorrection factor considerably (see Nättilä et al. 2015, and Appendix A for more indepth discussion about this burst). In the recent analysis of 4U 1608−52 by Poutanen et al. (2014) they also used only the hardstate 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 Xray burster; Muno et al. 2002). This means that for an exact and reliable M − R analysis, a rotationmodified cooling tail method should be used (Suleimanov et al., in prep. ^{7}; 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 highdensity 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 largemass and smaller radius for smallmass 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 heavyion experiments of neutron matter (Lattimer & Lim 2013). As well as the EoS constraints of the superdense 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 Xray bursting sources, quiescent LMXBs, thermally emitting isolated NSs, neutron star seismology results, pulse profiles in Xray pulsars, moment of inertia, and crust thickness measurements.
One possible interpretation is that in the soft state the accretion disk continues all the way down to the NS surface forming a spreading/boundary layer. A combination of emission from a partly visible NS and from the spreading/boundary layer itself can then create timeevolving spectra that appear to have almost constant colorcorrection factor (Suleimanov & Poutanen 2006).
In the soft state the inner disk may act as a mirror reflecting part of the burst emission, therefore boosting the observed flux (Lapidus & Sunyaev 1985).
Eddington temperature formulated using the F_{Edd} and A parameters is not strictly constant in the M − R plane because F_{Edd} has a log g dependency (because of the dependence of f_{c} on log g). This complication is introduced through the new models (Suleimanov et al. 2012) that formally have superEddington luminosities due to the KleinNishina reduction of the effective crosssection. We note that our new cooling tail formalism allows us to take this into account as we use M and R as our parameters (instead of F_{Edd} and A).
Quartic terms are ignored, see Steiner (2006).
More accurately speaking, the constraints originate from nuclear masses (Kortelainen et al. 2010), quantum Monte Carlo model (Gandolfi et al. 2012), chiral interactions (Tews et al. 2013), and from isobaric analog states (Danielewicz & Lee 2014).
Our definition of the second A parameter differs from the one used in the work by Özel & Psaltis (2015) where . Conclusions, however, remain the same.
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/481. 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
 Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 [Google Scholar]
 Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V. eds. G. H. Jacoby, & J. Barnes (San Francisco: ASP), ASP Conf. Ser. 101 [Google Scholar]
 Bauböck, M., Berti, E., Psaltis, D., & Özel, F. 2013, ApJ, 777, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Bauböck, M., Özel, F., Psaltis, D., & Morsink, S. M. 2015, ApJ, 799, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Damen, E., Magnier, E., Lewin, W. H. G., et al. 1990, A&A, 237, 103 [NASA ADS] [Google Scholar]
 Danielewicz, P., & Lee, J. 2014, Nucl. Phys. A, 922, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Gandolfi, S., Carlson, J., & Reddy, S. 2012, Phys. Rev. C, 85, 032801 [NASA ADS] [CrossRef] [Google Scholar]
 Grinstead, C., & Snell, J. 1997, Introduction to Probability, 2nd edn. (Providence, RI: American Mathematical Society) [Google Scholar]
 Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars 1: Equation of State and Structure (New York: Springer), Astrophys. Space Sci. Libr., 326 [Google Scholar]
 Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Kajava, J. J. E., Nättilä, J., Latvala, O.M., et al. 2014, MNRAS, 445, 4218 [NASA ADS] [CrossRef] [Google Scholar]
 Klüpfel, P., Reinhard, P.G., Bürvenich, T. J., & Maruhn, J. A. 2009, Phys. Rev. C, 79, 034310 [NASA ADS] [CrossRef] [Google Scholar]
 Kortelainen, M., Lesinski, T., Moré, J., et al. 2010, Phys. Rev. C, 82, 024313 [Google Scholar]
 Kuulkers, E., Homan, J., van der Klis, M., Lewin, W. H. G., & Méndez, M. 2002, A&A, 382, 947 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuulkers, E., den Hartog, P. R., in’t Zand, J. J. M., et al. 2003, A&A, 399, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lapidus, I. I., & Sunyaev, R. A. 1985, MNRAS, 217, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Lapidus, I. I., Syunyaev, R. A., & Titarchuk, L. G. 1986, Sov. Astron. Lett., 12, 383 [NASA ADS] [Google Scholar]
 Lattimer, J. M. 2012, Ann. Rev. Nucl. Part. Sci., 62, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Lattimer, J. M., & Lim, Y. 2013, ApJ, 771, 51 [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Lattimer, J. M., & Steiner, A. W. 2014, Eur. Phys. J. A, 50, 40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, Space Sci. Rev., 62, 223 [NASA ADS] [CrossRef] [Google Scholar]
 London, R. A., Taam, R. E., & Howard, W. M. 1986, ApJ, 306, 170 [Google Scholar]
 Malone, C. M., Nonaka, A., Almgren, A. S., Bell, J. B., & Zingale, M. 2011, ApJ, 728, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Malone, C. M., Zingale, M., Nonaka, A., Almgren, A. S., & Bell, J. B. 2014, 788, 115 [Google Scholar]
 Miller, M. C. 2013, ArXiv eprints [arXiv:1312.0029] [Google Scholar]
 Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244 [NASA ADS] [CrossRef] [Google Scholar]
 Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ, 580, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Natalucci, L., Bazzano, A., Cocchi, M., et al. 2000, ApJ, 536, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Nättilä, J., Suleimanov, V. F., Kajava, J. J. E., & Poutanen, J. 2015, A&A, 581, A83 [NASA ADS] [EDP Sciences] [Google Scholar]
 Oppenheimer, J. R., & Volkoff, G. M. 1939, Phys. Rev., 55, 374 [NASA ADS] [CrossRef] [Google Scholar]
 Ortolani, S., Bica, E., & Barbuy, B. 1997, A&A, 326, 614 [NASA ADS] [Google Scholar]
 Özel, F. 2013, Rep. Prog. Phys., 76, 1 [Google Scholar]
 Özel, F., & Psaltis, D. 2015, ApJ, 810, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlov, G. G., Shibanov, I. A., & Zavlin, V. E. 1991, MNRAS, 253, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Penninx, W., Damen, E., van Paradijs, J., Tan, J., & Lewin, W. H. G. 1989, A&A, 208, 146 [NASA ADS] [Google Scholar]
 Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., Nättilä, J., Kajava, J. J. E., et al. 2014, MNRAS, 442, 3777 [NASA ADS] [CrossRef] [Google Scholar]
 Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Phys. Rev. D, 79, 124032 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W. 2006, Phys. Rev. C, 74, 045808 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W. 2014a, bamr: Bayesian analysis of mass and radius observations, Astrophysics Source Code Library [record record ascl:1408.020] [Google Scholar]
 Steiner, A. W. 2014b, O2scl: Objectoriented scientific computing library, Astrophysics Source Code Library [record record ascl:1408.019] [Google Scholar]
 Steiner, A. W., & Gandolfi, S. 2012, Phys. Rev. Lett., 108, 081102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, ApJ, 765, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W., Gandolfi, S., Fattoyev, F. J., & Newton, W. G. 2015, Phys. Rev. C, 91, 015804 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., Revnivtsev, M., & Werner, K. 2011a, ApJ, 742, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2011b, A&A, 527, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., Klochkov, D, & Werner, K. 2016, EPJ A, 52, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Tews, I., Krüger, T., Hebeler, K., & Schwenk, A. 2013, Phys. Rev. Lett., 110, 032504 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tolman, R. C. 1939, Phys. Rev., 55, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Valenti, E., Ferraro, F. R., & Origlia, L. 2010, MNRAS, 402, 1729 [NASA ADS] [CrossRef] [Google Scholar]
 van Paradijs, J., Dotani, T., Tanaka, Y., & Tsuru, T. 1990, PASJ, 42, 633 [NASA ADS] [Google Scholar]
 Weinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018 [NASA ADS] [CrossRef] [Google Scholar]
 Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Zamfir, M., Cumming, A., & Galloway, D. K. 2012, ApJ, 749, 69 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Two hardstate bursts from 4U 1724–307
The cooling tail method has been previously applied to another hardstate burst (obsid: 100900101021) 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 hardstate 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 colorcorrection 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 χ^{2} values at the low luminosity tail (below about <0.5F_{td}) originating from an additional powerlawlike distribution of highenergy 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 colorcorrection 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^{4} law (see righthand 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 hardstate burst from the source that is consistent with passive cooling.
Fig. A.1 Lefthand panel: cooling tail in the F ∝ L/L_{Edd} vs. K^{− 1 / 4} ∝ f_{c} plane with 1σ errors presented by black crosses for a hardstate PRE burst from 4U 1724−307 analyzed by Suleimanov et al. (2011b). Bestfit 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. Righthand panel: temperature evolution of the Suleimanov et al. (2011b) burst. Blackbody temperature T_{bb} is shown for the cooling tail with black crosses. Blue crosses show the colorcorrected temperatures T_{eff}(1 + z)^{1}. The straight inclined lines show a powerlaw with an index of 4 corresponding to the F ∝ T^{4} 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 bestfit 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^{8} to M and R since the determinant of the Jacobian is (B.1)From here we can see that we end up ignoring the socalled critical line where R = 4GM/c^{2}, 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 jointfit 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^{2} line because of the distance constraints.
Fig. B.1 Massradius constraints for the sources from the hard state PRE bursts assuming uniform F_{Edd} and A priors. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upperleft 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 F_{Edd} and A priors. Lefthand 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. Righthand 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.
Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + model A.
Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + model C.
Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + model A.
Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + model C.
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
Most probable values for masses and radii for NSs constrained to lie on one mass versus radius curve.
Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + model A.
Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + model C.
Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + model A.
Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + model C.
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 hardstate PRE bursts. The black line shows the estimated bolometric flux (lefthand yaxis) in units of 10^{7} erg cm^{2} s^{1}. The blue ribbon shows the 1σ limits of the blackbody normalization (outer righthand yaxis) in (km / 10 kpc)^{2}. Similarly, the dashed orange line shows the colorcorrected angular size (R_{∞}/D_{10})^{2} using the same axis. The red ribbon show the 1σ errors for blackbody color temperature (inner righthand yaxis) 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/L_{Edd} vs. K^{− 1 / 4} ∝ f_{c} plane with 1σ error limits presented by crosses. Gray crosses show the burst evolution before the touchdown. Bestfit theoretical atmosphere models are shown by blue (SolA001) or red (He) curves. Right panel: temperature evolution of the bursts. Blackbody color temperature T_{c} is shown for each cooling tail with black crosses. Red (He) or blue (SolA001) crosses show the colorcorrected temperatures T_{eff}(1 + z)^{1}. Dashed lines show a powerlaw with an index of 4 corresponding to the F ∝ T^{4} relation. Highlighted gray area marks the region of the cooling tail used in the fitting procedures. 

In the text 
Fig. 4 Massradius 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 upperleft 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 highdensity 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 n_{1} ≈ 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 lefthand 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. Righthand 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 Lefthand panel: cooling tail in the F ∝ L/L_{Edd} vs. K^{− 1 / 4} ∝ f_{c} plane with 1σ errors presented by black crosses for a hardstate PRE burst from 4U 1724−307 analyzed by Suleimanov et al. (2011b). Bestfit 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. Righthand panel: temperature evolution of the Suleimanov et al. (2011b) burst. Blackbody temperature T_{bb} is shown for the cooling tail with black crosses. Blue crosses show the colorcorrected temperatures T_{eff}(1 + z)^{1}. The straight inclined lines show a powerlaw with an index of 4 corresponding to the F ∝ T^{4} relation. 

In the text 
Fig. B.1 Massradius constraints for the sources from the hard state PRE bursts assuming uniform F_{Edd} and A priors. Constraints are shown by 68% (dotted line) and 95% (solid line) confidence level contours. The upperleft 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 F_{Edd} and A priors. Lefthand 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. Righthand 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.