The direct cooling tail method for X-ray burst analysis to constrain neutron star masses and radii

Abstract

Determining neutron star (NS) radii and masses can help to understand the properties of matter at supra-nuclear densities. Thermal emission during thermonuclear X-ray bursts from NSs in low-mass X-ray binaries provides a unique opportunity to study NS parameters, because of the high fluxes, large luminosity variations and the related changes in the spectral properties. The standard cooling tail method uses hot NS atmosphere models to convert the observed spectral evolution during cooling stages of X-ray bursts to the Eddington flux F_Edd and the stellar angular size Ω. These are then translated to the constraints on the NS mass M and radius R. Here we present the improved, direct cooling tail method that generalizes the standard approach. First, we adjust the cooling tail method to account for the bolometric correction to the flux. Then, we fit the observed dependence of the blackbody normalization on flux with a theoretical model directly on the M–R plane by interpolating theoretical dependences to a given gravity, hence ensuring only weakly informative priors for M and R instead of F_Edd and Ω. The direct cooling method is demonstrated using a photospheric radius expansion burst from SAX J1810.8–2609, which has happened when the system was in the hard state. Comparing to the standard cooling tail method, the confidence regions are shifted by 1σ towards larger radii, giving R = 11.5–13.0 km at M = 1.3–1.8 M_⊙ for this NS.

stars: neutron, X-rays: bursts, X-rays: individual: (SAX J1810.8–2609), X-rays: stars

1 INTRODUCTION

The internal structure of neutron stars (NSs) is not exactly known because of uncertainties in the matter equation of state in NS cores, where the matter density is larger than nuclear density (Haensel, Potekhin & Yakovlev 2007). There are many theoretical models of the equation of state, which predict various mass–radius dependences for NSs (Lattimer & Prakash 2007). Therefore, the possible dense matter equation of states could be constrained using the measurements of the NS radii R and masses M (see, e.g. Lattimer & Steiner 2014).

NS masses are accurately measured in close binary systems with radio pulsars. These works establish a relatively narrow range of NS masses, between 1.2 and 2 M_⊙ (see, e.g. Thorsett & Chakrabarty 1999; Özel et al. 2012; Antoniadis et al. 2013; Kiziltan et al. 2013). Within the next 20 yr, radius measurement with ∼5 per cent accuracy in such systems will be possible only for the double pulsar, PSR J0737-3039, via its moment of inertia (Kramer & Wex 2009). Spectral studies of the thermal emission from NSs provide a framework for NS radii determination. Although the observed X-ray spectra are thermal, a simple fit with the blackbody gives too small NS radii (see, e.g. Pavlov & Luna 2009; Klochkov et al. 2013). This is explained by the fact that the spectrum emergent from the gaseous NS atmospheres is significantly harder than the blackbody with equivalent effective temperature T_eff and the corresponding flux level in the Rayleigh–Jeans part of the spectrum is lower (see e.g. Zavlin, Pavlov & Shibanov 1996; Suleimanov & Werner 2007). Therefore, to achieve the same bolometric flux and to fit the observed spectra with the atmospheric models, a larger NS apparent radius is required (see e.g. Zavlin, Pavlov & Trumper 1998; Ho & Heinke 2009; Klochkov et al. 2013). The atmosphere models fits give the ratio of the apparent NS radius at infinity R_∞ = R(1 + z) (here z is the surface redshift) to the source distance D. The known distance to the NSs in globular clusters or in supernova remnants allows then to get R_∞ and constrain R (Heinke et al. 2006, 2014; Guillot et al. 2013; Klochkov et al. 2013, 2015). An important part of this method is the availability and reliability of the NS model atmospheres and their emergent spectra.

One of the possible ways to obtain simultaneous M–R constraints is via fast X-ray timing (Watts et al. 2016) using waveforms of oscillations during thermonuclear type I X-ray bursts (Lo et al. 2013; Miller & Lamb 2015) or of the persistent emission of accreting millisecond pulsars (Poutanen & Gierliński 2003). Alternatively, we can use spectral information from X-ray bursting NSs in low-mass X-ray binaries (see the reviews by Lewin, van Paradijs & Taam 1993; Miller 2013; Özel 2013; Miller & Lamb 2016; Suleimanov et al. 2016). Both the observed X-ray burst spectra (Galloway et al. 2008; Güver, Psaltis & Özel 2012; Worpel, Galloway & Price 2015) and the theoretical model atmosphere spectra are reasonably well fitted by a diluted blackbody (Suleimanov, Poutanen & Werner 2011a, 2012). The flux escaping the atmosphere can be represented as F_E ≈ πw B_E(T_c = f_cT_eff) (where T_c is the colour temperature and the dilution factor w is related to the colour correction factor f_c as |$w\approx f_{\rm c}^{-4}$|⁠). If f_c is known, we can estimate R_∞ from the measurements of the blackbody radius R_bb (for given distance) through |$R_\infty =R_{\rm bb}f_{\rm c}^2$|⁠. Of particular interest are the so-called photospheric radius expansion (PRE) bursts, whose luminosity exceeds the Eddington limit L_Edd. These bursts allow us to obtain an additional limitation on NS M and R either via the ratio of the fluxes at the PRE phase and the touchdown (when the photosphere touches the NS surface) or via estimation of the Eddington flux when the atmosphere reaches the Eddington limit (Ebisuzaki 1987; Damen et al. 1990; van Paradijs et al. 1990; Lewin et al. 1993).

