Thermonuclear burst observations for model comparisons: a reference sample

TL;DR

This paper compiles and analyzes observations of various thermonuclear X-ray bursts, providing detailed data sets to serve as benchmarks for testing and improving theoretical ignition models.

Contribution

It presents a curated sample of burst observations with detailed properties, facilitating direct comparison with numerical ignition simulations and advancing understanding of burst mechanisms.

Findings

Identification of four distinct ignition cases

Measured recurrence times and burst profiles

Provided data sets for model validation

Abstract

We present a sample of observations of thermonuclear (type-I) X-ray bursts, selected for comparison with numerical models. Provided are examples of four distinct cases of thermonuclear ignition: He-ignition in mixed H/He fuel (case 1 of Fujimoto et al. 1981); He-ignition in pure He fuel, following exhaustion of accreted H by steady burning (case 2); ignition in (almost) pure He accumulated from an evolved donor in an ultracompact system; and an example of a superburst, thought to arise from ignition of a layer of carbon fuel produced as a by-product of more frequent bursts. For regular bursts, we measured the recurrence time and calculated averaged burst profiles from RXTE observations. We have also estimated the recurrence time for pairs of bursts, including those observed during a transient outburst modelled using a numerical ignition code. For each pair of bursts we list the burst properties including recurrence time, fluence and peak flux, the persistent flux level (and inferred accretion rate) as well as the ratio of persistent flux to fluence. In the accompanying material we provide a bolometric lightcurve for each burst, determined from time-resolved spectral analysis. Along with the inferred or adopted parameters for each burst system, including distance, surface gravity, and redshift, these data are suggested as a suitable test cases for ignition models.

Tables

Table 1. Table 1: Target thermonuclear burst source properties

	Dist.	Accreted fuel			$g$	$R$
Source	(kpc)	$X_0$	$Z_{\mathrm{CNO}}$	$1+z$	(${10}^{14}\mathrm{cm}{s}^{-2}$)	(km)	Ref.
GS 1826$-$24	6.1	0.7	0.02	1.23	2.34	12.1	[1,2]
SAX J1808.4$-$3658	$3.4\pm 0.1$	${0.48}_{-0.08}^{+0.12}$	${0.017}_{-0.005}^{+0.007}$	1.26	1.86	11.2	[2]
4U 1820$-$303	$7.6\pm 0.4$	$\lesssim 0.1$	0.02	1.409	2.96	$11.1\pm 1.8$	[3,4,5]
4U 1636$-$536	$5.6\pm 0.4$	0.7	0.02	1.26	1.86	11.2

Table 2. Table 2: Properties of thermonuclear bursts observed from target sources by the Rossi X-ray Timing Explorer

	$N_{\mathrm{burst}}$	burst	$\Delta t$	$F_{\mathrm{per}}^b$		$\dot{m}$	$E_b$	$F_{\mathrm{pk}}$
