A Bayesian Approach to Matching Thermonuclear X-ray Burst Observations with Models

TL;DR

This paper introduces a Bayesian method to match thermonuclear X-ray burst observations with models, enabling detailed parameter inference for neutron star systems like SAX J1808.4--3658.

Contribution

It develops a Bayesian analysis framework to infer system parameters from X-ray burst data, accounting for inclination and other factors, applicable to any accreting source with Type I bursts.

Findings

Inferred neutron star radius of approximately 11.8 km.

Estimated fuel composition with hydrogen mass fraction around 0.57.

Determined system inclination angle to be about 69 degrees.

Abstract

We present a new method of matching observations of Type I (thermonuclear) X-ray bursts with models, comparing the predictions of a semi-analytic ignition model with X-ray observations of the accretion-powered millisecond pulsar SAX J1808.4--3658 in outburst. We used a Bayesian analysis approach to marginalise over the parameters of interest and determine parameters such as fuel composition, distance/anisotropy factors, neutron star mass and neutron star radius. Our study includes a treatment of the system inclination effects, inferring that the rotation axis of the system is inclined $\left(69^{+4}_{-2}\right)^\circ$ from the observers line of sight, assuming a flat disc model. This method can be applied to any accreting source that exhibits Type I X-ray bursts. We find a hydrogen mass fraction of $0.57^{+0.13}_{-0.14}$ and CNO metallicity of $0.013^{+0.006}_{-0.004}$ for the accreted fuel is required by the model to match the observed burst energies, for a distance to the source of $3.3^{+0.3}_{-0.2}\,\mathrm{kpc}$. We infer a neutron star mass of $1.5^{+0.6}_{-0.3}\,\mathrm{M}_{\odot}$ and radius of $11.8^{+1.3}_{-0.9}\,\mathrm{km}$ for a surface gravity of $1.9^{+0.7}_{-0.4}\times10^{14}\,\mathrm{cm}\,\mathrm{s}^{-2}$ for SAX J1808.4--3658.

Tables

Table 1. Table 1: Observations of the 2002 October outburst of SAX J1808.4--3658

No.	Burst start time	Peak Flux	Fluence	$\alpha$
	(MJD)	(10${}_{}{}^{-9}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$)	(${10}^{-6}\mathrm{erg}{\mathrm{cm}}^{-2}$)
1	52562.4136	163$\pm$3	2.62$\pm$0.02	–
2	52564.3051	177$\pm$3	2.65$\pm$0.01	107$\pm$2
3	52565.1843	179$\pm$3	2.99$\pm$0.01	118$\pm$2
4	52566.4268	177$\pm$3	3.46$\pm$0.02	128$\pm$2

Table 2. Table 2: SAX J1808.4--3658 derived neutron star parameters

Parameter	Value	Units
X	${0.58}_{-0.14}^{+0.13}$
Z	${0.013}_{-0.004}^{+0.006}$
$Q_b$	${0.4}_{-0.2}^{+0.3}$	$\mathrm{MeV}/\mathrm{nucleon}$
$M$	${1.5}_{-0.3}^{+0.6}$	$M_{\odot }$
$R$	${11.8}_{-1.0}^{+1.3}$	$\mathrm{km}$
${\dot{m}}_{\max }$	${0.037}_{-0.002}^{+0.002}$	${\dot{m}}_{\mathrm{Edd}}$
$g$	${1.9}_{-0.4}^{+0.7}$	${10}^{14}\mathrm{cm}{s}^{-2}$
$1+z$	${1.27}_{-0.05}^{+0.13}$
$d$	${3.3}_{-0.2}^{+0.3}$	$\mathrm{kpc}$
${\xi }_b$	${0.74}_{-0.10}^{+0.10}$
${\xi }_p$	${0.87}_{-0.10}^{+0.12}$
$\cos i$	${0.36}_{-0.04}^{+0.07}$

Table 3. Table 3: SAX J1808.4--3658 predicted burst parameters

Burst	$\overline{X}$	$\alpha$	$\Delta t$ (hours)
2	${0.20}_{-0.06}^{+0.06}$	${108}_{-9}^{+12}$	$14.7$
3	${0.20}_{-0.06}^{+0.06}$	${107}_{-9}^{+12}$	$14.4$
4*	${0.19}_{-0.06}^{+0.05}$	110${}_{-9}{}^{}{}_{}{}^{+11}$	16.1
5*	${0.16}_{-0.04}^{+0.04}$	${117}_{-9}^{+10}$	$19.4$
6*	${0.12}_{-0.03}^{+0.03}$	${126}_{-7}^{+9}$	$31.2$
7	${0.07}_{-0.02}^{+0.02}$	${139}_{-5}^{+6}$	$56$

Table 4. Table 4: SAX J1808.4--3658 derived neutron star parameters for a fixed mass and radius

Parameter	Value
$X$	${0.57}_{-0.15}^{+0.13}$
$Z$	${0.012}_{-0.003}^{+0.004}$
$Q_b$ ($\mathrm{MeV}/\mathrm{nucleon}$)	${0.6}_{-0.1}^{+0.1}$
$M$ ($M_{\odot }$)	$1.4$
$R$ ($\mathrm{km}$)	$11.2$
${\dot{m}}_{\max }$ (${\dot{m}}_{\mathrm{Edd}}$)	${0.035}_{-0.001}^{+0.002}$
$g$ (${10}^{14}\mathrm{cm}{s}^{-2}$)	$1.86$
$1+z$	$1.26$
$d$ ($\mathrm{kpc}$)	${2.9}_{-0.1}^{+0.1}$
${\xi }_b$	${0.9}_{-0.1}^{+0.1}$
${\xi }_p$	${1.04}_{-0.04}^{+0.05}$
$\cos i$	${0.35}_{-0.05}^{+0.08}$

Table 5. Table 5: SAX J1808.4--3658 derived neutron star parameters

Parameter	Value
$X$	${0.58}_{-0.15}^{+0.13}$
$Z$	${0.013}_{-0.004}^{+0.006}$
$Q_b$ ($\mathrm{MeV}/\mathrm{nucleon}$)	${0.4}_{-0.2}^{+0.2}$
$M$ ($M_{\odot }$)	${1.6}_{-0.3}^{+0.5}$
$R$ ($\mathrm{km}$)	${11.9}_{-1.2}^{+1.4}$
${\dot{m}}_{\max }$ (${\dot{m}}_{\mathrm{Edd}}$)	${0.037}_{-0.002}^{+0.002}$
$g$ (${10}^{14}\mathrm{cm}{s}^{-2}$)	${1.95}_{-0.4}^{+0.6}$
$1+z$	${1.29}_{-0.06}^{+0.1}$
$d$ ($\mathrm{kpc}$)	${3.3}_{-0.2}^{+0.2}$
${\xi }_b$	${0.74}_{-0.10}^{+0.1}$
${\xi }_p$	${0.87}_{-0.10}^{+0.1}$
$\cos i$	${0.37}_{-0.05}^{+0.07}$

Equations

$$\alpha =\frac{F_{\mathrm{p}}{c}_{\mathrm{bol}}\Delta t}{E_{\mathrm{b}}}$$