Ebisuzaki (1987) proposed to fit the observed evolution of the colour temperature with luminosity in the cooling tail with his theoretical atmosphere models to find the Eddington temperature T_{Edd, ∞}, which is the redshifted effective temperature corresponding to the Eddington limit at the NS surface. Measurements of T_{Edd, ∞} gives a constraint on M and R (see also Suleimanov et al. 2011b). van Paradijs et al. (1990) fitted the evolution of the ratio of the observed colour temperature to F^1/4 (where F is the observed bolometric flux) by theoretical models of the evolution of f_c with luminosity (in units of the Eddington one) finding f_c at different points of the cooling tail and using them to determine R_∞. This method also allowed us to check the consistency of R_∞ measured from different points of the cooling tail.

More recently a simplified approach was taken by Özel (2006) and Özel, Güver & Psaltis (2009). They proposed the ‘touchdown method’, where the Eddington flux is identified with the touchdown flux and the average normalization of the blackbody in the burst cooling tail was used to estimate R_∞ with a fixed value of f_c corresponding to a low luminosity. For the analysis, they selected bursts which occur in the soft state of the source at high accretion rate and show minimum evolution of R_bb in the cooling tail. Such an approach was criticized for a number of reasons (see Steiner, Lattimer & Brown 2010; Suleimanov et al. 2011b; Miller 2013; Poutanen et al. 2014; Miller & Lamb 2016). First, it is not clear whether the Eddington limit is actually reached at the touchdown. Secondly, the estimation of the Eddington luminosity ignores the difference between actual Compton scattering cross-section and the Thomson one. Thirdly, assumption of the (nearly) constant f_c as well as the selection of the bursts based on the constancy of observed R_bb strongly contradicts the atmosphere models. Fourthly, the constraints on NS parameters are heavily dependent on the details of the assumed distance distribution. Furthermore, it was shown that the touchdown method applied to the bursts happening in the soft state with the assumed cuts on distance, assumption about chemical composition and f_c produce the results that have very small probability, i.e. the method is internally inconsistent. The most recent version of the touchdown method (Özel et al. 2016) accounts for the deviations of the cross-section from the Thomson one using the approximation from Suleimanov et al. (2012), but introduced other errors, such as considering systematic errors the same way as the statistical ones (see Miller & Lamb 2016).

Recently, we have suggested the cooling tail method for determination of NS masses and radii from the spectral evolution data during the cooling tail of the burst (Suleimanov et al. 2011b; see also appendix A of Poutanen et al. 2014). In this method, we fit the evolution of the observed blackbody normalization with flux by the theoretical model of |$f_{\rm c}^{-4}$| dependence on the relative luminosity ℓ = L/L_Edd. However, we relax the assumption often made that the Eddington limit is reached at the touchdown (Ebisuzaki 1987; van Paradijs et al. 1990; Özel 2006; Özel et al. 2016) and we use more data from the cooling tail to determine the Eddington flux. The method implicitly assumes that X-ray burst spectra as well as the model atmosphere spectra are well approximated by the blackbody emission and that the best-fitting blackbody flux well represents the bolometric flux of the source and the model.

For comparison with the data, we have computed an extensive set of hot NS atmosphere models (Suleimanov et al. 2011a, 2012) at a large grid of ℓ, surface gravity log g and chemical composition. Models with increased abundances of heavy elements, which mimic atmospheres polluted with nuclear burning ashes, were also computed (Nättilä et al. 2015). Similarly to the touchdown method, the fitting procedure provides two parameters: the observed Eddington flux F_Edd (defined for Thomson opacity) and the ratio R_∞/D. These parameters are then converted to the M and R constraints assuming some probability distribution for the distance. The cooling tail method was implemented to obtain NS mass–radius limitations for two X-ray bursting sources, 4U 1724–307 (Suleimanov et al. 2011b) and 4U 1608–52 (Poutanen et al. 2014) that obtained a lower limit of R > 13 km for the mass range 1.2 < M/M_⊙ < 2. However, transformation from a pair (F_Edd, R_∞/D) to (M, R) is not trouble-free: because of a divergence in the Jacobian, the resulting distribution of M and R can be biased away from the ‘2-Schwarzschild-radius’ line, R = 4GM/c², as was shown by Özel & Psaltis (2015). This problem can be easily avoided by performing a fit to the cooling tail data directly on a grid in the M–R plane by interpolating the computed theoretical f_c–ℓ relations to the corresponding value of log g. Such an approach within a Bayesian framework was adopted by Nättilä et al. (2016), allowing us to find mass–radius constraints for three bursters 4U 1702–429, 4U 1724–307 and SAX J1810.8–260.

