Issue 
A&A
Volume 608, December 2017



Article Number  A31  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201731082  
Published online  01 December 2017 
Neutron star mass and radius measurements from atmospheric model fits to Xray burst cooling tail spectra
^{1} Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland
email: joonas.a.nattila@utu.fi
^{2} Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{3} Department of Astronomy and Joint SpaceScience Institute, University of Maryland, College Park, MD 207422421, USA
^{4} Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA
^{5} Finnish Centre for Astronomy with ESO (FINCA), University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland
^{6} European Space Astronomy Centre (ESA/ESAC), Science Operations Department, 28691 Villanueva de la Cañada, Madrid, Spain
^{7} Institut für Astronomie und Astrophysik, Kepler Centre for Astro and Particle Physics, Universität Tübingen, Sand 1, 72076 Tübingen, Germany
^{8} Astronomy Department, Kazan (Volga region) Federal University, Kremlyovskaya str. 18, 420008 Kazan, Russia
^{9} Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA
Received: 2 May 2017
Accepted: 26 September 2017
Observations of thermonuclear Xray bursts from accreting neutron stars (NSs) in lowmass Xray binary systems can be used to constrain NS masses and radii. Most previous work of this type has set these constraints using Planck function fits as a proxy: the models and the data are both fit with diluted blackbody functions to yield normalizations and temperatures that are then compared with each other. For the first time, we here fit atmosphere models of Xray bursting NSs directly to the observed spectra. We present a hierarchical Bayesian fitting framework that uses current Xray bursting NS atmosphere models with realistic opacities and relativistic exact Compton scattering kernels as a model for the surface emission. We test our approach against synthetic data and find that for data that are well described by our model, we can obtain robust radius, mass, distance, and composition measurements. We then apply our technique to Rossi Xray Timing Explorer observations of five hardstate Xray bursts from 4U 1702−429. Our joint fit to all five bursts shows that the theoretical atmosphere models describe the data well, but there are still some unmodeled features in the spectrum corresponding to a relative error of 1–5% of the energy flux. After marginalizing over this intrinsic scatter, we find that at 68% credibility, the circumferential radius of the NS in 4U 1702−429 is R = 12.4±0.4 km, the gravitational mass is M = 1.9±0.3 M_{⊙}, the distance is 5.1 < D/ kpc < 6.2, and the hydrogen mass fraction is X < 0.09.
Key words: dense matter / stars: neutron / Xrays: binaries / Xrays: bursts
© ESO, 2017
1. Introduction
The masses and radii of neutron stars (NSs) encode valuable information about the properties of the matter in their cores (Lattimer 2012; Lattimer & Steiner 2014), which reaches several times nuclear saturation density and has strong isospin asymmetry, and which therefore cannot be analyzed in terrestrial laboratories. Hence, detailed measurements of NS masses and radii are invaluable in the study of cold dense matter, and in particular, in the equation of state (EoS) of the matter, that is, the relation between thermodynamic quantities such as the pressure and the energy density (see Miller 2013; Özel 2013; Miller & Lamb 2016 for recent discussions). The reliability of such constraints depends on the degree to which systematic errors can be controlled (in many current analyses, these errors are significantly larger than the formal statistical uncertainties; see Miller 2013; Miller & Lamb 2016), as well as on the precision of the astrophysical models that are applied to the data.
One type of source that has attracted considerable attention in this context are lowmass Xray binaries (LMXBs) that exhibit frequent thermonuclear Xray bursts (for reviews, see Lewin et al. 1993; Strohmayer & Bildsten 2006). By collecting observations from these bursts and modeling how they cool down, we can set constraints on the size of the emitting area (for early work in this field, see, e.g., Ebisuzaki 1987; Damen et al. 1990; van Paradijs et al. 1990; Lewin et al. 1993). It was, however, only the Rossi Xray Timing Explorer (RXTE) that was able to produce a large catalog of observations to study (see, e.g., Galloway et al. 2008). Since then, a large number of bursts from different sources have been put to use (for recent reviews, see Miller & Lamb 2016; Suleimanov et al. 2016). Two principal methods are currently used to infer the gravitational mass M and the circumferential radius R from burst cooling tails, both of which stem from the earlier work: the touchdown method (e.g., Özel 2006; Özel et al. 2009, 2016; Güver et al. 2010), and the cooling tail method (e.g., Suleimanov et al. 2011a, 2017; Poutanen et al. 2014; Nättilä et al. 2016). Both methods fit the observed emission using the Planck function and then compare the evolution of the observed temperature and normalization to the model predictions (see, however, Kuśmierek et al. 2011, for an early attempt to circumvent the usage of Planck function fits alone). These fits significantly simplify model comparison because the observed spectra are relatively well described by thermal emission. However, because neither model atmosphere spectra (e.g., Suleimanov et al. 2012) nor the most accurately measured observed spectra (Miller et al. 2011) are exactly Planckian, using Planck fits as proxies throws away information and could even introduce biases.
Here we present for the first time simultaneous direct atmosphere model fits to a set of Xray burst observations. We begin by studying the constraints that can be obtained from synthetic data, for which our model is a good description. This allows us to assess the accuracy of our method and to explore possible biases in the results. We then apply our method to data from five hardstate bursts of 4U 1702−429. We obtain interesting constraints on the mass and radius of this star and also study some of the previously neglected physical assumptions present in the fitting procedures. The LMXB system 4U 1702−429 is a particularly good testbed for the fitting as the cooling tail has previously been modeled for the five hardstate bursts from this source (Nättilä et al. 2016). Our initial analysis suggests that direct fitting of detailed atmosphere models to data is a promising avenue to extract neutron star masses and radii from Xray burst data.
In Sect. 2 we present the theoretical basis for our analysis. This section is split into two parts: in Sect. 2.1 we model the emerging radiation and couple it to the actual observations, and in Sect. 2.2 we formulate the Bayesian framework and present our hierarchical fitting model. We apply our model to synthetic data in Sect. 3.1 and then to real Xray burst observations from 4U 1702−429 in Sect. 3.2. In Sect. 3.2 we also present our new improved mass, radius, distance, and composition constraints for the source. In Sect. 4 we discuss our results. Finally, we present our summary in Sect. 5.
2. Methods
2.1. Model for the emerging radiation
We assume that radiation from a point on the surface of an NS is initially emitted with a local specific intensity at energy E′, as measured in the local frame of the emission. Assuming that the radiation propagates through vacuum to a distant observer, this observer will detect this radiation at energy E, where the energies are related by (1)where z takes into account the rotationinduced Doppler shifts and the gravitational redshift. In the limit of low spin frequency, ν → 0, the external spacetime is Schwarzschild and there are no Doppler shifts, and therefore the net redshift approaches (2)where G is the gravitational constant, c is the speed of light, and M and R are the gravitational mass and the circumferential radius of the NS, respectively. A distant observer will measure a specific (monochromatic) intensity I_{E} that is related to the original specific intensity (Liouville’s theorem for photons, see, e.g., Misner et al. 1973; Rybicki & Lightman 1979) by (3)The total observed monochromatic flux from the star, as seen by a distant observer, is then (see, e.g., Nättilä & Pihajoki 2017) (4)where dΩ is the solid angle that the surface element occupies on the observer’s sky. In this paper we consider a uniformly emitting slowly rotating NS. In this limit, the observed flux is related in a simple way to the flux F′ emitted at the NS surface: (5)where R_{∞} = R(1 + z) is the apparent NS radius, D is the distance to the source, and , where μ is the cosine of the angle between the local normal direction and the direction of emission of radiation.
In general, a burster is not expected to emit uniformly, and rotation rates of known bursters extend up to 620 Hz (Muno et al. 2002; Watts 2012). Rotation introduces Doppler shifts that vary over the surface of the star and therefore smear sharp spectral features such as emission lines. These shifts also broaden the continuum spectra, but such broadening can usually be neglected (Nättilä & Pihajoki 2017). Moreover, the assumption of uniform emission combined with slow rotation means that the observed flux depends on the surface flux and distance, but not on the angular dependence of the spesific intensity. This is not true in more general situations. We employ these approximations because they allow us to simplify the general Eq. (4)and avoid the usage of computationally costly ray tracing to combine the flux from different parts of the star. They also allow us to neglect a few potentially important but often unknown parameters, such as the NS rotation frequency and the observer’s inclination angle and the unknown latitude dependence of the flux.
The gravitational acceleration at the NS surface is given by (6)The shape of the emerging spectrum has a weak dependence on g. The composition of the atmosphere also affects the spectrum via the energy dependence of the opacity κ, which includes contributions from true absorption and scattering. For example, for an atmosphere with a hydrogen mass fraction X, the Thomson scattering opacity is (7)When we assume Thomson opacity and a spherically symmetric flux, the outward radiative acceleration balances the inward gravitational acceleration at the stellar surface when the stellar luminosity reaches the Eddington luminosity L_{Edd}, which is defined by (8)The actual critical luminosity is reached when the radiative acceleration (9)equals the surface gravitational acceleration g. Here κ_{R} is the flux mean opacity (equal in our case to the Rosseland mean opacity), (10)is the bolometric surface flux, T_{eff} is the effective temperature of radiation, and σ_{SB} is the StefanBoltzmann constant. At high temperatures close to the critical luminosity, the opacity is dominated by Compton scattering and is lower than κ_{T} because of the KleinNishina effect (Poutanen 2017), resulting in a critical luminosity exceeding L_{Edd} by 5–10% (Suleimanov et al. 2012). We use the ratio g_{rad}/g to measure the escaping flux from the star.
Using Eq. (8), we can also introduce the bolometric Eddington flux (11)which is characterized by the corresponding Eddington temperature T_{Edd}. The corresponding observed Eddington flux then can be obtained by integrating Eq. (5) over energies: (12)where A = (R_{∞}/D)^{2} is related to the apparent angular size of the star and T_{Edd,∞} = T_{Edd}/ (1 + z) is the redshifted Eddington temperature.
For the flux escaping from the NS surface F′(E′), we used the spectra from the atmosphere models computed in Suleimanov et al. (2011b, 2012) and Nättilä et al. (2015). These calculations were based on the stellar modeling program atlas (Kurucz 1970, 1993), but were modified to deal with high temperatures (Ibragimov et al. 2003; Suleimanov & Poutanen 2006; Suleimanov & Werner 2007) and to take into account Compton scattering (Suleimanov et al. 2012) using an exact relativistic redistribution function (see, e.g., Poutanen & Svensson 1996). The models were computed in hydrostatic equilibrium, using local thermodynamic equilibrium, and assuming a planeparallel atmosphere structure. Because of these approximations, our model is only valid when the atmospheric scale height is much smaller than the stellar radius, which means that g_{rad}/g must be lower than and not too close to unity.
Here we limited our model spectra to the range g_{rad}/g = 0.2−0.98 to avoid any physical complications occurring at low or high temperatures: at high temperatures, the scale height is too large, and at low temperatures (relevant late in the burst tails), it is likely that ongoing accretion breaks the assumption that the observed radiation emerges only from the passively cooling neutron star surface. In practice, limiting g_{rad}/g to this range implies focusing on the first approximately10 s of the cooling tail. The compositions computed by Suleimanov et al. (2012) were pure hydrogen (X = 1), pure helium (X = 0, Y = 1), and solar hydrogentohelium ratio X = 0.738 with different metallicities (0.01, 0.1, 0.3, and 1.0 of solar). For simplicity, we here considered atmospheres with a metallicity that is 0.01 of solar. This selection is possible and does not introduce a considerable error because metals will only affect the spectra at the very late stage of the burst, when the atmosphere has a sufficiently low effective temperature (Suleimanov et al. 2011b, 2012). On the other hand, we did not consider observations on a stage so late that this effect would play a role. We also note that an exact selection of the metallicity does not play a key role here because we did not consider cold atmosphere models (we used g_{rad}/g > 0.2), for which photoionization edges start to dominate the spectral features. In addition, we considered surface gravities of log _{10}g = 14.0, 14.3, and 14.6 (with g in cgs units).
We used these models to obtain spectra with any given g_{rad}/g, log g, and X by linearly interpolating (or in the case of log g also linearly extrapolating) the logarithm of the monochromatic fluxes on the model photon energy grid. The model parameter limits were g_{rad}/g = 0.2−0.98, log g = 13.7−14.9, and X = 0−1. We checked the accuracy of our interpolations by comparing our results against actual model spectra that were computed between the original grid points, and found that the relative accuracy of the spectral energy flux was better than 1%.
Of course, what we observe is not energy flux but rather photon counts in energy channels. We therefore, converted our model spectra into photon counts by convolving them with a response function R(I,E) of a detector: (13)where M(E) is the photon number flux of the model at energy E and is the observing time. Here the response function R is proportional to the probability that an incoming photon of energy E will be detected in channel i and is a discrete function (i.e., a response matrix) such that (14)for an energy range E_{j−1} to E_{j}. In addition, we must take into account that the data might have a nonzero background. In this case, we fit the observed background with some spectral model F_{bkg}(E), so that our total model photon flux at energy E is(15)The background flux is often estimated by observations be fore or after the burst, but observational (Yu et al. 1999; van Paradijs & Lewin 1986; Kuulkers et al. 2003; Chen et al. 2011; in’t Zand et al. 2011; Serino et al. 2012; Degenaar et al. 2013; Worpel et al. 2013; Peille et al. 2014; Worpel et al. 2015; Degenaar et al. 2016; Koljonen et al. 2016; Kajava et al. 2017c) and theoretical work (Walker 1992; Miller & Lamb 1996; Ballantyne & Strohmayer 2004; Ballantyne & Everett 2005) suggests that the burst can increase or decrease the background rate and even change its spectrum. Thus background estimates from times near the burst are unreliable and use of them could introduce bias. This is one reason why we focused on bursts that occur during the hard spectral state: for such bursts, the persistent background emission is very weak (≲ 1% of the peak flux). Therefore, although in practice we estimate the background using a 16 s observation before the burst, which we find is described well by blackbody plus powerlaw components (bbodyrad and powerlaw in xspec), the background modeling is unimportant because the emission is dominated by the burst radiation. Results supporting this kind of static persistent emission treatment were also presented in Kajava et al. (2017c), where it was concluded that even if the background emission varies during the burst, it is unlikely to contribute more than 1% of the burst flux in the hard state. In the soft state, on the other hand, one of the persistent emission components can brighten more than tenfold during the bursts (Kajava et al. 2017a). Finally, we note that we multiplied both the background and the theoretical burst spectra by an interstellar absorption model (similar to phabs in xspec) to account for the nonzero neutral hydrogen column depth.
2.2. Hierarchical fitting model for M and R constraints
Next, we construct a framework for comparing the emission models to the actual observations of Xray bursts. In order to do this, we formulate a hierarchical model for an NS that has been observed to have N^{B} bursts B_{k}, where k = 1,...,N^{B}. We denote the set of all bursts as . Each of these bursts B_{k} has spectra, which we label as S_{jk}, one for each time bin j. The set of all spectra in B_{k} is similarly denoted as . Each spectrum S_{jk} consists of a set with measurements of counts C_{𝒟,ijk}, measured in the detector channel i.
For a single channel we can define the likelihood function as so that is the probability of the data given the model ℳ and a set of assumptions ℋ. When the source of experimental noise is due to the number of events arriving at the detector, the counting statistics are Poisson distributed. Hence, to estimate the model goodnessoffit for one element in some arbitrary burst B_{k} and spectrum S_{jk}, we compute the likelihood for a Poisson distributed data as (16)The joint likelihood for a single spectrum is then (17)where i ranges from the first detector channel i_{min} to the last detector channel i_{max} used in the analysis. We note that because in practice likelihoods can be extremely large or small, we instead used log likelihoods in our analysis. The joint likelihood for a burst is , and the total joint likelihood for all bursts is .
In the limit of high count rates, the Poisson distribution is well approximated by a Gaussian distribution. In this case, the likelihood is proportional to exp(−χ^{2}/ 2), where (18)This is known as Pearson’s weighting when the statistical error in the denominator is taken from the model. Similarly, we could describe the error with the help of the data counts, as is done with Neyman’s weighting. Both of them are biased (in the opposite directions) estimators of the model parameters (see Humphrey et al. 2009). Because χ^{2} is proportional to the log likelihood, the joint log likelihood for the spectra in a burst, and the joint log likelihood for all bursts, is the sum of the individual χ^{2} values. We caution that particularly in the higherenergy channels, there may not be enough counts for the Poisson distribution to be well approximated by a Gaussian. In this case, χ^{2} is not a good approximation to the log likelihood, and this could contribute to the formally poor χ^{2} we find in Sect. 3.2 when we analyze data from 4U 1702−429.
It is also possible that we have underestimated the uncertainties in our data. In this case it is typical, in a Bayesian analysis, to introduce an intrinsic scatter σ_{int} into the system. Physically the intrinsic scatter can be understood to originate either from the instrument calibration error or from the uncertainty in the actual model used, which as we recall is interpolated from tabulated points in the space of composition, surface gravity, and temperature. The addition of intrinsic scatter is similar to the error expansions in frequentist methods where data errors are increased until the total χ^{2}/ d.o.f. is around unity. The underlying idea is that intrinsic scatter acts to quantify and penalize our ignorance of the model: by increasing σ_{int}, the possible credible regions for other parameters also inflate to take into account that the data are not fully described by the model. Mathematically this is done by convolving the original Gaussian distribution, (where σ can be taken to be in relation to the χ^{2} formulation), with another normal distribution with undefined error σ_{int}. It is standard to assume, given no other knowledge, that the errors can be added in quadrature: . The likelihood in Eq. (16)can then be replaced with (19)where theasterisk marks the convolution of the two functions. Taking the logarithm simplifies the latter expression to (20)From this we see that the expression reduces to χ^{2} when σ_{int} → 0 and , as the first two terms in the likelihood expression can then be ignored as constants. It is important to note that σ_{int} cannot be increased infinitely to obtain a better likelihood because the ln term in the loglikelihood expression compensates for the last term. Hence, there exists some balance between the two terms and σ_{int} can only grow to some finite value where the previously unexplained scatter in the observations is then explained by the model.
Parameters of hierarchical fitting models.
The actual model spaces are then constructed hierarchically in addition to different clusters of nested data groups. Such a model with nested hierarchy can be physically motivated by considering our problem at hand: we have a neutron star that has some model parameters that can characterize it (such as size and distance). This neutron star will then exhibit bursts that could also have some model parameters (such as composition and ignition depth). The bursts, however, all share the parameters that the neutron star has and, hence, in the model parameter hierarchy the burst parameters appear lower. The bursts, on the other hand, are constructed of a series of snapshots in time that we call energy spectra. Again, one individual spectrum could have parameters dedicated only to that one particular spectrum or share some parameters among the other spectra in the burst. This nesting of parameters we then call a hierarchical nested model in this paper.
We here consider four hierarchical models, A–D, presented in Table 1. As an example, we next describe the models in more detail. At the top level of model A, we define three shared global parameters: NS mass M, radius R, and distance D to the star. The combination of M and R then gives us the surface gravity g and the redshift 1 + z. These can be combined with the distance to give the quantity A, which is proportional to the solid angle occupied by an NS on the sky. In addition to these basic parameters, we can set a global composition of the accreted matter via the hydrogen mass fraction X, as is done in Model B. The next level in the model involves the different bursts for which we do not introduce any clusterspecific subparameters in this work. Proceeding in the hierarchy tree, each spectrum of the burst has always at least one individual parameter to sample: the effective temperature of the emerging spectrum as expressed through the parameter g_{rad}/g. In Model C, we also introduce the fraction S_{f}, which emits (and we assume that the emitting portion emits uniformly). S_{f} is a free parameter that enters the flux equation by modifying the apparent angular size A′ = S_{f}A. As we show below, the real data are not fully described by the model. To accommodate this deviation, in the end, we expanded model B by introducing a free intrinsic scatter to the system on the global scale. We labeled this model D. In contrast to models A to C, for D we chose a nonuniform phenomenological distance prior; this choice was motivated by the synthetic data results. Intrinsic scatter was always, when present, sampled as , as it is a scale parameter in the model.
When a parameter is not free, but has some fixed value, our model is said to have a set of physical assumptions ℋ that implicitly enter the likelihood calculations. The strictest set of assumptions ℋ was imposed for model A, which assumes a constant (uniform) emitting area S_{f} = 1 and a known nonvarying chemical composition (X = 0 or 0.73 in this work). In Model B, the assumption of the chemical composition was relaxed. Similarly, in Model C, we tested the validity of the constant emitting area assumption.
On a purely theoretical basis, we would expect a model in which every parameter is defined as free in the lowest hierarchy level to be the least informative: by allowing the hydrogen fraction X and the surface fraction S_{f} both to evolve freely in time for each spectrum, we could check the assumption of constant uniform emitting area and constancy of the chemical composition. In practice, such freedom in the model is not possible, however, as X and S_{f} are strongly correlated because they both affect the normalization of the flux. While composition does have a slight effect on the shape of the spectral energy distribution, the current data do not allow any meaningful constraints without the additional normalization dependency. These freedoms could be slightly limited by making the composition vary only from burst to burst. This would then allow us to study the time evolution of the composition on much longer timescales from burst to burst. Another possibility would be to introduce a burstspecific S_{f} term into the model to allow variations between bursts, for example, by a changing inner radius of the accretion disk.
Finally, we sampled the parameter model space using Bayesian inference. We introduced uniform prior distributions for M and R in the range 1.0−2.2 M_{⊙} and 8−16 km, respectively. For the distance, a uniform prior was taken in between D = 2−10 kpc For model D, we chose to sample not D, but D^{3/2}, corresponding to a weakly informative prior of for the distance, which slightly favors higher values. This selection seems to remove the otherwise strong preference for lower masses^{1}. We discuss this selection further in Sect. 3.1. When the hydrogen mass fraction was not fixed, we assumed a flat prior ranging from 0 (pure helium) to 1 (pure hydrogen). For the spectrumspecific nuisance parameters, we took similarly uniform limits so that g_{rad}/g = 0.2−0.98 and S_{f} = 0.5−1.5. We note that values of S_{f} > 1 were also allowed, as it might be possible that the apparent emitting area exceeds the one inferred from the star’s angular size A because of reflection from the accretion disk (Lapidus & Sunyaev 1985). We also did not impose any time relation in any of our models between neighboring time bins via adjacent g_{rad}/g values: every spectrum was free to attain any value in the prior range regardless of the adjacent spectra. We then constrained the model parameters by sampling from the posterior distributions using Markov chain Monte Carlo (MCMC) sampling. To explore the parameter spaces, we implemented an affineinvariant ensemble sampler, as discussed in Goodman & Weare 2010 (see also Guillot et al. 2013). The actual implementation is heavily based on the bamr code (Steiner 2014a), which in turn relies on the o2scl library (Steiner 2014b). The ensemble sampler is similar to a normal MetropolisHastings algorithm, but evolves not one, but many parallel sample values, called walkers, together. The random step for each walker is then performed using a socalled stretch move algorithm, where each walker makes a small step in the parameter space in relation to the whole ensemble. Acceptance is still performed using a MetropolisHastings scheme. With correlated distributions, this will improve the autocorrelation times of the chains tremendously, allowing us to sample the parameter space more thoroughly in a shorter time.
3. Analysis
3.1. Synthetic data
We begin our analysis by applying the methods described above to synthetic data. The mock data were created to resemble the observations from NASA’s RXTE Proportional Counter Array (PCA; Jahoda et al. 2006). We produced data using R = 12 km, M = 1.5 M_{⊙}, D = 6.0 kpc, and X = 0. These values are similar to those inferred from the five bursts of 4U 1702−429 that we analyze later (see Nättilä et al. 2016). We also studied similar NS configurations, but with X = 0.73, which thus corresponds to solar composition, to determine how the composition affects the results. The mock observations were created by computing the actual model spectra using the atmosphere models described in Sect. 2.1. In this process, 20 spectra for each burst were created for five bursts in total so that g_{rad}/g was linearly spaced between values of 0.2 and 0.95. For each spectrum we convolved the model with an actual RXTE/PCA response matrix, computed the number of observed counts in each energy channel, and then drew the observed number of counts from a Poisson distribution centered around the real value. For the background spectra, we used the real background files from 4U 1702−429 hardstate bursts. We increased the exposure time from 0.25 s to 0.5 s after ten time bins; this procedure is commonly used for real data to keep the signaltonoise ratio level per time bin approximately constant despite the decreasing flux. We fixed the neutral hydrogen column depth to N_{H} = 1.87 × 10^{22} cm^{2}, which again is similar to that of 4U 1702−429 (Worpel et al. 2013). Figure 1 shows a set of spectra from one such synthetic burst.
Fig. 1 Synthetic spectra (crosses) with corresponding bestfit atmosphere models (solid lines) for A. Different colors show spectra for varying g_{rad}/g. Individual spectra are shifted by factors of 2 in the ydirection for clarity. 
Next we fit these synthetic data with different models to assess how well we expect to constrain each model parameter. This also gives us information about the possible biases in the method, given that we know the input values of the parameters. From here on, when we discuss credible regions, we always refer to the highest posterior density regions. First we study how well we can obtain information about nuisance parameters such as the normalized gravity g_{rad}/g and the fraction S_{f} of the surface that emits. Figure 2 shows the evolution of these parameters. We also note that our other parameters, such as M, R, and D, were allowed to vary freely. When the surface fraction was fixed, we obtained the correct g_{rad}/g with a precision of about 0.02 (in units of g_{rad}), as can be seen from the width of the 68% posterior distributions for models A and B. If S_{f} was taken to be free (model C), the uncertainty in g_{rad}/g increased by an order of magnitude, although the input value was within the uncertainty region. S_{f} is determined with a precision of about 10%. The precision for S_{f} is lower because only the spectral shape, rather than the amplitude, was used to match the model to the data.
Fig. 2 Evolution of the normalized luminosity g_{rad}/g for burst 1 of the synthetic data. Constraints on the surface fraction S_{f} are also shown for model C, which is the only one of our models in which this parameter is free. The dotted black line with the scale on the right axis shows the corresponding standard deviation of the obtained parameter distributions. The blue dashed line shows the input value used to create the data. The intensity of the red coloring is proportional to the probability density. This figure shows that when the fitting model is consistent with the model used to produce the data, we obtain parameter values that are accurate and precise. 
Fig. 3 Posterior distributions for the MCMC run with synthetic data for models A and B. The top panels of the triangles show the marginalized parameter posteriors (in arbitrary units); here the dark and light orange shadings give the 68% and 95% credible regions. The lower panels show the projected twodimensional parameter posteriors against each other. For these panels the solid line encloses the 95% credible regions. The blue stripes show the original values of the synthetic data that were used to create the data: R = 12 km, M = 1.5 M_{⊙}, X = 0, and D = 6 kpc. If the composition is known, the radius is precisely recovered. When X is a free parameter, the prior limit of X > 0 and the correlation of M, X, and D with each other leads to asymmetric posteriors around the true value, and so the inferred radius is a lower limit. 
Fig. 4 Mass and radius posteriors for synthetic data created for R = 12 km and M = 1.5 M_{⊙}, which are shown here with cyan crosses. The left panel shows a spectral fit with fixed emitting area S_{f} = 1 and hydrogen mass fraction X = 0 (model A). The right panel shows a spectral fit with a free hydrogen fraction X (model B). In both panels, the dotted line encapsulates the 68%, the dashed line the 95%, and the solid line the 99.7% credible regions. The dashed blue lines show values of constant redshift: 1 + z = 1.2 and 1.3. The dotted orange lines are the contours of constant surface gravity: log g = 14.0, 14.3 and 14.6. The solid green line shows the critical radius R = 4GM/c^{2}. The dark orange region in the top left corner marks the region of parameter space forbidden by the requirement of causality (Haensel et al. 2007; Lattimer & Prakash 2007). The input (M,R) point is within the 95% credible regions in both cases, even when the hydrogen mass fraction is a free parameter. 
Fig. 5 Posterior distributions for the MCMC run with synthetic data for model D with helium and solar compositions. The red solid line in the D panel shows the prior distribution () that we used. Other symbols and legends are the same as in Fig. 3. With the inclusion of the weak distance prior, the fit is better at recovering the original values. When the hydrogen mass fraction is not exactly zero, the true R, M, and X are recovered when a correct family of solutions is considered. In many cases, the incorrect highmass smallradius family of solutions is easy to discard on a physical basis as it is close to or inside the region ruled out by causality. 
Fig. 6 Mass and radius posteriors for synthetic data created for R = 12 km and M = 1.5 M_{⊙}, which are shown here with cyan crosses. The left panel shows the results for a pure helium composition (X = 0), whereas the right panel shows the results for a solar composition (X = 0.73). The symbols and legends are the same as in Fig. 4. Even with a free hydrogen mass fraction X, the true M−R values are now more consistent with the constraints coming from the fits. 
As a next step, we consider the parameters of greatest interest. Figure 3 shows the marginalized parameter distributions for M, R, D, and X (when sampled) along with the twodimensional posterior space projections. Figure 4 shows the MR projection in more detail. Model A is able to recover the correct radius R = 12 km with an accuracy of about 0.7 km after marginalization over all other parameters. The mass, on the other hand, is always underestimated and shows a bimodal structure that is familiar from cooling tail fits (Suleimanov et al. 2011a; Poutanen et al. 2014; Nättilä et al. 2016; Suleimanov et al. 2017). The input (M,R) value is just at the boundary of the 95% credible region. We determined the distance with a 68% scatter of about 0.6 kpc and no bias. The sharp spike in the marginalized distance distribution near the maximum value originates from the solutions near the critical line where R = 4GM/c^{2}. Mass and radius values on this line correspond to the solutions that are close to the maximum distance attainable for the system (see, e.g., Appendix A of Poutanen et al. 2014, for more discussion).
When the hydrogen fraction X is a free parameter and the data are analyzed using model B, the radius and the mass are both underestimated. This effect originates from the asymmetric X distribution, which arises because the lower limit of X = 0 is set by physical assumptions. For larger X we obtain lower values of M and R than what we obtain using X = 0. In the M−R plane, our proper solution is now only inside the 95% credible regions of the posteriors because of the strong bias toward lower masses. For similar reasons, the distance in this case is underestimated. The hydrogen fraction is constrained to be X < 0.2 with 95% credibility.
When X = 0.73, there is a similar underestimation of the values of the parameters. This is again caused by the connection between X, D, and M. In this case, the posterior for X is not symmetric around the true value because the distance has a maximum set by the observed flux level. This causes X to favor higher values (i.e., X > 0.73) and so the constraints on M and D are similar to the results from the analysis of model B when X = 0.
Ideally, we would like the method to be free of any bias in M,R,D, or X. These parameters are not our observables, however. Quantities closer to the observations include the redshift given in Eq. (2), which depends on M and R, and the surface gravity g defined by Eq. (6), which is also a function of M and R. The distance D and the hydrogen fraction X enter the system of equations via the flux (Eq. (5)), which in turn is limited by the Eddington flux (12). What we observe directly is the number of counts, which is related to the photon (number) flux of the source by Eq. (13). Because of this, all of our parameters are interconnected in a complicated fashion. In Bayesian inference, it is typical to study such a system by defining some information criterion and then to minimize its value given the fit parameters. This would then give us a description of the least informative, often multidimensional, priors. We elected instead to impose simple unidimensional priors for the system based on experience with our analysis of synthetic data. We found that M,X, and D are the most tightly connected parameters in the system. Hence, imposing a prior distribution for one of them will strongly affect the rest. We here decided to optimize our results for M and X at the cost of introducing a nonflat prior for D. Based on our different test runs, we concluded that a prior of P(D) ∝ D^{1/2} (which therefore slightly favors higher values of D) produces the least biased constraints for M and X.
This leads us to propose a fourth and final model, model D, which is an extension of model B. In model D the hydrogen fraction X was a free parameter, but we also incorporated intrinsic scatter σ_{int} into the analysis and chose P(D) ∝ D^{1/2}. Using these assumptions, we recovered the input parameters for the synthetic data, without any significant bias when X ≠ 0. This is evident in Figs. 5 and 6, where we used model D with synthetic data that have X = 0 (lefthand panels of both figures) and X = 0.73 (righthand panels of both figures). We find that when X is not exactly 0, we are able to reproduce the input radius without bias. Additionally, if the second cluster of solutions at high masses is neglected, we also reproduce the input mass, hydrogen fraction, and distance. In practice, this can be done by imposing a mass cutoff of M < 2.0 M_{⊙} or by selecting only the lowmass solutions below the critical radius R = 4GM/c^{2}. For a pure helium atmosphere, the sharp boundary at X = 0 leads to asymmetric posteriors around the real value, and so the estimates are always biased towards lower or higher values. An additional check is that the intrinsic scatter is driven toward low values, which it must be because we created the data without any additional scatter.
For pure helium (X = 0, model A) we constrained the radius to be R = 11.3 ± 0.4 (0.7) km at 68% (95%) credibility. The mass was similarly constrained to be . Thus the input values were inside the 95% credible intervals. The distance is found to be . The hydrogen mass fraction is constrained to be X < 0.05 (0.09), which is consistent with the input value X = 0. Similarly, for the synthetic solar data (X = 0.73), the radius, distance, and hydrogen mass fraction results are , , , and . Most importantly, we see that the correct radius is obtained without any bias. When the second highmass cluster of solutions is omitted by asserting an additional M < 2 M_{⊙} prior, we obtain , , , and . In this case, both the radius and the mass are correctly recovered with 0.9 km and 0.4 M_{⊙} precision, respectively. Even accounting for the mass imprecision when the composition is pure helium, model D is still the best of our models in producing precise and unbiased estimates of the parameters.
3.2. 4U 1702–429
We now study the RXTE data from 4U 1702−429. The bursts we use are from obsid 50025010100, 80033010108, 80033011904, 80033012002, and 80033012100, starting at MJDs of 51 781.333039, 52 957.629763, 53 211.964665, 53 212.794286, and 53 311.806086, respectively. The data were reduced in a way similar to the reduction in Galloway et al. (2008; see also Kajava et al. 2014; Poutanen et al. 2014; Nättilä et al. 2016). As we described in Sect. 3.1, we binned the data in time: each time the count rate decreased by a factor of approximately , we doubled the exposure time so that the number of counts in each bin remained relatively high. The RXTE data were also deadtime corrected (see, for example, Nättilä et al. 2016), which was of course not necessary for the synthetic data. We also note that unlike the synthetic data, the real observations have a varying “quality” because of the varying number of PCUs between the five bursts (ranging from 2 to 5 active PCUs). Some sample spectra from one of the bursts are shown in Fig. 7. We also show the evolution of each individual spectral parameter of each burst in Fig. 8. Just as with the synthetic data, we see that the normalized luminosity g_{rad}/g is well constrained, and the evolution of this parameter is strikingly similar to its evolution in the mock data. In the model C fits we see that the surface emitting fraction is constrained to be very close to unity for the entirety of this burst, which provides some evidence that the full surface emits close to uniformly in this case.
Fig. 7 Spectra for one Xray burst from 4U 1702−429 (crosses) with corresponding bestfit atmosphere models (solid lines) for A. Different colors show spectra for varying g_{rad}/g. Individual spectra are shifted in powers of 2 in the ydirection for clarity. 
Fig. 8 Evolution of the normalized luminosity g_{rad}/g for burst 1 from 4U 1702−429. The surface fraction S_{f} constraints are also shown for model C, which is the only one of our models in which this parameter is free. The dotted black line with the scale on the right axis shows the corresponding standard deviation of the obtained parameter distributions. The intensity of the red coloring is proportional to the probability density. The constraints we find using these real RXTE data bear a clear similarity to the constraints from our synthetic data fits seen in Fig. 2. 
Fig. 9 Upper panel: ratio of the data to the bestfit model for the timeresolved spectra of five bursts from 4U 1702−429 for Model D. Lower panel: deviation Δχ of the data from the model. Only energy bins where the number of counts exceeds 50 are shown. 
χ^{2} values for the atmosphere model bestfits for model D.
As a final test of our atmosphere model goodnessoffit, we can study possible systematic deviations from the model spectra in an individual energy channel level. In Fig. 9 we show the ratio of the data and the bestfit model flux (from the model D run) together with a channelspecific Δχ value (defined as , that is, the difference between data and model ℳ in units of standard error σ) for each of the five analyzed bursts. Additionally, Table 2 shows the sum of χ^{2} values for the bestfit models. There do not appear to be any significant persistent structures in the residuals. There is a slight deficit of flux around 10−14 keV in the burst tails, which corresponds to a ~ 5% difference between the model and the data (Δχ ~ 0.4). No such features are seen in the synthetic data fits. It is therefore possible that this deficit has a physical origin. We discuss these deviations further in Sect. 4. Our total χ^{2} (for model D) is 3103.5 with 2438 degrees of freedom (reduced of 1.27) for the set of all five bursts, when all spectral energy bins with fewer than 20 counts are ignored. The quality of the fit is, however, decreased drastically by the lowcount channels, and when we instead select 50 as our channel count cutof,f we obtain χ^{2} = 2541.0 with 2122 degrees of freedom ().
For our actual parameter constraints we considered only models A and D, because model B produced somewhat biased constraints with synthetic data and model C has too much freedom in its parameters (particularly with the inclusion of the surface emitting fraction S_{f} as a free parameter). Figure 10 shows the full posterior distributions, and Fig. 11 shows the twodimensional M−R posterior distributions in more detail for each model. The precision of constraints is similar to what it was for the synthetic data. We find that R = 12.4 ± 0.4 (0.6) km and M = 1.4 ± 0.2 (0.4) M_{⊙} for model A (which has fixed chemical composition X = 0), where we list the 68% (and 95% in parentheses) error regions. In model D the constraints are similar: and . We also find that X < 0.09 (0.16) at 68% (95%) credibility. The most probable value is X = 0.06 instead of X = 0. We find a distance of D = 5.5 ± 0.4 (0.7) kpc with model A or with model D. We note that in the model D fits, the second nonphysical highmass smallradius family of solutions is now mostly located inside the causality region, so it is naturally ruled out by physical considerations. Some small group of solutions, however, remains at M = 2.2 M_{⊙} and R = 10 km, which then shifts the lower limit of the 95% radius credibility interval down to 10 km instead of 11.4 km, which would be obtained by omitting it entirely. In contrast to the synthetic data, the real bursts also have a nonzero intrinsic scatter of . This corresponds to about 30 counts per second per energy channel. In reality, the intrinsic scatter accommodates more local deviations such as those between 10−14 keV, as seen in Fig. 9. Such an error in the observed counts reflects a ~ 1−5% deviation from the model flux, depending on the g_{rad}/g value (higher g_{rad}/g corresponds to higher temperature, i.e., higher count rate, for which the deviation is closer to the 1% level, whereas the opposite is true for a small g_{rad}/g).
4. Discussion
The measurements we presented here result from the use of full atmospheric spectral models of thermonuclear Xray burst cooling, not from the usual use of diluted blackbody fits. This gives us access to additional information via the surface redshift and the surface gravity. Our new method also allows us to validate many of the assumptions that underlie previous work. For example, we find that the Eddington limit is reached (and exceeded) near the beginning of the cooling tail (see Fig. 8, which shows the fit results for the first burst in our sample; the remaining four bursts are almost identical). Hence, this is the most direct validation yet (assuming that the atmosphere models are correct) that at least the bursts that we analyzed here are photospheric radius expansion bursts. Moreover, the dependence of the spectral shape evolution on the atmosphere composition allows us for the first time to set reliable limits on the hydrogen mass fraction in the photosphere. This is made possible by the fact that the temperature evolution of the atmosphere is dependent on the composition.
Fig. 10 Posterior distributions for the MCMC runs with real data for five PRE bursts from 4U 1702−429. The panels and symbols are the same as in Figs. 3 and 5. The red solid line in the Dpanel shows the prior distribution () that we used. Both models are seen to produce posterior shapes that are similar to what we found in our synthetic data fits (Fig. 5). Both models also produce consistent estimates for the radius and distance. 
Fig. 11 Mass and radius posteriors for 4U 1702−429. The left panel shows the results for model A, which has fixed S_{f} = 1 and X = 0. The right panel shows the results for model D, which has free X and fixed S_{f} = 1. The symbols and legends are the same as in Fig. 4. We recall that for our synthetic data fits (shown in Fig. 6), model A underestimates the mass slightly. For X ≠ 0 the correct mass is obtained using model D. This suggests that the true mass of 4U 1702−429 is in the range M = 1.4−2.0 M_{⊙}. Models A and D are both capable of recovering the radius used to construct the synthetic data, which suggests that the radius of 4U 1702−429 is R = 12.0−12.9 km at 68% credibility. 
4.1. Uncertainties and systematic errors
By fitting the atmosphere models directly to the data, we can also assess the degree to which the models represent the data. Although we reiterate our caveats about the use of χ^{2}, particularly for model comparisons (see also Andrae et al. 2010), we note that the model D fit to 4U 1702−429 has χ^{2}/ d.o.f. = 2541.0/2122 (see also Table 2), whereas for a simple blackbody fit χ^{2}/ d.o.f. = 2716.3/2010. The blackbody fits were obtained using xspec version 12.9.1 (Arnaud 1996) with bbodyrad model. Thus the model atmosphere fit has an additional 112 degrees of freedom, but its χ^{2} is 175.5 smaller. This kind of comparison is, however, not strictly fair because individual blackbody best fits minimize the quoted χ^{2} value for channels with more than 50 counts, whereas the results from the hierarchical modeling are obtained from the MCMC chain that considers full Poisson or Gaussian likelihoods. The bestfit values are also very susceptible to the exact energy range used for the fit, and so the χ^{2} values reported here are only indicative. When an extra 0.5% calibration error is introduced, as is advised for RXTE spectral analysis (Shaposhnikov et al. 2012), the χ^{2} value also decreases, in both cases, by about 30. The main difference here that we wish to emphasize is that blackbody fits involve two parameters per spectrum (temperature and normalization), whereas direct spectral fitting only has g_{rad}/g as a parameter for each spectrum (this is true for models A, B and D; for model C, the surface emitting fraction is also a free parameter for each spectrum). In addition to the individual spectrum parameters, there are between three to five global parameters, and thus the total number of model parameters is 1 × 116 + [3,4] (mass, radius, distance, and hydrogen mass fraction in some cases). This means that when compared with the blackbody model, the atmosphere model is able to reduce the number of parameters needed from ~ 230 to roughly 120 while retaining comparable accuracy.
Nonetheless, the formal statistical fit is not good, and this is reflected in the nonzero value of σ_{int} in our model D fits. Hence it appears that there are unmodeled effects in the data. An obvious candidate for such complications is the rotation of the star. However, because 4U 1702−429 has a relatively low rotational frequency of 329 Hz (Markwardt et al. 1999), the effects are unlikely to be large for this star. Our preliminary studies show that the effect of rotational broadening of the spectrum is strongest at low and high energies, and hence broadening might account for some of the deviations seen at E ≲ 3 keV and E ≳ 12 keV in Fig. 9. The impact on the fit quality (i.e., χ^{2}), on the other hand, is small because rotational broadening tends to smooth out the spectra without producing any sharp features (Nättilä & Pihajoki 2017). Hence, distinguishing the effects of rotation from the atmosphere model fits based on the fit quality alone is hard.
As is the case for rotational smearing, we cannot detect whether there is nonuniform surface emission in the sense that the temperature varies across the surface. This is because our model C, which has a free surfaceemitting fraction S_{f}, can only capture effects where some part of the star is partially covered. There have been detections of burst oscillations from 4U 1702−429 that imply a nonuniform surface temperature (Markwardt et al. 1999; Galloway et al. 2008; Ootes et al. 2017). However, these have been detected only during the softstate bursts from this particular source (Ootes et al. 2017), and so it could be that the effect, if it exists, is small because it cannot be detected using RXTE data.
Another possible source of error is the treatment of heavy elements in the atmosphere. Nättilä et al. (2015) showed that heavy elements can have a significant impact on burst spectra. Most importantly, heavy elements produce photoionization edges around E ~ 9−14 keV. The metals are likely to originate from nuclear burning during the burst and might be brought to the surface layers of the star by convection (Weinberg et al. 2006; Malone et al. 2011, 2014). The spectral deviations are expected to appear mainly when the metals start to recombine at lower temperatures after the photosphere has cooled down. The detection of such features in other sources with longer more energetic bursts (in’t Zand & Weinberg 2010; Kajava et al. 2017b) implies that they might also play some role in the shorter bursts analyzed here.
A third possible source of deviations involves the persistent (nonburst) emission, which we currently assume to be constant during the fit. Some recent studies indicate that this might not be the case even in the hard state (Ji et al. 2015; Degenaar et al. 2016; Kajava et al. 2017c). However, the initial level of persistent emission for the five bursts from 4U 1702−42 is very low and we expect this effect to not have a very significant impact on the observed radiation. A crude estimate can be obtained by varying the background emission, not with the full hierarchical model, but with individual blackbody fits. In this case, the constraints come from the Eddington limit and from the normalization K_{tail} = (R_{bb} [km] /D_{10})^{2} in the burst tail. Here R_{bb} is the blackbody radius (∝ R) and D_{10} = D/ 10 kpc. By varying the background emission with factors ranging from 0.5 to 2, we obtain constraints where the Eddington flux is still accurately recovered, but the apparent normalization in the tail is decreased or increased, respectively. The value of the normalization in the tail, in this case, is within 2% of the original value. This is consistent with the fact that the persistent emission during the hard state is about 1% of the Eddington flux, whereas in the tail F ~ 0.5F_{Edd} and so a 2% scatter in the normalization K is expected as F = KT^{4}. The measured radius R scales roughly (for the fixed compactness) as (see, e.g., Eq. (A9) in Poutanen et al. 2014) so that such a deviation in the normalization results in a ~ 2% scatter in the radius. This means that the uncertainty related to the persistent emission would then translate into about 250 m absolute error in our measured radius.
The final possible source of error is the neutral hydrogen column density. It only affects the lowenergy channels of RXTE, which can have an impact on the parameters deduced late in the burst tail. Near the Eddington limit, the radiation peaks at E ~ 10 keV, but when the NS cools down, the bulk of the thermal radiation moves to lower energies. Hence, the effect is similar to the aforementioned persistent emission where the main effect is on the normalization in the tail of the burst. When the value of N_{H} is decreased, the modeled lowenergy radiation is affected less by the absorption, and so the inferred value of the normalization in the tail also decreases because the model flux is now higher. Unfortunately, it is also in the RayleighJeans tail that the surface gravity of the atmosphere models has the greatest effect on the emergent spectra. Similar considerations as in the persistent emission case show that varying N_{H} by a factor of 2 leads to an error in K_{tail} of 5% that then translates into similar relative error in radius. We do, however, note that the measured N_{H} is usually obtained by other instruments that operate at lower Xray energies where it is easier to measure the neutral hydrogen column density. Hence, an uncertainty of a factor of 2 certainly overestimates the error related to the value.
Our consideration of error sources leads us to propose that the emission above F ≳ 0.5F_{Edd} probably is the cleanest option for M−R measurements. However, the high flux near F ~ F_{Edd} is not free of problems: early in the cooling tail, the count rates from the source are highest and so the detector is affected by the deadtime correction the most. Deadtime correction near the peak can be as high as 5%, which would directly translate into error in the measured F_{Edd}. This would again translate into a similar uncertainty in the radius. In reality, of course, the deadtime correction scheme proposed by the instrument calibration team should be quite effective at covering this effect, and so errors as large as 5% originating from this are not expected. To be safe, the fluxes between (0.5−0.95)F_{Edd} should give the most stringent constraints. However, this decreases the amount of available data even more, and for example, here in this work, we pushed the aforementioned limits to cover g_{rad}/g = (0.2−0.98) that we still think are viable. Another option would be to try to model the varying background emission and also marginalize over some plausible hydrogen column density range to capture all of the known error sources. It would be useful to perform such an analysis in the future. All in all, this shows that we are approaching the absolute measuring accuracy of the RXTE satellite.
The bestfit results are also robust against any systematic calibration error in the flux normalization. Because the constraints for R mainly originate from F_{Edd} and K_{tail} (in contrast to the redshift 1 + z and surface gravity g, which have a much weaker effect), the radius is mainly constrained by the temperature evolution of the burst alone. For an unknown systematic energyindependent shift ζ affecting the observed spectra we still obtain , where T_{Edd,∞} is the Eddington temperature (see Eq. (A9) in Poutanen et al. 2014 and Eq. (12)). This is a distanceindependent quantity that makes the derived radius independent of any normalization factor ζ affecting the observed flux.
4.2. Comparison and robustness of the constraints
It is also interesting to compare our analysis of the 4U 1702−429 bursts to previous constraints that were obtained using the cooling tail method. By applying the cooling tail method to the same set of hardstate bursts that we analyze here, Nättilä et al. (2016) measured the NS radius to be R ≈ 13 km for M = 1.5 M_{⊙}. However, their M−R posteriors have a complicated bananalike shape (see Fig. 4 in Nättilä et al. 2016), and thus the inferred radius depends strongly on the assumed mass. Nättilä et al. (2016) found that introducing priors on the EoS led to better constraints on the mass. That is, the assumption that all of the sources analyzed in Nättilä et al. (2016; in addition to 4U 1702−429, they used 4U 1724−307 and SAX J1810.8−2609) originate from the same underlying EoS helps pin down the mass. They find that at 68% probability, M = 1.8 ± 0.3 M_{⊙} and R = 11.9 ± 0.6 km, assuming no phase transitions (QMC+model A in their paper). These constraints are in a good agreement with the values derived here, which from model D are M ≈ 1.9 M_{⊙} and R ≈ 12.4 km. The compositions are also in good agreement: Nättilä et al. (2016) assumed a pure helium composition (X = 0), whereas here the fit itself shows that X < 0.09 (68%). The distance constraints also agree well: D = 5.6 ± 0.9 kpc versus 5.5 ± 0.4 kpc, for the cooling tail and the direct spectral fitting methods (model A), respectively.
We note that our M−R results are located away from the critical radius R = 4GM/c^{2} (see Özel & Psaltis 2015, for a discussion). If the constraints from F_{Edd} and A (∝ K_{tail}) are not consistent with each other^{2}, then the M−R solution is forced to obey this relation because no real solution exists. This could happen, for example, if the model is applied to data that it does not describe, such as data from softstate bursts where the behavior of the cooling tail might not be entirely determined by the NS surface alone (see Steiner et al. 2010; Poutanen et al. 2014; Kajava et al. 2014; Nättilä et al. 2016).
If the true values of M and R are close to the R = 4GM/c^{2} relation, then it is very hard to distinguish this correct solution from an incorrect solution that is forced upon the system by a model that is inconsistent with the data. Hence, if all of the M−R solutions from multiple sources are located only at this line (see Özel et al. 2016, for such a situation), it either means that the model is applied inconsistently to data that it does not describe, or that all NSs happen to have the same compactness M/R = c^{2}/ 4G.
Another interesting aspect of our method is its ability to constrain the composition of the atmosphere. As can be seen from the synthetic data fits, it is possible to set limits for X with about 10% precision at 68% credibility. This opens up a whole new window to the study of accretion physics because we can correlate the burst behavior against the composition of the accreted matter. We should, however, note that the composition we probe here is the composition during the burst, and so in theory, the nuclear reprocessing might change the true composition during the measurement.
Although X = 0 is still consistent with the data, we can ask whether other aspects of the 4U 1702−429 bursts are consistent with there being some hydrogen in the atmosphere. One such consistency check involves the ratio of the fluence of the persistent emission between bursts to the burst fluence itself. This ratio, which is usually called α, is a measure of the ratio of the gravitational specific energy release to the thermonuclear energy release; because hydrogen fusion releases much more energy than helium fusion, we expect α to be larger when there is less hydrogen present. For 4U 1702−429, α ≈ 75 (Galloway et al. 2008), whereas for 4U 1820−303 α is in the range 125–155 (Haberl et al. 1987). The neutron star in 4U 1820−303 is usually assumed to have a nearly pure helium atmosphere (Cumming 2003), so the lower value of α in 4U 1702−429 is consistent with the presence of some hydrogen in the latter source. We note, however, that the values of α quoted here are the minimum values and might change from burst to burst. Bursts from 4U 1820−303 also exhibit a fast rise and have short timescales, both of which are believed to be consequences of fast helium burning, whereas the 4U 1702−429 bursts have longer durations and also longer rise times (Galloway et al. 2008). Neither of these findings is conclusive, but they do point into the same direction: there should be some traces of hydrogen in the 4U 1702−429 atmosphere, as is suggested by our analysis.
4.3. Future prospects
It is interesting to consider additional possibilities that are suggested by our new and detailed analysis. One obvious extension is to include the PRE phase in our fitting. To do this, however, we will need accurate atmosphere models of extended NS photospheres. The advantage would be to increase the available data for analysis, and it will probably also significantly improve the measurements of M and R because the expansion must be heavily dependent on the redshift z and surface gravity g. Preliminary work into this direction has already been reported in Medin et al. (2016). Their work also allows an independent validation of our atmosphere models. Our results agree well with theirs in the range g_{rad}/g ≈ 0.2−1, implying that at least in the context of the mutual assumptions set by both computations, the results are reproducible.
Another important, but computationally very expensive future prospect is to fit all possible Xray bursts to obtain M−R constraints. This would help to set groundbreaking constraints on the EoS of the dense matter. Preliminary studies already validate the previous results that the atmosphere models are not applicable to the softstate bursts (Poutanen et al. 2014; Kajava et al. 2014). Last, it is also important to understand why the models do not agree completely with the data and the physical origin of these deviations. For this, more work is needed in order to understand the physics and environments of the bursts better.
Despite the uncertainty about the mass, we have significantly improved the constraints on the compactness of the neutron star. Our mass and radius measurements are encouragingly consistent with recent theoretical analyses of the EoS of cold dense matter (Lattimer & Prakash 2016). This result gives hope that we may use astrophysical neutron star measurements to better constrain the behavior of the ultradense matter.
5. Summary
We have presented the first direct atmosphere model spectral fits to thermonuclear Xray burst cooling tails. Our method is a generalization of previous work, in which blackbody parameters were used as a proxy to trace the evolution of the energy spectrum. By fitting the atmosphere models directly to the data, we were able to extract more information from the data and also to test some of the physical assumptions made in previous analyses.
We find that fits to synthetic data, which are generated using the same model that we employed for our analyses, as expected reveal a lack of bias and also show the prospects for precise measurements of the mass and radius. When we applied our fitting procedure to RXTE data from five hardstate typeI Xray bursts from 4U 1702−429 in Sect. 3.2, the resulting posteriors were clearly similar to the synthetic data, although the formal quality of the fits is worse than in the ideal case. When we artificially added intrinsic noise in our analysis of the 4U 1702−429 data to produce a formally good fit, we found that the radius was constrained to be R = 12.4 ± 0.4 km at 68% credibility for the two models we employed. The source distance was constrained to be between 5.1 kpc <D < 6.2 kpc (68% combined credibility limits from models A and D). We found that the hydrogen mass fraction X for 4U 1702−429 could be constrained to X < 0.09 at 68% credibility. The highestprobability value was X = 0.06 rather than X = 0.
The mass seems to be the hardest parameter to constrain. When we applied our two models to synthetic data, model A typically underestimated M by ~ 0.1 M_{⊙}, whereas model D showed an even stronger underestimation for atmospheres with no hydrogen in them. When X > 0, model D was seen to reproduce the mass when applied to synthetic data. Our analysis of the RXTE data for 4U 1702−429 yielded similar results: model A gave M = 1.4 ± 0.2 M_{⊙}, whereas model D gave M = 1.9 ± 0.3 M_{⊙}. If a bias similar to what we find when analyzing synthetic data applies to the 4U 1702−429 analysis, then the real mass is expected to lie closer to the model D constraints. We therefore suggest that the 95% credible interval of model D is a trustworthy limit, and in this limit we find that the mass for 4U 1702−429 lies in the range 1.4 <M/M_{⊙} < 2.2.
In our case, asserting an informative distance prior leads to a flat posterior in mass. However, this will most likely also affect physical observables such as the flux. Thus, even though we appear unbiased in mass, we are now biased in the flux. Because of this, the aforementioned distance prior was only imposed for this one model to study the possible effects it might have.
Acknowledgments
We appreciate the detailed comments provided by the referee that helped to improve and clarify the paper. This research was supported by the University of Turku Graduate School in Physical and Chemical Sciences (JN). M.C.M. was supported in part by NASA NICER grant NNX16AD90G. A.W.S. was supported by grant NSF PHY 1554876. J.J.E.K. acknowledges support from the ESA research fellowship programme and the Academy of Finland grants 268740 and 295114. J.N. and J.J.E.K. acknowledge support from the Faculty of the European Space Astronomy Centre (ESAC). V.F.S. was supported by the German Research Foundation (DFG) grant WE 1312/481 and by the Russian Government Program of Competitive Growth of Kazan Federal University. J.P. thanks the Foundations’ Professor Pool, the Finnish Cultural Foundation and the National Science Foundation grant PHY1125915 for support. The computer resources of the Finnish IT Center for Science (CSC) and the FGCI project (Finland) are acknowledged.
References
 Andrae, R., SchulzeHartung, T., & Melchior, P. 2010, ArXiv eprints [arXiv:1012.3754] [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, 17 [Google Scholar]
 Ballantyne, D. R., & Everett, J. E. 2005, ApJ, 626, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Ballantyne, D. R., & Strohmayer, T. E. 2004, ApJ, 602, L105 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Y.P., Zhang, S., Torres, D. F., et al. 2011, A&A, 534, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cumming, A. 2003, ApJ, 595, 1077 [NASA ADS] [CrossRef] [Google Scholar]
 Damen, E., Magnier, E., Lewin, W. H. G., et al. 1990, A&A, 237, 103 [NASA ADS] [Google Scholar]
 Degenaar, N., Miller, J. M., Wijnands, R., Altamirano, D., & Fabian, A. C. 2013, ApJ, 767, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Degenaar, N., Koljonen, K. I. I., Chakrabarty, D., et al. 2016, MNRAS, 456, 4256 [NASA ADS] [CrossRef] [Google Scholar]
 Ebisuzaki, T. 1987, PASJ, 39, 287 [NASA ADS] [Google Scholar]
 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Guillot, S., Servillat, M., Webb, N. A., & Rutledge, R. E. 2013, ApJ, 772, 7 [Google Scholar]
 Güver, T., Özel, F., CabreraLavers, A., & Wroblewski, P. 2010, ApJ, 712, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Haberl, F., Stella, L., White, N. E., Gottwald, M., & Priedhorsky, W. C. 1987, ApJ, 314, 266 [NASA ADS] [CrossRef] [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. Lib., 326 [Google Scholar]
 Humphrey, P. J., Liu, W., & Buote, D. A. 2009, ApJ, 693, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Ibragimov, A. A., Suleimanov, V. F., Vikhlinin, A., & Sakhibullin, N. A. 2003, Astron. Rep., 47, 186 [NASA ADS] [CrossRef] [Google Scholar]
 in’t Zand, J. J. M., & Weinberg, N. N. 2010, A&A, 520, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 in’t Zand, J. J. M., Galloway, D. K., & Ballantyne, D. R. 2011, A&A, 525, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Ji, L., Zhang, S., Chen, Y., et al. 2015, ApJ, 806, 89 [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]
 Kajava, J. J. E., Koljonen, K. I. I., Nättilä, J., Suleimanov, V., & Poutanen, J. 2017a, MNRAS, 472, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Kajava, J. J. E., Nättilä, J., Poutanen, J., et al. 2017b, MNRAS, 464, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Kajava, J. J. E., SánchezFernández, C., Kuulkers, E., & Poutanen, J. 2017c, A&A, 599, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koljonen, K. I. I., Kajava, J. J. E., & Kuulkers, E. 2016, ApJ, 829, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. 1993, CDROMs (Cambridge, Mass.: Smithsonian Astrophysical Observatory), No. 11 [Google Scholar]
 Kurucz, R. L. 1970, SAO Special Report, 309 [Google Scholar]
 Kuśmierek, K., Madej, J., & Kuulkers, E. 2011, MNRAS, 415, 3344 [NASA ADS] [CrossRef] [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]
 Lattimer, J. M. 2012, Annu. Rev. Nucl. Part. Sci., 62, 485 [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2016, Phys. Rep., 621, 127 [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]
 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, ApJ, 788, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Markwardt, C. B., Strohmayer, T. E., & Swank, J. H. 1999, ApJ, 512, L125 [NASA ADS] [CrossRef] [Google Scholar]
 Medin, Z., von Steinkirch, M., Calder, A. C., et al. 2016, ApJ, 832, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C. 2013, ArXiv eprints [arXiv:1312.0029] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1996, ApJ, 470, 1033 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 2016, Eur. Phys. J. A, 52, 63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miller, M. C., Boutloukos, S., Lo, K. H., & Lamb, F. K. 2011, in Fast Xray Timing and Spectroscopy at Extreme Count Rates, PoS (HTRS 2011)24 [Google Scholar]
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co.) [Google Scholar]
 Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ, 580, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Nättilä, J., & Pihajoki, P. 2017, A&A, submitted [arXiv:1709.07292] [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]
 Nättilä, J., Steiner, A. W., Kajava, J. J. E., Suleimanov, V. F., & Poutanen, J. 2016, A&A, 591, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ootes, L. S., Watts, A. L., Galloway, D. K., & Wijnands, R. 2017, ApJ, 834, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F. 2006, Nature, 441, 1115 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Özel, F. 2013, Rept. Progr. Phys., 76, 016901 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., & Psaltis, D. 2015, ApJ, 810, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., Güver, T., & Psaltis, D. 2009, ApJ, 693, 1775 [NASA ADS] [CrossRef] [Google Scholar]
 Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Peille, P., Olive, J.F., & Barret, D. 2014, A&A, 567, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poutanen, J. 2017, ApJ, 835, 119 [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]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: WileyInterscience) [Google Scholar]
 Serino, M., Mihara, T., Matsuoka, M., et al. 2012, PASJ, 64, 91 [NASA ADS] [Google Scholar]
 Shaposhnikov, N., Jahoda, K., Markwardt, C., Swank, J., & Strohmayer, T. 2012, ApJ, 757, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W. 2014a, Astrophysics Source Code Library [record ascl: 1408.020] [Google Scholar]
 Steiner, A. W. 2014b, Astrophysics Source Code Library [record ascl: 1408.019] [Google Scholar]
 Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Strohmayer, T., & Bildsten, L. 2006, in Compact stellar Xray sources, eds. W. Lewin, & M. van der Klis (Cambridge: Cambridge University Press), Cambridge Astrophys. Ser., 39, 113 [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Werner, K. 2007, A&A, 466, 661 [NASA ADS] [CrossRef] [EDP Sciences] [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. F., Poutanen, J., Klochkov, D., & Werner, K. 2016, Eur. Phys. J. A, 52, 20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V. F., Poutanen, J., Nättilä, J., et al. 2017, MNRAS, 466, 906 [NASA ADS] [CrossRef] [Google Scholar]
 van Paradijs, J., & Lewin, H. G. 1986, A&A, 157, L10 [NASA ADS] [Google Scholar]
 van Paradijs, J., Dotani, T., Tanaka, Y., & Tsuru, T. 1990, PASJ, 42, 633 [NASA ADS] [Google Scholar]
 Walker, M. A. 1992, ApJ, 385, 642 [NASA ADS] [CrossRef] [Google Scholar]
 Watts, A. L. 2012, ARA&A, 50, 609 [NASA ADS] [CrossRef] [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]
 Worpel, H., Galloway, D. K., & Price, D. J. 2015, ApJ, 801, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, W., Li, T. P., Zhang, W., & Zhang, S. N. 1999, ApJ, 512, L35 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Synthetic spectra (crosses) with corresponding bestfit atmosphere models (solid lines) for A. Different colors show spectra for varying g_{rad}/g. Individual spectra are shifted by factors of 2 in the ydirection for clarity. 

In the text 
Fig. 2 Evolution of the normalized luminosity g_{rad}/g for burst 1 of the synthetic data. Constraints on the surface fraction S_{f} are also shown for model C, which is the only one of our models in which this parameter is free. The dotted black line with the scale on the right axis shows the corresponding standard deviation of the obtained parameter distributions. The blue dashed line shows the input value used to create the data. The intensity of the red coloring is proportional to the probability density. This figure shows that when the fitting model is consistent with the model used to produce the data, we obtain parameter values that are accurate and precise. 

In the text 
Fig. 3 Posterior distributions for the MCMC run with synthetic data for models A and B. The top panels of the triangles show the marginalized parameter posteriors (in arbitrary units); here the dark and light orange shadings give the 68% and 95% credible regions. The lower panels show the projected twodimensional parameter posteriors against each other. For these panels the solid line encloses the 95% credible regions. The blue stripes show the original values of the synthetic data that were used to create the data: R = 12 km, M = 1.5 M_{⊙}, X = 0, and D = 6 kpc. If the composition is known, the radius is precisely recovered. When X is a free parameter, the prior limit of X > 0 and the correlation of M, X, and D with each other leads to asymmetric posteriors around the true value, and so the inferred radius is a lower limit. 

In the text 
Fig. 4 Mass and radius posteriors for synthetic data created for R = 12 km and M = 1.5 M_{⊙}, which are shown here with cyan crosses. The left panel shows a spectral fit with fixed emitting area S_{f} = 1 and hydrogen mass fraction X = 0 (model A). The right panel shows a spectral fit with a free hydrogen fraction X (model B). In both panels, the dotted line encapsulates the 68%, the dashed line the 95%, and the solid line the 99.7% credible regions. The dashed blue lines show values of constant redshift: 1 + z = 1.2 and 1.3. The dotted orange lines are the contours of constant surface gravity: log g = 14.0, 14.3 and 14.6. The solid green line shows the critical radius R = 4GM/c^{2}. The dark orange region in the top left corner marks the region of parameter space forbidden by the requirement of causality (Haensel et al. 2007; Lattimer & Prakash 2007). The input (M,R) point is within the 95% credible regions in both cases, even when the hydrogen mass fraction is a free parameter. 

In the text 
Fig. 5 Posterior distributions for the MCMC run with synthetic data for model D with helium and solar compositions. The red solid line in the D panel shows the prior distribution () that we used. Other symbols and legends are the same as in Fig. 3. With the inclusion of the weak distance prior, the fit is better at recovering the original values. When the hydrogen mass fraction is not exactly zero, the true R, M, and X are recovered when a correct family of solutions is considered. In many cases, the incorrect highmass smallradius family of solutions is easy to discard on a physical basis as it is close to or inside the region ruled out by causality. 

In the text 
Fig. 6 Mass and radius posteriors for synthetic data created for R = 12 km and M = 1.5 M_{⊙}, which are shown here with cyan crosses. The left panel shows the results for a pure helium composition (X = 0), whereas the right panel shows the results for a solar composition (X = 0.73). The symbols and legends are the same as in Fig. 4. Even with a free hydrogen mass fraction X, the true M−R values are now more consistent with the constraints coming from the fits. 

In the text 
Fig. 7 Spectra for one Xray burst from 4U 1702−429 (crosses) with corresponding bestfit atmosphere models (solid lines) for A. Different colors show spectra for varying g_{rad}/g. Individual spectra are shifted in powers of 2 in the ydirection for clarity. 

In the text 
Fig. 8 Evolution of the normalized luminosity g_{rad}/g for burst 1 from 4U 1702−429. The surface fraction S_{f} constraints are also shown for model C, which is the only one of our models in which this parameter is free. The dotted black line with the scale on the right axis shows the corresponding standard deviation of the obtained parameter distributions. The intensity of the red coloring is proportional to the probability density. The constraints we find using these real RXTE data bear a clear similarity to the constraints from our synthetic data fits seen in Fig. 2. 

In the text 
Fig. 9 Upper panel: ratio of the data to the bestfit model for the timeresolved spectra of five bursts from 4U 1702−429 for Model D. Lower panel: deviation Δχ of the data from the model. Only energy bins where the number of counts exceeds 50 are shown. 

In the text 
Fig. 10 Posterior distributions for the MCMC runs with real data for five PRE bursts from 4U 1702−429. The panels and symbols are the same as in Figs. 3 and 5. The red solid line in the Dpanel shows the prior distribution () that we used. Both models are seen to produce posterior shapes that are similar to what we found in our synthetic data fits (Fig. 5). Both models also produce consistent estimates for the radius and distance. 

In the text 
Fig. 11 Mass and radius posteriors for 4U 1702−429. The left panel shows the results for model A, which has fixed S_{f} = 1 and X = 0. The right panel shows the results for model D, which has free X and fixed S_{f} = 1. The symbols and legends are the same as in Fig. 4. We recall that for our synthetic data fits (shown in Fig. 6), model A underestimates the mass slightly. For X ≠ 0 the correct mass is obtained using model D. This suggests that the true mass of 4U 1702−429 is in the range M = 1.4−2.0 M_{⊙}. Models A and D are both capable of recovering the radius used to construct the synthetic data, which suggests that the radius of 4U 1702−429 is R = 12.0−12.9 km at 68% credibility. 

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.