Neutron star mass and radius measurements from atmospheric model fits to X-ray burst cooling tail spectra

J. Nättilä; M. C. Miller; A. W. Steiner; J. J. E. Kajava; V. F. Suleimanov; J. Poutanen

doi:10.1051/0004-6361/201731082

Home

All issues

Volume 608 (December 2017)

A&A, 608 (2017) A31

Full HTML

Free Access

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/0004-6361/201731082
Published online		01 December 2017

A&A 608, A31 (2017)

Neutron star mass and radius measurements from atmospheric model fits to X-ray burst cooling tail spectra

J. Nättilä¹^,2, M. C. Miller³, A. W. Steiner⁴, J. J. E. Kajava⁵^,1^,6, V. F. Suleimanov⁷^,8 and J. Poutanen¹^,2^,9

¹ Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland
e-mail: joonas.a.nattila@utu.fi
² Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
³ Department of Astronomy and Joint Space-Science Institute, University of Maryland, College Park, MD 20742-2421, USA
⁴ Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA
⁵ Finnish Centre for Astronomy with ESO (FINCA), University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland
⁶ European Space Astronomy Centre (ESA/ESAC), Science Operations Department, 28691 Villanueva de la Cañada, Madrid, Spain
⁷ Institut für Astronomie und Astrophysik, Kepler Centre for Astro and Particle Physics, Universität Tübingen, Sand 1, 72076 Tübingen, Germany
⁸ Astronomy Department, Kazan (Volga region) Federal University, Kremlyovskaya str. 18, 420008 Kazan, Russia
⁹ Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA

Received: 2 May 2017
Accepted: 26 September 2017

Abstract