Here we present the improved, direct cooling tail method that generalizes the standard approach, allows deviations of w from |$f_{\rm c}^{-4}$| and is free from the problems with the aforementioned transformation. In the new method, for every (M, R) pair, we first get the gravity log g and find the relation between model dilution factor w and the colour correction factor f_c on ℓ for the given log g by interpolating these quantities on a set of pre-computed atmosphere models covering an extended range of log g. We then fit this relation to the dependence of the blackbody normalization on blackbody flux, K − F_BB, observed in the cooling tail, using the distance as the only fitting parameter. The resulting χ² map constrains the NS M and R. We demonstrate the method using a PRE burst that occurred while SAX J1810.8–2609 was in the hard spectral state.

We note here that the reliable constraints on NS parameters can be obtained only from the bursts that show evolution of the blackbody with flux consistent with the predictions of the atmosphere models of passively cooling NS. It turned out that only some bursts happening during the low/hard state of the sources satisfy this requirement (Kajava et al. 2014; Poutanen et al. 2014). Interestingly, some long, very hard-state bursts show too strong spectral evolution that can be explained by metal enrichment of the atmosphere (Suleimanov et al. 2011b; Nättilä et al. 2016; Kajava et al. 2017), but making their usage for (M, R) determination difficult. On the other hand, none of the bursts happening in the soft, high-accretion rate state shows an evolution consistent with the atmosphere model predictions (Kajava et al. 2014). This is likely because of the strong influence of the accreting matter on the atmosphere (Suleimanov et al. 2011b; Poutanen et al. 2014; Koljonen, Kajava & Kuulkers 2016) and/or eclipsing of the NS by the accretion disc, making all constraints coming from the soft state bursts and still using the passively cooling atmosphere models (see e.g. Özel et al. 2016) questionable.

We also note that constraints on NS masses and radii can be obtained by direct fitting of the observed spectra with the atmosphere models (without going through the intermediate step of fitting the blackbodies to the data and the models). For example, Kuśmierek, Madej & Kuulkers (2011) fitted two separate spectra of a strong X-ray burst in the ultracompact binary 4U 1820−30 with their own helium-rich model atmosphere spectra and obtained rather weak constraints on M and R. They did not, however, use any information about the spectral evolution of the burst. On the other hand, we could fit simultaneously all the burst spectra getting improved (M, R) constraints, but this is a much more demanding exercise (J. Nättilä et al., in preparation).

The present paper is constructed as follows. In Section 2, we present the basic formulae, describe the standard cooling tail method and the correction to the method that accounts for the difference between the dilution factor w and |$f_{\rm c}^{-4}$|⁠. In Section 3, we present the new direct cooling tail method that is implemented on the M–R plane and does not suffer from the problems related to the divergence of the Jacobian. We then apply the new method to a PRE burst of SAX J1810.8–2609, show how the constraints coming from the direct cooling tail method differ from those obtained by the standard method and discuss possible systematic uncertainties. We conclude in Section 4.

2 COOLING TAIL METHOD

2.1 Basic relations

Let us consider a non-rotating NS of gravitational mass M and circumferential radius R, emitting luminosity L homogeneously over its surface. The NS atmosphere models implicitly assume that there is no any external influence to the NS surface that might affect either the emission or its propagation to the observer (i.e. no eclipse by the accretion disc, no boundary/spreading layer). In addition to the chemical composition, the NS emission depends on the surface gravity

\begin{equation} g = \frac{GM}{R^2}(1+z) \end{equation}

(1)

and the surface effective temperature given by

\begin{equation} \sigma _{\rm SB}\,T_{\rm eff}^4 = \frac{L}{4\pi R^2}. \end{equation}

(2)

Here z is the surface gravitational redshift given by relation

\begin{equation} 1+z = \left(1-\frac{R_{\rm S}}{R}\right)^{-1/2}, \end{equation}

(3)

σ_SB is the Stefan–Boltzmann constant and R_S = 2GM/c² is the Schwarzschild radius. It is also useful to introduce the critical Eddington luminosity and the corresponding effective temperature T_Edd:

\begin{equation} L_{\rm Edd} = \frac{4\pi \,GM\,c}{\kappa _{\rm T}}(1+z) = 4\pi R^2\,\sigma _{\rm SB}T_{\rm Edd}^4, \end{equation}

(4)

where κ_T = 0.2(1 + X) cm² g⁻¹ is the Thomson electron scattering opacity and X is the hydrogen mass fraction of the atmospheric plasma.

