CONSTRAINTS ON NEUTRON STAR MASS AND RADIUS IN GS 1826−24 FROM SUB-EDDINGTON X-RAY BURSTS

Thermonuclear flashes from accreting neutron stars, observed as type I X-ray bursts, in principle provide a way to constrain neutron star masses and radii (for a review see Section 4 of Lewin et al. 1993). The large observational catalogs of type I X-ray bursts now available (Galloway et al. 2008a, hereafter G08) and new spectral models (Madej et al. 2004; Majczyna et al. 2005; Suleimanov et al. 2011b) have motivated fresh attempts to do this using photospheric radius expansion (PRE) bursts (Özel 2006; Özel et al. 2009, 2011; Güver et al. 2010a, 2010b; Steiner et al. 2010; Suleimanov et al. 2011a, 2011b). In this approach, the peak luminosity of the burst (specifically at the "touchdown" point when the photosphere returns to the neutron star surface) is related to the Eddington luminosity, and the normalization of the burst spectrum is related to the emitting area. If some information about the distance to the source is available, constraints on the neutron star mass and radius can be derived. These works have highlighted and spurred debate about some of the systematic errors that must be taken into account, such as uncertainty in identifying the moment at which the photosphere touches down (Galloway et al. 2008b; Steiner et al. 2010; Suleimanov et al. 2011a, 2011b; Güver et al. 2012a), and differences in derived radii when using bursts at different accretion rates from the same source (Suleimanov et al. 2011a, 2011b; Güver et al. 2012b).