$$\dot{m}=c_{\mathrm{bol}}{F}_{\mathrm{p}}{r}_1$$

$$\text{centering}g=\frac{G{M}_{\mathrm{NS}}}{R_{\mathrm{NS}}^2}\cdot (1+z)\text{@add@centering}$$

$$\frac{\mathrm{d}{Y}_{{}^1\mathrm{H}}}{\mathrm{d}t}=\frac{4\ln 2Y_{\mathrm{CNO}}}{t_{\frac{1}{2},\mathrm{CNO}}},$$

$$Y_{\mathrm{CNO}}=\frac{Z_{\mathrm{CNO},\beta }}{{\bar{A}}_{\mathrm{CNO},\beta }}$$

$${\bar{A}}_{\mathrm{CNO},\beta }=\frac{t_{\frac{1}{2}{,}^{15}\mathrm{O}}{A}_{{}^{15}\mathrm{O}}+t_{\frac{1}{2}{,}^{14}\mathrm{O}}{A}_{{}^{14}\mathrm{O}}}{t_{\frac{1}{2}{,}^{15}\mathrm{O}}+t_{\frac{1}{2}{,}^{14}\mathrm{O}}},$$

$${ϵ}_{\mathrm{H}}=\frac{\mathrm{d}}{\mathrm{d}t}{Y}_{{}^1\mathrm{H}}{E}_{\mathrm{H}},$$

$${ϵ}_{\mathrm{H}}=5.94\times 10^{13}\frac{Z_{\mathrm{CNO},\beta }}{0.01}\mathrm{erg}{\mathrm{g}}^{-1}{\mathrm{s}}^{-1},$$

$${ϵ}_{\mathrm{H}}=5.93\times 10^{13}\frac{Z_{\mathrm{CNO},\mathrm{scaled}\odot }}{0.01}\mathrm{erg}{\mathrm{g}}^{-1}{\mathrm{s}}^{-1}$$

$$t_{\mathrm{H}}=\frac{Y_{{}^1\mathrm{H}}}{Y_{\mathrm{CNO}}}\frac{t_{\frac{1}{2}{,}^{15}\mathrm{O}}+t_{\frac{1}{2}{,}^{14}\mathrm{O}}}{4\ln 2}=20.45\left(\frac{X_0}{0.7}\right)\left(\frac{0.01}{Z_{\mathrm{CNO}}}\right)\mathrm{h}$$

$$y_{\mathrm{d}}=6.9\times 10^8\left(\frac{\dot{m}}{0.01{\dot{m}}_{Edd}}\right)\left(\frac{0.01}{Z_{\mathrm{CNO}}}\right)\left(\frac{X_0}{0.7}\right)\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}$$

$$\text{centering}L(x_i|\theta ,\mathrm{model})=\frac{1}{{\sigma }_i{f}_i\sqrt{2\pi }}\exp \left[\frac{-(x_i-\mathrm{model}{)}^2}{2f_i{\sigma }_i}\right]\text{@add@centering}$$

$$r_1=\frac{{\dot{m}}_i}{c_{\mathrm{bol}}{F}_{\mathrm{p},i}{\xi }_{\mathrm{p}}}$$

$$r_2=\frac{\alpha }{{\alpha }_{\mathrm{pred}}}$$

$$r_3=\frac{E_{\mathrm{b}}/10^{-9}}{L_{\mathrm{b},\mathrm{pred}}{\xi }_{\mathrm{b}}}$$

$$[Fe/H]\sim \varphi (3+\frac{3}{4})-3$$

$$r_2=0.74816\left(\frac{M_{\mathrm{NS}}}{1.4{\mathrm{M}}_{\odot }}\right){\left(\frac{R_{\mathrm{NS}}}{11.2\mathrm{km}}\right)}^{-1}{\left(\frac{Q_{\mathrm{nuc}}}{Q_{\mathrm{nuc},\mathrm{pred}}}\right)}^{-1}\left(\frac{{\xi }_{\mathrm{p}}}{{\xi }_{\mathrm{b}}}\right)$$

$$r_3=63.23{\left(\frac{d}{10\mathrm{kpc}}\right)}^{-2}{\left(\frac{R_{\mathrm{NS}}}{11.2\mathrm{km}}\right)}^2{\left(\frac{1+z}{1.258}\right)}^{-1}\left(\frac{Q_{\mathrm{nuc}}}{Q_{\mathrm{nuc},\mathrm{pred}}}\right){\xi }_{\mathrm{b}}$$

$${\xi }_{\mathrm{p}}=r_1{\left(\frac{10}{d}\right)}^2$$

$${\xi }_{\mathrm{b}}=\frac{r_2{\xi }_{\mathrm{p}}}{0.74816}$$

$$d=10\sqrt{\frac{63.23{\xi }_{\mathrm{b}}}{r_3}}\mathrm{kpc}$$

$$\frac{{\xi }_{\mathrm{p}}}{{\xi }_{\mathrm{b}}}=\frac{0.5+|\cos i|}{2|\cos i|}$$

$$f_x=\frac{(M_{\mathrm{c}}\sin i{)}^3}{(M_{\mathrm{NS}}+M_{\mathrm{c}}{)}^2}=\frac{4{\pi }^2(a_x\sin i{)}^3}{G{P}_{\mathrm{orb}}^2}$$

Full text

††pubyear: 2019††pagerange: A Bayesian Approach to Matching Thermonuclear X-ray Burst Observations with Models–A.2

1 Introduction

Accretion-powered millisecond pulsars (AMSPs) are neutron stars with weak ($\sim 10^{8}$ G) magnetic fields that are accreting matter from a low mass companion star (e.g., Lewin et al., 1993; Chakrabarty & Morgan, 1998; Wijnands & van der Klis, 1998; Patruno & Watts, 2012; Galloway & Keek, 2017). These systems are typically transient X-ray sources, and exhibit coherent X-ray pulsations as well as Type I thermonuclear X-ray bursts periodically (e.g., Wijnands, 2004). A source goes into outburst when the accretion disc around the neutron star, formed by accretion of matter onto the companion star, reaches a critical state that causes matter to transfer directly onto the neutron star and ignites thermonuclear bursts (e.g., Strohmayer & Bildsten, 2006; Galloway & Keek, 2017). Type I X-ray bursts are the most frequently observed thermonuclear powered outbursts in nature and provide a unique insight into the nuclear reactions that can occur in extreme environments. They are characterised by a sudden increase in the X-ray flux of a source in outburst (e.g., Galloway et al., 2008). They can be observed from neutron stars in binary systems, and are thought to be caused by unstable ignition of hydrogen or helium on the surface of an accreting neutron star (e.g., Strohmayer & Bildsten, 2006).

The challenge inherent in understanding low-mass neutron star binary systems is estimating system parameters, including the composition of the fuel that drives the thermonuclear burst; and the distance/inclination of the source, which cannot be directly observed as they are kilometre-sized objects at kilo-parsec distances. These parameters are an essential ingredient in modelling and understanding the dense matter of neutron stars and the physical processes behind thermonuclear bursts. Applications range from stellar evolution (with fuel composition of the bursts enabling constraints on the evolution of the companion star and population synthesis) to constraining proton- and carbon-rich nuclear reactions that are difficult to replicate in Earth-bound laboratories.