The observed values change due to the general relativistic effects

\begin{equation} R_{\infty } = R\,(1+z), \quad T_{\infty } = \frac{T_{\rm eff}}{1+z}, \quad L_\infty = \frac{L}{(1+z)^2}. \end{equation}

(5)

The observed Eddington luminosity and temperature are also reduced

\begin{eqnarray} L_{\rm Edd,\infty } &= &\frac{4\pi \,GM\,c}{\kappa _{\rm T}} \frac{1}{1+z}, \end{eqnarray}

(6)

\begin{eqnarray} T_{\rm Edd,\infty } &= &\left(\frac{g\,c}{\sigma _{\rm SB}\,\kappa _{\rm T}}\right)^{1/4} \frac{1}{1+z}. \end{eqnarray}

(7)

It is well known that the majority of the observed X-ray burst spectra can be fitted by the Planck function with two parameters, the observed colour temperature T_BB and normalization K (Galloway et al. 2008; Güver et al. 2012; Worpel et al. 2015):

\begin{equation} F_{\rm E} \approx \pi B_{\rm E}(T_{\rm BB})\,K = \pi B_{\rm E}(T_{\rm BB})\,\frac{R_{\rm BB}^2}{D^2}. \end{equation}

(8)

The observed bolometric flux in this approximation is

\begin{equation} F_{\rm BB} = \sigma _{\rm SB}\,T_{\rm BB}^4\,K. \end{equation}

(9)

On the other hand, the observed bolometric flux is determined as

\begin{equation} F_\infty = \sigma _{\rm SB}\,T_\infty ^4\, \Omega , \end{equation}

(10)

where

\begin{equation} \Omega = \frac{R_{\infty } ^2}{D^2} \end{equation}

(11)

is the angular dilution factor proportional to the solid angle occupied by the NS on the sky. Equations (9) and (10) can be combined to find a relation between the observed effective and colour temperatures, blackbody normalization and Ω:

\begin{equation} T_{\rm BB}^4\, K = T_\infty ^4 \,\Omega . \end{equation}

(12)

2.2 Standard cooling tail method

According to the NS atmosphere models (Suleimanov et al. 2011a, 2012), the emergent burst spectra are not real blackbodies, but they have blackbody-like shapes due to strong interactions between electrons and photons by Compton scattering in the surface atmospheric layers. The model spectra can be fitted in the observed X-ray energy range by a diluted blackbody:

\begin{equation} \mathcal {F}_{\rm E} \approx \pi \frac{1}{f_{\rm c}^4} B_{\rm E} (T_{\rm c}), \end{equation}

(13)

with the colour temperature that is higher than the model atmosphere effective temperature by a colour corrector factor f_c:

\begin{equation} T_{\rm c} = f_{\rm c}T_{\rm eff}. \end{equation}

(14)

Approximation (13) conserves the bolometric flux

\begin{equation} \int _0^\infty {f_{\rm c}^{-4}} \pi B_{\rm E}(T_{\rm c}) \,{\rm d}E = f_{\rm c}^{-4}\sigma _{\rm SB}\,T_{\rm c}^4 = \sigma _{\rm SB}\,T_{\rm eff}^4 \end{equation}

(15)

and allows us to connect the intrinsic NS effective temperature T_eff to the observed blackbody (colour) temperature:

\begin{equation} T_{\rm BB} = \frac{T_{\rm c}}{1+z} = \frac{f_{\rm c}T_{\rm eff}}{1+z} = f_{\rm c}T_\infty . \end{equation}

(16)

Combining this with equation (12), we get

\begin{equation} K = \frac{1}{f_{\rm c}^4} \Omega , \end{equation}

(17)

or

\begin{equation} K^{-1/4}= f_{\rm c}A = f_{\rm c}\left(\frac{R_\infty }{D}\right)^{-1/2}. \end{equation}

(18)

The computed f_c vary in the range 1.0–1.8, depending on the fundamental parameters of the NS atmosphere: T_eff (or the relative NS luminosity ℓ = L/L_Edd), the surface gravity log g and the chemical composition (Suleimanov et al. 2011a, 2012; Nättilä et al. 2015). Equation (17) implies that the observed normalization K has to depend on the observed bolometric flux F_BB exactly the same way as |$f_{\rm c}^{-4}$| depends on ℓ. This behaviour is expected if our assumption of a passively cooling NS during the cooling tail of X-ray burst is correct. Indeed, X-ray bursts happening in a low/hard persistent spectral state show the predicted K − F_BB behaviour (Suleimanov et al. 2011b; Kajava et al. 2014; Poutanen et al. 2014). Therefore, we can fit the observed dependence K^−1/4 − F_BB by the theoretical one f_c − ℓ, and obtain two fitting parameters A, defined through equation (18) and the observed Eddington flux