Epoch	($N_{\mathrm{av}}$)	IDsa	(hr)	(${10}^{-9}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	$c_{\mathrm{bol}}$	(${\dot{m}}_{\mathrm{Edd}}$)	(${10}^{-6}\mathrm{erg}{\mathrm{cm}}^{-2}$)	(${10}^{-9}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	$\alpha$
GS 1826$-$24
1998 Jun	6(1)	7	$5.14\pm 0.07$	$1.167\pm 0.006$	$1.806\pm 0.009$	0.0513	$1.102\pm 0.011$	$30.9\pm 1.0$	$34.2\pm 0.5$
2000 Sep	11(7)	14–20	$4.177\pm 0.010$	$1.593\pm 0.017$	$1.787\pm 0.003$	0.0692	$1.126\pm 0.016$	$29.1\pm 0.5$	$38.6\pm 0.3$
2007 Mar	10(7)	58–65	$3.530\pm 0.004$	$1.87\pm 0.02$	$1.751\pm 0.003$	0.0796	$1.18\pm 0.04$	$28.4\pm 0.4$	$35.3\pm 1.0$
SAX J1808.4$-$3658
2002 Oct	1(1)	2	$16.55\pm 0.06$c	2.541–2.298	$2.085\pm 0.019$	0.0472d	$2.649\pm 0.018$	$229\pm 4$	$114.4\pm 1.9$
	1(1)	3	21.10	2.298–1.946	$2.13\pm 0.04$	0.0432d	$2.990\pm 0.017$	$232\pm 4$	$118.2\pm 1.9$
	1(1)	4	29.82	1.946–1.826	$2.157\pm 0.002$	0.0384d	$3.46\pm 0.02$	$232\pm 4$	$128.2\pm 2.1$
4U 1820$-$303
1997 May 4	3(1)	1	$2.681\pm 0.007$	$3.72\pm 0.18$	$1.45\pm 0.09$	0.144	$0.381\pm 0.003$	$61\pm 2$	$138.5\pm 1.4$
2009 Jun 12	2(2)		${1.892}^e$	$5.70\pm 0.04$	1.4981	0.226	$0.371\pm 0.010$	$56.6\pm 1.4$	$160.1\pm 1.8$
4U 1636$-$536
2001 Feb	1(1)		$\le 1.53\times {10}^4$ f	$4.73\pm 0.03$	1.5346	0.167	$110\pm 9$ g	$21.9\pm 0.6$ g

Equations

$$t_{\mathrm{CNO}}=9.8\left(\frac{X_0}{0.7}\right){\left(\frac{Z_{\mathrm{cno}}}{0.02}\right)}^{-1}\mathrm{hr}$$

$$\alpha =\frac{Q_{\mathrm{grav}}}{Q_{\mathrm{nuc}}}(1+z)\approx \frac{GM}{R}{Q}_{\mathrm{nuc}}^{-1}(1+z)$$

$${\alpha }_{\mathrm{obs}}=\frac{F_p{c}_{\mathrm{bol}}\Delta t}{E_b}$$

$$\dot{m}=4\pi d^2{F}_p{c}_{\mathrm{bol}}/{\dot{m}}_{\mathrm{Edd}}$$

$$F_X=\frac{L_X}{4\pi d^2{\xi }_b(1+z)}$$

$$F_b=\frac{L_b}{4\pi d^2{\xi }_b}$$

Full text

1 INTRODUCTION

Thermonuclear (type-I) bursts arise from unstable ignition of accreted fuel (typically mixed hydrogen and helium) on the surface of accreting neutron stars in low-mass binary systems. Such events are of high priority for observers, as they provide information about the fuel composition, accretion rate, and even neutron star spin, mass and radius (e.g. Heger et al. 2007; Strohmayer & Bildsten 2006; Steiner et al. 2013). A key component contributing to our understanding of the burning physics is numerical modelling of the complex series of nuclear reactions which trigger and power the bursts (e.g. Fujimoto et al. 1981). To date, the degree to which predictions of numerical models have been compared in detail to observations is limited; this paper is part of a wider effort to address this situation.

The physical conditions and processes which broadly influence the burst behaviour are relatively well understood. The primary determinant is the local accretion rate $\dot{m}$, frequently expressed as a fraction of the Eddington rate $\dot{m}_{\rm Edd}$111Throughout this paper, we adopt a value of $\dot{m}_{\rm Edd}=8.8\times 10^{4}\ {\rm g\,cm^{-2}\,s^{-1}}$, at which the outwards force due to radiation pressure (assuming spherical symmetry) equals the surface gravity. The Eddington rate is also the predicted threshold at which the burning is expected to stabilise (e.g. Narayan & Heyl 2003), although observationally bursting behaviour appears to cease for many sources at a substantially lower level Cornelisse et al. (2003); Galloway et al. (2008).

For systems accreting a mix of hydrogen and helium at typical accretion rates ($\sim 0.1\dot{m}_{\rm Edd}$), the fuel layer is hot enough that H will burn stably via the hot-CNO cycle. This process is limited by $\beta$-decay, and will exhaust the H at the base of the layer in a time

[TABLE]

(e.g. Lampe et al. 2016) which depends only on the accreted H-fraction $X_{0}$ and the metallicity $Z_{\rm cno}$. Thus, the composition of the fuel layer at ignition may vary, and a number of categories have been identified by numerical studies. If the recurrence time is short, the burst will ignite in a mixed H/He environment, leading to a “case 1” Fujimoto et al. (1981) or “prompt mixed” Narayan & Heyl (2003) burst. At lower accretion rates, the recurrence time may exceed $t_{\rm CNO}$, leading to ignition in a pure-He layer, i.e. “case 2” or “delayed helium” bursts. At even lower rates, it is thought that H-burning too will become unstable, leading to ignition via that mechanism; but no unambiguous examples of such bursts have been observed.

An additional ignition case is provided by neutron stars in ultracompact systems, which likely accrete (almost?) pure He. Here the ignition will always be in a H-poor environment, largely independently of the accretion rate.

The nuclear energy generation during the burst is related to the average H-fraction in the fuel $\left<X\right>$, with $Q_{\rm nuc}=1.6+4\left<X\right>\ {\rm MeV\,nucleon^{-1}}$. The composition of the fuel may be inferred from the $\alpha$ parameter, the ratio of burst energy to accretion energy over the burst recurrence time. For a neutron star of mass $M$, radius $R$,

[TABLE]

where $z$ is the gravitational redshift at the neutron star surface, given by $1+z=(1-2GM/Rc^{2})^{-1/2}$. Provided the accretion rate is proportional to the persistent flux $F_{p}$, $\alpha$ may be measured (up to a factor corresponding to the ratio of the anisotropy of burst and persistent emission) as

[TABLE]

where $c_{\rm bol}$ is the bolometric correction to the (band-limited) flux $F_{p}$, $\Delta t$ is the (regular) burst recurrence time, and $E_{b}$ the burst bolometric fluence (total energy).

The composition of the fuel may also be deduced from the shape of the lightcurve. He burns during the burst via the triple-$\alpha$ reaction, which proceeds on a much faster timescale than the hot-CNO, rp- and $(\alpha,p)$ reactions which burn H. The higher the He mass fraction, the faster the burning will proceed during the burst, reflecting in the burst rise and duration, as well as the timescale $\tau$ (the ratio of the fluence to the peak flux).

The numerical codes that have been used to compare to observations fall into three broad classes. The first class determines the ignition conditions given the accretion rate and fuel composition, and predicts the burst recurrence time and energetics. While this model has been developed to compare with the inferred atmospheric expansion during a burst as suggested by measurements of burst oscillations Cumming & Bildsten (2000), it has also been used to compare with observations of bursts at low accretion rates Galloway & Cumming (2006). The primary disadvantage of such models is that they do not follow the time-dependent compositional structure of the atmosphere, which is significantly modified by the thermonuclear burning.

The next class of models are one-zone time-dependent codes that simulate some fraction of the nuclear reaction network. Such codes have been used to demonstrate the extent of the rp-process that powers mixed H/He bursts Schatz et al. (2001), as well as serving to probe the sensitivity of burst lightcurves to individual reaction rates Cyburt et al. (2016).

The current (practical) state of the art for burst modelling is in 1-D multi-zone models that track the full extent of the nuclear reaction networks. Codes in this class include kepler Woosley et al. (2004); the general relativistic hydrodynamics code AGILE, coupled with a nuclear reaction network Fisker et al. (2006), and MESA Paxton et al. (2015).

1-D multi-zone models have been used to make the most detailed comparisons of observations to models, most notably for GS 1826$-$24 Heger et al. (2007); Zamfir et al. (2012). However, these comparisons are generally limited in the extent of the sources and the observational data utilised. These limitations, despite the extensive modeling capabilities (and accumulated observational data) motivate the present study.

Here we present a sample of observations of thermonuclear bursts, intended for comparison with numerical models. The objectives of the sample assembly are twofold: first, as test cases to quantify variations between different codes (and hence intrinsic model uncertainties); and second, as examples that may be used to refine the system parameters by direct comparison with individual codes.

We deliberately take a pragmatic approach which incorporates constraints on system parameters obtained by comparison with particular numerical models (of varying degrees of fidelity). While this approach introduces dependence of the system parameters on the specific models chosen for these analyses, this dependence is unimportant for either of the two objectives described above. For the code comparisons, the specific choice of input parameters does not matter, provided they are plausible; and the parameters we have adopted are the best current parameter estimates for the systems studied. As for improving the system parameters via more detailed comparisons with individual models, the adopted system parameters provided here should be viewed as a suggested starting point.

This paper is presented as follows. In section §2 we describe the sources selected to make up the sample. In §3 we describe the observational data from which the sample lightcurves and other parameters are drawn. In section §4 we briefly describe the inferred properties of the bursting sources, and compare the lightcurves in the sample. Finally, in §5 we suggest how the assembled data may be applied and used to test numerical models.

2 SOURCE SELECTION

We selected four well-studied burst sources, which span the range of theoretical ignition cases identified observationally (Table 1). We also selected a fourth source to serve as an example of a superburst, 4U 1636$-$536, which is notable as one of only two such events observed at high sensitivity with RXTE (e.g. Keek et al. 2014).

We list in Table 1 those quantities required to simulate thermonuclear bursts and compare to the observations, including the distance, mass fraction of hydrogen $X_{0}$ and of CNO nuclei, $Z_{\rm CNO}$, and the adopted redshift $1+z$ and surface gravity $g$. We caution that in most cases these parameters are not precisely known, however for the purposes of deciding upon suitable values for simulation tests, we adopt values which are generally consistent with the observed burst properties and/or the neutron star population.

The surface gravity and fuel composition parameters are key input to the models, while the redshift is used to transform the burst lightcurves (and recurrence times) from the neutron-star frame to the observer’s frame (see §5). The source distance $d$

is required to convert model predictions to observed quantities, and is subject to specific uncertainties depending upon the method by which it is estimated, as described further in appendix A.

Below we describe the basic observational properties of each source.

2.1 GS 1826$-$24

This source, discovered as a new transient by Ginga in 1988 September Tanaka (1989), has been X-ray bright and exhibiting bursts ever since. It is known as the “clocked” or “textbook” burster, due to its unusually regular and consistent bursts Ubertini et al. (1999); Galloway et al. (2004). The inferred accretion rate varies little on short timescales, and is within the range 5–13% $\dot{m}_{\rm Edd}$ Galloway et al. (2008); Chenevez et al. (2016). This pattern of regular bursting behaviour has been broken only recently, with the transition to weaker, irregular bursts accompanied by a soft spectral state, in 2014 June Chenevez et al. (2016). By virtue of the typically consistent bursting behaviour, GS 1826$-$24 is also one of the few sources for which detailed efforts have been made to deduce the composition of the accreted fuel. However, these efforts have not yet resulted in unambiguous results.

The measured $\alpha$-values for the source, at $\approx 40$, strongly suggest a high fraction of H in the burst fuel, and so we adopt the solar value, i.e. $X_{0}=0.7$, for the accreted H fraction. The lack of variation in burst fluence (i.e. total burst energy) over a range in accretion rate was taken as evidence by Galloway et al. (2004) that the fuel must be deficient in CNO nuclei. Otherwise, steady (hot-CNO) H-burning between the bursts would signficantly change the fuel composition at ignition, as the burst recurrence time decreased from around 6 hr to 4 hr. A subsequent comparison of burst properties, including the detailed shape of the lightcurve, with time-dependent 1-D model predictions by kepler Woosley et al. (2004) led instead to the conclusion that the accreted fuel must be approximately solar in composition, i.e. $Z_{\rm CNO}\approx 0.02$, at odds with the earlier analysis Heger et al. (2007). Since kepler includes substantially more physics than the code used in the earlier analysis, we adopt solar composition following the most recent study.

Despite the uncertainty in the fuel composition, the properties of the bursts from GS 1826$-$24 strongly suggest burning of mixed H/He fuel by unstable He ignition, i.e. “case 1” of Fujimoto et al. (1981) (or alternatively, “prompt mixed” bursts in regime 3 of Narayan & Heyl 2003).

2.2 SAX J1808.4$-$3658

This source was discovered as a new transient in 1996 with BeppoSAX in ’t Zand et al. (1998), and has been active every 2–3 years ever since. Observations made during the 1998 outburst by RXTE revealed 401 Hz X-ray pulsations, making it the first accretion-powered millisecond pulsar Wijnands & van der Klis (1998). Bright, radius-expansion bursts have been observed in almost every outburst, and the detection of burst oscillations also at 401 Hz in bursts during the outburst of 2002 October confirmed the link between oscillation frequency and neutron-star spin Chakrabarty et al. (2003).

The transient outbursts typically reach a maximum inferred accretion rate of only $0.05\ \dot{m}_{\rm Edd}$, and the bursts have fast ($\lesssim 0.5$ s) rises indicative of He-rich fuel. The wide (2.01 hr) orbit can accommodate a (H-rich) companion Chakrabarty & Morgan (1998), and it is thought that the accreted H is burned away by $\beta$-limited hot CNO burning prior to ignition. The $\approx 20$–30 hr recurrence times, given the low accretion rates, are sufficient provided the metallicity, $Z_{\rm cno}$ is not far below the solar value. This type of burst corresponds to “case 2” of Fujimoto et al. (1981), or “delayed helium” (regime 4) of Narayan & Heyl (2003).

Comparisons of the burst properties with the predictions of numerical ignition models indicate a range of possible compositions, with $Z_{\rm cno}\propto X$ Galloway & Cumming (2006). The lowest plausible H-fraction in the accreted fuel was $X_{0}\approx 0.35$.

2.3 3A 1820$-$303

Located in the globular cluster NGC 6624, 4U 1820$-$303 (also known as Sgr X-4) was first detected as a bursting source with the Astronomical Netherlands Satellite (ANS; Grindlay et al. 1976). When in the bursting state, the source exhibits quasi-regular radius-expansion bursts with recurrence times in the range 2–4 hr (e.g. Cumming 2003). The source intensity is modulated on a $\approx 176$-d periodicity Priedhorsky & Terrell (1984), and it only exhibits bursts in the (relatively short) low-intensity phase of the cycle, when the inferred accretion rate is $\approx 0.2$–0.95 $\dot{M}_{\rm Edd}$ (e.g. Galloway et al. 2008).

The orbital period of 685 s is known from periodic modulation of both persistent X-rays and the optical counterpart King & Watson (1986); Stella et al. (1987), and indicates an “ultracompact” binary which is too close for a Roche-lobe filling main-sequence star. The mass donor is thus assumed to be a very low-mass He white dwarf, and thus accreting pure He (i.e. $X_{0}=0$). Some evolutionary models predict a small ($X\sim 0.1$) amount of H in the surface layers, and this cannot be ruled out from comparing the observed bursts with ignition models Cumming (2003).

The mass and radius of the neutron star in this system have been estimated by equating the blackbody normalisation in the cooling tail and the peak flux of PRE bursts, and adopting the inferred cluster distance Özel et al. (2016). Although unresolved systematic errors may yet affect such measurements (e.g. Steiner et al. 2013), these values represent the best current estimates for the source, and are adopted in Table 1.

2.4 4U 1636$-$536

4U 1636$-$536 is a persistently-accreting burst source in a 3.8 hr orbit van Paradijs et al. (1990). The system is one of the most prolific type-I (thermonuclear) bursters known, and has been studied extensively with most major observatories (e.g. Lewin et al. 1993; Galloway et al. 2008). The neutron star spin has been measured from burst oscillations and transient pulsations during a superburst, at 579.3 Hz Zhang et al. (1997); Strohmayer & Markwardt (2002). The relatively wide orbit, spectral features from hydrogen Augusteijn et al. (1998), and the (at times) long profiles of the bursts strongly suggest the mass donor is a main-sequence star which accretes hydrogen-rich fuel.

The system has exhibited long-term variations in its persistent intensity. Between 1996 and 2001, the X-ray flux was in the range 3–$6\times 10^{-9}\ {\rm erg\,cm^{-2}\,s^{-1}}$, indicating an accretion rate of 11–21% $\dot{m}_{\rm Edd}$ (cf. with Galloway et al. 2008). Since then, the flux has been consistently approximately 50% lower. At these accretion rates, it would be expected that bursts would consistently exhibit long profiles arising from ignition of mixed H/He fuel, although this is not always the case. 4U 1636$-$536 is also well-known for short-recurrence time bursts occuring in trains of up to 4 Keek et al. (2010).

The first superbursts were detected from RXTE/ASM data Wijnands (2001), and two additional events were subsequently reported, one also detected by the PCA Kuulkers et al. (2004); Kuulkers (2009). All the superbursts were detected while the system was in it’s higher flux range, prior to 2001. The 2001 event, on February 22 (MJD 51962.7), was observed with the PCA instrument, and has offered the highest signal-to-noise data during such an event, comparable only to the event seen from 4U 1820$-$303 Strohmayer & Brown (2002). Preliminary spectral fitting of the 4U 1636$-$536 burst indicated that the spectrum could not be adequately fitted with blackbody models Kuulkers et al. (2004), and subsequently it was shown that the burst spectrum included substantial contributions arising from reflection of the burst emission from the accretion disk Keek et al. (2014).

3 OBSERVATIONS & ANALYSIS

We used preliminary results from the Multi-INstrument Burst ARchive (MINBAR222http://burst.sci.monash.edu/minbar), which includes bursts observed by the Rossi X-ray Timing Explorer Proportional Counter Array (RXTE/PCA), BeppoSAX Wide-Field Camera (WFC; Jager et al. 1997; in ’t Zand et al. 2004) and INTEGRAL Joint European X-Ray Monitor (JEM-X; Lund et al. 2003).

We calculated burst lightcurves from time-resolved spectroscopic analysis of selected bursts observed by RXTE/PCA Jahoda et al. (1996). The MINBAR analysis follows the approach of Galloway et al. (2008), incorporating the latest PCA responses and the effects of deadtime, as for the measured peak PRE burst fluxes (see appendix A). Earlier analyses of samples including these data have demonstrated evidence for an increased contribution of persistent flux during bursts, possibly arising from the effects of Poynting-Robertson drag on the inner accretion disk Worpel et al. (2013, 2015). While the fractional increase in the persistent flux can be $\sim 10$, this contribution is a small fraction of the burst flux, and so any correction required to the burst flux is also small (of order a few per cent, similar to the intrinsic uncertainty of the time-resolved flux measurements). Thus, we neglect this effect for the purposes of determining our mean profiles.

3.1 Selecting burst sequences

For GS 1826$-$24 and 4U 1820$-$303, we identified burst sequences in MINBAR for which we could reliably infer the burst recurrence time, even when the observations were interrupted. Such interruptions arise for instruments in low-Earth orbit, due to Earth occultations and passages through the South Atlantic Anomaly. An example burst train is shown in Figure 1.

We then selected and averaged bursts observed by the RXTE/PCA within each burst train. MINBAR also includes bursts observed by BeppoSAX/WFC and INTEGRAL/JEM-X, which are useful for establishing the regularity of the bursts, but the sensitivity of these instruments is much lower than RXTE, and so are unsuitable for producing high-quality burst lightcurves.

For GS 1826$-$24, we augmented our burst sample with optical bursts observed during RXTE observations in 1998 June, as also analysed by Thompson et al. (2008).

We measured the average recurrence time of each train by assigning a trial integer value to each burst, and then performing a linear least-squares fit to the times. The trial values were determined by dividing the time since the first burst, by the minimum burst separation of the train; we subsequently adjusted the sequence numbers to minimise the RMS value of the fit. We typically obtained an RMS value of $\approx 20$ min or better.

From the set of bursts in each train $N_{\rm burst}$, we then selected the RXTE bursts ($N_{\rm av}$, Table 2), and created average lightcurves from each train. We also calculated averaged burst fluences, peak fluxes and $\alpha=\Delta tF_{p}c_{\rm bol}/E_{b}$, where $F_{p}$ is the persistent (accretion) flux, $c_{\rm bol}$ the bolometric correction (see §3.2), and $E_{b}$ the burst fluence.

We estimated the accretion rate $\dot{m}$ as a fraction of the Eddington rate $\dot{m}_{\rm Edd}$, as

[TABLE]

where $F_{p}$ is the average persistent flux over the burst train, and $d$ the estimated source distance from Table 1. This estimate neglects the possible anisotropy of the persistent emission, so must be taken as a guide only.

The identification of suitable burst trains for 4U 1820$-$303 was made difficult by the relatively small number of bursts available in the MINBAR sample. Just 67 bursts were detected by BeppoSAX/WFC, INTEGRAL/JEM-X or RXTE/PCA, with only 16 events observed with RXTE. We ultimately identified eight burst trains with three or more bursts and evidence for regular bursting behaviour; the range of inferred recurrence time was 2.7–4.4 hr.

Only two of these sequences included bursts detected with RXTE, each with just three bursts. For each sequence, we observed a pair of bursts within a few hours of each other, and a third burst some days earlier or later. We ultimately rejected one of these sequences, in 2009 May (MJD 54980; see also in’t Zand et al. 2012) because there were significant variations in the persistent flux (and spectral shape) over the interval, and it was not possible to be confident about the recurrence time for the close pair of bursts, due to data gaps between them. This was not the case for the remaining train, in 1997 May (MJD 50572), and so we retained this set of bursts. With only one burst observed by RXTE in this sequence, we simply provide the time-resolved spectroscopy results from that burst, as was done for SAX J1808.4$-$3658 (see §3.4)

Due to the sparse data for 4U 1820$-$303, we augmented the sample with a pair of bursts observed by RXTE/PCA on 2009 June 12, separated by only 1.89 hr. The nearest observations in time preceeding (following) that observation were 13 d before (3 d after), so any possible test for periodicity over that time range would likely be uninformative given the likely phase drift of the burst sequence. Assuming that the separation of 1.89 hr does reflect a regular recurrence time, this pair would be the most frequent bursts ever observed from this system (cf. with Clark et al. 1977).

The constraints on recurrence times for the superbursts from 4U 1636$-$536 is poor due to the low duty cycle of observations between the detected events. The superburst in our sample was preceded by another $1.75$ years prior. However, the inferred accretion column is approximately twice the ignition column, suggesting that the actual recurrence time is closer to $0.9$ years (Keek et al. 2015). In addition to the properties of the burst itself reported in Table 2, we constrained the “quench time” corresponding to the period following the superburst during which normal thermonuclear burning is suppressed. The first thermonuclear (H/He) event detected from 4U 1636$-$536 following the superburst was on MJD 51985.5 (BeppoSAX observation 10898), giving $t_{\rm quench}\leq 22.8$ d. The observation during which that event was detected commenced 17.5 hr before the burst was detected, and had a duty cycle of $\approx 70$%, with interruptions due to the low-Earth orbit. The coverages suggests a $\approx 30$% chance that an earlier burst could have been missed, so that the quench time could have been shorter by up to 17.5 hr. However, no intervening observations were made to rule out additional bursts missed earlier than that observation.

3.2 Bolometric corrections

We measured the persistent flux in the 3–25 keV band from RXTE and BeppoSAX observations by fitting commonly-used phenomenological models, and converted the integrated flux to a bolometric value by applying a bolometric correction for each observation. The correction was estimated by fitting to RXTE PCA and HEXTE data a broad-band model typically comprising an absorbed Comptonisation component (compTT in xspec; Titarchuk 1994), usually with a Gaussian component representing Fe K$\alpha$ emission in the range 6.4–6.7 keV. We measured the 3–25 keV flux from each fit, and then created an ideal response (using dummyrsp in xspec) covering the energy range 0.1–1000 keV, and integrated the flux over this range, excluding the effects of neutral absorption. We estimated the bolometric correction as the ratio of the unabsorbed 0.1–1000 keV flux to the absorbed 3–25 keV flux; for the observations here, the correction was typically in the range 1.4–2.

This approach is reliable provided that the contribution to the bolometric flux outside the range to which PCA and HEXTE are sensitive (typically 3–100 keV) is small (or at least consistent). However, it is known that for GS 1826$-$24, that additional low-energy contributions may arise at times, so that the persistent flux measured by RXTE in the 3–25 keV band is not always a reliable estimator of the accretion rate (e.g. Thompson et al. 2008). In some observations, there is evidence for unusually low $\alpha$-values, suggesting that a significant fraction of the X-ray flux is emitted outside the RXTE band; contemporaneous Chandra or XMM-Newton observations suggest that this flux may be emitted as a low-temperature ($\lesssim 1$ keV) thermal component. Thus, we avoided from our data selection bursts from such epochs identified by Thompson et al. (2008), including those from 2003 April, and also the observations from 1998 June and 1997 November. An additional epoch of unusually low $\alpha\approx 26$ is from observations in 2006 August, and was also excluded. The remaining set of three burst trains have broadly consistent $\alpha$-values, but we caution that there may yet be additional systematic errors contributing to the determination of the accretion rate in this (and possibly other) systems.

3.3 Lightcurve fitting

Here we describe the approach by which the values of the neutron star radius, redshift, and surface gravity for GS 1826$-$24 were chosen. We replicated the analysis of Heger et al. (2007) with updated analysis results and an improved method to determine the best-fit system parameters. We selected the train of bursts observed in 2000 September, with a recurrence time of 4.177 hr, similar to that used previously. The 7 bursts observed with RXTE are augmented by the observation of 25 additional bursts with BeppoSAX over the same interval, with the complete train spanning 12 d. However, the entire sequence of bursts is not consistent with a steady recurrence time, likely due to small variations in accretion rate. Thus, we restricted the sequence to a shorter interval covering the RXTE bursts, beginning on MJD 51810 and spanning almost 4 d (as plotted in Fig. 1). We then compared the averaged lightcurve (calculated from the RXTE bursts only) with examples predicted by the kepler code, tabulated in Lampe et al. (2016).

Heger et al. (2007) found excellent agreement between the observed lightcurve and a model referred to as “A3”, with solar composition, an accretion rate of $1.58\times 10^{-9}\ M_{\odot}\,{\rm yr}^{-1}$. In the sample of Lampe et al. (2016), the accretion rate is quantified in terms of the Eddington value, defined as $1.75\times 10^{-8}\ M_{\odot}\,{\rm yr}^{-1}$. Thus, for model A3, the accretion rate was $\dot{M}=0.09\dot{M}_{\rm Edd}$, and this model corresponds to “a05d” in the Lampe et al. (2016) sample.

In our revised comparison, we simultaneously matched the lightcurve and the recurrence time, using a Markov-Chain Monte Carlo (MCMC) code emcee Foreman-Mackey et al. (2013) to marginalise over the three parameters of interest: the redshift $1+z$ (by which the lightcurve will be stretched when transforming into the observer’s frame), the relative intensity scale, including the distance and redshift:

[TABLE]

where $L_{X}$ is the burst luminosity predicted by kepler, $d$ is the distance, $\xi_{b}$ takes into account the possible anisotropy of the burst emission (see also appendix A, and $F_{X}$ the observed burst flux. A third parameter is related to the offset in time between the model and observation, and has no physical meaning.

Our initial runs with this code did not match the recurrence time, and to ensure that the model-predicted burst interval (when scaled by $1+z$) matched that observed, we weighted this parameter in the likelihood function by a factor of 50. This ensured in the subsequent comparisons that the redshifted model recurrence time matched that observed to within a few percent.

The comparison with model a05d, as used by Heger et al. (2007), gave good agreement with the observed lightcurves, although requiring a redshift of $1+z=1.392$. Following the approach of Lampe et al. (2016), we transformed the model results (calculated in a Newtonian frame) to the observer’s frame assuming that the gravity, and the Newtonian and the general-relativistic neutron star masses are identical. This implies that the radius required for the Newtonian and general-relativistic gravity will be different, and the corresponding radius (for a neutron star mass of $1.4\ M_{\odot}$) implied by our comparison of model a05d was 8.5 km, well below the lower limit of 10.4 km inferred for a wide range of neutron stars (e.g. Steiner et al. 2013).

Thus, we made a limited exploration of alternative models, focussing on those with solar composition. The comparison with model a028 gave a best-fit redshift of $1.234$, implying a radius of 12.1 km, well within the expected range. The accretion rate (in units of the Eddington rate) for this model is 0.082, which is somewhat higher than that inferred for the 2000 Sep epoch, at 0.0692; but may be consistent taking into account the persistent emission anistropy. Heger et al. (2007) estimated the ratio of persistent to burst anisotropies at $\xi_{p}/\xi_{b}=1.55$, implying $\xi_{p}>1$ (e.g. Fujimoto 1988); in that case, the luminosity (and hence accretion rate) inferred from the persistent flux would be higher, bringing the observation estimate roughly in line with the value assumed for run a028. The corresponding distance is $d\xi_{b}^{1/2}=6.1$ kpc. We do not quote errors on these parameters, as the lightcurve comparison (Fig. 2) is not of sufficient quality to have confidence in the posterior distributions (the estimated probability densities for the model parameters). However, the comparison is of comparable quality as in Heger et al. (2007), and so for the present purposes, we adopt the values for the redshift and distance.

3.4 Outburst fitting

The approach adopted for SAX J1808.4$-$3658 was necessarily different, as this source has not exhibited regular bursts in previous observations. The persistent flux (and hence accretion rate) varied significantly over the time at which the bursts were observed in 2002 October, and so the burst separations (and hence, presumably, the recurrence times) vary significantly. Furthermore, the burst separations are all longer than 12 hr, so that many data gaps are present between each burst, introduced by the interruptions in the visibility of the source due to the observing satellite’s $\approx 90$-min low-Earth orbit.

Galloway & Cumming (2006) compared the observed burst times with the predictions of a numerical ignition model, to constrain the system parameters. For the purposes of this paper, we have revised the analysis of those authors, to reconcile the changes in the RXTE/PCA calibration since that work was completed, to correctly treat the persistent and burst anisotropy factors, and to use a more rigorous parameter space exploration to determine the best-fit system parameters and uncertainties.

We matched the burst times, fluences and $\alpha$-measurements to the settle code developed by Cumming & Bildsten (2000) using MCMC and marginalising over the accreted hydrogen fraction $X_{0}$, CNO metallicity $Z_{\rm CNO}$, base flux $Q_{b}$, and scaling factors incorporating the source distance and emission anisotropy. The settle model predicts burst parameters in the neutron star frame, given an inferred accretion rate. Necessarily, we adopt the average inferred accretion rate (based on linear interpolation between each observation) for each burst interval, as settle does not simulate varying $\dot{m}$. The redshift $1+z$ and (Newtonian) surface gravity $g$ are fixed in the settle code, at the values listed in Table 1 so the results are implicitly for those values.

The best-fit parameters we report are dependent on the fidelity of the model, which does not take into account all of the underlying physics, such as the effects of thermal and compositional inertia. More detailed models, such as kepler Woosley et al. (2004) can provide more rigorous modelling of the underlying physics, however are computationally expensive and so impractical for the Monte Carlo approach adopted here. However, kepler simulations to reproduce the bursts observed from SAX J1808.4$-$3658 during the 2002 outburst have been performed recently, demonstrating that the assumption of a constant accretion rate between bursts results in the recurrence time being underestimated, when compared to the scenario in which the accretion rate decreases between bursts (Johnston et al. 2017, in preparation). To account for this bias, we corrected the recurrence time predictions of settle using correction factors determined by kepler. The correction factor depends on the $\dot{m}$ gradient, and thus is different for each burst. We adopt correction factors for each of the seven predicted bursts of SAX J1808.4$-$3658, which vary from $\approx$ 1 for the first burst to $\approx$ 1.2 for the final burst.

We infer the persistent and burst anisotropy factors ($\xi_{p}$ and $\xi_{b}$ respectively) from the scaling factors and find $\xi_{p}$ = 1.05$\pm$0.01 and $\xi_{b}$ = 0.91${}^{+0.04}_{-0.03}$. This implies an inclination of 62$\pm$1$\degree$. Note that the error on the inclination is statistical only, and does not include the contribution from uncertainties in the mass and radius, due to the fixed values adopted in the settle model.

As with the analysis of Galloway & Cumming (2006), our maximum-likelihood solution predicts an additional pair of bursts between the first pair of observed bursts, on MJD 52562.41363 and 52564.30515 (#1 and 2 from Galloway et al. 2008). These events would have fallen in data gaps, which explains why they were not observed. We thus report the burst separations for the last two bursts observed during the 2002 October outburst as the recurrence time; for the second burst, we estimate the recurrence time inferred from the burst train modelling based on the simulated time of the preceding event. We note that the inferred recurrence time of 16.55 hr for the second burst observed with RXTE is comfortably within the overall range for the source, with the minimum separation of 12.6 hr set by a pair of bursts observed by BeppoSAX in ’t Zand et al. (1998).

4 RESULTS

The derived (or adopted) properties for the sources of interest are listed in Table 1. The properties of the bursts selected for our sample are summarised in Table 2. We plot the burst lightcurves for each of the sequences in Figure 3. The accompanying material333The lightcurves are intended to be made available in the PASA data store upon publication of the paper, but in the meantime can be found at http://burst.sci.monash.edu/reference. for this paper includes lightcurves (averaged in cases where there are more than one RXTE burst in the sequence) for comparison with numerical model results.

The burst lightcurves show clear differences related to the inferred composition and accretion rate. For 4U 1820$-$303, the burst lightcurves are short, with durations (defined as the interval over which the luminosity exceeds 10% of the maximum) of 15 s or so. Each of these bursts exhibits strong photospheric radius expansion, and the quality of the spectral fits during the radius expansion episodes are relatively poor (reduced-$\chi^{2}$ values in the range 2–8). As $\dot{m}$ increases, the burst timescale increases slightly ($\tau=6.24$ to $6.55$ s), which qualitatively supports the presence of a small amount of H in the accreted fuel. With smaller recurrence time, there is less time to burn the accreted H, and the mass fraction at ignition will be larger.

The bursts from SAX J1808.4$-$3658, although also thought to be powered primarily by He, have a notably different morphology, with an extended radius-expansion phase at roughly constant luminosity. The timescale variation for SAX J1808.4$-$3658 is in the opposite direction, with the burst timescale and duration becoming shorter at higher accretion rates. Here the recurrence time is long enough that any accreted H has already been exhausted at the base, so that lower accretion rate allows a larger pile of He to accumulate. The resulting increased fluence thus requires a longer interval emitting at the Eddington luminosity, with some of the additional energy likely exported as kinetic energy of outflowing material.

The bursts from GS 1826$-$24 are different again, showing the much slower ($\approx 5$ s) rises and long rp-process powered tails characteristic of these mixed H/He bursts. The peak luminosity is significantly smaller than for the He-rich bursts, and durations are approximately a factor of 5 longer. The burst profiles show only very subtle variations from epoch to epoch, with the peak flux decreasing slightly, while the timescale increases, as the accretion rate increases.

The superburst profile is principally distinguished from the shorter, H or H/He events by the orders-of-magnitude longer timescale, with a duration of more than 2 hr.

5 DISCUSSION

We have assembled a sample of burst observations intended for comparisons with, and between, numerical models. Comparisons of different model codes for the same test cases will allow estimates of the intrinsic uncertainty that affects burst simulations. These uncertainties may also depend on the particular type of burst being simulated; that is, the variation in model predictions for some types of bursts may be larger than for others. For example, Galloway & Cumming (2006) suggest that the simple analytic model adopted for the simulations of bursts from SAX J1808.4$-$3658 offer comparable fidelity to more sophisticated models, including kepler. This is less true for simulations of H-rich bursts from GS 1826$-$24 observed at higher accretion rates, as the simple models do not include the effects of compositional or thermal inertia.

It is also anticipated that detailed comparisons of the burst measurements here may allow more stringent constraints to be placed on the neutron star parameters, as has been demonstrated as a proof-of-principle in sections §3.3 and §3.4.

The anticipated procedure to replicate each of the burst measurements in Table 2 is as follows

Set the surface gravity $g$ and neutron star radius $R$ to the values appropriate for the source of interest in Table 1 2. 2.

Set the composition in the accreted material, $X_{0}$ and $Z_{\rm CNO}$ 3. 3.

Perform a model run with a constant accretion rate chosen from one of the entries for that source from Table 2, and extract the burst luminosity as a function of time. 4. 4.

Re-scale the predicted luminosity appropriately for the redshift parameter $1+z$ for this source, listed in Table 2, and also the burst anisotropy parameter $\xi_{b}$, where known. For the second and third bursts from SAX J1808.4$-$3658, the predicted recurrence time and burst fluence may be adjusted by the factor described in §3.4 to account for the effects of the declining accretion rate during each burst 5. 5.

Calculate from the simulated burst train the recurrence time $\Delta t_{\rm pred}$, and from the burst lightcurve, the burst fluence, peak flux, $\alpha$-value, and timescale 6. 6.

Compare these predictions with the observed lightcurve and the measured values in Table 2.

In the event that the model recurrence time is significantly longer (shorter) than the observation, the cause may be that the true accretion rate is different from the estimate due to the unknown degree of anisotropy $\xi_{p}$ of the persistent emission. In that case, the adopted accretion rate might be increased (decreased), by roughly the ratio of the predicted and observed recurrence times, and the procedure repeated from step 3 onwards.

One additional parameter that is commonly included in burst simulations is the heat flux from below the simulation region, commonly parameterised as a “base flux” $Q_{b}$. This quantity is not well constrained by observations, and furthermore depends on the accretion rate (e.g. Fig. 18 of Cumming et al. 2006); most studies to date have adopted values in the range 0.1–0.15 MeV nucleon*-1*. For the purposes of model-to-model comparisons we recommend adopting a value of 0.1, but suggest that for the purposes of matching the observations, this parameter can be considered free within the range 0.1–0.8 MeV nucleon*-1*, depending upon the choice of core neutrino emissivity model.

The approach to replicate the superburst observation is necessarily different. Modeling a superburst by accreting hydrogen/helium comes at great computational expense, because hundreds of hydrogen/helium flashes must be simulated to create the carbon fuel. Furthermore, it remains challenging to reproduce the observationally inferred carbon mass fraction (e.g., Stevens et al. 2014; Keek & Heger 2016). An alternative approach is to directly accrete a mixture of carbon and heavier elements (e.g., Fe or Ru; Cumming & Bildsten 2001; Schatz et al. 2003). For the 2001 superburst from 4U 1636$-$536, we recommend a fuel composition of $26\%$ 12C and $74\%$ 56Fe (Keek et al. 2015). Such simulations are able to describe the observed superburst light curves in detail, but generally predict larger ignition depths than observed (Keek & Heger 2011).

The ultimate goal of this study is to constrain the rates of particular nuclear reactions, or masses of particular isotopes, based on the burst observations. It has been shown that several reactions are likely to affect the shape of the burst lightcurve (e.g. Cyburt et al. 2016). Similarly, nuclear masses have also been shown to have an effect Schatz & Ong (2016). It is hoped that the carefully-assembled and calibrated data presented here offer the first step on the road towards precision nuclear physics from thermonuclear bursts.

Acknowledgements.

This work benefited from discussions at the International Symposium on Neutron Stars in the Multi-Messenger Era in May 2016, and more broadly by support from the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements). The authors are grateful for support received as part of the International Team on Nuclear Reactions in Superdense Matter by the International Space Science Institute in Bern, Switzerland. The data analysed here are derived from preliminary analysis results of the Multi-INstrument Burst ARchive (MINBAR), which is supported under the Australian Academy of Science’s Scientific Visits to Europe program, and the Australian Research Council’s Discovery Projects and Future Fellowship funding schemes. LK is supported by NASA under award number NNG06EO90A.

Appendix A Source distances

An important source of systematic uncertainty in the determination of persistent and burst luminosities is the source distance $d$, and in combination, the degree of anisotropy of the burst and persistent emission, parameterised as $\xi_{b}$ and $\xi_{p}$, respectively. Anisotropy of the burst (and persistent) flux is expected given the strongly non-spherically symmetric geometry introduced by the accretion disk and binary companion. While the distances to burst sources have been estimated with moderate precision from observations of radius-expansion bursts (e.g. Galloway et al. 2008), these measurements implicitly also depend on the burst anisotropy. Strictly speaking, distances estimated in such a manner correspond to $d\xi_{b}^{1/2}$, and in general it is not possible to separate these quantities to derive the physical distance $d$, independently.

The definition of the anisotropy factor $\xi_{b}$ follows that of Fujimoto (1988), such that for a total burst luminosity $L_{b}$, the measured burst flux would be given by

[TABLE]

From the above relation, we understand that $\xi_{b}=1$ corresponds to isotropic emission in all directions, while $\xi_{b}<1$ indicates that the burst flux toward the observer is in excess of the isotropic value, i.e. preferential beaming in our direction. As pointed out by He & Keek (2016), the values of $\xi_{b}$ and $\xi_{p}$ depend upon the system inclination and the structure of the accretion disk, neither of which are well constrained.

For 4U 1636$-$536, we estimated the distance based on measurements of the burst emission from photospheric radius-expansion (PRE) bursts in MINBAR, assuming that they reach the Eddington limit. This follows the approach of Galloway et al. (2008), but the burst flux measurements have been updated with the latest PCA response matrices, and the effects of deadtime in the burst spectra have been included. The overall effect on the distances (compared to the previous analyses) is a reduction of 5–10%. We also incorporated the assumed gravitational redshift at the surface, as listed in Table 1.

The anisotropy of burst emission for 4U 1636$-$536 can be estimated from the inclination of the binary system with respect to the observer’s line of sight (e.g., Fujimoto 1988). Optical observations constrain the inclination angle to the range $36^{\circ}-74^{\circ}$ (Casares et al. 2006), which implies $0.7\lesssim\xi_{\mathrm{b}}^{-1}\lesssim 1.3$ for a thin disk (larger values for smaller angles; He & Keek 2016). The observation of a broad fluorescent iron line favors a large inclination angle (Pandel et al. 2008), but this depends strongly on the interpretation of the reflection spectrum. Similarly, different constraints can be derived from the reflection signal during the superburst from this source (Keek et al. 2015). Therefore, an uncertainty of several tens of percents must be taken into account for the burst luminosity of 4U 1636$-$536.

For GS 1826$-$24 and SAX J1808.4$-$3658, we instead adopted distances derived from comparing the burst measurements with models, as described in sections §3.3 and §3.4.

Zamfir et al. (2012) derived an upper limit on the distance (combined with the burst anisotropy) for GS 1826$-$24 of $d\xi_{b}^{1/2}<5.5$ kpc. While our lightcurve comparisons for this source do not take into account the blackbody normalisation, as was analysed by Zamfir et al. (2012), our derived system constraints are otherwise broadly consistent with those authors. However, our best-fit value based on the comparison with model a028 of Lampe et al. (2016) is 6.1 kpc, slightly higher than the inferred limit.

One additional consistency check that we can make is to compare the expected Eddington flux for the adopted system parameters, with that measured from the first photospheric radius-expansion burst from the source, observed in 2014 June with NuSTAR Chenevez et al. (2016). Assuming that the degree of burst anisotropy $\xi_{b}$ does not change during radius expansion, and that the peak flux is achieved at touchdown (when the expanding photosphere has returned to the neutron star surface), the expected flux would be (e.g. Galloway et al. 2008) $3.7\times 10^{-8}\ {\rm erg\,cm^{-2}\,s^{-1}}$. This is well within the confidence interval of $(4.0\pm 0.3)\times 10^{-8}\ {\rm erg\,cm^{-2}\,s^{-1}}$ derived for the radius-expansion burst observed by Chenevez et al. (2016). In the absence of a more detailed comparison with kepler simulations, which would allow us to incorporate the results from all three burst trains presented in Table 2, we adopt the current value, but note that the true value may be lower.

For 4U 1820$-$303, we adopted the distance to the host globular cluster, NGC 6624 Kuulkers et al. (2003).