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. Reference Heger, Cumming, Galloway and Woosley2007; Strohmayer & Bildsten Reference Strohmayer, Bildsten, Lewin and van der Klis2006; Steiner, Lattimer, & Brown Reference Steiner, Lattimer and Brown2013). 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. Reference Fujimoto, Hanawa and Miyaji1981). 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}$
Footnote
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 Reference Narayan and Heyl2003), although observationally bursting behaviour appears to cease for many sources at a substantially lower level (Cornelisse et al. Reference Cornelisse2003; Galloway et al. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008).
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 β-decay, and will exhaust the H at the base of the layer in a time
$$\begin{equation}
t_{\rm CNO} = 9.8 \left(\frac{X_0}{0.7}\right)\left(\frac{Z_{\rm cno}}{0.02}\right)^{-1}\ {\rm h}
\end{equation}$$
(e.g. Lampe, Heger, & Galloway Reference Lampe, Heger and Galloway2016) which depends only on the accreted H-fraction X 0 and the metallicity Z 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. Reference Fujimoto, Hanawa and Miyaji1981) or ‘prompt mixed’ (Narayan & Heyl Reference Narayan and Heyl2003) burst. At lower accretion rates, the recurrence time may exceed t 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 < X > , with Q nuc = 1.6 + 4 < X >MeVnucleon−1. The composition of the fuel may be inferred from the α parameter, the ratio of burst energy to accretion energy over the burst recurrence time. For a neutron star of mass M, radius R,
$$\begin{equation}
\alpha = \frac{Q_{\rm grav}}{Q_{\rm nuc}}(1+z) \approx \frac{GM}{R} Q_{\rm nuc}^{-1} (1+z),
\end{equation}$$
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, α may be measured (up to a factor corresponding to the ratio of the anisotropy of burst and persistent emission) as
$$\begin{equation}
\alpha _{\rm obs} = \frac{F_{\rm p} c_{\rm bol} \Delta t}{E_{\rm b}},
\end{equation}$$
where c bol is the bolometric correction to the (band-limited) flux F p, Δ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-α reaction, which proceeds on a much faster timescale than the hot-CNO, rp, and (α, 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 τ (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 Reference Cumming and Bildsten2000), it has also been used to compare with observations of bursts at low accretion rates (Galloway & Cumming Reference Galloway and Cumming2006). 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. Reference Schatz2001), as well as serving to probe the sensitivity of burst lightcurves to individual reaction rates (Cyburt et al. Reference Cyburt, Amthor, Heger, Johnson, Keek, Meisel, Schatz and Smith2016).
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. Reference Woosley2004); the general relativistic hydrodynamics code AGILE, coupled with a nuclear reaction network (Fisker et al. Reference Fisker, Görres, Wiescher and Davids2006), and MESA (Paxton et al. Reference Paxton2015).
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. Reference Heger, Cumming, Galloway and Woosley2007; Zamfir, Cumming, & Galloway Reference Zamfir, Cumming and Galloway2012). However, these comparisons are generally limited in the extent of the sources and the observational data utilised. These limitations, despite the extensive modelling 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 two-fold: 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 Section 3, we describe the observational data from which the sample lightcurves and other parameters are drawn. In section Section 4, we briefly describe the inferred properties of the bursting sources, and compare the lightcurves in the sample. Finally, in Section 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. Reference Keek, Ballantyne, Kuulkers and Strohmayer2014).
Table 1. Target thermonuclear burst source properties.
Values in italics indicate there are no constraints specific for that source.
References: (1) Heger et al. (Reference Heger, Cumming, Galloway and Woosley2007); (2) this work; (3) Kuulkers et al. (Reference Kuulkers, den Hartog, in ’t Zand, Verbunt, Harris and Cocchi2003); (4) Cumming (Reference Cumming2003); (5) Özel et al. (Reference Özel, Psaltis, Güver, Baym, Heinke and Guillot2016).
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 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 Section 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 Reference Tanaka, Hunt and Battrick1989), 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. Reference Ubertini, Bazzano, Cocchi, Natalucci, Heise, Muller and Zand1999; Galloway et al. Reference Galloway, Cumming, Kuulkers, Bildsten, Chakrabarty and Rothschild2004). The inferred accretion rate varies little on short timescales, and is within the range 5–13%
$\dot{m}_{\rm Edd}$
(Galloway et al. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008; Chenevez et al. Reference Chenevez2016). 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. Reference Chenevez2016). 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 α-values for the source, at ≈ 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. (Reference Galloway, Cumming, Kuulkers, Bildsten, Chakrabarty and Rothschild2004) that the fuel must be deficient in CNO nuclei. Otherwise, steady (hot-CNO) H-burning between the bursts would significantly change the fuel composition at ignition, as the burst recurrence time decreased from around 6 to 4 h. 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. Reference Woosley2004) led instead to the conclusion that the accreted fuel must be approximately solar in composition, i.e. Z CNO ≈ 0.02, at odds with the earlier analysis (Heger et al. Reference Heger, Cumming, Galloway and Woosley2007). 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. (Reference Fujimoto, Hanawa and Miyaji1981) [or alternatively, ‘prompt mixed’ bursts in regime 3 of Narayan & Heyl (Reference Narayan and Heyl2003)].
2.2. SAX J1808.4 − 3658
This source was discovered as a new transient in 1996 with BeppoSAX (in ’t Zand et al. Reference in ’t Zand, Heise, Muller, Bazzano, Cocchi, Natalucci and Ubertini1998), and has been active every 2–3 yrs 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 Reference Wijnands and van der Klis1998). 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. Reference Chakrabarty, Morgan, Muno, Galloway, Wijnands, van der Klis and Markwardt2003).
The transient outbursts typically reach a maximum inferred accretion rate of only
$0.05\ \dot{m}_{\rm Edd}$
, and the bursts have fast ( ≲ 0.5 s) rises indicative of He-rich fuel. The wide (2.01 h) orbit can accommodate a (H-rich) companion (Chakrabarty & Morgan Reference Chakrabarty and Morgan1998), and it is thought that the accreted H is burned away by β-limited hot CNO burning prior to ignition. The ≈ 20–30 h recurrence times, given the low accretion rates, are sufficient provided the metallicity, Z
cno is not far below the solar value. This type of burst corresponds to ‘case 2’ of Fujimoto et al. (Reference Fujimoto, Hanawa and Miyaji1981), or ‘delayed helium’ (regime 4) of Narayan & Heyl (Reference Narayan and Heyl2003).
Comparisons of the burst properties with the predictions of numerical ignition models indicate a range of possible compositions, with Z cno∝X (Galloway & Cumming Reference Galloway and Cumming2006). The lowest plausible H-fraction in the accreted fuel was X 0 ≈ 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. Reference Grindlay, Gursky, Schnopper, Parsignault, Heise, Brinkman and Schrijver1976). When in the bursting state, the source exhibits quasi-regular radius-expansion bursts with recurrence times in the range 2–4 h (e.g. Cumming Reference Cumming2003). The source intensity is modulated on a ≈ 176-d periodicity (Priedhorsky & Terrell Reference Priedhorsky and Terrell1984), and it only exhibits bursts in the (relatively short) low-intensity phase of the cycle, when the inferred accretion rate is ≈ 0.2–0.95
$\dot{M}_{\rm Edd}$
(e.g. Galloway et al. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008).
The orbital period of 685 s is known from periodic modulation of both persistent X-rays and the optical counterpart (King & Watson Reference King and Watson1986; Stella, White, & Priedhorsky Reference Stella, White and Priedhorsky1987), 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 ~ 0.1) amount of H in the surface layers, and this cannot be ruled out from comparing the observed bursts with ignition models (Cumming Reference Cumming2003).
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. Reference Özel, Psaltis, Güver, Baym, Heinke and Guillot2016). Although unresolved systematic errors may yet affect such measurements (e.g. Steiner et al. Reference Steiner, Lattimer and Brown2013), 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 h orbit (van Paradijs et al. Reference van Paradijs1990). 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, van Paradijs, & Taam Reference Lewin, van Paradijs and Taam1993; Galloway et al. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008). The neutron star spin has been measured from burst oscillations and transient pulsations during a superburst at 579.3 Hz (Zhang et al. Reference Zhang, Lapidus, Swank, White and Titarchuk1997; Strohmayer & Markwardt Reference Strohmayer and Markwardt2002). The relatively wide orbit, spectral features from hydrogen (Augusteijn et al. Reference Augusteijn, van der Hooft, de Jong, van Kerkwijk and van Paradijs1998), 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 × 10−9ergcm−2s−1, indicating an accretion rate of 11–21%
$\dot{m}_{\rm Edd}$
(cf. with Galloway et al. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008). 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 occurring in trains of up to 4 (Keek et al. Reference Keek, Galloway, Zand and Heger2010).
The first superbursts were detected from RXTE/ASM data (Wijnands Reference Wijnands2001), and two additional events were subsequently reported, one also detected by the PCA (Kuulkers et al. Reference Kuulkers, in ’t Zand, Homan, van Straaten, Altamirano, van der Klis, Kaaret, Lamb and Swank2004; Kuulkers Reference Kuulkers2009). All the superbursts were detected while the system was in its 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 Reference Strohmayer and Brown2002). Preliminary spectral fitting of the 4U 1636 − 536 burst indicated that the spectrum could not be adequately fitted with blackbody models (Kuulkers et al. Reference Kuulkers, in ’t Zand, Homan, van Straaten, Altamirano, van der Klis, Kaaret, Lamb and Swank2004), 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. Reference Keek, Ballantyne, Kuulkers and Strohmayer2014).
3 OBSERVATIONS AND ANALYSIS
We used preliminary results from the Multi-INstrument Burst ARchive (MINBARFootnote 2 ), which includes bursts observed by the Rossi X-ray Timing Explorer Proportional Counter Array (RXTE/PCA), BeppoSAX Wide-Field Camera (WFC; Jager et al. Reference Jager1997; in ’t Zand et al. Reference in ’t Zand, van den Heuvel, Wijers and Zand2004) and INTEGRAL Joint European X-Ray Monitor (JEM-X; Lund et al. Reference Lund2003).
We calculated burst lightcurves from time-resolved spectroscopic analysis of selected bursts observed by RXTE/PCA (Jahoda et al. Reference Jahoda, Swank, Giles, Stark, Strohmayer, Zhang and Morgan1996). The MINBAR analysis follows the approach of Galloway et al. (Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008), 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, Galloway, & Price Reference Worpel, Galloway and Price2013, Reference Worpel, Galloway and Price2015). While the fractional increase in the persistent flux can be ~ 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.
Figure 1. Example lightcurve showing bursts observed in 2000 September by BeppoSAX/WFC and the Rossi X-ray Timing Explorer/PCA from GS 1826 − 24. The top panel shows the lightcurves from the two instruments, rebinned to 10-s resolution, with the detected bursts marked with open symbols. The WFC lightcurve is shown in black, with open diamonds indicating the bursts; for PCA, the lightcurve is red, with open squares marking the bursts. The vertical lines above the symbols indicate the predicted times of bursts according to a constant recurrence time model; note that the unobserved bursts in the model train consistently fall in the data gaps. In the lower panel, we show the residual to the constant recurrence time model. Here, the best-fit recurrence time is 4.177 ± 0.010 h, with an RMS error of 9.53 min.
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. (Reference Thompson, Galloway, Rothschild and Homer2008).
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 ≈ 20 min or better.
From the set of bursts in each train N burst, we then selected the RXTE bursts (N av, Table 2), and created average lightcurves from each train. We also calculated averaged burst fluences, peak fluxes and α = ΔtF p c bol/E b, where F p is the persistent (accretion) flux, c bol the bolometric correction (see Section 3.2), and E b the burst fluence.
Table 2. Properties of thermonuclear bursts observed from target sources by the Rossi X-ray Timing Explorer.
Burst index number in the catalogue of Galloway et al. (Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008).
b Persistent flux in the energy range 3–25 keV.
c Inferred from ignition model comparisons.
d Calculated based on linear interpolation between observations falling within the burst interval.
e Only two bursts detected, so may not have been part of a regular train.
f Recurrence time upper limit of 1.75 yr based on the next earliest event, detected by the RXTE/ASM (Kuulkers et al. Reference Kuulkers, in ’t Zand, Homan, van Straaten, Altamirano, van der Klis, Kaaret, Lamb and Swank2004).
g Values taken from the ‘optimal’ fits of Keek et al. (Reference Keek, Ballantyne, Kuulkers and Strohmayer2014).
We estimated the accretion rate
$\dot{m}$
as a fraction of the Eddington rate
$\dot{m}_{\rm Edd}$
as
$$\begin{equation}
\dot{m}=4\pi d^2 F_{\rm p}c_{\rm bol}/\dot{m}_{\rm Edd},
\end{equation}$$
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 h.
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. Reference in’t Zand, Homan, Keek and Palmer2012) 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 Section 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 h. The nearest observations in time preceding (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 h does reflect a regular recurrence time, this pair would be the most frequent bursts ever observed from this system (cf. with Clark et al. Reference Clark, Li, Canizares, Hayakawa, Jernigan and Lewin1977).
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 yrs prior. However, the inferred accretion column is approximately twice the ignition column, suggesting that the actual recurrence time is closer to 0.9 yrs (Keek et al. Reference Keek, Cumming, Wolf, Ballantyne, Suleimanov, Kuulkers and Strohmayer2015). 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 quench ⩽ 22.8 d. The observation during which that event was detected commenced 17.5 h before the burst was detected, and had a duty cycle of ≈ 70%, with interruptions due to the low-Earth orbit. The coverages suggests a ≈ 30% chance that an earlier burst could have been missed, so that the quench time could have been shorter by up to 17.5 h. 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 broadband model typically comprising an absorbed Comptonisation component (compTT in xspec; Titarchuk Reference Titarchuk1994), usually with a Gaussian component representing Fe Kα 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–1 000 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–1 000 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. Reference Thompson, Galloway, Rothschild and Homer2008). In some observations, there is evidence for unusually low α-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 ( ≲ 1 keV) thermal component. Thus, we avoided from our data selection bursts from such epochs identified by Thompson et al. (Reference Thompson, Galloway, Rothschild and Homer2008), including those from 2003 April, and also the observations from 1998 June and 1997 November. An additional epoch of unusually low α ≈ 26 is from observations in 2006 August, and was also excluded. The remaining set of three burst trains have broadly consistent α-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. (Reference Heger, Cumming, Galloway and Woosley2007) 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 h, similar to that used previously. The seven 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 Figure 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. (Reference Lampe, Heger and Galloway2016).
Heger et al. (Reference Heger, Cumming, Galloway and Woosley2007) found excellent agreement between the observed lightcurve and a model referred to as ‘A3’, with solar composition, an accretion rate of 1.58 × 10−9M⊙yr−1. In the sample of Lampe et al. (Reference Lampe, Heger and Galloway2016), the accretion rate is quantified in terms of the Eddington value, defined as 1.75 × 10−8M⊙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. (Reference Lampe, Heger and Galloway2016) 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. Reference Foreman-Mackey, Hogg, Lang and Goodman2013) 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:
$$\begin{equation}
F_X = \frac{L_X}{4\pi d^2\xi _{\rm b}(1+z)},
\end{equation}$$
where LX is the burst luminosity predicted by kepler, d is the distance, ξb takes into account the possible anisotropy of the burst emission (see also appendix A, and FX 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. (Reference Heger, Cumming, Galloway and Woosley2007), gave good agreement with the observed lightcurves, although requiring a redshift of 1 + z = 1.392. Following the approach of Lampe et al. (Reference Lampe, Heger and Galloway2016), 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.4M⊙) 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. Reference Steiner, Lattimer and Brown2013).
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. (Reference Heger, Cumming, Galloway and Woosley2007) estimated the ratio of persistent to burst anisotropies at ξp/ξb = 1.55, implying ξp > 1 (e.g. Fujimoto Reference Fujimoto1988); 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ξ1/2 b = 6.1 kpc. We do not quote errors on these parameters, as the lightcurve comparison (Figure 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. (Reference Heger, Cumming, Galloway and Woosley2007), and so for the present purposes, we adopt the values for the redshift and distance.
Figure 2. Best-fitting comparison of averaged, observed burst lightcurves from GS 1826 − 24, with kepler model a028. The symbols with errors show the average lightcurve observed by RXTE in 2000 September, at which time the recurrence time was 4.177 h The model curve, rescaled on both axis based on the best-fit distance and redshift, is the red solid curve. The best-fitting distance and redshift, determined from an MCMC analysis, were dξ1/2 = 6.1 kpc and 1 + z = 1.23, respectively.
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 h, 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 ≈ 90-min low-Earth orbit.
Galloway & Cumming (Reference Galloway and Cumming2006) 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 α-measurements to the settle code developed by Cumming & Bildsten (Reference Cumming and Bildsten2000) using MCMC and marginalising over the accreted hydrogen fraction X
0, CNO metallicity Z
CNO, base flux Qb
, 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. Reference Woosley2004) 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., 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 ≈ 1 for the first burst to ≈ 1.2 for the final burst.
We infer the persistent and burst anisotropy factors (ξp and ξb, respectively) from the scaling factors and find ξp = 1.05 ± 0.01 and ξb = 0.91+0.04 −0.03. This implies an inclination of 62 ±1○. 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 (Reference Galloway and Cumming2006), 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. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008). 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 h for the second burst observed with RXTE is comfortably within the overall range for the source, with the minimum separation of 12.6 h set by a pair of bursts observed by BeppoSAX (in ’t Zand et al. Reference in ’t Zand, Heise, Muller, Bazzano, Cocchi, Natalucci and Ubertini1998).
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 material 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.
Figure 3. Representative burst lightcurves from 4U 1820 − 303, SAX J1808.4 − 3658, GS 1826 − 24, and 4U 1636 − 536 observed by the Rossi X-ray Timing Explorer. In each panel, we plot the lightcurve at the indicated recurrence time, with flux converted to luminosity at the distances indicated in Table 1. The bursts contributing to each profile are listed in Table 2; where more than one profile is observed by RXTE in the train, we show the average profile. Note the change in x- and y-axis between the top panels (4U 1820 − 303 and SAX J1808.4 − 3658) and the lower two panels.
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-χ2 values in the range 2–8). As
$\dot{m}$
increases, the burst timescale increases slightly (τ = 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 ( ≈ 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 h.
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 (Reference Galloway and Cumming2006) 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:
-
1. Set the surface gravity g and neutron star radius R to the values appropriate for the source of interest in Table 1.
-
2. Set the composition in the accreted material, X 0 and Z CNO.
-
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. 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 ξ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 Section 3.4 to account for the effects of the declining accretion rate during each burst.
-
5. Calculate from the simulated burst train the recurrence time Δt pred, and from the burst lightcurve, the burst fluence, peak flux, α-value, and timescale.
-
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 ξ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. Figure 18 of Cumming et al. Reference Cumming, Macbeth, Zand and Page2006); 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. Modelling 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. Reference Stevens, Brown, Cumming, Cyburt and Schatz2014; Keek & Heger Reference Keek and Heger2016). An alternative approach is to directly accrete a mixture of carbon and heavier elements (e.g. Fe or Ru; Cumming & Bildsten Reference Cumming and Bildsten2001; Schatz, Bildsten, & Cumming Reference Schatz, Bildsten and Cumming2003). For the 2001 superburst from 4U 1636 − 536, we recommend a fuel composition of 26% 12C and 74% 56Fe (Keek et al. Reference Keek, Cumming, Wolf, Ballantyne, Suleimanov, Kuulkers and Strohmayer2015). Such simulations are able to describe the observed superburst light curves in detail, but generally predict larger ignition depths than observed (Keek & Heger Reference Keek and Heger2011).
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. Reference Cyburt, Amthor, Heger, Johnson, Keek, Meisel, Schatz and Smith2016). Similarly, nuclear masses have also been shown to have an effect (Schatz & Ong Reference Schatz and Ong2016). 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 programme, and the Australian Research Council’s Discovery Projects and Future Fellowship funding schemes. LK is supported by NASA under award number NNG06EO90A.
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 ξb and ξ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. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008), these measurements implicitly also depend on the burst anisotropy. Strictly speaking, distances estimated in such a manner correspond to dξ1/2 b, and in general, it is not possible to separate these quantities to derive the physical distance d, independently.
The definition of the anisotropy factor ξb follows that of Fujimoto (Reference Fujimoto1988), such that for a total burst luminosity L b, the measured burst flux would be given by
$$\begin{equation}
F_{\rm b} = \frac{L_b}{4\pi d^2\xi _{\rm b}}.
\end{equation}$$
From the above relation, we understand that ξb = 1 corresponds to isotropic emission in all directions, while ξb < 1 indicates that the burst flux towards the observer is in excess of the isotropic value, i.e. preferential beaming in our direction. As pointed out by He & Keek (Reference He and Keek2016), the values of ξb and ξ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. (Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008), 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 Reference Fujimoto1988). Optical observations constrain the inclination angle to the range 36°–74° (Casares et al. Reference Casares, Cornelisse, Steeghs, Charles, Hynes, O’Brien and Strohmayer2006), which implies 0.7 ≲ ξ−1 b ≲ 1.3 for a thin disk (larger values for smaller angles; He & Keek Reference He and Keek2016). The observation of a broad fluorescent iron line favours a large inclination angle (Pandel, Kaaret, & Corbel Reference Pandel, Kaaret and Corbel2008), 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. Reference Keek, Cumming, Wolf, Ballantyne, Suleimanov, Kuulkers and Strohmayer2015). 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 Sections 3.3 and 3.4.
Zamfir et al. (Reference Zamfir, Cumming and Galloway2012) derived an upper limit on the distance (combined with the burst anisotropy) for GS 1826 − 24 of dξ1/2 b < 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. (Reference Zamfir, Cumming and Galloway2012), 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. (Reference Lampe, Heger and Galloway2016) 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. Reference Chenevez2016). Assuming that the degree of burst anisotropy ξ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. Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008) 3.7 × 10−8ergcm−2s−1. This is well within the confidence interval of (4.0 ± 0.3) × 10−8ergcm−2s−1 derived for the radius-expansion burst observed by Chenevez et al. (Reference Chenevez2016). 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. Reference Kuulkers, den Hartog, in ’t Zand, Verbunt, Harris and Cocchi2003).