\begin{equation} F_{\rm Edd} = \frac{L_{\rm Edd, \infty }}{4\pi D^2}=\frac{GM\,c}{\kappa _{\rm T}\,D^2} \frac{1}{1+z}. \end{equation}

(19)

They can be combined to the observed Eddington temperature, which is independent of D

\begin{equation} T_{\rm Edd, \infty } = AF_{\rm Edd}^{1/4} = 9.81\,A^{\prime }F_{\rm Edd, -7}^{1/4}\,{\rm keV}, \end{equation}

(20)

where F_{Edd, -7} = F_Edd/10⁻⁷ erg s⁻¹ cm⁻² and the normalized fitting parameter A΄ = (R_∞[km]/D₁₀)^−1/2 with D₁₀ = D/10 kpc. The T_{Edd, ∞} found from observation determines a specific curve on the M–R plane, which allows the NS radius to be evaluated, as NS masses are in a finite range 1.2–2 M_⊙. This approach was successfully applied to some bursts from NSs in low-mass X-ray binaries 4U 1724–307 (Suleimanov et al. 2011b), 4U 1608–52 (Poutanen et al. 2014) and 4U 1702–429, 4U 1724–307 and SAX J1810.8–260 (Nättilä et al. 2016). We note that there is a difference between the NS radius in 4U 1724–307 obtained by Suleimanov et al. (2011b) and Nättilä et al. (2016) because they considered different hard-state bursts. The burst discussed by Suleimanov et al. (2011b) has likely a metal-enriched atmosphere (Nättilä et al. 2016), which was not accounted for by Suleimanov et al. (2011b).

2.3 Importance of the spectral dilution factor

The actual spectra from the NS atmosphere models have been fitted in the energy band 3–20 keV (corresponding to the largest effective area of the Proportional Counter Array onboard of the Rossi X-ray Timing Explorer) with the two-parameter model (Suleimanov et al. 2011a, 2012; Nättilä et al. 2015)

\begin{equation} \mathcal {F}_{\rm E} \approx w \pi B_{\rm E}(f_{\rm c}T_{\rm eff}), \end{equation}

(21)

allowing the spectral dilution factor w to be different from |$f_{\rm c}^{-4}$|⁠. The values |$wf_{\rm c}^4$| were also tabulated. For the standard method, we used the presentation of the normalization (17), because the values |$wf_{\rm c}^4$| are close to unity for almost all ℓ and the deviations become significant only at low luminosities, ℓ ≲ 0.1, which are unimportant for the cooling tail method (see Fig. 1).

Figure 1.

Dependence of the correction factor |$wf_{\rm c}^4$| on the relative luminosity ℓ for various log g and two chemical compositions.

Let us now adjust the cooling tail method accounting for the correct spectral dilution factor. We note that the two parameter fitting (21) does not conserve the surface bolometric flux

\begin{equation} \int _0^\infty w \pi B_{\rm E}(f_{\rm c}T_{\rm eff}) \,{\rm d}E = wf_{\rm c}^4\,\sigma _{\rm SB}T_{\rm eff}^4, \end{equation}

(22)

and we need a bolometric correction factor |$(wf_{\rm c}^4)^{-1}$| to convert the observed blackbody flux to the actual bolometric flux

\begin{equation} F_\infty = \left(wf_{\rm c}^4\right)^{-1}\, F_{\rm BB}. \end{equation}

(23)

Thus, we can relate the observed blackbody flux to the relative luminosity as

\begin{equation} F_{\rm BB} = \left(wf_{\rm c}^4\right)\ l \ F_{\rm Edd}. \end{equation}

(24)

Therefore, the relation between the observed blackbody normalization K and the NS radius, distance and the spectral fitting parameters also changes:

\begin{equation} \left(wf_{\rm c}^4\right)^{-1}\, T_{\rm BB}^4\, K = T_\infty ^4 \,\Omega , \end{equation}

(25)

or

\begin{equation} K = w\,\Omega . \end{equation}

(26)

Therefore, formally we have to fit the observed dependences K^−1/4 − F_BB with the computed relations |$w^{-1/4} - wf_{\rm c}^4\,\ell$|⁠, or simply just fitting the observed K − F_BB dependence by |$w - wf_{\rm c}^4\,\ell$| (see equations 24 and 26 and Fig. 2). From here on, in order to simplify the presentation, we use the latter dependence.

Figure 2.

Dependences of the dilution factor w on the model corrected relative luminosity |$wf_{\rm c}^4 \ell$| for fixed log g = 14.3 and various chemical compositions. The dashed curve corresponds to the solar H/He mix and reduced metal abundance (Z = 0.01 Z_⊙).

There are two fitting parameters for this procedure, F_Edd and Ω. They can be combined to get the observed Eddington temperature T_{Edd, ∞} given by equation (7):

\begin{equation} T_{\rm Edd, \infty } = 9.81\, \left( F_{\rm Edd, -7}/\Omega {[(km/10 kpc)}^{2}] \right) ^{1/4}\,{\rm keV}. \end{equation}