Neutron stars are the densest directly observable stellar objects known. The study of neutron stars through various electromagnetic bands has seen much progress in recent years, furthering the understanding of high-density cold matter and the dense equation of state (e.g. Miller et al., 2019; Degenaar & Suleimanov, 2018; Özel & Freire, 2016) The recent observation of gravitational waves from merging neutron stars (Abbott et al., 2017) has also seen an increase in efforts to constrain neutron star properties, such as the neutron star mass, radius, spin, and equation of state (e.g., Abbott et al., 2018; Annala et al., 2018; Radice et al., 2018). Constraining these properties is fundamental in the understanding of ultra-dense matter, such as the material that makes up a neutron star, as well as providing additional constraints in the search for more gravitational waves.

Observations of AMSPs that exhibit Type I X-ray bursts are becoming increasingly more accessible, with the Multi-INstrument Burst Archive (MINBAR)111https://burst.sci.monash.edu/minbar/ providing an extensive collection of processed and anlaysed X-ray burst data. Observations are limited to the X-ray activity of these sources, providing measurements of the accretion flux, burst flux, and recurrence time of bursts.

Efforts to model burst sources have evolved from early models (e.g., Fujimoto et al., 1987; Taam, 1980; van Paradijs et al., 1988) that focused on the correlation between burst energies and recurrence times, to varying success, to more recent models (e.g., Cumming, 2003; Galloway et al., 2004; Woosley et al., 2004) which have a heavier focus on the nuclear physics driving the bursts, predicting fuel compositions and accretion rates by implementing more detailed models of the nuclear reaction rates that produce the observed energy generation rates.

These observations and models motivate an effort to combine the data obtained from both, to obtain a deeper understanding of the thermonuclear processes that occur during an X-ray burst, as well as identifying where current models are failing to predict observed properties (e.g., Galloway et al., 2017). The most well-studied transients, such as the accretion-powered millisecond pulsar SAX J1808.4--3658, provide crucial test cases for numerical models of thermonuclear bursts and can be used as standard cases in modelling studies.

Galloway & Cumming (2006) showed that the bursts from SAX J1808.4--3658 were consistent with ignition in an almost pure He environment. The long recurrence times of $\approx$1 day allow enough time for hydrogen to be depleted by thermally-stable hot CNO burning between bursts. By modelling the burst properties, they were able to constrain the hydrogen fraction in the accreted layer, and showed that the observed recurrence times imply a relation between the hydrogen fraction and fraction of CNO elements. As such, SAX J1808.4--3658 provides the best example of pure helium ignited bursts, and validates the theory in which hydrogen depletes before ignition of the X-ray burst in a helium environment.

In this study we develop a new method for matching burst observations with a semi-analytic ignition model, extending the previous work of Galloway & Cumming (2006), to determine system parameters such as fuel composition, neutron star mass and radius, distance, and inclination. We use a Markov Chain Monte Carlo approach to find the best-fitting system parameters, with uncertainties. In Section 2 we describe the semi-analytic model used in this study, compare it to detailed time-dependent simulations, and describe the MCMC implementation. In Section 3 we report the parameter limits predicted by the model, and describe how well our model matches the observations. Finally, in Section 4, we summarise the results and give our conclusions.

2 Methods

We used the ignition model of Cumming & Bildsten (2000) (hereafter called settle) to generate burst sequences, which were then compared to Rossi X-ray Timing Explorer (RXTE) Proportional Counter Array (PCA) telescope observations of the 2002 outburst of SAX J1808.4--3658 using a Markov Chain Monte Carlo (MCMC) algorithm. We obtained the ‘best’ model run that matched the observations, providing predictions for parameters such as fuel composition, distance, anisotropies, neutron star mass and radius, that cannot be measured otherwise. We used the python implementation of MCMC, emcee (Foreman-Mackey et al., 2013). RXTE observations were obtained from archival data available in MINBAR and shown in Table 1, updated since the initial analysis of this outburst by Galloway & Cumming (2006). Observed properties of the bursts include times, persistent flux, burst fluence and $\alpha$ values. Here, $\alpha$ is defined as the ratio of integrated persistent flux to the fluence, given by

[TABLE]

where $E_{\mathrm{b}}$ is the burst fluence, $\Delta t$ is the recurrence time, $F_{\mathrm{p}}$ is the persistent flux, and $c_{\mathrm{bol}}$ is the bolometric correction. For $c_{\mathrm{bol}}$ we adopted a constant value of $2.12$, as in Galloway & Cumming (2006).

In order to determine the accretion rate at the time of a burst, we used persistent flux measurements of the outburst from the PCA and linearly interpolated the values between measurements either side of the burst. This enabled us to obtain a value at the time of each burst, from which an accretion rate was inferred using

[TABLE]

where $F_{\mathrm{p}}$ is the persistent flux, $\dot{m}$ is the accretion rate, and $r_{\mathrm{1}}$ is the scaling factor between the accretion rate and the persistent flux, described in Section 2.3.

2.1 settle

settle is a semi-analytic ignition model developed by Cumming & Bildsten (2000) that integrates the thermal profile of a settling atmosphere to find ignition conditions for a single burst, from which the burst train can be predicted. It applies a one-zone ignition criterion to simple models of the accumulating layer, which allows a survey of parameter space whilst also enabling estimation of the ignition depth. Researchers have applied settle to observations of regular Type I bursters, such as 4U 1820--30 (Cumming, 2003), GS 1826--24 (Cumming, 2003), and SAX J1808.4--3658 (Galloway & Cumming, 2006). The settle model assumes the accumulating fuel layer is heated by hot CNO hydrogen burning within the layer, and a constant flux from deeper layers. It does not include heating from energy released in a previous burst, or the possibility of ignition of unburnt fuel within the ashes of a previous burst (commonly referred to as thermal and compositional inertia; Woosley et al. 2004; see Section 2.2 for corrections we applied to settle to account for some of these simple assumptions). Assuming spherical symmetry, settle does not address complex physics that is currently poorly understood, such as how the burning front spreads across the star after ignition occurs, what determines the number of hotspots on the surface of the star, or the cause of asymmetry near the end of bursts as the material is cooling (e.g., Cavecchi & Spitkovsky, 2019).

settle adopts a neutron star mass and radius (e.g., $M_{\mathrm{NS}}=1.4\,\mathrm{M}_{\odot}$ and $R_{\mathrm{NS}}=11.2\,\mathrm{km}$) which then sets the gravitational redshift and surface gravity, where the surface gravity, $g$, in the X-ray burst layer is given by

[TABLE]

where G is the Gravitational constant, $z$ is the gravitational redshift, $M_{\mathrm{NS}}$ is the gravitational neutron star mass and $R_{\mathrm{NS}}$ is the neutron star radius in the observer frame.