Observations of thermonuclear X-ray bursts from accreting neutron stars (NSs) in low-mass X-ray 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 X-ray bursting NSs directly to the observed spectra. We present a hierarchical Bayesian fitting framework that uses current X-ray 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 X-ray Timing Explorer observations of five hard-state X-ray 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 / X-rays: binaries / X-rays: 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 low-mass X-ray binaries (LMXBs) that exhibit frequent thermonuclear X-ray 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 X-ray 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 X-ray 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 hard-state 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 hard-state 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 X-ray 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 X-ray 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 $I_{E^{'}}^{'}$ $\hbox{$I'_{E'}$}$ 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 $\frac{E}{E^{'}} = \frac{1}{1 + z},$ $\begin{equation} \frac{E}{E'} = \frac{1}{1+z}, \end{equation}$ (1)where z takes into account the rotation-induced 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 $\lim_{ν \to 0} (1 + z) = {(1 - \frac{2 GM}{c^{2} R})}^{- 1 / 2},$ $\begin{equation} \lim_{\nu \rightarrow 0} (1+z) = \left(1 - \frac{2 G M}{c^2 R} \right)^{-1/2}, \label{eq:redshift} \end{equation}$ (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 $I_{E^{'}}^{'}$ $\hbox{$I'_{E'}$}$ (Liouville’s theorem for photons, see, e.g., Misner et al. 1973; Rybicki & Lightman 1979) by $I_{E} = {(\frac{E}{E^{'}})}^{3} I_{E^{'}}^{'} .$ $\begin{equation} I_{E} = \left( \frac{E}{E'} \right)^3 I'_{E'}. \end{equation}$ (3)The total observed monochromatic flux from the star, as seen by a distant observer, is then (see, e.g., Nättilä & Pihajoki 2017) $F_{obs} (E) = \int I_{E} dΩ,$ $\begin{equation} F_{\mathrm{obs}}(E) = \int I_E {\rm d}\Omega, \label{eq:fluxint} \end{equation}$ (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: $F_{obs} (E) = \frac{F^{'} (E^{'})}{(1 + z)^{3}} {(\frac{R_{\infty}}{D})}^{2} = \frac{F^{'} (E^{'})}{1 + z} {(\frac{R}{D})}^{2},$ $\begin{equation} F_{\mathrm{obs}}(E) = \frac{F'(E')}{(1+z)^3} \left( \frac{R_{\infty}}{D} \right)^2 = \frac{F'(E')}{1+z} \left( \frac{R}{D} \right)^2, \label{eq:FobsE} \end{equation}$ (5)where R_∞ = R(1 + z) is the apparent NS radius, D is the distance to the source, and $F^{'} (E^{'}) = 2 π {^{\int}}_{0}^{1} I_{E^{'}}^{'} (μ) μ d μ$ $\hbox{$F'(E')=2\pi \int_0^1 I'_{E'} (\mu)\mu \,{\rm d}\mu$}$ , 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 $g = \frac{GM}{R^{2}} (1 + z) .$ $\begin{equation} g = \frac{G M}{R^2} (1+z). \label{eq:gravity} \end{equation}$ (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 $κ_{T} \approx 0.2 (1 + X) {cm}^{2} g^{-1} .$ $\begin{equation} \kappa_{\mathrm{T}} \approx 0.2(1 + X)~\mathrm{cm}^2\,\mathrm{g}^{-1}. \end{equation}$ (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 $L_{Edd} = \frac{4 πGMc}{κ_{T}} (1 + z) .$ $\begin{equation} L_{\mathrm{Edd}} = \frac{4\pi G M c}{\kappa_{\mathrm{T}}} (1+z). \label{eq:Ledd} \end{equation}$ (8)The actual critical luminosity is reached when the radiative acceleration $g_{rad} = \frac{κ_{R}}{c} ℱ$ $\begin{equation} g_{\mathrm{rad}} = \frac{\kappa_{\rm R}}{c} \mathcal{F} \end{equation}$ (9)equals the surface gravitational acceleration g. Here κ_R is the flux mean opacity (equal in our case to the Rosseland mean opacity), $ℱ = \int F^{'} (E^{'}) d E^{'} = σ_{SB} T_{eff}^{4}$ $\begin{equation} \mathcal{F} = \int F'(E') {\rm d}E' = \sigma_{\mathrm{SB}} T_{\mathrm{eff}}^4 \end{equation}$ (10)is the bolometric surface flux, T_eff is the effective temperature of radiation, and σ_SB is the Stefan-Boltzmann constant. At high temperatures close to the critical luminosity, the opacity is dominated by Compton scattering and is lower than κ_T because of the Klein-Nishina 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 $ℱ_{Edd} = \frac{L_{Edd}}{4 π R^{2}} = σ_{SB} T_{Edd}^{4} = \frac{gc}{κ_{T}},$ $\begin{equation} \mathcal{F}_{\mathrm{Edd}} = \frac{L_{\mathrm{Edd}}}{4 \pi R^2} = \sigma_{\mathrm{SB}} T_{\mathrm{Edd}}^4 = \frac{g c}{\kappa_{\mathrm{T}}}, \end{equation}$ (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: $F_{Edd} = A \frac{ℱ_{Edd}}{(1 + z)^{4}} = A σ_{SB} T_{Edd, \infty}^{4},$ $\begin{equation} F_{\mathrm{Edd}} = A\ \frac{\mathcal{F}_{\mathrm{Edd}}}{(1+z)^4} = A\ \sigma_{\mathrm{SB}} T_{\mathrm{Edd,\infty}}^4, \label{eq:Feddobs} \end{equation}$ (12)where A = (R_∞/D)² 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 plane-parallel 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 hydrogen-to-helium 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 ₁₀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: $C_{ℳ,i} = t_{𝒟} \int_{0}^{\infty} M (E) R (I,E) d E,$ $\begin{equation} C_{\mathcal{M},i} = t_{\mathcal{D}} \int_0^{\infty} M(E) R(I, E)~{\rm d}E, \label{eq:counts} \end{equation}$ (13)where M(E) is the photon number flux of the model at energy E and $\hbox{$t_{\mathcal{D}}$}$ 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 $R_{𝒟} (i,j) = \frac{\int_{E_{j - 1}}^{E_{j}} R (i,E) d E}{E_{j} - E_{j - 1}},$ $\begin{equation} R_{\mathcal{D}}(i,j) = \frac{\int_{E_{j-1}}^{E_j} R(i,E)~{\rm d}E}{E_j - E_{j-1}}, \end{equation}$ (14)for an energy range E_j−1 to E_j. In addition, we must take into account that the data might have a non-zero 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 $M (E) = \frac{F_{bkg} (E)}{E} + \frac{F_{obs} (E)}{E} \cdot$ $\begin{equation} M(E) = \frac{F_{\mathrm{bkg}}(E)}{E} + \frac{F_{\mathrm{obs}}(E)}{E}\cdot \end{equation}$ (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 power-law 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 non-zero 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 X-ray 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 $𝒢^{B} \equiv {B_{k}} \begin{matrix} N^{B} \\ k = 1 \end{matrix}$ $\hbox{$\clusterB \equiv \left\{ {\clusterBelem{k}} \right\}_{k=1}^{\NCB}$}$ . Each of these bursts B_k has $N_{k}^{S}$ $\hbox{$\NCS{k}$}$ 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 $𝒢_{k}^{S} \equiv {}^{{}{S_{jk}}^{}}_{j = 1}^{N_{k}^{S}}$ $\hbox{$\clusterS{k} \equiv \left\{ {\clusterSelem{j}{k}} \right\}_{j=1}^{\NCS{k}}$}$ . Each spectrum S_jk consists of a set $𝒢_{jk}^{C} \equiv {}^{{}{C_{𝒟, ijk}}^{}}_{i = 1}^{N_{jk}^{C}}$ $\hbox{$\clusterC{j}{k} \equiv \left\{ {\clusterCelem{i}{j}{k}} \right\}_{i=1}^{\NCC{j}{k}}$}$ with $N_{jk}^{C}$ $\hbox{$\NCC{j}{k}$}$ measurements of counts C_𝒟,ijk, measured in the detector channel i.

For a single channel we can define the likelihood function as $\hbox{$L(\mathcal{M})_{i} = P(\mathcal{D} | \mathcal{M}, \mathcal{H})$}$ so that $\hbox{$P(\mathcal{D} | \mathcal{M}, \mathcal{H})$}$ is the probability of the data $\hbox{$\mathcal{D}$}$ 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 goodness-of-fit for one element $\hbox{$C_{\mathcal{D},i}={\clusterCelem{i}{j}{k}}$}$ in some arbitrary burst B_k and spectrum S_jk, we compute the likelihood for a Poisson distributed data as $L_{i} = \frac{(C_{ℳ,i})^{C_{𝒟,i}} e^{- C_{ℳ,i}}}{(C_{𝒟,i})!} \cdot$ $\begin{equation} L_{i} = \frac{(C_{\mathcal{M},i} )^{C_{\mathcal{D},i}} ~e^{-C_{\mathcal{M},i}} }{(C_{\mathcal{D},i})!}\cdot \label{eq:ind_prob} \end{equation}$ (16)The joint likelihood for a single spectrum is then $L_{S} (ℳ) = P ({C_{𝒟,i}} \begin{matrix} i_{\max} \\ i = i_{\min} \end{matrix} | ℳ, ℋ) = \underset{i}{􏽙} L_{i},$ $\begin{equation} L_{\mathrm{S}}(\mathcal{M}) = P(\left\{ {C_{\mathcal{D},i}} \right\}_{i=i_{\mathrm{min}}}^{i_{\mathrm{max}}} | \mathcal{M}, \mathcal{H}) = \prod_{i} L_{i}, \end{equation}$ (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 $\hbox{$L_{\mathrm{B}} = \prod_{S} L_{\mathrm{S}}$}$ , and the total joint likelihood for all bursts is $\hbox{$L = \prod_{B} L_{\mathrm{B}}$}$ .

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), where $χ^{2} = \sum_{i} \frac{(C_{𝒟,i} - C_{ℳ,i})^{2}}{C_{ℳ,i}} \cdot$ $\begin{equation} \chi^2 = \sum_{i} \frac{( C_{\mathcal{D},i} - C_{\mathcal{M},i} )^2}{C_{\mathcal{M},i} }\cdot \label{eq:chi} \end{equation}$ (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 χ² 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 χ² values. We caution that particularly in the higher-energy channels, there may not be enough counts for the Poisson distribution to be well approximated by a Gaussian. In this case, χ² is not a good approximation to the log likelihood, and this could contribute to the formally poor χ² 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 χ²/ 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, $\hbox{$\mathcal{N}_{\sigma}(x)$}$ (where σ can be taken to be $\sqrt{C_{ℳ}}$ $\hbox{$\sqrt{C_{\mathcal{M}} }$}$ in relation to the χ² 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: $σ_{tot}^{2} = σ^{2} + σ_{int}^{2}$ $\hbox{$\sigma_{\rm tot}^2=\sigma^2+\intscat^2$}$ . The likelihood in Eq. (16)can then be replaced with $\begin{matrix} \begin{matrix} \begin{matrix} \end{matrix} L_{i} = 𝒩_{σ} * 𝒩_{σ_{int}} (C_{𝒟,i} - C_{ℳ,i}) \begin{matrix} \end{matrix} \\ \begin{matrix} \end{matrix} = \frac{1}{\sqrt{2 π (σ^{2} + σ_{int}^{2})}} \exp (- \frac{1}{2} \frac{(C_{𝒟,i} - C_{ℳ,i})^{2}}{σ^{2} + σ_{int}^{2}}), \end{matrix} \end{matrix}$ $\begin{eqnarray} \begin{split} L_i &= \mathcal{N}_{\sigma} \ast \mathcal{N}_{\intscat}(C_{\mathcal{D},i} - C_{\mathcal{M},i} ) \\ & = \frac{1}{\sqrt{2\pi(\sigma^2 + \intscat^2)}} \mathrm{exp}\left( - \frac{1}{2} \frac{(C_{\mathcal{D},i} - C_{\mathcal{M},i} )^2}{\sigma^2 + \intscat^2} \right), \end{split} \end{eqnarray}$ (19)where theasterisk marks the convolution of the two functions. Taking the logarithm simplifies the latter expression to $\ln L_{i} = const - \frac{1}{2} \ln (σ^{2} + σ_{int}^{2}) - \frac{1}{2} \frac{(C_{𝒟,i} - C_{ℳ,i})^{2}}{σ^{2} + σ_{int}^{2}} \cdot$ $\begin{equation} \ln L_i = \mathrm{const} - \frac{1}{2} \ln{(\sigma^2 + \intscat^2)} - \frac{1}{2} \frac{(C_{\mathcal{D},i}- C_{\mathcal{M},i})^2}{\sigma^2 + \intscat^2}\cdot \end{equation}$ (20)From this we see that the expression reduces to χ² when σ_int → 0 and $σ = \sqrt{C_{ℳ,i}}$ $\hbox{$\sigma = \sqrt{C_{\mathcal{M},i}}$}$ , 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 log-likelihood 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.

Table 1

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 $\hbox{$\clusterBelem{k} \in \clusterB$}$ for which we do not introduce any cluster-specific subparameters in this work. Proceeding in the hierarchy tree, each spectrum of the burst $S_{jk} \in 𝒢_{k}^{S} \in B_{k}$ $\hbox{$\clusterSelem{j}{k} \in \clusterS{k} \in \clusterBelem{k}$}$ 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_fA. 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 non-uniform phenomenological distance prior; this choice was motivated by the synthetic data results. Intrinsic scatter was always, when present, sampled as $\ln σ_{int}^{2}$ $\hbox{$\lognat \intscat^2$}$ , 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 non-varying 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 burst-specific 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 $\sqrt{D}$ $\hbox{$\sqrt{D}$}$ for the distance, which slightly favors higher values. This selection seems to remove the otherwise strong preference for lower masses¹. 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 spectrum-specific 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 affine-invariant 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 Metropolis-Hastings algorithm, but evolves not one, but many parallel sample values, called walkers, together. The random step for each walker is then performed using a so-called 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 Metropolis-Hastings 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 hard-state 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 signal-to-noise ratio level per time bin approximately constant despite the decreasing flux. We fixed the neutral hydrogen column depth to N_H = 1.87 × 10²² 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 best-fit 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 y-direction 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 two-dimensional 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². 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 ( $\sqrt{D}$ $\hbox{$\!\sqrt{D}$}$ ) 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 high-mass small-radius 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 two-dimensional posterior space projections. Figure 4 shows the M-R 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². 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 non-flat 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 (left-hand panels of both figures) and X = 0.73 (right-hand 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 low-mass solutions below the critical radius R = 4GM/c². 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 $M = 1. 2_{- 0.2 (0.2)}^{+ 0.2 (0.4)} M_{⊙}$ $\hbox{$M= 1.2^{+ 0.2 ~(0.4)}_{-0.2 ~(0.2)}~\Msun$}$ . Thus the input values were inside the 95% credible intervals. The distance is found to be $D = 5. 4_{- 0.4 (0.5)}^{+ 0.3 (0.7)} kpc$ $\hbox{$D= 5.4^{+ 0.3 ~(0.7)}_{-0.4 ~(0.5)}~\mathrm{kpc}$}$ . 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 $R = 12. 0_{- 0.9 (2.1)}^{+ 0.7 (1.1)} km$ $\hbox{$R= 12.0 ^{+ 0.7 ~( 1.1)}_{ -0.9 ~( 2.1)}~\mathrm{km}$}$ , $M = 2. 2_{- 0.7 (1.0)}^{+ 0.1 (0.1)} M_{⊙}$ $\hbox{$M= 2.2 ^{+ 0.1 ~( 0.1)}_{-0.7 ~( 1.0)}~\Msun$}$ , $D = 6. 9_{- 0.8 (1.6)}^{+ 0.3 (0.4)} kpc$ $\hbox{$D= 6.9^{+ 0.3 ~( 0.4)}_{-0.8 ~( 1.6)}~\mathrm{kpc}$}$ , and $X = 0.5 9_{- 0.04 (0.06)}^{+ 0.16 (0.33)}$ $\hbox{$X= 0.59^{+ 0.16 ~( 0.33)}_{-0.04 ~(0.06)}$}$ . Most importantly, we see that the correct radius is obtained without any bias. When the second high-mass cluster of solutions is omitted by asserting an additional M < 2 M_⊙ prior, we obtain $R = 12. 0_{- 0.9 (2.2)}^{+ 0.6 (1.0)} km$ $\hbox{$R= 12.0 ^{+ 0.6 ~( 1.0)}_{-0.9 ~(2.2)}~\mathrm{km}$}$ , $M = 1. 5_{- 0.3 (0.4)}^{+ 0.4 (0.4)} M_{⊙}$ $\hbox{$M= 1.5^{+ 0.4 ~( 0.4)}_{-0.3 ~(0.4)}~\Msun$}$ , $D = 6. 9_{- 1.0 (1.6)}^{+ 0.4 (0.4)} kpc$ $\hbox{$D= 6.9^{+ 0.4 ~(0.4)}_{-1.0 ~(1.6)}~\mathrm{kpc}$}$ , and $X = 0.5 9_{- 0.05 (0.07)}^{+ 0.15 (0.31)}$ $\hbox{$X= 0.59^{+ 0.15 ~( 0.31)}_{-0.05 ~( 0.07)}$}$ . 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 50025-01-01-00, 80033-01-01-08, 80033-01-19-04, 80033-01-20-02, and 80033-01-21-00, 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 $\sqrt{2}$ $\hbox{$\sqrt{2}$}$ , 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 X-ray burst from 4U 1702−429 (crosses) with corresponding best-fit 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 y-direction 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 best-fit model for the time-resolved 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.

Table 2

χ² values for the atmosphere model best-fits for model D.

As a final test of our atmosphere model goodness-of-fit, 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 best-fit model flux (from the model D run) together with a channel-specific Δχ value (defined as $\hbox{$(\mathcal{D} - \mathcal{M})/\sigma$}$ , that is, the difference between data $\hbox{$\mathcal{D}$}$ and model ℳ in units of standard error σ) for each of the five analyzed bursts. Additionally, Table 2 shows the sum of χ² values for the best-fit 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 χ² (for model D) is 3103.5 with 2438 degrees of freedom (reduced $χ_{r}^{2}$ $\hbox{$\chi^2_{\mathrm{r}}$}$ 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 low-count channels, and when we instead select 50 as our channel count cutof,f we obtain χ² = 2541.0 with 2122 degrees of freedom ( $χ_{r}^{2} = 1.20$ $\hbox{$\chi^2_{\mathrm{r}} = 1.20$}$ ).

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 two-dimensional 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: $R = 12. 4_{- 0.4 (2.6)}^{+ 0.3 (0.6)} km$ $\hbox{$R = 12.4^{+ 0.3 ~( 0.6)}_{-0.4 ~( 2.6)}~\mathrm{km}$}$ and $M = 1. 9_{- 0.3 (0.5)}^{+ 0.1 (0.3)} M_{⊙}$ $\hbox{$M= 1.9 ^{+ 0.1 ~( 0.3)}_{-0.3 ~( 0.5)}~\Msun$}$ . 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 $D = 5. 9_{- 0.3 (0.8)}^{+ 0.2 (0.3)} kpc$ $\hbox{$D= 5.9 ^{+ 0.2 ~( 0.3)}_{-0.3 ~( 0.8)}~\mathrm{kpc}$}$ with model D. We note that in the model D fits, the second non-physical high-mass small-radius 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 non-zero intrinsic scatter of $\ln σ_{int}^{2} = 3.2$ $\hbox{$\lognat \intscat^2 = 3.2$}$ . 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 X-ray 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 ( $\sqrt{D}$ $\hbox{$\!\sqrt{D}$}$ ) 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 χ², particularly for model comparisons (see also Andrae et al. 2010), we note that the model D fit to 4U 1702−429 has χ²/ d.o.f. = 2541.0/2122 (see also Table 2), whereas for a simple blackbody fit χ²/ 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 χ² is 175.5 smaller. This kind of comparison is, however, not strictly fair because individual blackbody best fits minimize the quoted χ² 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 best-fit values are also very susceptible to the exact energy range used for the fit, and so the χ² 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 χ² 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 non-zero 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., χ²), 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 non-uniform surface emission in the sense that the temperature varies across the surface. This is because our model C, which has a free surface-emitting 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 non-uniform surface temperature (Markwardt et al. 1999; Galloway et al. 2008; Ootes et al. 2017). However, these have been detected only during the soft-state 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 (non-burst) 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₁₀)² in the burst tail. Here R_bb is the blackbody radius (∝ R) and D₁₀ = 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⁴. The measured radius R scales roughly (for the fixed compactness) as $R \propto T_{Edd, \infty}^{-4} \propto K_{tail} / F_{Edd}$ $\hbox{$R \propto T_{\rm Edd,\infty}^{-4} \propto K_{\mathrm{tail}}/F_{\mathrm{Edd}}$}$ (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 low-energy 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 low-energy 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 Rayleigh-Jeans 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 X-ray 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 best-fit 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 energy-independent shift ζ affecting the observed spectra we still obtain $R \propto ζ K_{tail} / F_{Edd} \propto ζ K_{tail} / ζK T_{Edd, \infty}^{4} \propto 1 / T_{Edd, \infty}^{4}$ $\hbox{$R \propto \zeta K_{\mathrm{tail}}/F_{\mathrm{Edd}} \propto \zeta K_{\mathrm{tail}}/ \zeta K T_{\mathrm{Edd},\infty}^4 \propto 1/T_{\mathrm{Edd},\infty}^4$}$ , where T_Edd,∞ is the Eddington temperature (see Eq. (A9) in Poutanen et al. 2014 and Eq. (12)). This is a distance-independent 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 hard-state 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 banana-like 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² (see Özel & Psaltis 2015, for a discussion). If the constraints from F_Edd and A (∝ K_tail) are not consistent with each other², 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 soft-state 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² 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²/ 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 X-ray 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 soft-state 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 ultra-dense matter.

5. Summary

We have presented the first direct atmosphere model spectral fits to thermonuclear X-ray 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 hard-state type-I X-ray 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 highest-probability 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.

²

Together, F_Edd and A set the normalization of the model spectra because both depend on the distance. They both rely on the assumption that it is only the NS surface that is emitting. If this assumption is invalid, then the observed values might not coincide with the theoretical values.

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/48-1 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 PHY-1125915 for support. The computer resources of the Finnish IT Center for Science (CSC) and the FGCI project (Finland) are acknowledged.

References

Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, ArXiv e-prints [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., Cabrera-Lavers, 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ánchez-Ferná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, CD-ROMs (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 e-prints [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 X-ray 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: Wiley-Interscience) [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 X-ray 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

Table 1

Parameters of hierarchical fitting models.

In the text

Table 2

χ² values for the atmosphere model best-fits for model D.

In the text

All Figures

	Fig. 1 Synthetic spectra (crosses) with corresponding best-fit 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 y-direction 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 two-dimensional 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². 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 ( $\sqrt{D}$ $\hbox{$\!\sqrt{D}$}$ ) 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 high-mass small-radius 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 X-ray burst from 4U 1702−429 (crosses) with corresponding best-fit 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 y-direction 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 best-fit model for the time-resolved 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 ( $\sqrt{D}$ $\hbox{$\!\sqrt{D}$}$ ) 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, ArXiv e-prints [arXiv:1012.3754] [Google Scholar]

[2] 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]

[3] Ballantyne, D. R., & Everett, J. E. 2005, ApJ, 626, 364 [NASA ADS] [CrossRef] [Google Scholar]

[4] Ballantyne, D. R., & Strohmayer, T. E. 2004, ApJ, 602, L105 [NASA ADS] [CrossRef] [Google Scholar]

[5] Chen, Y.-P., Zhang, S., Torres, D. F., et al. 2011, A&A, 534, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[6] Cumming, A. 2003, ApJ, 595, 1077 [NASA ADS] [CrossRef] [Google Scholar]

[7] Damen, E., Magnier, E., Lewin, W. H. G., et al. 1990, A&A, 237, 103 [NASA ADS] [Google Scholar]

[8] Degenaar, N., Miller, J. M., Wijnands, R., Altamirano, D., & Fabian, A. C. 2013, ApJ, 767, L37 [NASA ADS] [CrossRef] [Google Scholar]

[9] Degenaar, N., Koljonen, K. I. I., Chakrabarty, D., et al. 2016, MNRAS, 456, 4256 [NASA ADS] [CrossRef] [Google Scholar]

[10] Ebisuzaki, T. 1987, PASJ, 39, 287 [NASA ADS] [Google Scholar]

[11] Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]

[12] Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]

[13] Guillot, S., Servillat, M., Webb, N. A., & Rutledge, R. E. 2013, ApJ, 772, 7 [Google Scholar]

[14] Güver, T., Özel, F., Cabrera-Lavers, A., & Wroblewski, P. 2010, ApJ, 712, 964 [NASA ADS] [CrossRef] [Google Scholar]

[15] Haberl, F., Stella, L., White, N. E., Gottwald, M., & Priedhorsky, W. C. 1987, ApJ, 314, 266 [NASA ADS] [CrossRef] [Google Scholar]

[16] 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]

[17] Humphrey, P. J., Liu, W., & Buote, D. A. 2009, ApJ, 693, 822 [NASA ADS] [CrossRef] [Google Scholar]

[18] Ibragimov, A. A., Suleimanov, V. F., Vikhlinin, A., & Sakhibullin, N. A. 2003, Astron. Rep., 47, 186 [NASA ADS] [CrossRef] [Google Scholar]

[19] in’t Zand, J. J. M., & Weinberg, N. N. 2010, A&A, 520, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[20] in’t Zand, J. J. M., Galloway, D. K., & Ballantyne, D. R. 2011, A&A, 525, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[21] Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401 [NASA ADS] [CrossRef] [Google Scholar]

[22] Ji, L., Zhang, S., Chen, Y., et al. 2015, ApJ, 806, 89 [NASA ADS] [CrossRef] [Google Scholar]

[23] Kajava, J. J. E., Nättilä, J., Latvala, O.-M., et al. 2014, MNRAS, 445, 4218 [NASA ADS] [CrossRef] [Google Scholar]

[24] Kajava, J. J. E., Koljonen, K. I. I., Nättilä, J., Suleimanov, V., & Poutanen, J. 2017a, MNRAS, 472, 78 [NASA ADS] [CrossRef] [Google Scholar]

[25] Kajava, J. J. E., Nättilä, J., Poutanen, J., et al. 2017b, MNRAS, 464, L6 [NASA ADS] [CrossRef] [Google Scholar]

[26] Kajava, J. J. E., Sánchez-Fernández, C., Kuulkers, E., & Poutanen, J. 2017c, A&A, 599, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[27] Koljonen, K. I. I., Kajava, J. J. E., & Kuulkers, E. 2016, ApJ, 829, 91 [NASA ADS] [CrossRef] [Google Scholar]

[28] Kurucz, R. 1993, CD-ROMs (Cambridge, Mass.: Smithsonian Astrophysical Observatory), No. 11 [Google Scholar]

[29] Kurucz, R. L. 1970, SAO Special Report, 309 [Google Scholar]

[30] Kuśmierek, K., Madej, J., & Kuulkers, E. 2011, MNRAS, 415, 3344 [NASA ADS] [CrossRef] [Google Scholar]

[31] 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]

[32] Lapidus, I. I., & Sunyaev, R. A. 1985, MNRAS, 217, 291 [NASA ADS] [CrossRef] [Google Scholar]

[33] Lattimer, J. M. 2012, Annu. Rev. Nucl. Part. Sci., 62, 485 [Google Scholar]

[34] Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109 [NASA ADS] [CrossRef] [Google Scholar]

[35] Lattimer, J. M., & Prakash, M. 2016, Phys. Rep., 621, 127 [NASA ADS] [CrossRef] [Google Scholar]

[36] Lattimer, J. M., & Steiner, A. W. 2014, Eur. Phys. J. A, 50, 40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[37] Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, Space Sci. Rev., 62, 223 [NASA ADS] [CrossRef] [Google Scholar]

[38] Malone, C. M., Nonaka, A., Almgren, A. S., Bell, J. B., & Zingale, M. 2011, ApJ, 728, 118 [NASA ADS] [CrossRef] [Google Scholar]

[39] Malone, C. M., Zingale, M., Nonaka, A., Almgren, A. S., & Bell, J. B. 2014, ApJ, 788, 115 [NASA ADS] [CrossRef] [Google Scholar]

[40] Markwardt, C. B., Strohmayer, T. E., & Swank, J. H. 1999, ApJ, 512, L125 [NASA ADS] [CrossRef] [Google Scholar]

[41] Medin, Z., von Steinkirch, M., Calder, A. C., et al. 2016, ApJ, 832, 102 [NASA ADS] [CrossRef] [Google Scholar]

[42] Miller, M. C. 2013, ArXiv e-prints [arXiv:1312.0029] [Google Scholar]

[43] Miller, M. C., & Lamb, F. K. 1996, ApJ, 470, 1033 [NASA ADS] [CrossRef] [Google Scholar]

[44] Miller, M. C., & Lamb, F. K. 2016, Eur. Phys. J. A, 52, 63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[45] Miller, M. C., Boutloukos, S., Lo, K. H., & Lamb, F. K. 2011, in Fast X-ray Timing and Spectroscopy at Extreme Count Rates, PoS (HTRS 2011)24 [Google Scholar]

[46] Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co.) [Google Scholar]

[47] Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ, 580, 1048 [NASA ADS] [CrossRef] [Google Scholar]

[48] Nättilä, J., & Pihajoki, P. 2017, A&A, submitted [arXiv:1709.07292] [Google Scholar]

[49] Nättilä, J., Suleimanov, V. F., Kajava, J. J. E., & Poutanen, J. 2015, A&A, 581, A83 [NASA ADS] [EDP Sciences] [Google Scholar]

[50] 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]

[51] Ootes, L. S., Watts, A. L., Galloway, D. K., & Wijnands, R. 2017, ApJ, 834, 21 [NASA ADS] [CrossRef] [Google Scholar]

[52] Özel, F. 2006, Nature, 441, 1115 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]

[53] Özel, F. 2013, Rept. Progr. Phys., 76, 016901 [NASA ADS] [CrossRef] [Google Scholar]

[54] Özel, F., & Psaltis, D. 2015, ApJ, 810, 135 [NASA ADS] [CrossRef] [Google Scholar]

[55] Özel, F., Güver, T., & Psaltis, D. 2009, ApJ, 693, 1775 [NASA ADS] [CrossRef] [Google Scholar]

[56] Özel, F., Psaltis, D., Güver, T., et al. 2016, ApJ, 820, 28 [NASA ADS] [CrossRef] [Google Scholar]

[57] Peille, P., Olive, J.-F., & Barret, D. 2014, A&A, 567, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[58] Poutanen, J. 2017, ApJ, 835, 119 [NASA ADS] [CrossRef] [Google Scholar]

[59] Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [NASA ADS] [CrossRef] [Google Scholar]

[60] Poutanen, J., Nättilä, J., Kajava, J. J. E., et al. 2014, MNRAS, 442, 3777 [NASA ADS] [CrossRef] [Google Scholar]

[61] Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley-Interscience) [Google Scholar]

[62] Serino, M., Mihara, T., Matsuoka, M., et al. 2012, PASJ, 64, 91 [NASA ADS] [Google Scholar]

[63] Shaposhnikov, N., Jahoda, K., Markwardt, C., Swank, J., & Strohmayer, T. 2012, ApJ, 757, 159 [NASA ADS] [CrossRef] [Google Scholar]

[64] Steiner, A. W. 2014a, Astrophysics Source Code Library [record ascl: 1408.020] [Google Scholar]

[65] Steiner, A. W. 2014b, Astrophysics Source Code Library [record ascl: 1408.019] [Google Scholar]

[66] Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2010, ApJ, 722, 33 [NASA ADS] [CrossRef] [Google Scholar]

[67] Strohmayer, T., & Bildsten, L. 2006, in Compact stellar X-ray sources, eds. W. Lewin, & M. van der Klis (Cambridge: Cambridge University Press), Cambridge Astrophys. Ser., 39, 113 [Google Scholar]

[68] Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]

[69] Suleimanov, V., & Werner, K. 2007, A&A, 466, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[70] Suleimanov, V., Poutanen, J., Revnivtsev, M., & Werner, K. 2011a, ApJ, 742, 122 [NASA ADS] [CrossRef] [Google Scholar]

[71] Suleimanov, V., Poutanen, J., & Werner, K. 2011b, A&A, 527, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[72] Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[73] Suleimanov, V. F., Poutanen, J., Klochkov, D., & Werner, K. 2016, Eur. Phys. J. A, 52, 20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[74] Suleimanov, V. F., Poutanen, J., Nättilä, J., et al. 2017, MNRAS, 466, 906 [NASA ADS] [CrossRef] [Google Scholar]

[75] van Paradijs, J., & Lewin, H. G. 1986, A&A, 157, L10 [NASA ADS] [Google Scholar]

[76] van Paradijs, J., Dotani, T., Tanaka, Y., & Tsuru, T. 1990, PASJ, 42, 633 [NASA ADS] [Google Scholar]

[77] Walker, M. A. 1992, ApJ, 385, 642 [NASA ADS] [CrossRef] [Google Scholar]

[78] Watts, A. L. 2012, ARA&A, 50, 609 [NASA ADS] [CrossRef] [Google Scholar]

[79] Weinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018 [NASA ADS] [CrossRef] [Google Scholar]

[80] Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94 [NASA ADS] [CrossRef] [Google Scholar]

[81] Worpel, H., Galloway, D. K., & Price, D. J. 2015, ApJ, 801, 60 [NASA ADS] [CrossRef] [Google Scholar]

[82] Yu, W., Li, T. P., Zhang, W., & Zhang, S. N. 1999, ApJ, 512, L35 [NASA ADS] [CrossRef] [Google Scholar]