(27)

The curves of equal T_{Edd, ∞} on the M–R plane allow us to limit the NS mass and radius. The second parameter, Ω, could be useful if the distance is known, otherwise the distance can be estimated using its value.

To illustrate the method and demonstrate the differences between the new and the old approaches, we use the observational data of a PRE burst from SAX J1810.8–2609 analysed by Nättilä et al. (2016). We leave the analysis of other sources for the future. To allow an easy comparison with previous results, we select the data points after the touchdown down to the minimum flux F_min corresponding to 20 per cent of the touchdown flux F_td (see Fig. 3). To evaluate the effect of the data selection, we also considered F_min = 0.4F_td. The best-fitting solution was found with the simple Deming (2011) regression method by minimizing the function

\begin{equation} \chi ^2 \! = \! \sum _{i=1}^{N_{\rm obs}}\left[ \frac{(w\,\Omega - K_i)^2}{\Delta K_i^2}+ \frac{\left(wf_{\rm c}^4\,\ell \,F_{\rm Edd,-7} - F_{\rm BB,i}\right)^2}{\Delta F_{\rm BB,i}^2}\right], \end{equation}

(28)

which takes into account errors on both axes. The term in the square brackets is the square of the minimal dimensionless distance from the given observed point to the model curve |$w -wf_{\rm c}^4\,\ell$| with the given fitting parameters Ω and F_{Edd, -7}. The curve of constant T_{Edd, ∞} corresponding to the best-fitting parameters for log g =14.3 is shown in Fig. 4 by the solid red curve, while similar curve obtained with the standard cooling tail method (see table 1 in Nättilä et al. 2016) is shown with the dotted red curve. It is clear that the new method gives NS radii larger by about 0.5 km. The reason is clear from equations (24) and (26). The inverse bolometric correction |$wf_{\rm c}^4$| stretches the model curve, shifting it to higher fluxes at ℓ ≳ 0.9 and to lower fluxes for ℓ ≲ 0.9 (see Fig. 1). In order to fit the data, we would need to decrease F_Edd and at the same time to increase Ω. Both of these changes lead to smaller T_{Edd, ∞} and shift the corresponding solution to larger NS radii on the M–R plane.

Figure 3.

The observed K − F_BB dependence (circles with error bars) for the cooling tail (after touchdown) of the PRE burst of SAX J1810.8–2609 (see Nättilä et al. 2016). The red solid and the blue dashed curves are the best-fitting theoretical dependences |$w-wf_{\rm c}^4 \ell$| for solar H/He mix (Z = 0.01 Z_⊙) and log g = 14.3 corresponding to the data above F_min = 0.2F_td (black circles) and F_min = 0.4F_td (points to the right of the vertical dashed line), respectively. The best-fitting parameters are given in Table 1.

Figure 4.

Curves of constant T_{Edd, ∞} on the M–R plane, obtained from the data with F_min = 0.2F_td and three different log g =14.6 (dashed blue curve), 14.3 (solid red curve), and 14.0 (dot–dashed brown curve). The best-fitting parameters are given in Table 1. The red thin dotted curve corresponds to the best-fitting T_{Edd, ∞} from Nättilä et al. (2016) for the same data. The red thin dashed curve corresponds to F_min = 0.4F_td and log g = 14.3. The constant log g curves are shown by the thin dot–dashed lines.

The computed curves f_c − ℓ (as well as the curves |$w -wf_{\rm c}^4\,\ell$|⁠) are slightly different for different log g (Suleimanov et al. 2011a, 2012; Nättilä et al. 2015). Therefore, if we fit the observed curve K − F_BB with the theoretical dependences computed for different log g, we will obtain slightly different fitting parameters Ω and F_{Edd, -7}. As a result, the observed Eddington temperature T_{Edd, ∞} will be also different together with corresponding curves at the M–R plane (see Fig. 4 and Table 1). In fact, the constant T_{Edd, ∞} curves corresponding to different log g give correct values of NS M and R only at the crossing points with the corresponding constant log g curves. Fortunately, all constant T_{Edd, ∞} curves are close to each other within the statistical uncertainties.

Table 1.

Parameters of the fits of the K − F_BB dependence with the corrected cooling tail models for various log g.

log g	FEdd	Ω	TEdd, ∞	χ2/dof
	10−7 erg cm−2 s−1	(km/10 kpc)2	(keV)
Fmin  = 0.2Ftd
14.6	0.771	1222	1.555	58.6/33
14.3	0.776	1261	1.545	45.8/33
14.0	0.776	1315	1.529	44.5/33
Fmin  = 0.4Ftd
14.6	0.753	1293	1.524	30.6/17
14.3	0.765	1308	1.525	28.8/17
14.0	0.771	1338	1.520	27.8/17