An X-ray burst is triggered when helium burning becomes unstable at the base of the accumulated layer, at a temperature of $\approx 2\times 10^{8}$ K and a density range of $\sim 10^{5}$–$10^{6}\,\mathrm{g}\,\mathrm{cm}^{-3}$. The temperature profile of the accumulating layer of hydrogen and helium is calculated given the accretion rate, and the accumulating layer is allowed to build up until the one-zone condition for a thermal runaway is satisfied at the base of the layer. This condition is met when $\mathrm{d}\epsilon_{3\alpha}/\mathrm{d}T>\mathrm{d}\epsilon_{\mathrm{cool}}/\mathrm{d}T$, where $\epsilon_{3\alpha}$ is the energy production rate of helium burning (via the triple alpha reaction) and $\epsilon_{\mathrm{cool}}$ is the effective local cooling rate. At this point, a thermal instability occurs and the conditions for thermonuclear runaway are met. The temperature profile at ignition is determined such that $\mathrm{d}\epsilon_{3\alpha}/\mathrm{d}T=\mathrm{d}\epsilon_{\mathrm{cool}}/\mathrm{d}T$ at the base of the ignition column.

The composition of the ignition column at ignition depends on how much hydrogen has burned before the ignition conditions were met, during the accumulation of the fuel. The amount of hydrogen burned is determined by the local accretion rate, since hydrogen is assumed to burn steadily at a constant rate via the hot CNO cycle. The rate at which hydrogen burns is thus given by

[TABLE]

where $Y_{\mathrm{CNO}}$ is the parts per nucleon of CNO elements and $t_{\mathrm{\frac{1}{2},CNO}}$ is the cycle half life, i.e., the sum of the half lives of the isoptopes limiting the cycle. In settle, we assume the CNO elements are the $\beta$-limited CNO cycle isotopes, 14O and ${}^{15}\mathrm{O}$, in their cycle equilibrium abundance values and $Y_{\mathrm{CNO}}$ can be computed from

[TABLE]

where $\bar{A}_{\mathrm{CNO,\beta}}$ is the average mass number of the CNO isotopes in the $\beta$-limited CNO cycle, given by

[TABLE]

where $A_{i}$ is the mass number of species $i$, $t_{\frac{1}{2},^{15}\mathrm{O}}=122\,\mathrm{s}$ and $t_{\frac{1}{2},^{14}\mathrm{O}}=71\,\mathrm{s}$.

The ratio of the parts per nucleon of 14O to 15O is given by the ratio of their half lives. Thus the mass fractions $Z_{\mathrm{14}}=0.352\,Z_{\mathrm{CNO}}$ and $Z_{\mathrm{15}}=0.648\,Z_{\mathrm{CNO}}$, giving $\bar{A}_{\mathrm{CNO},\beta}=14.653\,\mathrm{g}\,\mathrm{mol}^{-1}$.

The energy production rate of hydrogen burning via the CNO cycle, $\epsilon_{\mathrm{H}}$, is given by

[TABLE]

where $E_{\mathrm{H}}$ is the energy produced from burning $4\,^{1}\mathrm{H}$ to ${}^{4}\mathrm{He}$. $E_{\mathrm{H}}$ is given by the change in mass excess less the energy carried away by neutrinos, $E_{\mathrm{H}}=6.03\times 10^{18}\,\mathrm{erg}\,\mathrm{g}^{-1}$. Thus, $\epsilon_{\mathrm{H}}$ is given by

[TABLE]

where $Z_{\mathrm{CNO,\beta}}$ is the mass fraction of CNO nuclei in the $\beta$-limited CNO cycle, computed as outlined above.

We also compare to the Sun as a reference point with its own, different, distribution of the CNO isotopes, with $\bar{A}_{\mathrm{CNO},\odot}=14.688$, such that a given CNO metallicity translates into a slightly different value of $Y_{\mathrm{CNO}}$. Using the solar abundance ratios for CNO isotopes from Lodders et al. (2009), one obtains:

[TABLE]

which is remarkably similar to the value using the $\beta$-limited $Z$ value reference, though only by coincidence, not by design.

Throughout the remainder of this paper all values of metallicity refer exclusively to the $\beta$-limited CNO metallicity, which is the relevant value for settle. The reference value of $0.01$ used here is indeed very close to the solar mass fraction of CNO isotopes given by Lodders et al. (2009) as $0.01001$.

Summarising, the time to burn all of the hydrogen ($t_{\mathrm{H}}$) is dependent only on the metallicity and initial hydrogen abundance. The time to burn hydrogen is thus

[TABLE]

where $X_{0}$ is the initial hydrogen mass fraction of the accreted fuel. This sets the minimum recurrence time required for pure helium bursts. Note that we have updated this equation from the original work of Cumming & Bildsten (2000), and it is the time to burn hydrogen locally on the neutron star surface, for the observer this timescale is longer by a redshift factor $1+z$. For higher accretion rates ($\dot{m}\gtrsim 0.04\,\dot{m}_{\mathrm{Edd}}$222$\dot{m}_{\mathrm{Edd}}$ is the Eddington accretion rate, given by $\dot{m}_{\mathrm{Edd}}=1.75\frac{0.7}{1+X}\times 10^{-8}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, where $X$ is the H mass fraction) helium burning becomes unstable before sufficient time has passed for all of the hydrogen to be burned, resulting in ignition occurring in a hydrogen rich environment. For lower ($\dot{m}\lesssim 0.04\,\dot{m}_{\mathrm{Edd}}$) accretion rates, there is sufficient time for all of the hydrogen to be burned before ignition, allowing a pure helium layer to accumulate before ignition, and a helium rich environment for the burst.

The column depth at which hydrogen runs out ($y_{\mathrm{d}}$) is given by:

[TABLE]

If helium ignites at a column depth $y<y_{\mathrm{d}}$, a mixed H/He burst occurs; and if helium ignites at a column depth $y>y_{\mathrm{d}}$, a pure helium layer has time to accumulate and ignition occurs in a pure helium environment.

Once ignition conditions are met, the burst parameters are calculated based on the total mass in the ignition column and its composition. The burst energy is calculated by assuming complete burning of the H/He layer, and the accreted material is assumed to cover the whole surface of the neutron star. The thickness of the accumulated layer at ignition determines the recurrence time of the bursts.

2.1.1 Updates to settle

When hydrogen remains at the time of ignition, and bursts ignite in a mixed H/He environment, there is extra energy released by the triple-alpha process due to proton capture by 12C. To avoid a discontinuity in accretion rate–recurrence time parameter space at the transition zone between mixed H/He bursts and pure He bursts, we gradually switched on the factor that multiplies the triple alpha reaction energy generation rate to account for the ${}^{12}\mathrm{C}(\mathrm{p},\gamma)$ reaction when protons begin running out at the critical H fraction of $X=1/7$ (since the ratio of protons to alpha particles is $2:3$). The gradual switch on of this multiplicative factor was not a feature of Cumming & Bildsten (2000)’s original code.