GS 1826−24 is a unique X-ray burster that shows remarkable agreement with theoretical models of recurrence times, energetics, and light curves (Galloway et al. 2004; Heger et al. 2007; in't Zand et al. 2009). The observed recurrence time (typically 3–5 hr) in a given epoch is the same from burst to burst to within a few minutes (Cocchi et al. 2001), and the burst light curves in a given epoch are very uniform (Galloway et al. 2004), implying the same conditions on the neutron star surface from burst to burst, and a regular limit cycle. Heger et al. (2007) compared the observed light curves to the theoretical models of Woosley et al. (2004). By choosing a model with approximately the same recurrence time as the data, and by varying the distance only (which scales the observed peak flux to match the model peak luminosity), the theoretical light curve fit most of the observed light curve well, except for deviations during the burst rise and at late times in the cooling tail. This is of great interest because the long ∼100 s tails of these bursts are powered by the rp-process (Wallace & Woosley 1981), and offer a way to test the nuclear physics input, such as masses of proton-rich heavy nuclei and their reaction rates, both of which are usually highly uncertain (e.g., Schatz et al. 1998; Schatz 2006). Even more remarkably, in't Zand et al. (2009) extracted the observed light curve out to more than 1000 s by combining multiple bursts. The late time cooling observed matched the late time cooling in the theoretical model of Heger et al. (2007), which arises from heat initially conducted inward to deeper layers that then emerges on long timescales.

In this paper, we take the comparison between observations and theory for GS 1826−24 one step further. We first determine the constraints on neutron star mass and radius that can be derived from the light curve comparison carried out by Heger et al. (2007). The basic idea here is that even though the X-ray bursts from GS 1826−24 do not reach the Eddington limit, we can still determine the intrinsic luminosity of the bursts by comparing with the light curve models. This replaces the touchdown measurement used for PRE bursts with a different condition, which we use to constrain M and R for GS 1826−24. In our analysis, we take care to include the possible anisotropy in burst and persistent emission and show how that could affect the mass and radius determination. We then compare the spectral evolution during the tail of the burst with spectral models. The good understanding of bursts from GS 1826−24 suggests that they could be a good testing ground for spectral models. We show that in the initial cooling phase following peak luminosity the spectral evolution agrees well with the models of Suleimanov et al. (2011b), and we derive the associated constraints on M and R. In both cases, we look for constraints that are independent of distance and emission anisotropies since neither are well constrained for GS 1826−24.

The outline of the paper is as follows. The data analysis is described in Section 2. In Section 3, we discuss the possible anisotropy of the burst and persistent emission and review calculations of the expected degree of anisotropy in the literature. In Section 4, we use the model light curve from Heger et al. (2007) to set the luminosity scale of the observed bursts and show that this gives a distance-independent relation between the redshift and color-correction factor f_c. Suleimanov et al. (2011a) argued that rather than using a single measurement of touchdown flux, the entire cooling track of the burst should be fit to spectral models. We do this in Section 5, and show that even though the peak flux is below Eddington, the fits provide a constraint on the value of F_Edd as well as the normalization of the spectrum. These two measurements translate into a distance-independent upper limit on R_∞. We compare these two different constraints and discuss their implications in Section 6.

We used data taken with the Proportional Counter Array (PCA; Jahoda et al. 1996) on board the Rossi X-ray Timing Explorer from the catalog of bursts detected over the mission lifetime (G08). Where not explicitly stated, the data analysis procedures are as in G08. Time-resolved spectra in the range 2–60 keV covering the burst duration were extracted on intervals as short as 0.25 s during the burst rise and peak, with the bin size increasing stepwise into the burst tail to maintain roughly the same signal-to-noise level. A spectrum taken from a 16 s interval prior to the burst was adopted as the background.

We re-fit the spectra over the energy range 2.5–20 keV using the revised PCA response matrices, v11.7,³ and adopted the recommended systematic error of 0.5%. The fitting was undertaken using XSPEC version 12. In order to accommodate spectral bins with low count rates, we adopted Churazov weighting. We modeled the effects of interstellar absorption using a multiplicative model component (wabs in XSpec), with the column density N_H frozen at 4 × 10²¹ cm⁻² (e.g., in't Zand et al. 1999). In the original analysis carried out by G08, the neutral absorption was determined separately for each burst from the mean value obtained for spectral fits carried out with the N_H value free to vary. This has a negligible effect on the fluxes, but can introduce spurious burst-to-burst variations in the blackbody normalization.

The burst data used here have been corrected for "dead time," a short period of inactivity in the detectors following the detection of an X-ray photon. There are, however, concerns regarding the absolute flux calibration of the PCA associated with variations in the flux from the Crab nebula and the effective area of the PCA. We will show that such absolute uncertainties will not influence our derived constraints.

The possibility that the burst or persistent emission is not isotropic has been long discussed (e.g., Lapidus et al. 1985), but has not always been included in recent work using X-ray bursts to constrain neutron star mass and radius. For example, in Özel (2006), Steiner et al. (2010), and Suleimanov et al. (2011a) the burst emission is assumed to be isotropic. Here, we review the expected size of the anisotropy. We follow Fujimoto (1988) and define an anisotropy parameter ξ by the relation 4πd²Fξ = L between the observed flux F and the luminosity of the source L over the whole sky, where d is the distance to the source. When ξ < 1 (>1), the radiation is beamed toward (away from) the observer. We write the anisotropy factor for the burst and persistent emission as ξ_b and ξ_p, respectively.

Lapidus et al. (1985) showed that if the accretion disk extends to the neutron star surface during the flash, it will intercept ≈1/4 of the radiation from the burst, reflecting it preferentially along the disk axis. They provide the approximate expression

$$\begin{equation} \xi _{\rm b}^{-1}={1\over 2}+\left|\cos i\right|, \end{equation} \tag{ 1 }$$

where i is the inclination angle (i = 0° means the system is viewed face on, looking down the disk axis), which closely fits their more detailed results derived from solving the radiative transfer equations for a disk geometry (they found a maximum value of 1.39 rather than 1.5). The range of ξ⁻¹_b is from 0.5 (edge on) to 1.5 (face on), implying an uncertainty of a factor of three depending on inclination angle.

The anisotropy factor for the persistent emission is perhaps even more uncertain than that for the burst flux, depending on the specific model for the inner accretion disk, boundary layer, and corona, etc. Lapidus et al. (1985) and Fujimoto (1988) derive opposite behaviors for the factor ξ_p as a function of inclination. The model presented in Fujimoto (1988) for the persistent emission assumes that radiation from the boundary layer, which encircles the neutron star in a "belt" about its equator, is largely screened by the inflated inner part of the accretion disk and scattered preferentially in a direction along the disk axis. In Lapidus et al. (1985), however, the inner part of the disk is assumed to be thin, and less than one half of the boundary layer radiation falls on the accretion disk and is re-scattered, again preferentially along the disk axis, while the remainder of the emission is beamed preferentially in the direction i = 90° (along the plane of the disk). This difference in their modeling of the inner accretion disk is made apparent by the fact that while Fujimoto (1988) predicts no radiation to be emitted in the i = 90° direction, Lapidus et al. (1985) find a substantial portion of the persistent emission will be beamed in that direction. The ratio ξ_p/ξ_b varies by up to a factor of ∼3 with inclination for both models, although while Fujimoto (1988) finds that the ratio is monotonically increasing with inclination, Lapidus et al. (1985) find the opposite trend. We note that these two models do not consider the effects of general relativity on the trajectories of photons near the neutron star surface. However, Lapidus et al. (1985) show that when considering this effect, a substantially larger proportion (∼28% for a ratio of the neutron star radius to the Schwarzschild radius R/r_s of 3) of radiation falls on the accretion disk, further enhancing beaming along the disk axis.

The definition of ξ is such that it always appears with distance d in the combination ξ^1/2d. Therefore, the uncertainty in the anisotropy factor acts in the same way as an additional uncertainty in the distance to the source. For GS 1826−24, Homer et al. (1998) suggest a limit i < 70° based on the low amplitude of the optical modulation at the orbital frequency, in which case Equation (1) gives 0.84 ≲ ξ⁻¹_b ≲ 1.5, or 0.85 < ξ^1/2_b < 1.1. Therefore, even if the distance to GS 1826−24 was perfectly known, anisotropy of the burst emission would represent an uncertainty of about ±15% in any quantity that depends on distance. For example, using spectral fits to determine R_∞ is subject to this uncertainty since the normalization of the spectrum depends on the solid angle R²_∞/d²ξ. Given these uncertainties, in this paper we look for constraints on M and R that are independent of distance and anisotropy.

We first ask what constraints on M and R arise from the comparison between the observed light curve and theoretical models of Heger et al. (2007). PRE bursts are often used in work to constrain neutron star properties from X-ray bursts because the peak luminosity of the burst can then be taken to be the Eddington luminosity. This cannot be done for GS 1826−24 because the bursts do not show PRE, implying that they have a peak luminosity below Eddington. Instead, here we pursue the idea that the model light curves that fit the observed light curves so well tell us the peak luminosity of the bursts.

Heger et al. (2007) selected from their models one that had a similar recurrence time to the observed bursts in 2000 (the model had t_recur = 3.9 hr as opposed to the observed t_recur = 4.07 hr). They showed that, when the distance to GS 1826−24 (actually ξ^1/2_bd) is chosen to make the predicted peak flux match the observed light curve, the theoretical and observed burst light curves show remarkable agreement. They considered fixed values of M and R, but the choice of those two parameters also changes the mapping between observed and model burst fluxes. Therefore, rather than vary distance alone, we find the value of the ratio of the observed flux to the model flux

$$\begin{equation} {F_{\rm obs}\over F_{\rm model}}=\xi _{\rm b}^{-1}\left({R\over d}\right)^2 {1\over \left(1+z\right)^2} \end{equation} \tag{ 2 }$$

that gives the best fit between the model and the data. Taking the peak values, F_{model, pk} = 1.29 × 10²⁵ erg cm⁻² s⁻¹ (computed from the redshifted peak luminosity quoted in Heger et al. 2007, with R = 11.2 km and z = 0.26) and F_{obs, pk} = 2.84 × 10⁻⁸ erg cm⁻² s⁻¹, we find F_obs/F_model = 2.20 × 10⁻³³. Substituting this value into Equation (2), we find

$$\begin{equation} {R\over \xi _{\rm b}^{1/2}d}=10.9\ {\rm km/6.0\ {\rm kpc}}\ {1+z\over 1.26} \left({F_{\rm obs}/F_{\rm model}\over 2.2\times 10^{-33}}\right)^{1/2}, \end{equation} \tag{ 3 }$$

where we use the redshift assumed by Heger et al. (2007). Note that the model light curve is likely to be insensitive to the model gravity, so that the ratio F_obs/F_model does not depend sensitively on the M and R used in the model. For example, the ignition column depth is weakly dependent on gravity in this burning regime (Bildsten 1998 derives y_ign∝g^−2/9). However, this is something that should be explored in further simulations. For now, we assume that F_obs/F_model is a constant, and take Equation (3) as a joint constraint on R and 1 + z.

The theoretical uncertainty in F_model is at present unknown. The predicted light curves depend on the input nuclear physics, and prescription for convection and other mixing processes, for example. These prescriptions vary from code to code, and currently only simulations from the Kepler code (Woosley et al. 2004) have been compared to the observations of GS 1826−24. Further simulations and comparisons are required to determine what range of predicted peak fluxes still produces light curves with the correct shape to fit the data. For now, in order to put an error bar on the prefactor in Equation (3), we assume that the theoretical uncertainty in F_obs/F_model is ±10%, and keep in mind the fact that this number is uncertain.

This raises the point that rather than use the peak flux only, we could also fit the entire light curve. In that case there is an extra parameter, the redshift 1 + z, which stretches the light curve in time. In principle, this provides a constraint on 1 + z. In practice, however, we find that the value of 1 + z obtained in the fit is sensitive to how much of the light curve is included in the fit. For example, fitting the entire light curve (until about 130 s after the peak) we find best-fit values 1 + z = 1.44, F_obs/F_model = 2.10 × 10⁻³³. If we fit the first 30 s only, which includes only the initial decline after the peak rather than the whole tail, we get a best fit of 1 + z = 1.32 and F_obs/F_model = 2.17 × 10⁻³³. We show in Figure 1 the separate fits to the entire light curve and the first 30 s, and we also include the model light curve fitted only by matching the peak fluxes, with the value for the redshift of 1 + z = 1.26, as assumed by Heger et al. (2007). We see that while the redshift is sensitive to the details of the fitting, the normalization F_obs/F_model is well determined. Therefore, here we use the normalization, but leave fits to the shape of the entire light curve to future work when a greater number of simulations are available.

Zoom In Zoom Out Reset image size

A second constraint comes from spectral fitting. Fitting the observed burst spectrum with a blackbody gives the total flux F_∞, color temperature T_{c, ∞}, and the blackbody normalization:

$$\begin{equation} K=\left({F_\infty \over \sigma T_{c,\infty }^4}\right)={R^2_\infty \over d^2 f_{\rm c}^4}\xi _{\rm b}^{-1}. \end{equation} \tag{ 4 }$$

The color-correction factor f_c = T_c/T_eff takes into account the hardening of the burst spectrum compared to a blackbody at the same effective temperature T_eff. We plot K as a function of time in Figure 2 for two average burst profiles for recurrence times 5.74 and 4.07 hr, respectively. Following the burst rise, which lasts for approximately 5 s, the normalization levels off until ≈60 s into the burst, when the normalization drops dramatically over 100 s to only about 25% of its original value. We discuss the variation of K in the tail, and the difference in K for the two different recurrence times in the next section. For our purposes here we find the mean value of K for the t_recur = 4.07 hr profile during the period following the peak where it is constant, giving K = 110 ± 2 (km/10 kpc)². If we take into account the dead-time correction near the burst peak (≈6%), the value we find is consistent with the more detailed analysis of blackbody normalization in these bursts carried out by Galloway & Lampe (2012).

Zoom In Zoom Out Reset image size

Dividing Equations (2) and (4), R/d and ξ_b drop out, giving

$$\begin{eqnarray} {f_{\rm c}\over 1+z} &=& K^{-1/4} \left({F_{\rm obs}\over F_{\rm model}}\right)^{1/4} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &=&1.17\ \left({K\over 110\ ({\rm km}/10\ {\rm kpc})^2}\right)^{-1/4}\left({F_{\rm obs}/F_{\rm model}\over 2.2\times 10^{-33}}\right)^{1/4}. \end{eqnarray} \tag{ 6 }$$

We show model calculations of f_c from Suleimanov et al. (2011b) in the next section that typically have f_c ≈ 1.4–1.5 during the phase where K is relatively constant. For f_c = 1.4 and 1.5 we get 1 + z = 1.19 and 1.28, respectively. Note that as well as being independent of d and ξ_b, the value of 1 + z determined in this way is not very sensitive to the values of K and F_model/F_obs (proportional to the 1/4 power of each). For example, introducing an uncertainty in F_model/F_obs of ±10% gives a prefactor in Equation (9) of 1.17 ± 0.03 or 1 + z = 1.28 ± 0.03 (f_c/1.5).

There is one more constraint, which comes from the agreement between the measured and model recurrence times, which effectively measures the local-mass accretion rate $\dot{m}$ onto the star. We define $\dot{m}$ to be the rest-mass accretion rate at the stellar surface. Then the accretion flux as observed at infinity is

$$\begin{equation} F_X=\dot{m} \left({R\over d}\right)^2 \xi _p^{-1} {c^2z\over (1+z)^2}. \end{equation} \tag{ 7 }$$

Dividing Equations (2) and (7) gives the observed quantity

$$\begin{equation} {F_X\over \dot{m} c^2}{F_{\rm model}\over F_{\rm obs}} = \left({\xi _{\rm b}\over \xi _p}\right) z, \end{equation} \tag{ 8 }$$

a direct measure of redshift, independent of distance, but dependent on the anisotropy parameter ratio ξ_p/ξ_b. The accretion rate in model A3 of Heger et al. (2007) was $\dot{m}=7980\ {\rm g\ cm^{-2}\ s^{-1}}$, and the measured persistent flux in 2000 was F_X = 2.91 ± 0.03 × 10⁻⁹ erg cm⁻² s⁻¹, giving

$$\begin{equation} z = 0.18\pm 0.02\ \left({\xi _p\over \xi _{\rm b}}\right). \end{equation} \tag{ 9 }$$

An alternative way to derive this result is to match the theoretical α value, the ratio of persistent fluence between bursts to burst fluence, to the observed value. The models of Heger et al. (2007) with recurrence time of 4 hr have a theoretical value α_model = ΔMc²z_model/E_burst(ξ_b/ξ_p) = 55 (ξ_b/ξ_p), where E_nuc is the burst energy, ΔM the ignition mass, and a redshift z_model = 0.26. The observed α at the same recurrence time is α_obs ≈ 37 (Figure 2 of Heger et al. 2007), giving z = z_model(α_obs/α_model) = 0.17(ξ_p/ξ_b), in agreement with the value in Equation (9).

If the anisotropy parameters were known, Equations (9), (2), and (4) uniquely determine the three quantities 1 + z, f_c, and R/d, respectively. However, as noted in Section 2, the anisotropy parameters are not well constrained. In particular, Equation (9) for the redshift is not constraining once the large uncertainty in ξ_p/ξ_b is taken into account. The main result of this section is therefore the relation between 1 + z and f_c in Equation (5), since it is independent of d and ξ_b.

In the next section, we constrain the Eddington flux by fitting the burst cooling tracks to spectral models. We can use the model light curve to say something about the expected value of F_Edd, the observed flux which corresponds to the Eddington flux at the surface of the star. In the model, the Eddington flux locally is F_Edd = cg/κ = 0.882 × 10²⁵g₁₄ erg cm⁻² s⁻¹ for X = 0.7, giving F_{model, pk}/F_{model, Edd} = 1.46/g₁₄. Scaling the observed peak flux, the Eddington flux as observed at infinity should be 1.95 × 10⁻⁸ g₁₄(1.7/1 + X) erg cm⁻² s⁻¹, or F_Edd/(10⁻⁸ erg cm⁻² s⁻¹) = 1.95, 3.89, 7.79 for log₁₀g = 14.0, 14.3, 14.6.

We now turn to fitting theoretical calculations of the color-correction factor f_c to the data. First, we describe the fitting procedure and results (Section 5.1) and then the constraints on neutron star parameters, in particular an upper limit on R_∞ (Section 5.2). In Section 5.3, we discuss the variation of K with accretion rate (Figure 2) in the context of the spectral models.

5.1. Fit for A and F_Edd

Suleimanov et al. (2011b) calculated f_c as a function of F/F_Edd for a range of surface gravities and atmospheric compositions, and discussed how these models could be applied to data. We follow their analysis and fit the theoretical f_c-F/F_Edd curves to the observed relation between K^−1/4 and flux F. The fitting parameters are F_Edd and A = K^−1/4/f_c. Comparing with Equation (4), we see that

$$\begin{equation} A=\left({R_\infty \over d}\right)^{-1/2}\xi _{\rm b}^{1/4}. \end{equation} \tag{ 10 }$$

If f_c was a constant independent of flux, and K was constant in the cooling tail of the burst, fitting for A would be equivalent to using the measured normalization K and a value of f_c to extract R_∞/d. Instead, here we are using the entire cooling track to obtain A. In addition, by fitting the shape of the cooling track we can obtain the overall flux scale F_Edd even though the burst itself does not reach Eddington luminosity.

We start by fitting the data from GS 1826−24 with recurrence time 5.74 hr to the different models from Suleimanov et al. (2011b). Below we use the fits to obtain an upper limit on R_∞, which motivates us to start with the bursts with the largest value of K and therefore larger R_∞ values. At the end of this section, we discuss whether it is possible to include the 4.07 hr recurrence time bursts that have smaller values of K (Figure 2) in a consistent picture. When fitting, for simplicity we calculate χ² based on comparing K^−1/4 and f_c, and do not include the errors in the flux measurement. This seems reasonable because the observational error in the overall flux scale is ${\approx} \delta F/\sqrt{N}$, where the individual flux error δF ≈ 10⁻⁹ erg cm⁻² s⁻¹, smaller than the overall uncertainty in the parameter F_Edd that we obtain from our fits.

Suleimanov et al. (2011b) calculate spectral models for pure H and pure He atmospheres, and solar H/He fractions with different metallicity. Based on the light curve comparison and energetics, Galloway et al. (2004) and Heger et al. (2007) concluded that the accreted layer has solar metallicity and a substantial amount of hydrogen (a solar H/He ratio in their models). Here we fit to the full range of models from Suleimanov et al. (2011b) to investigate how changing composition affects our derived limits on neutron star parameters. Also, the photospheric abundances could be different from the abundances near the base of the layer where the X-ray burst ignites. For example, the metallicity in the burning layer could be enhanced by partially burned fuel left over from a previous burst. The hydrogen fraction in the accreted material could be lower than solar, and so it is useful to consider the pure He limit as a limiting case when the hydrogen fraction at the photosphere is reduced.

Figure 3 shows example fits to the 5.74 hr recurrence time bursts. By varying A and F_Edd, we are able to obtain good fits for fluxes down to approximately 1/3 of the peak flux. At lower fluxes, the behavior of the model and observations is qualitatively similar, in that f_c rises rapidly at low fluxes, but the detailed behavior does not match the models. At late times or low fluxes, K^−1/4 rises more rapidly than predicted. We therefore confine the fit to the initial part of the cooling tail and use it to derive A and F_Edd. To do this, we fit Δt_fit seconds of data starting at the time of peak flux. The time Δt_fit is chosen so that as much data is included in the fit as possible while still giving an adequate fit (with the late time data excluded, we find reduced χ² values in the range 0.23–0.47 for the models listed in Table 1). For all except the solar metallicity models we take Δt_fit = 35 s, corresponding to fluxes down to ≈1/3 of the peak flux. The solar metallicity models begin to deviate from the data after Δt_fit = 10 s (about 1/2 of the peak flux) because of the dip in f_c at low fluxes.

Zoom In Zoom Out Reset image size

Table 1. Fits of Spectral Models to the Burst Cooling Tail

Composition	log10g	Δtfita	A	FEdd	χ2reduced	Upper limits
		(s)	(108)	(10−8 erg cm−2 s−1)	(dof)	dξ1/2b	R∞	Mmax	$R_{M>M_\odot }$
						(kpc)	(km)	(M☉)b	(km)c
Pure H	14.3	35	1.144 ± 0.002	5.15 ± 0.20	0.27 (59)	4.3	10.2	1.3	8.2
Pure He	14.3	35	1.186 ± 0.002	4.14 ± 0.09	0.26 (59)	9.7	21.5	2.8	19.8
0.01 Z☉d	14.0	35	1.138 ± 0.002	4.71 ± 0.14	0.24 (59)	5.5	13.2	1.7	11.3
0.01 Z☉d	14.3	35	1.154 ± 0.002	5.08 ± 0.15	0.23 (59)	5.0	11.6	1.5	9.7
0.01 Z☉d	14.6	35	1.167 ± 0.002	5.43 ± 0.19	0.26 (59)	4.6	10.5	1.4	8.4
0.1 Z☉d	14.0	35	1.170 ± 0.003	7.35 ± 0.66	0.47 (59)	4.0	9.0	1.2	6.8
0.1 Z☉d	14.3	35	1.168 ± 0.004	5.80 ± 0.41	0.45 (59)	4.5	10.2	1.3	8.1
0.1 Z☉d	14.6	35	1.174 ± 0.003	5.75 ± 0.31	0.36 (59)	4.4	9.9	1.3	7.8
Z☉d,e	14.0	10	1.178 ± 0.020	5.92 ± 1.22	0.36 (31)	5.3	12.4	1.6	10.6
Z☉d,e	14.3	10	1.159 ± 0.017	4.81 ± 1.08	0.29 (31)	5.6	13.2	1.7	11.3
Z☉d,e	14.6	10	1.164 ± 0.011	4.84 ± 0.52	0.25 (31)	5.4	12.6	1.6	10.7

Notes. ^aWe fit to data from the time of peak luminosity until Δt_fit seconds later. ^bThe maximum neutron star mass consistent with the upper limit on R_∞ = R(1 + z). ^cThe upper limit on radius assuming a lower limit on mass M > 1 M_☉. ^dThe composition is solar H/He abundance plus the indicated proportion of solar metallicity. ^eThese fits yielded more than one local χ² minima. Here, we report only the minima located at the lowest values of F_Edd and A. See the text for more details.

Download table as:  ASCIITypeset image

The results of the fits for different spectral models are listed in Table 1. We used Markov Chain Monte Carlo implemented with the Metropolis–Hastings algorithm (Gregory 2005, chap. 12) to sample the parameter space and find the distributions for the fitting parameters, A and F_Edd. Those distributions were then each fitted by a Gaussian profile in order to derive their respective central values and 1σ uncertainties. In certain cases, as noted in Table 1, the distributions were not well described by a single Gaussian profile due to the presence of more than one peak. In those cases, we fit a Gaussian profile to the peak at the lowest values of F_Edd and A. Note that these parameters are correlated since an increase (decrease) in A, which moves the model curves upward (downward) with respect to the observations, can be offset by a corresponding increase (decrease) in F_Edd, which moves the curves rightward (leftward). The range of F_Edd for all the fitted models is from 4.1 to 7.4 × 10⁻⁸ erg cm⁻² s⁻¹, which lies within the range of F_Edd from the Heger et al. (2007) models used in Section 4. Excluding the 0.1 Z_☉, log₁₀g = 14.0 model, the range of F_Edd is relatively narrow, 4.1–5.9 × 10⁻⁸ erg cm⁻² s⁻¹.

The behavior at low fluxes is shown in more detail in the second panel of Figure 3, which has a logarithmic flux axis. The slope of the increase in f_c with decreasing flux is steeper in the data than in the models for low-metallicity models. For solar metallicity models, the slopes are similar. This would enable a good fit of the whole data set to those models, particularly at low fluxes, but only if F_Edd is in the range 15–25 × 10⁻⁸ erg cm⁻² s⁻¹. In fact, even with the restriction of using only 10 s of data following the burst peak we found other adequate fits, as separate local χ² minima, in that range of F_Edd. This is much larger than expected and as can be seen from the relations derived below would give very small limits on R_∞, and so we do not consider these fits further.

5.2. Upper Limit on Distance and R_∞

A measurement of A and F_Edd typically translates into two values for M and R as follows. When the flux at the surface of the star is at the local Eddington flux cg/κ, the observed flux is

$$\begin{equation} F_{\rm Edd}={{\it GMc}\over \kappa d^2}{1\over 1+z}\xi _{\rm b}^{-1}, \end{equation} \tag{ 11 }$$

where we take the opacity to be κ = 0.2 cm² g⁻¹ (1 + X), as used in Suleimanov et al. (2011b). Following Steiner et al. (2010) we define the quantities

$$\begin{equation} \alpha \equiv {\kappa d\over c^3}F_{\rm Edd}A^2 \xi _{\rm b}^{1/2} = {u\over 2}(1-u) \end{equation} \tag{ 12 }$$

$$\begin{equation} \gamma \equiv {c^3\over \kappa }{1\over A^4F_{\rm Edd}} = {R\over (u/2)(1-u)^{3/2}}, \end{equation} \tag{ 13 }$$

where u = 2GM/Rc². These definitions differ slightly from those of Steiner et al. (2010) in that they include the anisotropy parameter ξ_b. Then

$$\begin{equation} u = {1\over 2}\pm {1\over 2}\left(1-8\alpha \right)^{1/2} \end{equation} \tag{ 14 }$$

$$\begin{equation} R=\alpha \gamma \left(1-u\right)^{1/2}. \end{equation} \tag{ 15 }$$

To calculate α and γ from A and F_Edd obtained from the fits, we require distance d and composition X. A given α and γ then give two solutions for M and R.

We treat the A and F_Edd values as given independently of the derived M and R. In fact the derived A and F_Edd values depend on the gravity assumed for the spectral models. The color correction f_c decreases with increasing gravity in the models of Suleimanov et al. (2011b; see their Figure 5). This implies an increasing A with gravity, which we find in our results for the Z = 0.01 Z_☉ models. A similar trend is not seen in the Z = 0.1 Z_☉, Z_☉ models. The color correction still decreases with gravity, but the shapes of the models change in such a way as to favor smaller values of F_Edd, which, given the correlation between the two fitting parameters, leads to values for A that are smaller than expected for the higher gravities. The value of F_Edd also increases with gravity for the Z = 0.01 Z_☉ models because the slope of f_c with F/F_Edd steepens with increasing gravity, requiring a larger value of F_Edd to agree with observed K^−1/4–F slope.

Overall, we found that the fits are sensitive to the detailed shape of the atmosphere models. This could be due to some small irregularities in the model slopes attributable to the coarseness of the flux grid on which the color corrections are evaluated, and the limited range of fluxes spanned by the bursts we are analyzing.

Table 1 shows that for a given metallicity, varying the surface gravity from log₁₀g = 14.0 to 14.6 changes A by up to ≈3% for low metallicity, and F_Edd by up to ≈20%. The resulting changes in the limit on R_∞ are ≈30% for low-metallicity models. This gives a measure of the error introduced by not carrying out a self-consistent fit in which the gravity of the spectral model used and derived M and R are consistent.

Equation (12) shows that real-valued solutions for u require α ⩽ 1/8 (Steiner et al. 2010). As emphasized by Suleimanov et al. (2011a, 2011b), this gives an upper limit on the distance for which solutions are possible,

$$\begin{eqnarray} \xi _{\rm b}^{1/2} d &\le & {1\over 8} {c^3\over \kappa }{1\over F_{\rm Edd}A^2} \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &=& 5.6\ {\rm kpc}\ \left({A_8\over 1.2}\right)^{-2}\left({F_{{\rm Edd},-8}\over 4.0}\right)^{-1}\left({1+X\over 1.7}\right)^{-1}, \end{eqnarray} \tag{ 17 }$$

where A₈ = A/10⁸ and F_{Edd, −8} = F_Edd/(10⁻⁸ erg cm⁻² s⁻¹). It also provides a limit on R_∞ = R(1 + z). To see this, note that the neutron star radius is

$$\begin{eqnarray} R_\infty &=&R(1+z)=\alpha \gamma \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &=&12.0\ {\rm km}\ \alpha \ \left({A_8\over 1.2}\right)^{-4}\left({F_{{\rm Edd},-8}\over 4.0}\right)^{-1}\left({1+X\over 1.7}\right)^{-1}. \end{eqnarray} \tag{ 19 }$$

where the definition of γ from Equation (13) was used and (1 − u)^1/2 was substituted with (1 + z)⁻¹. An upper limit on R_∞ is obtained by setting α = 1/8 in Equation (19).

The upper limits on ξ^1/2_bd and R_∞ are given in Table 1. To calculate them we use Equations (16) and (19) with 95% lower limits on the quantities A²F_Edd and A⁴F_Edd derived from our fits. A slightly different procedure is used for the cases where the fits yielded multiple χ² minima. To derive the most conservative upper limits on ξ^1/2_bd and R_∞, we consider only the χ² local minimum located at the lowest value of F_Edd and A, manifested as a distinct, Gaussian-like peak in the respective distributions for the quantities A²F_Edd and A⁴F_Edd. Considering only the part of the Gaussian-like distribution lying below the peak value, we find the 90% lower limits for A²F_Edd and A⁴F_Edd. This is equivalent to taking the 95% lower limit of the whole peak, but has the advantage of allowing us to isolate the χ² minimum of interest from the rest of the distribution. As a check, we applied this procedure to model fits showing a single χ² minimum, and found very small differences (<1%) in the derived upper limits when compared with those found by considering the entire distributions.

An upper limit on R_∞ implies an upper limit on the neutron star mass M_max = c²R_∞/3^3/2G (at that mass the radius is $R_\infty /\sqrt{3}$), also given in Table 1. Note that the upper limits on R_∞ and d are correlated. Since ξ^1/2_bd_lim = γA²/8 (compare Equations (13) and (16)), we can rewrite Equation (19) as

$$\begin{equation} R_\infty <12.0\ {\rm km}\ \left({\xi _{\rm b}^{1/2} d_{\rm lim}\over 5.6\ {\rm kpc}}\right)\left({A_8\over 1.2}\right)^{-2}, \end{equation} \tag{ 20 }$$

a larger distance limit allows larger radii.

For solar abundance of hydrogen at the photosphere, we find ξ^1/2d  ≲  4.0–5.6 km and R_∞  <  9.0–13.2 km. This represents quite stringent limits on the neutron star mass and radius. For this range of R_∞, the maximum neutron star mass is in the range 1.2–1.7 M_☉. If we consider a lower mass limit of 1 M_☉, the neutron star radius must be smaller than R(1 M_☉) = 6.8–11.3 km (the individual values for each model are given in Table 1).

5.3. Variation of K with Accretion Rate

Figure 2 shows that the bursts with recurrence times of 4.07 hr have significantly smaller values of K than the 5.74 hr bursts, by ≈20% (see Galloway & Lampe 2012 for a detailed discussion of the variation of K in the sample of bursts from GS 1826−24). Variations in K between bursts have been seen in other sources. For example, Damen et al. (1989) found that the blackbody temperature (evaluated at a fixed flux level) depended on burst duration. They suggested that variations in chemical composition at the photosphere and the resulting changes in color correction might explain the changing blackbody temperature (and therefore normalization).

We investigate two possible composition variations: changing metallicity with solar H/He abundance, and changing the hydrogen fraction. First, we consider solar H/He abundance and changing metallicity. Suleimanov et al. (2011b) show that f_c drops with increasing metallicity. Therefore, we fit the solar metallicity model to the 5.74 hr bursts to determine values of A and F_Edd (as given in Table 1). These values are then used to compare a low-metallicity model to the 4.07 hr data. This comparison is shown in the left panel of Figure 4. The low-metallicity model lies below the 4.07 hr data, showing that the difference in K cannot be explained by a decrease in metallicity from solar to a fraction of solar.

Zoom In Zoom Out Reset image size

Second, we consider a change in hydrogen fraction at the photosphere. The right panel of Figure 4 shows the pure He atmosphere fit for the 5.74 hr bursts (see Table 1), and a low-metallicity solar H/He abundance model for the 4.07 hr bursts in which we use the same value of A determined by the pure He atmosphere fit, but decrease the derived F_Edd by a factor of 1 + X = 1.7 to account for the difference in Eddington flux with composition. This plot shows that the change in f_c in going from pure He to solar H composition is enough to account for the variation in K observed. However, the solar composition model does not match the 4.07 hr data in terms of location on the F/F_Edd axis. Another way to say it is that if we fit the 4.07 hr data with a solar composition model, the required F_Edd would be larger than for the 5.74 hr data, instead of being a factor 1 + X times smaller, as is required for simultaneous fits. Furthermore, we see in the lower panel of Figure 4 that reducing the derived F_Edd by a factor of 1 + X = 1.7 for the 4.07 hr bursts implies that the peak flux for those bursts exceeds the Eddington limit, which is known not to be the case. Therefore, a consistent explanation of the variation in K in terms of changing H fraction at the photosphere is not possible.

We have compared light curve and spectral models with observations of type I X-ray bursts from GS 1826−24. Here we summarize the main conclusions and discuss our results further.

A general point is that anisotropy in the burst emission enters as an additional uncertainty in any derived quantity that depends on distance. Since it changes the relation between the source luminosity and observed flux, the anisotropy parameter ξ_b (defined in Section 3) always enters in combination with distance as ξ^1/2_bd. Even in cases where the distance to a source can be accurately determined, the anisotropy introduces an effective uncertainty of up to a factor of 20%–30%. Anisotropy could be a smaller effect for PRE bursts if the inner disk is disrupted during the burst and intercepts a smaller amount of radiation than a disk extending all the way to the stellar surface. Nonetheless, it remains a source of systematic error on derived neutron star radii that needs to be investigated further. For GS 1826−24, the limit i < 70° from Homer et al. (1998) gives ξ^−1/2_b = 0.9–1.2. Given this uncertainty and the fact that the distance to GS 1826−24 is not well constrained, we focused on deriving limits on M and R that are independent of distance and anisotropy.

The first of these constraints comes from using the model light curve from Heger et al. (2007) to fix the overall luminosity scale of the observed bursts. We showed that this leads to a distance- and anisotropy-independent relation between the redshift 1 + z and color-correction factor f_c (Equation (5)) that depends weakly on the measured normalization K and the ratio of observed and model peak fluxes. For a color correction between 1.4 and 1.5, which spans the range of values in Figure 6 of Suleimanov et al. (2011a), for example, the inferred redshift is between z = 0.19 and 0.28.

The second constraint comes from comparing the spectral evolution during the cooling tail with the spectral models of Suleimanov et al. (2011b), which determines the Eddington flux F_Edd and the quantity A = K^−1/4/f_c. As noted by Suleimanov et al. (2011b), for a given set of measured F_Edd, A parameters, there is an upper limit to the distance of the source beyond which there is no solution for M and R. We point out here that measuring A and F_Edd also places an upper limit on R_∞ = R(1 + z) (and therefore upper limits on M and R for a given source). This limit is independent of distance and anisotropy and depends only on the measured values of A and F_Edd and the surface hydrogen fraction. For GS 1826−24, atmospheric models with solar hydrogen fractions give R_∞ < 9.0–13.2 km (Table 1), which implies a neutron star mass M < 1.2–1.7 M_☉ and R < 6.8–11.3 km assuming a lower mass limit of 1 M_☉. The corresponding distance limits are d < 4.0–5.6 kpc ξ^−1/2_b.

Uncertainties associated with absolute flux calibration do not affect our results; they are equivalent to an incorrect measurement of the distance to the source, which our constraints are independent of.

The constraints on M and R are summarized in Figure 5. We show the upper limits on R_∞ from Table 1 for all the solar hydrogen composition models each plotted at the respective surface gravity and the pure helium model with log g = 14.3, and the redshift range 1 + z = 1.16–1.31 from Equation (5) with f_c = 1.4–1.5 and a 10% uncertainty in the ratio F_obs/F_model. The limits on radii for the solar hydrogen composition are comparable to but a little lower than current theoretical expectations based on dense matter calculations that have radii of 10–13 km for neutron star equations of state that reach a maximum mass >2 M_☉ (Hebeler et al. 2010; Gandolfi et al. 2011). The mass–radius relation found in Steiner et al. (2010), derived from a set of PRE X-ray bursts and hydrogen atmosphere fits for transiently accreting neutron stars in quiescence, also lies at slightly larger radii than our R_∞ limits for solar composition. It should be noted that Suleimanov et al. (2011a) call into question the results of Steiner et al. (2010) by suggesting that "short" PRE bursts should be excluded from analysis as they show smaller blackbody normalizations in the burst tail and also do not follow the theoretically expected spectral evolution. The implication is that the mass–radius relation derived in Steiner et al. (2010) would shift to higher radii as a result of using the more reliable "long" PRE bursts, and thus farther away from our derived upper limits.

Zoom In Zoom Out Reset image size

A smaller hydrogen fraction at the photosphere in GS 1826−24 would increase the R_∞ limits and make them consistent with the theoretical calculations and the mass–radius curve from Steiner et al. (2010). We can get an impression of what the upper limit on R_∞ would be for an atmosphere with a reduced hydrogen fraction by first looking at the extreme case of the pure helium atmosphere and its derived upper limit of 21.5 km for log₁₀(g) = 14.3 (see Figure 5). Such an upper limit is consistent with theoretical calculations and the results from Steiner et al. (2010). We can go one step further and estimate the hydrogen fraction we would require to have such a consistency with previous results using Equation (19). Assuming a surface gravity of log₁₀(g) = 14.3, an upper limit on R_∞ of ∼16 km or more would be required. Using values for F_Edd and A averaged across the solar H/He model fits with log₁₀(g) = 14.3, we estimate that a hydrogen fraction of X ≈ 0.5 or less would be needed. We are assuming that such a spectral model would not differ too greatly in shape from the solar H/He models. Galloway et al. (2004) find that the theoretical variations in burst properties with persistent flux between ignition models with X = 0.7 and X = 0.5 are largely indistinguishable. However, more burst light curve simulations would be necessary to establish whether agreement with observed light curves is still possible with a reduced accreted hydrogen fraction.

Given that the upper limit on R_∞ depends on the color correction as f⁴_c (via A in Equation (19)), even a small 5% increase in the value of f_c would yield a ∼22% increase in the upper limit on R_∞. Furthermore, Suleimanov et al. (2011b) discuss large color-correction factors of f_c = 1.6–1.8 (cf. also Suleimanov & Poutanen 2006) possibly arising from a spreading layer associated with accretion onto the neutron star equator. Perhaps in GS 1826−24 something similar is happening, although the increase in color correction required is not as large. It is worth noting that changing the visible area, for example, by blocking one hemisphere of the neutron star with the accretion disk during the burst, does not change the inferred limits on radius because the limit on R_∞ is independent of the anisotropy factor ξ_b (isotropic emission from only half the area is equivalent to setting ξ_b = 2).

There are several points to keep in mind when looking at our derived constraints on M and R. First, the constraints are only partly self-consistent in the sense that the light curve model used to fit the data does not have the same gravity as the derived M and R. Heger et al. (2007; and Woosley et al. 2004) used a specific choice of gravity in their X-ray burst simulations. As we argue in Section 4, the light curve probably does not depend too sensitively on gravity, but additional simulations are needed to check this, and to calculate the uncertainty in the predicted model flux which enters in Equation (5) relating f_c and 1 + z. On the other hand, for the comparison to spectral models, in Figure 5, the upper limits on R_∞, represented by the solid colored curves, are placed in such a way as to coincide with the appropriate curve of constant surface gravity, consistent with the atmosphere spectral models used to derive those upper limits.

A second issue is that the upper limit on R_∞ from the spectral models is based on fitting the initial part of the cooling tail only. We found that at first the slope of K^−1/4 with flux agrees well with the theoretical models of f_c. In the latter part of the burst, however, at lower fluxes, the agreement breaks down. For F/F_Edd ≲ 0.2–0.3, K^−1/4 increases with decreasing flux, but more rapidly than expected based on the predicted f_c values, particularly those given by the solar metallicity models. Some other explanation is required for the rapid increase in f_c and corresponding decrease in blackbody normalization in the tail of the burst. in't Zand et al. (2009) suggest that this decrease could be due to incorrect subtraction of the persistent emission, in particular, the subtraction of a thermal component that comes from the neutron star surface during accretion, which is no longer present during the burst. Van Paradijs & Lewin (1986) pointed out that this effect should become important during the tail of the burst when the burst flux becomes comparable to that of the accretion. Looking at Figure 3, the observations and the low-metallicity (solar metallicity) models begin to deviate at fluxes below ∼5 × 10⁻⁹ erg cm⁻² s⁻¹(≳ 10⁻⁸ erg cm⁻² s⁻¹) compared with the persistent flux of 2.1 × 10⁻⁹ erg cm⁻² s⁻¹.

The disagreement between the observations and models begins sooner following the burst peak for solar metallicity than for low-metallicity models because the former have a depression in f_c at low fluxes F/F_Edd ≲ 0.3 (see Figure 3), arising from absorption edges in partially ionized Fe (Suleimanov et al. 2011b). There is no sign of such a dip in the observations of GS 1826−24. This suggests a low metallicity in the photosphere, contrasting with the conclusions of Galloway et al. (2004) and Heger et al. (2007) who argued that the metallicity was solar, based on burst light curves and energetics. A way to reconcile these disparate results is to consider the possibility that Fe, whose presence has a significant influence on the spectrum but not on the burst energetics, may be absent from the atmosphere during bursts. If accretion halts during the burst, then Fe will rapidly sink through the atmosphere (Bildsten et al. 2003). On the other hand, since bursts from GS 1826−24 are all sub-Eddington, accretion may continue during the burst, resupplying Fe to the photosphere. Furthermore, while disk accretion only deposits mass near the equator, the accreted mass spreads faster latitudinally (≲ 0.1 s; Inogamov & Sunyaev 1999; Piro & Bildsten 2007) than the timescale for Fe to sink through the atmosphere of ∼1 s. Proton spallation could also destroy a substantial amount of the accreted Fe (Bildsten et al. 2003). It should be noted that X-ray bursts produce a wide range of elements in their ashes that could significantly alter the spectrum. However, mixing of burned material to the photosphere is thought not to occur due to the substantial entropy barrier (Joss 1977; Weinberg et al. 2006).

Another important unresolved issue is the changing spectral normalization K with accretion rate. The blackbody normalization K is ≈20% smaller for the 4.07 hr recurrence time bursts than the 5.74 hr recurrence time bursts. We cannot explain this difference by changing the composition at the photosphere and therefore changing f_c (see the discussion in Section 5.3). Also, it seems unlikely that a major change in composition would occur with only a ≈50% change in accretion rate and a smaller change in burst energy and light curves. As mentioned previously in Section 6, Suleimanov et al. (2011b) discuss large color-correction factors associated with accretion onto the neutron star equator, which they suggest accounts for the variations in measured K for the different spectral states of 4U 1724−307. Whether the 50% increase of accretion rate seen in GS 1826−24 could result in the amount of hardening of the spectrum observed needs to be investigated. If disk accretion onto the star significantly hardens the burst spectrum, it considerably complicates inference of the mass and radius from burst observations, and means that the range of f_c included when calculating errors on mass–radius determinations should allow for a larger range of values than that given by spectral models.

We have found that, even though they do not reach Eddington luminosity, the bursts from GS 1826−24 show enough dynamic range in flux as they cool to significantly constrain F_Edd by comparing with spectral models. A promising source to look at further is KS 1731−254 which shows both mixed H/He bursts with similar spectral evolution to GS 1826−24 (Galloway & Lampe 2012) and PRE bursts. Analysis of these different bursts, which occur at different persistent fluxes and involve different fuel compositions (based on their energetics, peak luminosities, and durations), would give a stringent test of the spectral models and help to constrain any additional spectral components.

This work was supported by the National Sciences and Engineering Research Council of Canada (NSERC) and the Canadian Institute for Advanced Research (CIFAR). D.K.G. is the recipient of an Australian Research Council Future Fellowship (project FT0991598). The authors are members of an International Team in Space Science on type I X-ray bursts sponsored by the International Space Science Institute (ISSI) in Bern, and we thank ISSI for hospitality during part of this work.