We note here that the solution depends slightly on the number of data points used in the fits and the method to find the best-fitting solution. For example, if we take the points above F_min = 0.4F_td, the constant T_{Edd, ∞} curve corresponding to the best-fitting χ² solution shifts to larger NS radii (dashed red curve in Fig. 4 and Table 1).

3 DIRECT COOLING TAIL METHOD

Formally, the cooling tail method gives a correct result only at the crossing point of the curves of the constant log g and the corresponding T_{Edd, ∞}. Therefore, we suggest here a direct cooling tail method with M and R as fitting parameters. This method is free from the weakness of the standard method (see also Özel & Psaltis 2015, for more discussion) and allows us to obtain solutions at the R = 2R_S line as well. Moreover, the new formalism gives a possibility to take into account different priors in the physical parameters M and R and, for example, to control the minimum (or maximum) allowed mass in the fit.

For every pair M and R, we can compute the gravity log g, surface gravitational redshift z, and (assuming some chemical composition) the observed Eddington temperature T_{Edd, ∞} (see equations 1, 3 and 27). We use interpolation in the sets of f_c − ℓ and |$wf_{\rm c}^4 -\ell$| curves, computed for different values of log g, for finding the theoretical curve |$w - wf_{\rm c}^4\,\ell$| corresponding to a given log g.

3.1 Method of interpolation

To improve the accuracy of the interpolation, we have extended the existing set of hot NS model atmospheres for log g = 14.0, 14.3 and 14.6 (Suleimanov et al. 2012) to values from 13.7 up to 14.9 with the step 0.15 (i.e. nine values). We use the same relative luminosity set ℓ as in the previous work of Suleimanov et al. (2012), but for ℓ ≥ 0.1 only. Such computations were performed for four chemical compositions, pure hydrogen, pure helium and solar H/He mix with solar (Z = Z_⊙) and sub-solar (Z = 0.01 Z_⊙) metal abundances. The model emergent spectra were fitted with a diluted blackbody (see the details in Suleimanov et al. 2011a) and corresponding values of f_c, |$wf_{\rm c}^{4}$| and of the radiative acceleration g_rad were obtained.

Spectral fit parameters f_c and w depend on the relative luminosity ℓ in a similar way at different log g, when ℓ is small enough (ℓ < 0.8). But significant differences arise at larger ℓ because real electron scattering opacity decreases when the plasma temperature increases (see Suleimanov et al. 2012). As a result, the maximum possible luminosity relative to the Eddington luminosity for Thomson opacity, ℓ, also increases. Therefore, it is natural to use dependences of the type f_c − g_rad/g for interpolation, which are similar to each other for all surface gravities (see fig. 8 in Suleimanov et al. 2012). At the first step, we find g_rad/g values corresponding to the ℓ grid. Then we introduce the fixed grid of g_rad/g and interpolate all the necessary parameters (T_eff, f_c and w), to the new grid g_rad/g at every log g. At the final step, we get the theoretical dependence |$w - wf_{\rm c}^4\,\ell$| for the log g computed for the investigated M and R pair by interpolating all the quantities on the grid of nine values of log g. In this work, we use the interpolation along weighting backward and forward parabolas (see the detail in Kurucz 1970).

3.2 Fitting procedure

For the corrected cooling tail method, we can use two fitting parameters, Ω and F_{Edd, -7}. However, for a given M and R pair, L_Edd is already determined by their values and both fitting parameters depend on only one parameter – the distance to the source D. This allows us to obtain an estimation of the distance to the source given the chemical composition of the NS atmosphere and to limit the final solution if there exist independent constraints on the distance.

In the direct cooling tail method, we minimize χ² at a grid in the M–R plane by fitting the observed K − F_BB curve with the theoretical |$w-wf_{\rm c}^4\,\ell$| curve interpolated to the current log g using D₁₀ as the only fitting parameter. We consider a grid of masses from 1 to 3 M_⊙ with the step 0.01 M_⊙ and a grid of radii from 9 to 17 km with the step 0.01 km. We ignore the pairs that do not satisfy the causality condition R > 1.45R_S (Lattimer & Prakash 2007). We thus obtain the map of best-fitting D and the χ² values at the M–R plane and the overall minimum of χ², which allows us to find the confidence regions (e.g. for two parameters, Δχ² = 2.3, 4.61, 9.21 give probabilities of 68, 90 and 99 per cent). We note here that the method used for finding the best-fitting solution, e.g. the χ²-method like here and in Suleimanov et al. (2011b), or the robust likelihood method from Poutanen et al. (2014), affects only slightly (within the 1σ error) the centroid of the solution. On the other hand, the resulting error contours on the M–R plane are affected more: the χ²-method generally underestimates the errors, because of the presence of the significant scattering in the data as well as because the best-fitting solution may be located at the boundary of the prior distributions. The robust estimator gives more correct presentation of the errors when benchmarked against a full Bayesian fit with intrinsic scatter in the system. For simplicity, we, however, show here only the results obtained by minimizing a non-robust χ² function given by equation (28).

