THE INFLUENCE OF ACCRETION RATE AND METALLICITY ON THERMONUCLEAR BURSTS: PREDICTIONS FROM KEPLER MODELS

We calculated parameters describing each burst in the train, which were averaged to give characteristic values for each set of initial conditions. We briefly cover the definition of each parameter in the following paragraphs.

The burst fluence, E_b, is determined by numerically integrating the burst luminosity with time

$$\begin{eqnarray}&&{E}_{{\rm{b}}}={\displaystyle \int }_{{t}_{\mathrm{start}}}^{{t}_{\mathrm{end}}}(L(t)-{L}_{\mathrm{th}})\ {dt},\end{eqnarray} \tag{ 1 }$$

where the persistent thermal emission is subtracted from the burst prior to integration. While ${t}_{\mathrm{start}}$ may be well before the burst rise, luminosities here contribute negligibly to the burst fluence, as they are close to the thermal level before the burst begins. We calculated the ratio of fluence to peak flux ($\tau ={E}_{{\rm{b}}}/{L}_{{\rm{p}}}$), as it is often used as a measure of burst morphologies (van Paradijs et al. 1988).

We also consider the ratio of persistent fluence between bursts to the burst fluence, commonly referred to as α (Gottwald et al. 1986). The alpha ratio is typically found by integrating the persistent emission; however, as accretion is constant within the models studied and the dominant source of persistent emission, we calculate α as follows:

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{E}_{{\rm{p}}}}{{E}_{{\rm{b}}}}=\displaystyle \frac{{L}_{\mathrm{acc}}+{L}_{\mathrm{th}}}{{E}_{{\rm{b}}}}{\rm{\Delta }}t,\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Delta }}t$ is the time since the previous burst and L_acc is found from the accretion rate. Note that this definition restricts us from calculating the value of α for the first burst in a train, as the recurrence time since the last burst is unknown. In the Newtonian frame L_acc is

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}=\displaystyle \frac{{GM}\dot{M}}{r}.\end{eqnarray} \tag{ 3 }$$

This neglects the gravitational redshift under which a correction for a distant observer is required (described in Appendix B). Note that the accretion luminosity exceeds the thermal persistent luminosity by at least an order of magnitude (typically two orders of magnitude).

To quantify the burst rise both its duration and morphology are measured. The duration is measured by the time taken to reach 10%, 25%, and 90% of the peak luminosity above the thermal level, given by t₁₀, t₂₅, and t₉₀, respectively (we additionally define the overall burst duration as the time from t₂₅ to the end of the burst). The rise times are then given in two ways (to enable comparison with existing measures), by ${t}_{10-90}$ and ${t}_{25-90}$, the times respectively taken to rise from 10% and 25% of the peak flux to 90% of the peak flux. The shape of the rise can be measured by the convexity (${ \mathcal C }$) of the burst. Maurer & Watts (2008) define the burst convexity as the deviation of the light curve from linear during the rise from t₁₀ to t₉₀. Observationally, variations in convexity can be correlated with the latitude of burst ignition on a neutron star. The convexity is calculated by rescaling the burst between t₁₀ and t₉₀ in both time and luminosity so each parameter varies from 0 to 10. Calling the rescaled luminosity l, and the rescaled time x, the line joining the two parameters is equal along the line l = x. The light curve is then a distance of l − x above the straight line joining the two boundary luminosity points. Integrating this difference, the convexity can be found as

$$\begin{eqnarray}&&{ \mathcal C }={\displaystyle \int }_{0}^{10}(l-x)\ {dx}.\end{eqnarray} \tag{ 4 }$$

We evaluate this integral numerically for each burst profile. Additionally, as the models do not have a uniform time resolution, we approximate t₁₀ and t₉₀ by linearly interpolating between the two points closest to $L=0.1{L}_{{\rm{p}}}$ and $L=0.9{L}_{{\rm{p}}}$.

Exponential fits are frequently made to the cooling tail of bursts (e.g., Galloway et al. 2004). Recently, in't Zand et al. (2014) suggested that the burst tail is better fit by a power law. We fitted single tailed exponential curves and power laws to each light curve. The exponential fit is given by

$$\begin{eqnarray}&&L(t)={L}_{0}\mathrm{exp}(-(t-{t}_{0})/{t}_{{\rm{b}}})+{L}_{\mathrm{bg}},\end{eqnarray} \tag{ 5 }$$

while the power law is fitted by

$$\begin{eqnarray}&&L(t)={L}_{0}{\left(\displaystyle \frac{t-{t}_{{\rm{s}}}}{{t}_{0}-{t}_{{\rm{s}}}}\right)}^{-\kappa }+{L}_{\mathrm{bg}}.\end{eqnarray} \tag{ 6 }$$

In both of these cases the persistent luminosity is measured as the luminosity 20 s before the burst peak, and t₀ is the time of the first fitted data point. In the case of the exponential fit, this leaves two free parameters, L₀ and t_b, while the power law has three free parameters, L₀, t_s, and κ. The most noticeable difference between the two curves is the divergence of the power law near t_s, a time loosely corresponding to the burst start. Due to the nature of the power law, t_s and κ are strongly coupled.

To determine the range of the fit, in't Zand et al. (2014) varied t₀, increasing this parameter's value until a local minimum in ${\chi }_{\nu }^{2}$ is reached. We instead fitted the burst tail, from the time where the luminosity drops to 90% of its peak value to the time where it falls to 10% of its peak value above the persistent thermal level. As each separated light curve does not have any error estimates, we fit via a least-squares minimization.