We updated some of the energy and rate values used in settle. We corrected the energy generation of hydrogen burning via the CNO cycle such that $\epsilon_{\mathrm{H}}=5.94\times 10^{15}\,\mathrm{erg}\,\mathrm{g}^{-1}\,\mathrm{s}^{-1}$ (used to be $5.8\times 10^{15}\,\mathrm{erg}\,\mathrm{g}^{-1}\,\mathrm{s}^{-1}$). We also used an updated version for the total nuclear energy released per unit mass in an X-ray burst as calibrated by Goodwin et al. (2018) using Kepler models, $Q_{\mathrm{nuc}}=1.35+6.05\,\bar{X}\,\mathrm{MeV}/\mathrm{nucleon}$, where $\bar{X}$ is the average hydrogen mass fraction of the ignition column. The ignition depth correction we applied due to differences between settle and Kepler is described in Section 2.2

Using Kepler models with fully time-dependent accretion rates, Johnston et al. (2018) found that using constant accretion rates (averaged between bursts) can systematically underestimate the burst recurrence times when the slope of the accretion rate is negative (see Figure 4 of that paper). Because settle implicitly uses averaged accretion rates, we accounted for this effect by multiplying each predicted recurrence time by a correction factor, as determined by the linear approximation from Johnston et al. (2018). These correction factors ranged from approximately $0.9$–$1.0$ for our modelled bursts. The correction factors we adopted are specific to SAX J1808.4–3658 and will either need to be recalculated or made more generally applicable if this method is applied to a different system. To rigorously test the influence of these correction factors we ran the code without the scaling factors and found that while the central value of the distance, hydrogen fraction, metallicity, and base flux shifted, the parameters still agreed within 1-$\sigma$ uncertainty of the model run with the scaling factors.

2.2 Code Verification

A simple code such as settle can be run in seconds, enabling multiple instances of the code to be run, and a marginalising method such as MCMC possible. Other codes, such as Kepler (Woosley et al., 2004), use bigger nuclear reaction networks and a more complex treatment of the neutron star atmosphere during outburst, at the cost of computational expense and runtime. Kepler solves the full time-dependent stellar structure equations, following the evolution of temperature and composition at each depth in the neutron star envelope over many cycles of bursts. To assess the limitations of settle, we ran a comparison with some Kepler models, and implemented corrections to settle to account for some of the differences we found. Woosley et al. (2004) provided a limited comparison between some settle and Kepler models, and found that, in general, the agreement was good, although with settle consistently overestimating the burst energies and recurrence times. The authors attribute the discrepancy between the settle and Kepler models to the lack of compositional inertia in the settle models.

We compared the predicted fluences ($E_{\mathrm{b}}$) and recurrence times ($\Delta t$) of the two codes for a range of hydrogen mass fractions ($X$), CNO mass fractions ($Z$), accretion rates ($\dot{m}$) and base fluxes ($Q_{\mathrm{b}}$). In order to allow meaningful comparison, we redshifted the Kepler recurrence times by multiplying by $1+z$, as Kepler calculates recurrence time in the reference frame of the neutron star, whereas settle calculates recurrence time in the reference frame of a distant observer. In both codes we fixed $M_{\mathrm{NS}}=1.4\,\mathrm{M}_{\odot}$ and $R_{\mathrm{NS}}=11.2\,\mathrm{km}$. We also increased the settle $Q_{\mathrm{b}}$ value by $0.1\,\mathrm{MeV}/\mathrm{nucleon}$ when comparing to Kepler to account for thermal inertia (Cumming & Bildsten, 2000).

On initial comparison we found the settle recurrence times and fluences were at least $\sim$ 30$\%$ higher than the Kepler predictions. We compared the temperature profiles of the ignition column for each code, shown in Figure 1.

We found the temperature profile of settle was overall cooler than Kepler, causing a longer recurrence time, deeper ignition, and thus brighter bursts. Since Kepler provides a more complex analysis of the underlying physics, we decided to correct the ignition depth of settle such that it matched the Kepler ignition depth. To do this we multiplied the settle ignition depth by a factor of $0.65$. Reducing the ignition depth reduced the recurrence time and fluences of the settle bursts and produced the comparison shown in Figure 2.