We apply the direct cooling tail method to the data described above and presented in Fig. 3. The minimum χ² is 45.4 for 32 dof. The confidence regions at the M–R plane are shown in Fig. 5. To compare with the results presented in Fig. 4, we also draw the curve of equal T_{Edd, ∞} obtained using the corrected cooling tail method for log g = 14.3 (red solid curve). There is good correspondence between this curve and the confidence regions given by the direct method. To compare with the previous result, the curve of equal T_{Edd, ∞} obtained with the standard cooling tail method (table 2 in Nättilä et al. 2016) is also shown (red dotted curve). We remind here that the method does not suffer from the problems with the Jacobian transformation from F_Edd, Ω to M, R (see the discussion in Özel & Psaltis 2015), because we fit directly on the M–R plane. We also see that instead of the two regions of allowed solutions (small R, large M and large R, smaller M; see e.g. Suleimanov et al. 2011b; Poutanen et al. 2014), the solution here prefers larger R and smaller M values. This is a direct consequence of the fact that the shape of the |$w-wf_{\rm c}^4\,\ell$| curves depends on log g and a smaller log g gives a better description of the data (see Table 1).

Figure 5.

The χ² confidence regions (68, 90 and 99 per cent probabilities) in the mass–radius plane for SAX J1810.8–2609 obtained using the direct cooling tail method. The data above F_min = 0.2F_td (see Fig. 3) are used. The red solid and dashed curves correspond to the best-fitting T_{Edd, ∞} obtained for log g = 14.3 for F_min = 0.2 and 0.4F_td, respectively, while the red dotted curve gives the results from Nättilä et al. (2016). The black dot–dashed curves correspond to the constant distance of 4.0, 4.5 and 5.0 kpc.

We conclude that the most probable radius of the NS in SAX J1810.8–2609 lies in the range between 11.5 and 13.0 km for the mass range of 1.3–1.8 M_⊙. However, the smaller radii of ∼10 km for a high-mass (∼2 M_⊙) NS cannot be formally excluded. The obtained constraints also limit the possible distance to the source that likely lies in the interval between 3.7 and 4.7 kpc.

4 SUMMARY

We have presented the direct cooling tail method that takes NS mass and radius together with the distance as the direct fitting parameters. This method implies fitting of an observed dependence between blackbody fitting parameters, the normalization and the observed blackbody flux, K − F_BB, with the corresponding theoretical dependence, |$w - wf_{\rm c}^4\,\ell$|⁠, which is found for every given M and R pair by interpolation in the extended grid (for nine values of log g) of the computed dependences f_c − ℓ and |$wf_{\rm c}^4 - \ell$|⁠. Such grids were computed for four chemical compositions: pure H and pure He, and solar H/He mix with the solar (Z = Z_⊙) and sub-solar (Z = 0.01 Z_⊙) metal abundances. The fitting procedure results in the χ² maps and corresponding confidence regions in the M–R plane. Other fitting methods (robust likelihood, Bayesian) give very similar results with slightly larger errors. Compared with the standard cooling tail method, the new method does not suffer from the problems with the Jacobian transformation from (F_Edd, Ω) to (M, R) allowing us to obtain solution anywhere on the M–R plane.

We applied the new method to a PRE burst that occurred in the low/hard state of SAX J1810.8–2609. The direct cooling tail method gives the Eddington temperature slightly smaller than the standard cooling tail method, resulting in a systematic shift of the radii by ∼1 σ to larger values (which is still within the statistical uncertainties). However, instead of the two separated (M, R) regions usually obtained with the standard method, our solution prefers larger R and smaller M values corresponding to smaller log g values, which is a direct consequence of the dependence of the theoretical cooling curves |$w - wf_{\rm c}^4\,\ell$| on gravity. We finally constrain the NS radius in SAX J1810.8–2609 to lie between 11.5 and 13.0 km for the assumed mass range of 1.3–1.8 M_⊙. The distance to the source likely lies between 3.7 and 4.7 kpc.

Acknowledgments

The work was supported by the German Research Foundation (DFG) grant WE 1312/48-1, the Magnus Ehrnrooth Foundation, the Russian Foundation for Basic Research grant 16-02-01145-a (VFS), the Foundations’ Professor Pool (JP), the Finnish Cultural Foundation (JP), the Academy of Finland grant 268740 (JP, JJEK), the University of Turku Graduate School in Physical and Chemical Sciences (JN), the Faculty of the European Space Astronomy Centre (JN, JJEK), the European Space Agency research fellowship programm (JJEK), and the Russian Science Foundation grant 14-12-01287 (MGR). We also acknowledge the support from the International Space Science Institute (Bern, Switzerland) and COST Action MP1304 NewCompStar.