The individual bursts are also averaged together to produce mean luminosity and radius profiles for each model. This provides a single representative burst for each run. Where three or more bursts occur, the first burst is ignored, as it is typically substantially different from following bursts, as it lacks an ashes layer from previous bursts (Woosley et al. 2004). In models that produce only two bursts the first two bursts are averaged. While this leads to models that have only two bursts showing a significantly higher variation in their parameters than those with three or more, on account of the different fuel column in the first burst, keeping this burst in the averaging process yields a measure of luminosity variation that would be otherwise lost. The actual luminosities that are averaged use the duration of the longest burst found for the train, rather than averaging bursts of different durations. During the averaging process, bursts are aligned so that the times of their peak luminosities coincide.

We present the averaged light curves and their data files online⁵ , to accompany the averaged burst parameters. The averaged light curves provide luminosities and the 1σ variation in luminosity. The averaged radius of the last zone in KEPLER is also provided for each model with its 1σ variation.

We also compare observations of GS 1826−24 to model predictions within this paper, using bursts observed by RXTE (Jahoda et al. 1996) in the MINBAR⁶ burst archive. The data analysis techniques used to produce MINBAR are similar to those found in Galloway et al. (2008). RXTE provides time resolved spectra across the energy range 2–60 $\mathrm{keV}$, which, during each burst, were binned in $0.25\ {\rm{s}}$ intervals during the burst rise and peak. During the tail, the time-bins were increased in order to preserve the S/N. The burst background was estimated from a $16\ {\rm{s}}$ interval prior to the burst.