In Figure 2 for higher hydrogen fractions ($X=0.6$ and $0.7$) there is a significant difference between the predicted fluences and recurrence times of the settle and Kepler bursts across all accretion rates. There is much better agreement at lower hydrogen fractions ($X<0.5$) and lower accretion rates ($\dot{m}<0.1\,\dot{m}_{\mathrm{Edd}}$ between the two codes. For the purpose of this paper we note that for the source we are investigating, SAX J1808.4--3658, observations indicate that is has a slightly depleted hydrogen fraction and low ($\approx 4\%$ Eddington) accretion rate, so it is within the parameter space that settle agrees well with Kepler.

Fundamental differences between the two models can be attributed primarily to the fact that settle has a one-zone ignition criterion which uses a local approximation for the cooling, whereas Kepler follows the detailed cooling as a function of depth. Other differences could arise due to the fact that settle assumes there is no triple-alpha burning during the accumulation of fuel, settle assumes that all of the accreted fuel is burnt in a single burst, and the burst train is generated without accounting for subsequent bursts occurring in the ashes of previous bursts. Finally, settles uses an approximate formula for the nuclear energy generation of the bursts, that has been calibrated to Kepler predictions, but still is not the same as calculating the full nuclear reaction network, with thermal inertia included. As a caveat we note that calibrating to kepler introduces a different kind of model bias. Since there have been no comprehensive studies done between kepler and other similar codes, we do not understand how accurately kepler models the true physics of these systems.

2.3 MCMC

We used a Python implementation of MCMC (emcee; Foreman-Mackey et al., 2013) to estimate posterior distributions of the properties predicted by settle, based on the known values of the observed properties of the burst train (time, $E_{\mathrm{b}}$ and $\alpha$). We defined a likelihood function, given by Equation 11, and prior ranges for each parameter. We used a simple Gaussian likelihood function where the variance was underestimated by some fractional amount, $f$, and the overall likelihood is the sum of the individual likelihoods for each parameter $i$,

[TABLE]

where $\sigma_{i}$ is the observed uncertainty, $x_{i}$ is the observed parameter value (time, fluence, and alpha), model is the parameter value predicted by the model, and $\theta$ are our parameters of interest. The parameters of interest we defined are hydrogen mass fraction $X$, metallicity $Z$, which is defined as the mass fraction of CNO elements (all in 14O and 15O), base flux $Q_{\mathrm{b}}$, neutron star mass $M$, neutron star radius $R$ and three scaling factors, $r_{\mathrm{1}}$, $r_{\mathrm{2}}$, and $r_{\mathrm{3}}$.

The 3 scaling factors, $r_{\mathrm{1}}$, $r_{\mathrm{2}}$, and $r_{\mathrm{3}}$ are defined as in Galloway & Cumming (2006), given in Equations 12, 13, and 14,

[TABLE]

where $\dot{m}_{i}$ is the accretion rate between burst $i-1$ and burst $i$, $\xi_{\mathrm{p}}$ is the persistent anisotropy (defined in Section 3.1, $F_{\mathrm{p},i}$ is the mean $2-25\,\mathrm{keV}$ persistent flux between burst $i-1$ and $i$ (in units of $10^{-9}\,\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$, from interpolation) and $c_{\mathrm{bol}}$ is the bolometric flux correction.

[TABLE]

where $\alpha$ is the observed value (Equation 1) and $\alpha_{\mathrm{pred}}=290/Q_{\mathrm{nuc,pred}}$.

[TABLE]

where $L_{\mathrm{b,pred}}$ is the predicted burst luminosity in units of $10^{39}\,\mathrm{erg}\,\mathrm{s}^{-1}$ and $E_{\mathrm{b}}$ is the observed burst fluence.

We initialised the parameters in a tiny ($\sim 10^{-3}$) Gaussian ball around the initial value. We ran the MCMC chains with 300 walkers for 2,000 steps, and discarded the first 200 steps in each chain to account for burn-in. See Section 3.4 for information on the tests we ran to ensure the MCMC chains were fully converged and independent of initial walker positions.

2.4 Priors

For each parameter except metallicity, mass and radius, we used a uniform prior, where we chose an acceptable range based on observations or theoretical assumptions. For metallicity ($Z_{\mathrm{CNO}}$) we chose a more informed prior, using a synthetic survey of the Milky Way in the location of SAX J1808.4--3658 at $2$–$4\,\mathrm{kpc}$ generated by the Galaxia code (Sharma et al., 2011) to determine a likely distribution of the metallicity, given by Equation 15 and shown in Figure 3.

[TABLE]

where $\phi\sim\mathcal{B}(10,3)$ is a Beta distribution with $\alpha=10$ and $\beta=3$, and the linear translation limits the metallicity prior to have support between $[\textrm{Fe/H}]=-3$ and $+0.75$ (i.e., $0.00001<Z_{\mathrm{CNO}}<0.056$).

For mass and radius we also imposed a more informed prior based on the models of Steiner et al. (2018), who used observations of eight quiescent low mass X-ray binaries in globular clusters to determine the neutron star mass–radius curve and the equation of state. Steiner et al. (2018) used a fixed mass grid and varied radius according to their chosen equation of state and observations of radii. We chose the simplest of their models (‘baseline’) and used the probability distribution given in Figure 4 for the relationship between mass and radius. We also constrained mass and radius to: $1.15\,\mathrm{M_{\odot}}<M<2.5\,\mathrm{M_{\odot}}$ and $9\,\mathrm{km}<R<17\,\mathrm{km}$. This prior effectively checks if the masses and radii predicted by the model are consistent with the expected mass-radius relation of low mass X-ray binaries in globular clusters. It provides a probability of the predicted radius given the predicted mass.

3 Results

We report the posterior limits for the predicted system parameters in Table 2, reporting the mean value and the 68th percentile limits of the distributions for the uncertainties. Additional model runs are shown in the Appendix, including one that replicates Galloway & Cumming (2006) with fixed mass and radius and one with a uniform prior in mass and radius.

The times and fluences of the 6 predicted bursts are plotted in Figure 5, as well as the 4 observed bursts. We predict 3 bursts that were not observed, but that cannot be ruled out by observations, since they occur during periods that the telescope was not observing. The two bursts predicted between the first and second observed bursts were also predicted by Galloway & Cumming (2006) and Johnston et al. (2018). The predicted burst properties agree reasonably well with the predictions from Johnston et al. (2018), who used Kepler to model this outburst without doing an exhaustive parameter space exploration such as we have done in this work. Interestingly, Johnston et al. (2018) also predict a very bright burst at late time, $\approx 7$ days after the first observed burst.

Burst properties for individual bursts inferred by the model are listed in Table 3. The recurrence time for each burst is predicted to fall in the range 14 – 56 hours. The time to exhaust hydrogen and ignite bursts in a pure helium environment is 13.0 hours at this accreted fuel composition (Equation 9). These predicted recurrence times and their respective accretion rates indicate that sufficient time passes between each burst for the accreted hydrogen to be completely burned by the hot CNO cycle, and we have ignition occurring in a pure helium environment, as has been previously concluded by Galloway & Cumming (2006); Johnston et al. (2018). The low $\bar{X}$ values for each burst confirms this, with $\bar{X}<0.2$ for all bursts.

The probability contours and posterior distributions for $X$, $Z$, $Q_{\mathrm{b}}$, $d$, $\xi_{\mathrm{b}}$, and $\xi_{\mathrm{p}}$ are plotted in Figure 6. There is a positive correlation between $X$ and $Z$, which was also found by Galloway & Cumming (2006). This correlation is due to the relationship between recurrence time and fuel composition, since provided all of the H is exhausted before ignition, there are multiple combinations of $X$ and $Z$ that reproduce the observed recurrence time.

3.1 Distance and Anisotropies

We used the scaling parameters $r_{\mathrm{1}}$, $r_{\mathrm{2}}$, and $r_{\mathrm{3}}$ to obtain posterior distributions for the distance and emission anisotropies of the source, and thus the inclination, shown in Figure 6. Relating the scaling parameters to their corresponding physical parameters gives:

[TABLE]

[TABLE]

Assuming $Q_{\mathrm{nuc}}=Q_{\mathrm{nuc,pred}}$, this gives the following set of equations which enable calculation of d, $\xi_{\mathrm{b}}$ and $\xi_{\mathrm{p}}$ using the three scaling factors.

[TABLE]

[TABLE]

[TABLE]

Using these equations, we obtained posteriors on distance ($d$), burst anisotropy ($\xi_{\mathrm{b}}$) and persistent anisotropy ($\xi_{\mathrm{p}}$), shown in Figure 6 and values reported in Table 2. He & Keek (2016) modelled the relationship between anisotropy factors and inclination in low mass binaries with different disc shapes. They considered 4 different disc geometries, motivated by an effort to match observations of high reflection fractions observed in a couple of superbursts. Since SAX J1808.4--3658 has never shown evidence of such high reflection fractions, we chose the simplest disc geometry described by He & Keek (2016), Model A, corresponding to a flat disc. He & Keek (2016) Model A is approximated by Equation 21, which we used to also obtain a posterior distribution for cos $i$, based on the inferred anisotropy factors.

[TABLE]

where $i$ is the inclination of the system, the angle between the rotation axis and the direction from which the system is viewed by an observer.

For this disc Model A, we infer a distance to SAX J1808.4--3658 of $3.3^{+0.3}_{-0.2}\,\mathrm{kpc}$, which agrees within uncertainty of the predictions of Galloway & Cumming (2006) of 3.1–3.8 kpc. We note that Galloway & Cumming (2006) did not account for the anisotropy of emission when inferring their distance estimate. Applying our predicted anisotropy factors to their distance estimate reduces their estimate to 2.7–3 kpc, still within our predicted distance range.

3.2 Neutron Star Mass and Radius

We infer a neutron star mass of $M_{\mathrm{NS}}=1.5^{+0.6}_{-0.3}\,\mathrm{M}_{\odot}$ based on the gravity and gravitational redshift required by the model that gives the best match to the observations. Since the neutron star in this system has been accreting matter periodically, we expect it to be more massive than a single non-binary neutron star (e.g. Kiziltan et al., 2013). Canonically, the neutron star mass is assumed to be $M_{\mathrm{NS}}=1.4\,\mathrm{M}_{\odot}$, and recent works constraining NS mass and radii have found the NS mass could be between $1.1$–$2.2\,\mathrm{M}_{\odot}$ (e.g., Antoniadis et al., 2016). There has been extensive work in trying to estimate the mass of SAX J1808.4--3658 with no consensus on the value, but the general finding that it is more massive than the canonical value of 1.4 $\mathrm{M}_{\odot}$ (e.g. Heinke et al., 2009; Elebert et al., 2009; Morsink & Leahy, 2011). If the NS mass is larger than the usually assumed NS mass, this could indicate that some of the accreted matter is not burned in the observed X-ray bursts, as was modelled by Goodwin et al. (2018), however, due to the uncertainty of our mass prediction this result is not significant.

The posterior distributions and probability contours for mass and radius are plotted in Figure 7, coloured by surface gravity. The prior we used for mass and radius is plotted in red.

From Figure 7 it appears that there is a preferred surface gravity loosely driving the correlation between mass and radius, with the preferred surface gravity being $1.9\times 10^{14}\,\mathrm{g}\,\mathrm{cm}^{-2}$. Interestingly, the radius posterior peaks at a radius lower than the peak radius of the prior distribution, indicating that the model fit prefers a smaller radius. The model is thus sensitive to changes in radius. We explore this further in Section A Bayesian Approach to Matching Thermonuclear X-ray Burst Observations with Models, in which we compare model runs with fixed mass and radius and flat mass and radius priors.

3.3 Binary System Parameters and Inclination

The mass function of the binary orbit is given by

[TABLE]

where $i$ is the binary inclination, $M_{\mathrm{c}}$ is the companion star mass and $M_{\mathrm{NS}}$ is the neutron star mass (e.g., Chakrabarty & Morgan, 1998; Bildsten & Chakrabarty, 2001).

Setting the projected semi-major axis, $a_{x}\,\sin\,i=62.809\,\mathrm{light-ms}$ (Chakrabarty & Morgan, 1998) and orbital period, $P_{\mathrm{orb}}=2.01\,\mathrm{h}$ gives $f_{x}=3.8\times 10^{-5}\,\mathrm{M}_{\odot}$. This leaves $M_{\mathrm{NS}}$, $i$, and $M_{\mathrm{c}}$. Our MCMC algorithm uses settle to find $M_{\mathrm{NS}}=1.5^{+0.6}_{-0.3}\,\mathrm{M}_{\odot}$. Given estimates of $M_{\mathrm{c}}=0.02-0.011\,\mathrm{M}_{\odot}$ (e.g., Bildsten & Chakrabarty, 2001) and using Equation 22, the mass function gives cos $i=0.62^{+0.09}_{-0.08}$ ($i=\left(64^{+37}_{-20}\right)^{\circ}$).

Lack of deep X-ray eclipse rules out $i>82^{\circ}$ ($\cos i<0.15$) (Chakrabarty & Morgan, 1998). There is a $2\%$ modulation in X-ray intensity at the orbital period, which suggests that $i$ must be large enough to allow partial X-ray blockage from some circumbinary material (Bildsten & Chakrabarty, 2001). Those authors rule out face on inclination due to the 2 hour single-peaked modulation in the optical intensity, confirmed by Homer et al. (2001). The flux minimum occurs when the neutron star is behind the companion, implying that the X-ray heating of the companion is the origin of the modulation, which would have no effect in a face on system. Our best estimate gives $\cos i=0.36^{+0.07}_{-0.04}$ ($i=\left(69^{+4}_{-2}\right)^{\circ}$), which agrees within $3\,\sigma$ of $\cos i=0.62^{+0.09}_{-0.08}$ found using the mass function, and is within the ranges predicted by Bildsten & Chakrabarty (2001) and Wang et al. (2001). More recently, Wang et al. (2013) and Cackett et al. (2009) have attempted to constrain the inclination through optical modelling and iron line modelling respectively. Wang et al. (2013) found $i=\left(50^{+6}_{-5}\right)^{\circ}$ and Cackett et al. (2009) found $i=\left(55^{+8}_{-4}\right)^{\circ}$ ($1\,\sigma$ uncertainties). The Cackett et al. (2009) estimate agrees within $2\,\sigma$ of our best estimate, and the Wang et al. (2013) estimate agrees within $3\,\sigma$.

Since the anisotropy factors are determined independently in this method, we checked the self consistency of these factors with the relationships predicted by the Fujimoto (1988) and He & Keek (2016) models in Figure 8.

SAX J1808.4--3658 is in the hard spectral state during its outbursts and due to its nature as an accreting pulsar, it is predicted to have a truncated disc (e.g., Psaltis & Chakrabarty, 1999). The Fujimoto et al. (1981) and He & Keek (2016) models do not model for a truncated disc, and thus we do not expect any of them to match our predictions perfectly. We could not decisively conclude which of the models matched our predictions best, as they all pass through the $1\,\sigma$ confidence area in Figure 8, motivating us to use the simplest model, He & Keek (2016) Model A.

3.3.1 Wind loss

Recent research has found that a significant amount ($\approx 30\%$) of the nuclear energy of a Type I X-ray burst could be used to unbind matter from the neutron star surface, ejected in a wind. Yu & Weinberg (2018) found that wind losses can cause up to 30$\%$ loss of nuclear energy in ejecting material from the system. We do not correct for wind loss in our models as upon testing it simply reduced the predicted distance by $\approx 7\%$, with mass loss causing the source to be predicted closer due to fainter bursts being predicted.

3.4 Chain Convergence and MCMC performance

The convergence of MCMC chains is notoriously difficult to assess. Goodman & Weare (2010) recommend using the integrated autocorrelation time ($\tau$) to quantify the effects of sampling error on the results. The autocorrelation time provides a measure of the large time error of the Monte Carlo estimator, and so is a reliable indicator of the accuracy of the sampler by measuring the asymptotic variance in the limit of long chains. This essentially measures the effective number of independent samples, and enables estimation of the number of samples required to reduce the relative error on the target integral to a few percent. Foreman-Mackey et al. (2013) recommend running the chains for at least $50\,\tau$ samples (or conservatively $100\,\tau$ samples).

In Figure 9 the autocorrelation time of each parameter we varied is plotted as a function of the number of samples. The autocorrelation time is not reliably estimated early on, as there are not enough samples. As more samples are taken, the autocorrelation time estimate becomes more accurate, and after $\approx 500$ samples, the autocorrelation time converged on a single value for all of the parameters. By 2,000 samples, all of the autocorrelation estimates are well below the $\tau=N/100$ line, indicating this is an accurate estimate of the autocorrelation time. The largest autocorrelation time is $<$10 so the chains should be converged according to this convergence test after 1,000 steps. We ran the chains for 2,000 steps and 300 walkers to ensure they would be fully converged.

The acceptance fraction of the samples is a widely used measure to check the performance of MCMC. The acceptance fraction is the number of proposed steps that are accepted, and will be small if too many steps are rejected and large if too many steps are accepted (Foreman-Mackey et al., 2013). The general consensus is that the acceptance fraction should be between $0.2$–$0.5$ (e.g., Gelman et al., 1996), however, it depends on the model and situation. The acceptance fraction of our final run was $0.14$, which is lower than $0.2$ but not unreasonably low. As a test, we ran the MCMC code returning only the prior (excluding the likelihood and thus model call) and found an acceptance fraction of $0.11$. We can thus conclude that the prior limits are constraining parameter space and causing a larger than usual amount of steps to be proposed outside the allowed domain and be rejected.

A final test we ran to ensure the end walker positions were independent of the initial positions was to change the size of the Gaussian ball the walkers were initialised around by a factor of 100. This did not affect the marginalised posterior distributions of each parameter.

4 Discussion and Conclusions

We have used a Bayesian approach to successfully match observed burst fluence, alpha, and recurrence times with the settle model to infer system parameters of an accreting millisecond pulsar in outburst. The properties of the bursts observed from SAX J1808.4--3658 indicate they arrive in helium rich fuel, given the relatively short length, low flux, and high alpha values (e.g., Galloway & Cumming, 2006). We inferred $\bar{X}<0.2$ for all bursts, as expected for helium rich bursts. Overall, we matched the observed burst fluences, alphas, and recurrence times to within the $1\,\sigma$ uncertainty of the model.

Based on the agreement between the observed burst properties, we then inferred and modelled composition, base flux, neutron star mass, neutron star radius, distance, and inclination. We set the prior for metallicity to the expected distribution of metallicities based on the location of SAX J1808.3--3658 in the Galaxy, and deduced the most likely CNO mass fraction (metallicity) of the source to be between 0.009–0.019. This range covers the expected value for solar CNO metallicity, which recent work puts at $0.01$ (Lodders et al., 2009). We inferred a slightly depleted hydrogen mass fraction of the accreted fuel, of $0.58^{+0.13}_{-0.14}$, indicating that the companion star in this system could be significantly evolved.

We inferred a base flux of $0.25$–$0.7\,\mathrm{MeV}/\mathrm{nucleon}$, which gives an indication of the base heating of the neutron star. We used the constraints on the relationship between NS mass and radius calculated by Steiner et al. (2018) for neutron stars in low mass X-ray binaries in globular clusters as an informed prior for the distribution of acceptable mass and radius combinations and inferred a NS mass between $0.9$–$1.8\,\mathrm{M}_{\odot}$ and a radius between $10.8$–$13.1\,\mathrm{km}$. The mass is constrained to within $40\%$ ($1\,\sigma$ limit) and the radius is constrained to within $11\%$ ($1\,\sigma$ limit), indicating that this method cannot constrain neutron star mass and radius to a high precision. This is most likely due to degeneracies in some of the other model parameters, such as base heating ($Q_{\mathrm{b}}$), that can adjust to produce the same burst energies and recurrence time when the mass, radius, and thus surface gravity and redshift are varied.

Finally, we used the scaling factors we defined and the posterior distributions of the parameters to infer a distance and inclination of SAX J1808.4--3658, based on He & Keek (2016) disc Model A. We found a distance range of $3.1$–$3.6\,\mathrm{kpc}$ for an inclination from the rotation axis of $67^{\circ}$–$73^{\circ}$. Galloway & Cumming (2006) used the fact that the observed bursts showed photospheric radius expansion to infer an upper limit of the distance to the source of $3.6\,\mathrm{kpc}$, just within our $1\,\sigma$ upper limit on the distance. An inclination of $69^{\circ}$, implies both the observed burst and persistent fluxes are artificially increased by reflection from the disc and by the disc itself, bringing the source closer when compared to distance estimates such as Galloway & Cumming (2006) that do not account for inclination.

We have reported the median value and 68th percentile limits of the predicted posterior distributions to give an indication of the parameter values. While this is a representation of the data, the true solutions are the posterior distributions that were calculated by the MCMC algorithm. Parameters estimated using this technique are limited to the accuracy of the settle model predictions. The benefit of using such a simple model is that it can be used in a MCMC algorithm to effectively sample parameter space.

We have demonstrated that we can infer system parameters of a bursting source using observations of an outburst and the settle model. Future work could apply this method of matching observed properties with models to the other known bursting AMSPs with observations of burst trains, or any accreting neutron star with observations of burst trains. This method enables initial exploration of parameter space, to determine an estimate for the range of parameters for each source. These ranges could then be used in more advanced, computationally expensive models such as Kepler to narrow down the parameter limits.

Acknowledgements

We thank Andrew Casey for his contribution of the metallicity prior we used. We thank Andrew Steiner for helpful discussions regarding the neutron star mass-radius prior. This paper utilizes preliminary analysis results from 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. AG and ZJ acknowledge support by an Australian Government Research Training (RTP) Scholarship. AH was supported by a grant from Science and Technology Commission of Shanghai Municipality (Grant No. 16DZ2260200) and National Natural Science Foundation of China (Grant No. 11655002). This work benefited from support by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements). We thank the anonymous referee for constructive and insightful comments that helped to improve the manuscript.

Appendix A Additional Model Runs

In this section we provide the results of additional model runs carried out in order to robustly test the MCMC code, as well as to replicate past works.

A.1 Fixed mass and radius

We held mass and radius constant at $M=1.4\,$M*⊙* and $R=11.2\,$km to replicate the analysis of Galloway & Cumming (2006). The parameter predictions are presented in Table 4 and the probability contours for X, Z, Qb, d, $\xi_{\mathrm{p}}$ and $\xi_{\mathrm{b}}$ are presented in Figure 10.

The predicted hydrogen mass fraction in Table 4 matches within uncertainty of those predicted by Galloway & Cumming (2006) ($X_{\mathrm{0}}=0.54$ for $Z=0.02$, or $X_{\mathrm{0}}=0.5$ for $Z=0.016$). Galloway & Cumming (2006) inferred a base flux of $Q_{\mathrm{b}}=0.325$ and a distance of 3.1–3.8 kpc which is just outside the lower limit of our 1-$\sigma$ range of predictions for these parameters. Galloway & Cumming (2006) used a chi-squared minimisation approach to determine the parameters of interest given a fixed grid of parameters in $X$, $Z$ and $Q_{\mathrm{b}}$ space to find the best fit burst fluence, recurrence times and alpha values. We used a more robust approach to determine the parameters of interest, and vary more parameters than just $X$, $Z$, and $Q_{\mathrm{b}}$ and so we do not expect our results to match exactly with those of Galloway & Cumming (2006). The composition constraints and distance estimate broadly agree with those found by Galloway & Cumming (2006) for a fixed NS mass and radius.

A.2 Flat priors on mass and radius

Here we show the parameter limits when assuming a flat prior for mass and radius, rather than the Steiner et al. (2018) probability distribution. The parameter limits are reported in Table 5 and posterior distributions are plotted in Figure 11.

When assuming a flat prior range in mass and radius compared to assuming a more informative prior, we find that the mass is equally well constrained in both cases but that the radius is better constrained for the more informative prior run. Otherwise, all parameters remain approximately the same, except for surface gravity, with a stronger surface gravity required by the flat prior run.