We accounted for the "deadtime" experienced between photons within the RXTE PCA detectors, which reduces the count rate below that experienced by the detector (by approximately 3% at an incident event rate of $\approx 400\ \mathrm{counts}\ {{\rm{s}}}^{-1}\ {\mathrm{PCA}}^{-1}$). Each spectrum was corrected for deadtime by taking an effective exposure that accounts for the deadtime fraction. The spectra were subsequently re-fitted with a blackbody model over the interval 2.5–20 $\mathrm{keV}$ using XSpec version 12 (Arnaud et al. 1996) with Churazov weighting to accomodate spectral bins with low count rates. The effects of interstellar absorption were modeled using a multiplicative component (wabs in XSpec) with the column density fixed at ${nH}=0.4$ (in't Zand et al. 1999).

The relationship between accretion rate and recurrence time has been studied for a number of sources, including GS 1826−24 (Galloway et al. 2004; Thompson et al. 2008) and MXB 1730−335, also known as the rapid burster (Bagnoli et al. 2013). Observationally, a power law relationship has frequently been found, relating the persistent flux, $({f}_{\mathrm{pers}})$ of the source to the burst recurrence time, with the form

$$\begin{eqnarray}&&{\rm{\Delta }}t\propto {f}_{\mathrm{pers}}^{-\eta }.\end{eqnarray} \tag{ 7 }$$

The persistent flux of a neutron star is primarily due to an energy deposition of ${Q}_{\mathrm{grav}}{{\rm{g}}}^{-1}$ from accretion onto the star. For a star of radius R at a distance d the flux (given an anisotropy of ${\xi }_{\mathrm{pers}}$) is given by

$$\begin{eqnarray}&&{f}_{\mathrm{pers}}=\displaystyle \frac{{L}_{\mathrm{pers}}}{4\pi {d}^{2}}=\displaystyle \frac{\dot{M}{Q}_{\mathrm{grav}}}{1+z}{\left(\displaystyle \frac{R}{d}\right)}^{2}{\xi }_{\mathrm{pers}}^{-1},\end{eqnarray} \tag{ 8 }$$

and hence ideally, the persistent flux is directly proportional to the mass accretion rate. It then follows that variations in recurrence time are due to changes in the accretion rate,

$$\begin{eqnarray}&&{\rm{\Delta }}t\propto {\dot{M}}^{-\eta }.\end{eqnarray} \tag{ 9 }$$

Naïvely, we can expect the power law gradient to be $\eta \sim 1$ if we assume that a constant critical mass M_crit is required for a burst, as ${\rm{\Delta }}t={M}_{\mathrm{crit}}/\dot{M}$. Deviations from this relation suggest that the mass necessary for an instability varies with accretion.

For each accretion composition in Table 2, we identified model runs that exhibit five or more bursts, and fitted power laws for the variation in recurrence time with accretion rate across the region of Case III burning. Recurrence times were redshifted by $z=0.26$, to be observationally equivalent to a $M=1.4\;{M}_{\odot }$, $r=11.2\;\mathrm{km}$ neutron star. We excluded models where ${\rm{\Delta }}t\lt 100\;{\rm{s}}$, though this only occurs for $Z\geqslant 0.10$. The best fit power law was found via a ${\chi }^{2}$-minimization, using the variation in recurrence time throughout the model run to weight each point. At the high accretion rate end of the power law, we stopped fitting when the recurrence time stopped decreasing. At the low accretion rate end, we stopped fitting when the points followed a noticeably different linear trend in log-log space (noting that the nature of the deviation changes with metallicity). For metallicities $z=0.020$ and above, this transition corresponds to a change in the burst behavior, from the long-tailed light curve morphology (characteristic of mixed H/He bursts) to an Eddington-limited light curve caused by pure He bursts. For these runs, the power law fit spans the entire range of Case III burning. The linear relationship becomes considerably steeper in the pure He bursting region. In the low metallicity cases ($z=0.000,z=0.001\;\mathrm{and}\ z=0.002$), mixed H/He burning still continues at accretion rates below where the power law stops being fitted, as the trend is translated vertically up for ${l}_{\mathrm{acc}}\lt 0.08$. The value of η and its associated 1σ uncertainty from the ${\chi }^{2}$ minimization are shown in Table 3, and the actual data and fits are shown in Figure 1, with the fitted points plotted in red.

Zoom In Zoom Out Reset image size

Table 3.  Fitted Power Law Index η for Case III Burning

Z	η	n	${l}_{\mathrm{acc},\mathrm{IV}}$	${l}_{\mathrm{acc},\mathrm{III}}$	${l}_{\mathrm{acc},\mathrm{II}}$
0.000	1.19 ± 0.03	10	⋯	0.75	0.92
0.001	1.17 ± 0.02	21	⋯	0.80	0.92
0.002	1.237 ± 0.014	10	⋯	0.80	0.97
0.004	1.17 ± 0.02	16	⋯	0.80	1.00
0.010	1.109 ± 0.006	32	⋯	1.10	1.11
0.020	1.106 ± 0.016	21	0.03	1.20	1.26
0.040	1.16 ± 0.02	19	0.07	1.55	1.71
0.100	1.18 ± 0.02	17	0.20	4.20	4.20
0.200	1.16 ± 0.02	7	0.70	6.50	7.50
N&H03			0.04	0.15	0.25
F81			0.08	⋯	1.0

Note. The power laws were found following a ${\chi }^{2}$ minimization of Equation (9) with respect to each data set. The number of fitted points is denoted by n. The accretion luminosities (in terms of l_Edd) where burning transitions occur are also shown for each metallicity. The transition times of Narayan & Heyl (2003) are reprinted from Table 1 alongside the transitions predicted by Fujimoto et al. (1981), who do not distinguish Case IV from steady burning.

Download table as:  ASCIITypeset image

We find that in the KEPLER models, the value of η depends upon the stellar composition and accretion rate. Using the recurrence times derived from our analysis of the model catalog we find that the variation with accretion rate is consistent with a power law across a wide range of accretion rates and metallicities.

There is no significant correlation between power law index and metallicity (Spearman $\rho =-0.48$, $p=0.187\approx 1.3\sigma$). Interestingly, all metallicities show $\eta \gtrsim 1.1$, while observations do not show such steep gradients. Galloway et al. (2004) find $\eta =1.05\pm 0.02$ for GS 1826−24 while Bagnoli et al. (2013) find $\eta =0.95\pm 0.03$ for MXB 1730−335. Linares et al. (2012) find that the X-ray pulsar IGR J17480−2446 shows η ≈ 1 when ${L}_{\mathrm{pers}}\gt 0.3{L}_{\mathrm{Edd}}$, and $\eta \approx 3$ when $0.2{L}_{\mathrm{Edd}}\lt {L}_{\mathrm{pers}}\lt 0.3{L}_{\mathrm{Edd}}$. The assumption that the burst anisotropy is independent of accretion rate is necessary to compare observational and theoretical power law indices. It is likely, however, that there are burst behaviors that violate this assumption.

We found that models with 1% and 2% metallicity come closest to the observed η values ($\eta \approx 1.11$). However, it would be overstating the available evidence to suggest this implies that the observational sources discussed in the preceeding paragraph are consequently metal-rich, as modeled power laws could be systematically high. It is nevertheless notable that metallicities around solar show a minimum η value in the parameter range tested, suggesting that if a large sample of power law gradients were compiled, the lowest gradients would correspond to metal-rich population I stars.

Using theoretical burst ignition models, Galloway et al. (2004) found $\eta \lt 1$ for solar metallicities, with η increasing as metallicity drops. The KEPLER models show noticeably higher power law indices than this, though the trend for steeper power laws to be found at lower metallicities is replicated. Part of the reason for our disagreement with these models is that the models used in Galloway et al. neglect additional heating of the atmosphere by accumulated CNO, leading to systematic differences from KEPLER.

We find also that the transition accretion rates, where the burning regime changes, increase with metallicity (Table 3). The transition accretion rates are defined where the burst behavior of the model changes in a qualitatively similar way to that defined by the burst cases presented by Narayan & Heyl (2003). We denote ${l}_{\mathrm{acc},\mathrm{IV}}={L}_{\mathrm{acc},\mathrm{IV}}/{L}_{\mathrm{Edd}}$ as the accretion luminosity (in terms of Eddington) as the lowest accretion luminosity at which we observe Case III bursts. ${l}_{\mathrm{acc},\mathrm{III}}$ denotes the accretion luminosity at which the shortest recurrence time is observed, and ${l}_{\mathrm{acc},\mathrm{II}}$ denotes the highest accretion rate for which we observe bursts. KEPLER simulations produce roughly similar transition accretion rates as stability analyses suggest, noting that both Fujimoto et al. (1981) and Narayan & Heyl (2003) consider solar accretion compositions.

Fujimoto et al. (1981) consider three neutron star potentials, all significantly different from the surface gravity adopted in KEPLER. Of the two closest to our considered surface gravity, the first overestimates the surface gravity by $\sim 55\%$ ($M=1.00\;{M}_{\odot }$, $r=7.49\;\mathrm{km}$), while the second underestimates it by $\sim 35\%$ ($M=0.476\;{M}_{\odot }$, $r=8.67\;\mathrm{km}$). These analyses predict ${l}_{\mathrm{acc},\mathrm{IV}}=0.08$ for the higher gravity, and ${l}_{\mathrm{acc},\mathrm{IV}}=0.03$ for the lighter gravity. We find ${l}_{\mathrm{acc},\mathrm{IV}}=0.03$ for $Z=0.02$, agreeing with the lower surface gravity case. Narayan & Heyl (2003) analyzed the stability of burning on a star with only a 20% ($M=1.40\;{M}_{\odot }$, $r=10.0\;\mathrm{km}$) higher surface gravity than KEPLER using a more rigorous approach, finding ${l}_{\mathrm{acc},\mathrm{IV}}=0.04$. While each study uses different surface gravities, our value for the transition rate from Case IV to Case III burning broadly agrees with earlier stability analyses.

Despite this, the transition from Case III to the Case II and I regimes occurs at noticeably higher accretion rates than predicted by Narayan & Heyl (2003). We instead find these changes occur near ${l}_{\mathrm{acc}}\sim 1.2$, rather than the ${l}_{\mathrm{acc}}\sim 0.25$ predicted by Narayan & Heyl (2003). This is in closer agreement to Fujimoto et al. (1981), who find that the transition to stable burning occurs at ${l}_{\mathrm{acc}}\sim 1$ regardless of surface gravity. Additionally, Case II burning occurs in general over a smaller range of accretion rates than Narayan & Heyl (2003) predict; however, our method of choosing when Case III burning transitions to Case II arguably misses part of the Case II region. This is most noticeable in the $Z=0.01$ runs as Case II burning shows oscillations in recurrence rate that increase the burst rate to faster than Case III burning. While in general the minimum recurrence time provides an approximate measurement of the Case III to Case II transition, burning compositions and burst morphologies would be more robust.

The general shape of the $\dot{M}\mbox{--}{\rm{\Delta }}t$ relationship outside of the region of Case III burning follows the predictions of Narayan & Heyl (2003), with a turnup where overstable burning commences, and a steeper power law index in the pure He bursting regime. In the overstable burning regime (Case II) the recurrence time increases with accretion rate, as a large proportion of the accreted fuel is being burnt during accretion. As both hydrogen and helium are burnt quickly in this regime a longer time is required to accrete sufficient fuel for a burst to occur. The turnup is relatively small within the KEPLER simulations, increasing the recurrence time by a fraction of an hour. This is in disagreement with the stability analyses of Narayan & Heyl (2003) who find that the recurrence time may increase by up to a day beyond the shortest recurrence times they observe. While these two results are quantitively different, the upturn we observe does support qualitatively the overstable behavior predicted for Case II bursting.

The behavior of the upturn is especially notable for the $Z=0.04$ case (Figure 2). While the upturn is reasonably brief in most compositions, it shows substantial structure for this metallicity. Physically, at these accretion regimes the increase in recurrence time is linked to helium in addition to hydrogen burning between bursts depleting the available fuel. In this regime luminosity between bursts tends not to fall quickly back to its persistent level. Interestingly, near ${l}_{\mathrm{acc}}\sim 1.65$ the uncertainty in ${\rm{\Delta }}t$ again becomes small, similar to the Case III uncertainties, and ${\rm{\Delta }}t$ slightly decreases, before resuming its upward trend. In general, the large uncertainties in ${\rm{\Delta }}t$ in this region are a consequence of metastable H/He burning, and the large impact small changes in the composition of the atmosphere can have on the stability of the accreted column. It is possible that in the intermediate area where ${l}_{\mathrm{acc}}\sim 1.65$ the column is slightly more stable; the exact reason for this will be investigated in more detail in the future.

Zoom In Zoom Out Reset image size

The steeper power law index in the region of pure He bursts is caused by the higher amount of accumulated helium required for a burst as compared to hydrogen, and the lack of heating from H burning. We could better constrain the power law behavior in this region with more simulations at low accretion rates; however, these typically require very long computation times per burst, as in the helium bursting regime more integration steps are needed per burst.

The transition between Case IV (pure He) burning and Case III (mixed H/He) burning is expected to cause a significant change in the recurrence time due to the higher temperatures required for runaway helium burning. This leads to a discontinuity in the ${\rm{\Delta }}t-\dot{M}$ relationship at accretion rates low enough for the accreted hydrogen to be depleted by hot CNO burning between bursts.The discontinuity occurs at a recurrence time similar to the time for the hot CNO cycle to deplete all hydrogen present. Fujimoto et al. (1981) express this relation as ${\rm{\Delta }}t=9.7\;\mathrm{hr}\;{({\text{}}Z/0.02)}^{-1}$ in the neutron star frame, while Galloway et al. (2004) suggest ${\rm{\Delta }}t=11\;\mathrm{hr}\;{({\rm{Z}}/0.02)}^{-1}({\rm{X}}/0.7)$. We attribute the discrepancy in the Galloway et al. expression to its adaptation of an earlier calculation in Bildsten (1998) that is based on energy arguments fail to incorporate neutrino energy losses in the hot CNO cycle, and was based on a different value of ${\bar{\mu }}_{\mathrm{CNO}}$ (see below) using CNO ratios from older solar abundance determinations. Given this difference, however, we provide an independent derivation of this relationship below, laying out the assumptions inherent within it.

The hot CNO cycle (Mathews & Dietrich 1984) occurs when temperatures are sufficiently high for proton capture by ${}^{13}{\rm{N}}$ to occur more frequently than ${}^{13}{\rm{N}}$ undergoes ${\beta }^{+}$ decay. This leads to the hot CNO reaction cycle

$$\begin{eqnarray}&&{}^{12}{\rm{C}}{(p,\gamma )}^{13}{\rm{N}}{(p,\gamma )}^{14}{\rm{O}}{({\beta }^{+}\nu )}^{14}{\rm{N}}(p,\gamma )\\ &&{}^{15}{\rm{O}}{({\beta }^{+}\nu )}^{15}{\rm{N}}{(p,\alpha )}^{12}{\rm{C}}.\end{eqnarray} \tag{ 10 }$$

The cycle is limited by the ${\beta }^{+}$ decay of ${}^{14}{\rm{O}}$ and ${}^{15}{\rm{O}}$. To find the time to burn all available hydrogen, we consider that each decay of ${}^{14}{\rm{O}}$ in the cycle requires four proton captures. This gives the burning rate of ${}^{1}{\rm{H}}$ as follows:

$$\begin{eqnarray}&&{\partial }_{t}{Y}_{{}^{1}{\rm{H}}}=-4\displaystyle \frac{\mathrm{ln}2}{{t}_{\frac{1}{2}{,}^{14}{\rm{O}}}}{Y}_{{}^{14}{\rm{O}}}.\end{eqnarray} \tag{ 11 }$$

Integrating this expression with the initial condition that ${Y}_{{}^{1}{\rm{H}}(t=0)}={Y}_{{}^{1}{\rm{H}}}^{\mathrm{init}}$ and finding when the hydrogen is depleted gives the time taken to burn an accumulated fuel layer by the hot CNO cycle,

$$\begin{eqnarray}&&{t}_{\mathrm{CNO}}=\displaystyle \frac{{Y}_{{}^{1}{\rm{H}}}^{\mathrm{init}}}{{Y}_{{}^{14}{\rm{O}}}}\cdot \displaystyle \frac{{t}_{\frac{1}{2}{,}^{14}{\rm{O}}}}{4\mathrm{ln}2}\end{eqnarray} \tag{ 12 }$$

In the hot CNO cycle, the reaction rate is restricted by the beta decay of ${}^{14}{\rm{O}}$ and ${}^{15}{\rm{O}}$, so at any point the majority of CNO is in these two isotopes, meaning ${Y}_{{}^{14}{\rm{O}}}+{Y}_{{}^{15}{\rm{O}}}\approx {Y}_{\mathrm{CNO}}$. Noting that $Y{(}^{14}{\rm{O}})/Y{(}^{15}{\rm{O}})={t}_{\frac{1}{2}{,}^{14}{\rm{O}}}/{t}_{\frac{1}{2}{,}^{15}{\rm{O}}}$ the number fraction of $Y{(}^{14}{\rm{O}})$ is expressible in terms of the CNO abundance:

$$\begin{eqnarray}&&\displaystyle \frac{{Y}_{{}^{14}{\rm{O}}}}{{Y}_{\mathrm{CNO}}}\approx \displaystyle \frac{{Y}_{{}^{14}{\rm{O}}}}{{Y}_{{}^{14}{\rm{O}}+{Y}_{{}^{15}{\rm{O}}}}}=\displaystyle \frac{{t}_{\frac{1}{2}{,}^{14}{\rm{O}}}}{{t}_{\frac{1}{2}{,}^{14}{\rm{O}}}+{t}_{\frac{1}{2}{,}^{15}{\rm{O}}}}.\end{eqnarray} \tag{ 13 }$$

Proceeding from Equation (12), the CNO burning time can now be expressed as

$$\begin{eqnarray}&&{t}_{\mathrm{CNO}}=\displaystyle \frac{{Y}_{{}^{1}{\rm{H}}}}{{Y}_{\mathrm{CNO}}}\cdot \displaystyle \frac{{t}_{\frac{1}{2}{,}^{14}{\rm{O}}}+{t}_{\frac{1}{2}{,}^{15}{\rm{O}}}}{4\mathrm{ln}2}.\end{eqnarray} \tag{ 14 }$$

Conversion of Y_CNO to a metallicity mass fraction requires assumptions about the composition in order to find the mean molecular mass. Assuming that all of the accreted metals are in CNO with their solar ratios (Asplund et al. 2009) the number fraction ${Y}_{\mathrm{CNO}}=Z/{\bar{\mu }}_{\mathrm{CNO}}$, with ${\bar{\mu }}_{\mathrm{CNO}}=14.5\;{\rm{g}}\;{\mathrm{mol}}^{-1}.$ Substituting the half-lives of both oxygen species (${t}_{\frac{1}{2}{,}^{14}{\rm{O}}}=71\;{\rm{s}}$ and ${t}_{\frac{1}{2}{,}^{15}{\rm{O}}}=122\;{\rm{s}}$) the hydrogen burning time is given by

$$\begin{eqnarray}&&{t}_{\mathrm{CNO}}=9.8\;\mathrm{hr}\;\left(\displaystyle \frac{{\rm{X}}}{0.7}\right){\left(\displaystyle \frac{Z}{0.02}\right)}^{-1},\end{eqnarray} \tag{ 15 }$$

noting that the prefactor is within the precision of Fujimoto et al. (1981), a difference that could be accounted for by a small change in the measurement of ${\bar{\mu }}_{\mathrm{CNO}}$ used. This prefactor can show some variability depending on what the metal component is assumed to be (noting that this variability comes from using an expression based on mass rather than number fraction). For example, assuming purely ${}^{16}{\rm{O}}$ (${}^{14}{\rm{N}}$, ${}^{12}{\rm{C}}$) leads to a prefactor of $10.8\;\mathrm{hr}$ ($9.5\;\mathrm{hr}$, $8.1\;\mathrm{hr}$). Additionally, as the CNO fraction does not typically comprise all the metallicity within a star, $9.8\;\mathrm{hr}$ is likely an underestimate for the CNO burning time. For solar composition (Asplund et al. 2009), only 66% of the metals (by mass) are CNO. Taking this into consideration the prefactor increases to $15\;\mathrm{hr}$ in the neutron star frame. We do not expect the KEPLER models to be best modeled by this last case, as the metallicity component accreted is entirely ${}^{14}{\rm{N}}$ rather than solar.

In runs at metallicities of $Z=0.02$ and higher, pure He bursts, indicated by PRE, are observed at low accretion rates (blue points at the low accretion end of the power law in Figure 1). For each composition, we expect the power law that describes Case III burning to have an upper bound in recurrence time for each composition to be described by the hot CNO cycle hydrogen burning time (Equation (12)). For $Z=0.02\ (0.04,0.10,0.20)$ this time is ${\rm{\Delta }}t=12\ (5.7,1.7,\ 0.4)\;\mathrm{hr}$ after redshifting by $1+z=1.26$. In each model run we observe Case III burning at recurrence times shorter than the hot CNO burning time and Case IV burning above these times, as expected.

This agreement is also noticed using the $11\;\mathrm{hr}$ prefactor suggested by Galloway et al. (2004). The modeled results in Galloway et al. (2004) have a redshift of $1+z=1.31$, leading to a rest frame prefactor of $8.4\;\mathrm{hr}$. Such a value could arise from accretion if most of the metallicity fraction is carbon. Most of the suggested prefactors thus far match our simulated data. Figure 1 features blue shading to highlight the range of times wherein the CNO cycle is expected to deplete available hydrogen, using the unredshifted prefactors of $8.1\;\mathrm{hr}$ and $15\;\mathrm{hr}$ as limits. Only the low end ($\lesssim 9\;\mathrm{hr}$) of this range is excluded by the KEPLER models.

Considering all this detail, we emphasize the approximate nature of Equation (15). Changes in accretion composition, gravitational redshift and the presence of ashes from previous bursts can foreseeably combine to greatly vary this prefactor. Increasing the resolution of simulations in accretion rate in order to better isolate the transition from Case IV to Case III burning may lead to a better determination of the prefactor (given the assumptions within KEPLER) and also allow some quantification of the deviation CNO ashes cause from the theoretical prefactor. However as observational systems have significantly greater deviations from theory than KEPLER, better constraints are unlikely to be very useful.

For models run at metallicities lower than $Z=0.02$, no runs show Case IV burning, as they do not span accretion rates low enough to deplete all the available hydrogen. Rather, a noticeable vertical translation of the power law index occurs for accretion rates with ${l}_{\mathrm{acc}}\lesssim 0.08$. To determine the significance of this deviation we considered the ${\chi }_{\nu }^{2}$ for fits to the recurrence time versus the accretion rate, across both the range of accretion rates from ${l}_{\mathrm{acc}}=0.08$ to the Case II transition, and for all accretion rates up to the Case II transition. For $Z=0.002$ ($Z=0.001$, $Z=0.0$) we found that using the larger fitting range drastically increased the ${\chi }_{\nu }^{2}$ from 0.98 (2.65,3.36) to 21.8 (33.0, 29.5). While these fits are phenomenological and a high ${\chi }^{2}$ may occur on account of this, the variation in these two fits suggests that at accretion rates ${l}_{\mathrm{acc}}\lesssim 0.08$ the power law that describes these bursts is different.

To this end, we examined the morphology of bursts at this low accretion rate, noticing that there is a change in the structure of the burst tail at accretion rates lower than where we cease fitting. At such low accretion rates, the burst light curve shows a two stage decay (Figure 3) rather than a single decay. The presence of a second decay, reached after the first decay flattens indicates an additional burning process, which could act to consume additional fuel, thereby extending the accretion time required to recuperate an unstable column. Further investigation into the nuclear networks active in KEPLER during these bursts, as well as the rapidity of the onset of the two stage decay, could present an interesting avenue for further study.

Zoom In Zoom Out Reset image size

Both equivalent burst duration (τ) and recurrence time (${\rm{\Delta }}t$) are commonly used observational measures in burst studies. Importantly, both measures are independent of distance. As such, the loci traced by simulated bursts in $\tau -{\rm{\Delta }}t$ space provide a series of curves which can be compared to observations (for an assumed gravitational redshift). This is considered in Figure 4, along with observed values of τ and ${\rm{\Delta }}t$ for GS 1826−24. The simulated parameters have been redshifted by $1+z=1.26$, corresponding to a $1.4\;{M}_{\odot }$ neutron star $11.2\;\mathrm{km}$ in radius.

Zoom In Zoom Out Reset image size

We compared these simulated parameters to observations of GS 1826−24 for which we had a reliable recurrence time (determined from finding two bursts close in time). The values of recurrence time and τ for GS 1826−24 appear between the locus traced by the $Z=0.01$ and $Z=0.02$ models. This suggests an accretion metallicity slightly lower than $Z=0.02$, in agreement with Heger et al. (2007a), although we note that their comparison also uses KEPLER models. The span of the observations corresponds to modeled accretion rates of ${l}_{\mathrm{acc}}=0.065$ to ${l}_{\mathrm{acc}}=0.123$. To illustrate the behavior of GS 1826−24 we have fitted a linear model ($\tau =A{\rm{\Delta }}t+B$) to the observations by ${\chi }^{2}$-minimization, finding $A=-2.8\pm 0.3\;{\rm{s}}\;{\mathrm{hr}}^{-1}$ and $B=50\pm 1\;{\rm{s}}$ (${\chi }_{\nu }^{2}=1.2$).

Interestingly, the gradient from the GS 1826−24 observations in Figure 4 shows a shallower drop in ${\rm{\Delta }}t$ with τ than the modeled points. As higher recurrence times are indicative of lower accretion rates, we see higher fluences than KEPLER predicts as accretion rate drops. The deviation from the locus of $Z=0.02$ is largely driven by the three data points at recurrence times $\gt 5\ \mathrm{hr}$. While it could be tempting to to eliminate these points as examples of non-standard behavior of the star, recent observations of GS 1826−24 have shown the star exhibiting an uncharacteristic soft spectral state (Chenevez et al. 2015), strengthening the idea that these points reflect the behavior of GS 1826−24, and that its behavior is thus not in accordance with that predicted by KEPLER for $Z=0.02$ accretion. Part of the reason for this may be that the accretion regime used by us, whereby only ${}^{14}{\rm{N}}$ is being accreted, does not provide adequate seed nuclei to reproduce the behaviors of naturally occurring stars where a more diverse quantity of metals are accreted. It woud be interesting in further studies to investigate how variations in the accretion composition affect these loci.

A recurring puzzle in the morphologies of modeled bursts is their variablility. Naïvely, models would be expected to show relatively uniform bursts (apart from the first few bursts in a simulation) as the composition of the neutron star atmosphere approaches a steady level. This is seen to some extent: the first burst in a train is disproportionately energetic, and the following bursts show variations of about 10% between the peak luminosity of the brightest and the dimmest.

This behavior implies that there is a natural variability in burst light curves that arises from nuclear physics, yet some observational sources, notably GS 1826−24, do not show as much variability within a burst epoch as simulated bursts show within a model run. We suggest that this is due to the number of convective cells on a bursting neutron star collectively removing the variation we see in the models. Each model effectively assumes the neutron star contains a single convective cell. Rather, assuming a cell diameter similar to the atmosphere height ($\sim 10\;\mathrm{cm}$), a neutron star has on the order of 10⁷ convective cells. This would produce a light curve equivalent to the average of 10⁷ modeled bursts and thereby suppress the variations in burst behavior that are caused by nuclear physics.

Variations in the shape of the burst rise can be expected as a function of the latitude on the neutron star at which the burst commences; nuclear processes; burst spreading; and also a contribution from scattering within the accretion disk. The KEPLER models simulate nuclear processes without simulating accretion disk interactions or ignition latitude, allowing the contribution of nuclear burning to burst convexity to be determined.

For our discussion of convexity we restrict the sample to bursts with a metallicity of $Z\leqslant 0.04$ that exhibit a smooth rise. This excludes many PRE bursts, as these often exhibit sudden increases in luminosity from convective zones rising to the stellar surface as burning commences. Bursts from runs with higher metallicities were excluded, as their convexities are skewed by the development of two peaks in the rise (see Section 4.5).

Figure 5 shows the variations in burst convexity with both accretion rate and metallicity. For metallicities less than $Z=0.002$, the convexity changes very little as the accretion luminosity changes from $0.1\;{L}_{\mathrm{Edd}}$ to $0.6\;{L}_{\mathrm{Edd}}$. Across this range the convexity rises with increasing metallicity, from near ${ \mathcal C }=1\%$ for $Z=0.000$ to ${ \mathcal C }=9\%$ for $z=0.020$. Above $0.6\;{L}_{\mathrm{Edd}}$ bursts with metallicities below solar converge towards ${ \mathcal C }=2\%$.

Zoom In Zoom Out Reset image size

Very few "concave" bursts (${ \mathcal C }\lt 0$) are observed in the KEPLER sample. Rather, $Z=0.0$ bursts have ${ \mathcal C }=0\%$ and convexity increases with metallicity. Variations in convexity within a run are small, with the mean uncertainty typically on the order of one percentage point. The burst-to-burst variation in convexity increases at high accretion rates, correlated with the onset of Case II burning, and also at low accretion rates, where a larger amount of the accreted hydrogen is depleted.

The KEPLER models show a far more constrained range of convexities within a run than those observed in single sources. Following Maurer & Watts (2008), we find the mean convexity of 4U 1636−536 to be ${{ \mathcal C }}_{\mathrm{avg}}=4.7\%$, with a $1\sigma$ spread of 10.8 percentage points, across a large range of persistent emissions. Even within a restricted range of persistent flux ($2\times {10}^{-9}\;\mathrm{erg}\ {\mathrm{cm}}^{-2}\;{\rm{s}}\lt {f}_{\mathrm{pers}}\lt 3\times {10}^{-9}\;\mathrm{erg}\ {\mathrm{cm}}^{-2}\;{\rm{s}}$) the spread in convexities is still notably larger than nuclear burning predicts. We also consider the variation in convexity shown by GS 1826−24, as it is known for showing consistent light curves dominated by nuclear processes (Heger et al. 2007a). GS 1826−24 shows a narrower distribution of convexities with ${{ \mathcal C }}_{\mathrm{avg}}=4.1\pm 2.9\%$, though the distribution of convexities is again wider than a single run would on average predict. Note that the error in observational convexity values here does not consider the error in convexity from each observation, as it is markedly smaller. We calculated this using the observed flux $\pm 1\sigma$, and found that for 4U 1636−536 and GS 1826−24, these varied on average by 0.7 and 0.3 percentage points from the convexity value calculated using the mean flux.

The larger spread in convexity for observations is in part due to the span of persistent fluxes (and thereby accretion rates) covered by observations, as convexity does show some variation with accretion rate. For both GS 1826−24 and 4U 1636−536, this is insufficient to explain the full range in convexities spanned by both sources. 4U 1636−536 shows convexities ranging from ${ \mathcal C }\approx -30\%$ to ${ \mathcal C }\approx 40\%$. No single composition shows this amount of variation in convexity across the range of accretion rates simulated, and nuclear burning does not produce noticeably concave bursts. Hence it is likely that physical processes (e.g., spreading) on the neutron star have a more significant impact upon the burst rise in 4U 1636−536 than nuclear processes.

GS 1826−24 shows a far more constrained range of convexities than 4U 1636−536; however, the observed dispersion in convexities is still more than typical within the models. This indicates that the physical processes that influence convexity are far more constrained in the GS 1826−24 system than in 4U 1636−536. In their comparison of GS 1826−24 to KEPLER models, Heger et al. (2007a) find that GS 1826−24 is likely accreting fuel with a solar metallicity at $\dot{M}=1.58\times {10}^{-9}\;{M}_{\odot }\ {\mathrm{yr}}^{-1}$. The predicted convexity from this comparison (${ \mathcal C }\approx 11\%$), however, is higher than is observed. This could indicate a physical process on GS 1826−24 that reduces the convexity in the rise, such as ignition at a certain latitude and burst spreading, or it could be due to systematic uncertainties in the KEPLER models.

We find α predicted by KEPLER models match the range of α seen observationally, with $\alpha \approx 40$ in the region $0.1\lesssim {l}_{\mathrm{acc}}\lesssim 1.0$. In general, α tends to grow with metallicity (Figure 6). At high accretion rates, variations in α between bursts obscure this trend, as fluence measurements become more variable when the recurrence time is of a similar order of magnitude to the burst duration (typically 10%–50% of t_b). The variation in α with accretion rate matches the expected trends for each metallicity, being approximately flat for no metallicity and increasing in curvature with metallicity following Galloway et al. (2004). For the data points with zero metallicity, the uncertainties at low accretion rates are sometimes very large; this is a manifestation of large variations in peak luminosity and recurrence time between bursts, which is accentuated in the cases where large uncertainties are seen.

Zoom In Zoom Out Reset image size

The tails of thermonuclear bursts can yield information about the cooling processes within the star, as well as provide a metric to enable comparisons with observations. We found that the power law index that best describes a burst tail increases dramatically at low accretion rates, as hydrogen is depleted and helium becomes more prevalent in a burst. This is quite noticeable at $Z=2\%$ and $Z=4\%$ (Figure 7). In PRE bursts, similiarly high power laws were found. If the heat in the burst is lost by radiation and the specific heat capacity is independent of temperature, the expected power law gradient is $\kappa =4/3$ (in't Zand et al. 2014), although an extended rp-process in the burst tail will obscure this cooling. This behavior is especially noticeable in the zero metallicity case, where artificially low power laws are caused by a two stage decay (see Figure 3). At low accretion rates where extended rp-process burning is minimal, power law gradients approach $\kappa =3$ and can be as high as $\kappa =6$ for some PRE bursts. Such high power laws can be caused by degenerate photons and electrons influencing the specific heat capacity of the atmosphere at high energies.

Zoom In Zoom Out Reset image size

Across most of the range of Case III burning, little variation in κ is observed with accretion rate (Figure 7), while the burst decay time increases. This is likely a consequence of the self-similarity of power laws with scaling. The trend we see suggests increased accretion leads to longer bursts with a similar decay power-law.

For metallicities of $Z=0.02$ and higher, the uncertainty in the power law constant grows considerably. This occurs as a result of the burst morphology showing large burst to burst variations as the burst duration grows similar to the recurrence time. In most morphological parameters, the measurements tend to show large burst to burst variations in the parameter value once the burst duration exceeds $\approx 10\%$ of the recurrence time.

Using the KEPLER simulations the nuclear origins of twin peak bursts can be investigated. As noted in simulations by Fisker et al. (2004), where two peaks occur, the first peak (which we refer to as the helium peak) arises from the rapid combustion of a helium layer at the bottom of the neutron star atmosphere, which convects rapidly to the surface, causing a large increase in surface brightness. The second peak (hereafter the primary peak) occurs due to a stalling of the rp-process at a nuclear waiting point. Notably, when there is insufficient helium accumulated deep in the neutron star atmosphere the peaks merge to create a typical mixed H/He burst.

In the KEPLER models, twin peaked structures are prominent in the higher metallicity models ($Z=0.10$) at low accretion rates (Figure 8), while the helium peak is suppressed at higher accretion rates, manifesting as a distinct two stage rise. Of the two peaks present at $Z=0.10$ the helium peak is larger than the primary peak at low $\dot{M}$, and also rises faster ($\sim 1\;{\rm{s}}$) than it does at higher accretion rates ($\sim 3\;{\rm{s}}$), bearing a slight resemblance to the quick rise of a helium flash. At solar metallicities the helium peak is not present as a local maximum but manifests rather as a kink in the burst rise (Figure 9) that is suppressed as accretion rate increases. In the intermediate case ($Z=0.04$) the helium peak is present at low accretion rates, though the models transition to pure He bursts before the helium peak grows to be equal in magnitude to the primary peak. For each case discussed, the burst train throughout the models consistently shows similar structures (twin-peaked or kinked) beyond the first few bursts in each train, which is markedly different from observations, where multi-peaked bursts rarely occur regularly. We have plotted the 11th burst for each simulation in Figures 8 and 9 to ensure that a consistent ashes layer has been found.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The presence of twin peaked structures in our models indicates that double peaked structures can have a nuclear origin; however, a larger than typical metallicity is required to produce these bursts consistently. The magnitude of the helium peak is largely dependent upon accretion rate. Lower accretion rates cause longer periods between bursts, depleting more hydrogen by hot CNO burning. The higher relative fraction of helium, especially at the base, leads to a more dominant initial helium powered rise.

The potential observability of a twin peaked burst, however, is not guaranteed, as spreading of the light curve in time due to the propogation of the nuclear burning front may blur the peaks. In observational sources that show repeated peaks, multiple peaks are not observed with regularity, possibly due to their obfuscation by the propagating burning front. Also possible is that most observational sources do not accrete sufficiently metallic material to display a pronounced helium peak; however, sometimes burning increases the local metallicity enough to cause occasional double peaked bursts as a consequence of compositional inertia.