NICER Discovers mHz Oscillations in the "Clocked" Burster GS 1826−238

Neutron stars in accreting X-ray binaries are known to exhibit phenomena associated with unstable nuclear burning on their surfaces. At present we know of 111 such systems¹⁶ (Galloway & Keek 2017). These objects produce hydrogen and helium-powered thermonuclear X-ray flashes, also known as Type I X-ray bursts (Strohmayer & Bildsten 2006), as well as the longer and rarer "superbursts," likely powered by carbon burning (Cumming & Bildsten 2001; Strohmayer & Brown 2002). A much smaller number of objects, five at present count, have shown low-frequency (≈7–10 mHz) oscillations very likely associated with the transition from stable to unstable burning, known as marginally stable nuclear burning (Heger et al. 2007b; Galloway & Keek 2017).

Revnivtsev et al. (2001) reported the first detections of mHz quasiperiodic oscillations (QPOs) in three accreting neutron star systems, 4U 1608−52, 4U 1636−536, and Aql X-1. They argued that the mHz oscillations were very likely associated with the neutron star surface and not the accretion flow. They also showed that they only occur within a narrow range of X-ray luminosity. Since then, mHz oscillations have also been reported from 4U 1323−619 (Strohmayer & Altamirano 2012) and the 11 Hz pulsar IGR J17480−2446 in the globular cluster Terzan 5 (Linares et al. 2012). We note that this object showed different behavior in its mHz oscillations compared to the others. During its outburst in 2011, as the flux increased, the recurrence time of bursts steadily decreased to about 3 minutes, finally appearing as oscillations in the light curve with this period. These oscillations also appeared at a much higher inferred accretion rate than in the other mHz sources (Linares et al. 2012; Galloway & Keek 2017).

Yu & van der Klis (2002) found that the frequency of kilohertz QPOs are anticorrelated with the brightness of the mHz oscillations in 4U 1608−52 and, because this relationship is opposite to the positive correlation of kHz QPO frequency with X-ray luminosity typically seen, argued that the flux variation associated with the mHz QPO must arise on the neutron star surface and not within the accretion flow. Altamirano et al. (2008) reported on extensive Rossi X-ray Timing Explorer (RXTE) observations of 4U 1636−536 that revealed a close connection between the mHz oscillations and the occurrence of X-ray bursts. They found that during intervals of mHz oscillations, as the frequency drifted down and dropped below about 7.5 mHz, the oscillations faded and X-ray bursting resumed, and similarly to Revnivtsev et al. (2001), they found no mHz oscillations immediately following the X-ray bursts.

More recently, Lyu et al. (2015) also studied episodes of mHz QPOs in 4U 1636−536 using XMM-Newton and RXTE data. Similarly to Altamirano et al. (2008), they found frequency drift of the mHz oscillations and suggested the drift timescale may be set by cooling of the deeper layers, as was previously argued by Keek et al. (2009) based on hydrodynamic calculations of helium burning with rotational mixing. Lyu et al. (2016) found that the bursts that occurred immediately after an episode of mHz oscillations preferentially showed a rising light curve shape with so-called positive "convexity." Such bursts show fast rise times, and this has been linked to a burst ignition location at or near the star's rotational equator (Cooper & Narayan 2007; Maurer & Watts 2008). Thus, Lyu et al. (2016) argued that these bursts, and the burning responsible for the mHz QPOs, may occur at the neutron star equator. Stiele et al. (2016) reported results from phase-resolved spectroscopy of the mHz oscillations in 4U 1636−536 obtained with XMM-Newton. They found an approximately constant color temperature with a pulse phase of ≈0.7 keV that they associated with the neutron star surface and argued that the oscillations result from a variation in the surface emitting area and not the temperature. They derived a surface emitting area consistent with that expected for a neutron star.

The physics relevant to marginally stable burning was recognized early on by Paczynski (1983) and has been further elucidated by Heger et al. (2007b). As the helium burning stability boundary is approached, the temperature dependence of the nuclear heating rate almost precisely balances that of the cooling rate. The result is a slow, quasiperiodic mode of burning, with the oscillation period given approximately by the geometric mean of the accretion and thermal timescales. For typical neutron star parameters and the relevant accretion rates, this gives an oscillation period close to 2 minutes, which is in good agreement with the range of oscillation frequencies reported for the mHz oscillations.

Theory predicts that for unstable burning associated with helium ignition, the boundary between stable and unstable burning should occur at local accretion rates close to the Eddington rate (Bildsten 1998). The fact that observational indications based on estimates of the X-ray luminosity suggest the onset of stability at accretion rates closer to ≈1/10 of the Eddington rate remains a major puzzle (van Paradijs et al. 1988; Cornelisse et al. 2003). It has been suggested that confinement of the accreted fuel and/or rotational mixing of the fuel to greater column depths could alleviate this apparent discrepancy (Keek et al. 2009). The fact that the mHz oscillations associated with marginal stability occur at a well-defined local accretion rate makes them, in principle, a particularly important probe of the accretion rate. Keek et al. (2014) found that the range of accretion rates in which marginally stable oscillations can occur was also sensitive to the nuclear reaction rates. Moreover, Heger et al. (2007b) showed that the oscillation frequency associated with marginal stability is sensitive to the hydrogen mass fraction in the accreted fuel, as well as the surface gravity. Thus, a more detailed understanding of these oscillations could lead to new probes of these quantities.

The accreting neutron star binary GS 1826−238 (hereafter GS 1826) is in many ways a prototype, having earned the moniker "clocked burster," for producing bursts with extreme regularity (Ubertini et al. 1999). The recurrence times of these regular bursts were shown to match well the predictions of theory for mixed H/He bursts (Galloway et al. 2004), and the light curve shapes were accurately modeled using theoretical calculations of rapid-proton (rp)-process burning of approximately solar composition fuel (Heger et al. 2007a), hence it has also been referred to as the "text-book" burster. Indeed, the bursts were so uniform in their properties that Zamfir et al. (2012) used them as "standard candles" to place constraints on the mass and radius of the neutron star.

The "clocked" bursts all occurred, however, while GS 1826 was in a "hard" spectral state. This is typically diagnosed with the use of X-ray colors (van der Klis 2006). Based on this, GS 1826 is classified as an "atoll" source, and the hard state is also referred to as the "island" state. Interestingly, Chenevez et al. (2016) recently reported on observations of GS 1826 in 2014 June with Swift and NuSTAR during a soft spectral state, including the detection of several X-ray bursts. They found that in this spectral state the X-ray burst recurrence times were no longer regular, and the bursts themselves also differed from the hard state bursts, being generally shorter in duration. These findings are indicative of a lower hydrogen fraction in the fuel and suggest that an additional source of stable burning occurs at the higher accretion rates in the soft state. They also detected the first photospheric radius expansion (PRE) burst from GS 1826 and used it to estimate the distance as 5.7 ± 0.2 kpc (assuming isotropic emission). Since then, the source has mostly remained in a soft state.

In this paper, we report the discovery of mHz oscillations from GS 1826 with the Neutron Star Interior Composition Explorer (NICER) that are very likely associated with marginally stable burning. The plan of the paper is as follows. We begin with a description of the observations and present the detection of mHz oscillations and an X-ray burst. Next, we explore the source accretion state by studying the X-ray colors, showing that the source was in a soft, "banana" state in the color–color diagram (Hasinger & van der Klis 1989). Next, we study the average mHz pulse profile and present the results of phase-resolved spectroscopy. We conclude with a brief summary and discussion of the implications of our findings for models of marginally stable burning and for constraining the neutron star properties of GS 1826.

NICER was installed on the International Space Station (ISS) in 2017 June and provides low background, high throughput (≈1900 cm² at 1.5 keV), fast timing observations across the 0.2–12 keV X-ray band (Gendreau et al. 2012). NICER achieves an absolute timing precision of ≈100 ns with the aid of a GPS receiver. NICER observed GS 1826 as part of a science team program with a major goal to detect and study X-ray bursts. Observations were obtained in the periods 2017 June 20–29, July 11–13, and September 9–15. Relevant "observation IDs" (obsids) include $00503101{nn}$ and $10503101{mm}$, where nn runs from 01 to 14 and mm from 01 to 16. We used HEASOFT version 6.23 and NICERDAS version 2018-03-01_V003 to process the data, and we applied the standard filtering criteria, which include time intervals with a pointing offset <0015,¹⁷ bright Earth limb angles >40°, dark Earth limb angles >30°, and outside the South Atlantic Anomaly. This resulted in total good exposures of 18, 6, and 10 ks for the June, July, and September epochs, respectively.

2.1. mHz QPOs and an X-Ray Burst

We found that a few data intervals showed enhanced backgrounds at energies above ≈8 keV, so to mitigate this we restricted the timing analysis to events with energies less than 7.5 keV. We note that because the observed count rate is decreasing with energy, only ≈0.5% of all observed events have energies greater than this limit. We computed light curves for each epoch using events in the range 0.4–7.5 keV. Figure 1 shows the resulting light curve, where the time axis is broken up into three panels, one each for the June, July, and September exposures, respectively. We identified a thermonuclear X-ray burst that occurred on 2017 June 22 at 22:31:42 UTC (labeled in Figure 1). The burst shows a hardness ratio evolution consistent with PRE, and also shows a shorter-duration light curve that differs from those of the hard state, "clocked" bursts (Galloway et al. 2004; Heger et al. 2007a), but is similar to bursts recently observed from GS 1826 in a high (soft) spectral state (Chenevez et al. 2016). We note that the publicly available data from the Monitor of All-sky X-ray Image (MAXI) and the Swift Burst Alert Telescope (BAT), are consistent with GS 1826 being predominantly in a soft spectral state from about 2016 December to the present. We will present a detailed spectral study of this burst in a subsequent paper.

Zoom In Zoom Out Reset image size

Visual inspection of the light curves from the first two dwells of the September data epoch (these dwells are the first in the right panel of Figure 1, and are also labeled "mHz QPO") suggested the presence of apparent pulsations with a few minute period. The light curves of these dwells in the 0.4–7.5 keV band, binned at 10 s, are shown in Figure 2. Pulses are rather clearly evident "by eye." To explore this apparent variability more quantitatively we computed power spectra using continuous intervals of 2200 s for each of these dwells. We again used events in the 0.4–7.5 keV band and sampled the light curves at 8192 Hz. The resulting power spectra, normalized as in Leahy et al. (1983), are shown in Figure 3, where the black and red curves correspond to the first and second dwells, respectively. Both spectra show strong, narrow peaks near 8 mHz, a frequency consistent with the few minute separation of pulses evident in Figure 2. These peaks are highly significant, for pure Poisson noise the chance probability to exceed a power value of 200 is 3.7 × 10⁻⁴⁴. From the power spectra shown in Figure 3 we estimate fractional amplitudes (rms) in the 0.4–7.5 keV band for the first and second dwells of 0.79% and 0.95%, respectively. These values for the frequency and amplitude are largely consistent with those previously reported for mHz oscillations in other sources (Revnivtsev et al. 2001).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mHz signal power from the first dwell is confined to a single Fourier bin, being effectively coherent over the length of the interval. The Fourier frequency of this bin is 8.18 mHz. In the second dwell, the peak power shifts down one Fourier bin, and significant power is spread over a few bins. This suggests a downward drift in the oscillation frequency from the first to second dwell. To quantify this drift we oversampled the power spectrum for each dwell by padding the light curve to 8000 s using the mean count rate determined from the good exposure (2200 s). We then found the mHz frequency associated with the peak power in each spectrum. From the frequency separation, Δν, of the peak power values we estimate a frequency drift rate of Δν/Δt ≈ −0.3 mHz hr⁻¹. Previous studies of mHz oscillations indicate that the frequency can both increase and decrease with time and that the frequency drift rate can vary (Altamirano et al. 2008; Lyu et al. 2014, 2015). The best-studied object in this regard is 4U 1636−536. The rate of decrease we estimate here for GS 1826 does not appear to be inconsistent with previous reports (Altamirano et al. 2008), that is, one can find time intervals reported in the literature for 4U 1636−536 that have approximately similar drift rates (see, for example, Figure 4 in Lyu et al. 2015).

After identifying the strong mHz pulsations described above, we also searched for similar signals in the remaining data. We note that many of the other individual dwells have exposures that are shorter than those of the 2017 September dwells with strong mHz oscillations. Therefore, not all observations are equally sensitive to such oscillations. Nevertheless, we searched for additional signals by computing power spectra of all dwells with exposures equal to or longer than 750 s (0.4–7.5 keV band). We detected mHz oscillation signals near 8 mHz in three other dwells, while for the remainder no signal was detected, and we obtained upper limits (90% confidence) to an oscillation in the 7–9 mHz frequency band that ranged from approximately 0.3%–0.5% (rms), with the precise value depending on the count rate and exposure of the particular dwell. Figure 4 shows power spectra from these three additional dwells with mHz oscillations. The spectra are labeled with the start date (MJD) of the dwell, and the number in parenthesis indicates the corresponding day in Figure 1. The successive spectra are displaced vertically by 100 for clarity. From figure bottom to top we estimate fractional amplitudes (rms) of 0.60%, 0.85%, and 0.65%. We also searched for mHz oscillations in the immediate aftermath of the X-ray burst detected on 2017 June 22, but no mHz signal was detected in a power spectrum computed from a 600 s interval (0.4–7.5 keV) following the burst. The upper limit on the amplitude is 0.5% (rms). Finally, we did not detect mHz oscillations in any of the dwells preceding the X-ray burst.

Zoom In Zoom Out Reset image size

2.2. Spectral State

To further explore the accretion state context of the mHz oscillations we calculated color–color and hardness–intensity diagrams for all of the data. We computed both a "soft" and "hard" color. For the former we determine the ratio of count rates in the bands 1.8–3.5 and 0.5–1.8 keV, while for the latter we use 5.2–6.8 and 3.5–5.2 keV. For the intensity we compute the count rate from 0.5–6.8 keV. We computed colors and intensities using 32 s intervals for all of the data. Figures 5 and 6 show the resulting color–color and hardness–intensity diagrams, respectively. In both figures, the blue (triangle) symbols represent data obtained from the two September epoch intervals when strong mHz oscillations were present, and the red (diamond) symbols were obtained from intervals immediately following the X-ray burst. A characteristic error bar is also shown in each plot to indicate the statistical precision.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further elucidate the relationship between the spectral colors and mHz oscillations we computed colors for all dwells in which oscillations were detected, as well as some of the longer dwells with no detections, and only upper limits. The colors measured for all of these dwells are shown in Figure 7, where the red diamond symbols represent those with upper limits, the blue square symbols represent those with mHz oscillation detections, and the size of the blue symbols is proportional to the oscillation amplitude (rms). The X-ray colors are defined in the same way as in Figure 5. These results suggest that the presence and amplitude of mHz oscillations are related to position in the color–color diagram, and therefore the mass accretion rate, as dwells with the highest amplitudes seem to occur preferentially at the higher values of soft and hard color. Such behavior is theoretically predicted for mHz oscillations associated with marginally stable burning (Heger et al. 2007b), as they are only evident close to the mass accretion rate at which burning stabilizes. It may be that we are seeing such a transition in Figure 7, since mHz oscillations are not detected at the lowest inferred accretion rates (the lower left portion of the color–color diagram). However, it is also possible that this behavior is similar to that seen in 4U 1636−536, where mHz oscillations can "come and go" while on the "banana" branch (Altamirano et al. 2008). The present data are too sparse to be conclusive, and more observations will be needed to distinguish between these two possibilities.

Zoom In Zoom Out Reset image size

Since the power spectra (Figure 3) from the two 2017 September dwells with strong mHz oscillations are consistent with rather coherent oscillations during these intervals, we computed phases for all events in each of these dwells separately using the frequencies corresponding to the peak Fourier power obtained from each dwell. To combine the two dwells we first found a constant phase offset for the second dwell that maximized the ${Z}_{1}^{2}$ signal power of the sum of both dwells (Buccheri et al. 1983). The resulting phase-folded pulse profile in the 0.4–7.5 keV band is shown in Figure 8. We fit a model of the form $A+B\sin (\phi -{\phi }_{0})$ and find an amplitude B/A = 1.12 ± 0.05%. The fit is good with a minimum χ² = 11.0 for 13 degrees of freedom, and the best-fitting model is also plotted in Figure 8 (red curve). To explore the energy dependence we also computed profiles for different energy bands and then fit them with the same model as above. Figure 9 shows the derived amplitudes, B/A, and phases, ϕ₀, versus energy. The amplitude shows a modest increase from about 0.7% at 1 keV to ≈2% at 7 keV. There is no strong evidence for a variation in ϕ₀ with energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To explore the pulse-phase spectrum we extracted two spectra, from regions around the pulse minimum and maximum. The phase ranges used for these extractions are indicated by the dashed and dashed–dotted vertical lines in Figure 8, respectively. In each case, we used a range of pulse phase equal to 0.25 cycles, so that each spectrum has an exposure of 0.25 × (2(2200)) = 1100 s. We used XSPEC version 12.9.1 with version 1.02 of the NICER instrument response, and we fit the spectra in the range 0.6–9 keV. We use a model comprising thermal comptonization (comptt in XSPEC) for the persistent accretion-driven emission and a blackbody (bbodyrad in XSPEC) to capture the neutron star thermal emission (Thompson et al. 2008; Chenevez et al. 2016). To model Galactic absorption we employ the tbabs model with abundances from Wilms et al. (2000), and we fix (freeze) the value of n_H = 0.4 × 10²² cm⁻² based on prior measurements, since the source is in the same (soft) state (Pinto et al. 2010; Chenevez et al. 2016). We employ a background spectrum extracted from NICER observations of RXTE background field number 5 (Jahoda et al. 2006), but we note that the spectrum is strongly source dominated across the fitted band. Finally, we add a 1% systematic error to accommodate current uncertainties in the NICER response.

As with many X-ray spectral modeling efforts, there is often no single, unique approach. The same is true in this case. In order to model either individual phase-resolved spectrum we find that both of the spectral components discussed above are required. However, the data do not have sufficient statistical quality to simultaneously fit for all parameters and determine unambiguously whether either of the components alone is associated with the mHz modulation. That is, one can either fix the persistent component (comptt) to be the same between the phase maximum and minimum spectra, or one can fix the blackbody (bbodyrad) component. Given the statistical quality of the present spectra one can achieve a similarly good fit using either approach. Below, we explore the approach that the persistent emission is not strongly affected by the mHz oscillation, and we assume that the difference between the two spectra is accounted for by the blackbody component. This appears reasonable given the modest 1%–2% amplitude of the modulation and the theoretical suggestion that such oscillations arise from the neutron star surface. However, we emphasize that this approach places constraints on the blackbody component under the assumption that it is the varying component, but it does not demonstrate unambiguously that it is the varying component.

Thus, we carry out a joint fit to both the minimum and maximum phase spectra, keeping the comptt parameters and the column depth, n_H, tied together, but allowing the blackbody parameters to vary independently between the two spectra.

We find this model provides a reasonable description of the observed spectra, with a reduced ${\chi }_{r}^{2}\approx 1.01$ (χ² = 1698 with 1676 degrees of freedom). The largest remaining residuals are centered near 2.3 keV, and this feature is associated with a known instrumental edge (see Ludlam et al. 2018 for further discussion of the NICER spectral response). Our spectral results are summarized in Table 1. Figure 10 shows the unfolded photon spectra, including the individual model components (top) and the fit residuals (bottom). The mean (unabsorbed) flux in the 0.6–9 keV band is 6.90 ± 0.01 × 10⁻⁹ erg cm⁻² s⁻¹. For the comptonization component parameters we find results that are qualitatively similar to Chenevez et al. (2016), but differ somewhat in the details. For the blackbody components we find that the fitted normalizations are consistent with remaining constant from pulse minimum to maximum, with a mean value of K_avg = 601.3. This normalization corresponds to an area, assuming a spherical emitting source and isotropic emission, with radius $R={({K}_{\mathrm{avg}}{d}_{10}^{2})}^{1/2}=13.98$ km, where d₁₀ is the distance in units of 10 kpc, and we have assumed d = 5.7 kpc (Chenevez et al. 2016). We note that this estimate ignores any atmospheric color corrections. Applying such a correction would increase the inferred radius by a factor of ${f}_{c}^{2}={({T}_{c}/{T}_{\mathrm{eff}})}^{2}$, where T_c and T_eff are the color (measured) and effective temperatures, respectively, and f_c is the so-called hardening factor. Based on model atmosphere calculations, f_c is likely to fall within the range ≈1.2–1.5 (Suleimanov et al. 2011; Medin et al. 2016), resulting in an increase in inferred radius by a factor of ≈1.4–2.2. However, given the other assumptions and uncertainties, the value is at least approximately consistent with expectations for surface emission from a neutron star.

Zoom In Zoom Out Reset image size

For the spectra at the maximum and minimum oscillation phases we find temperature values of kT = 0.701 ± 0.005 keV and kT = 0.684 ± 0.006 (both 1σ confidence), respectively. This suggests that the flux modulation associated with the mHz oscillation appears to be consistent with a change in the surface temperature, kT, of the neutron star thermal emission. However, we note that these values differ only at the ≈3σ level, and more data are needed to demonstrate this conclusively.

We report the detection with NICER of mHz oscillations in GS 1826 for the first time. The oscillation properties are generally consistent with those of the mHz QPOs observed in 4U 1608−52, 4U 1636−536, Aql X-1, and 4U 1323−619. The NICER data were obtained while GS 1826 was in a soft spectral state, as indicated by both MAXI (Sugizaki et al. 2011) and Swift-BAT (Krimm et al. 2013) long-term monitoring from about 2015 December to the present. Indeed, the current long-term trend appears to be a slow but steady increase in the 2–10 keV flux (based on the MAXI data) and very weak hard X-ray flux (from Swift-BAT). Our spectral color analysis appears consistent with this conclusion as well. The observed behavior (Figures 5–7) appears consistent with GS 1826 tracing out a portion of the so-called "banana" branch in the color–intensity diagram during the NICER observations. The mass accretion rate is generally inferred to increase from left to right along the arc of the "banana." To further test this we extracted spectra from several intervals that populate the lower left portion of the color–color diagram (see Figure 7), in order to estimate their fluxes and compare with the mean flux measured during the September epoch intervals with strong mHz oscillations (see Section 2.2). We find these spectra can be fit with the same model described above, and we find the unabsorbed flux (0.6–9 keV) in the range 5.0–5.2 × 10⁻⁹ erg cm⁻² s⁻¹, supporting the conclusion that the X-ray flux, and therefore mass accretion rate, increases from the lower left to the upper right in the color–color diagram (Figure 7).

Assuming the X-ray luminosity can be written as

$$\begin{eqnarray}&&{L}_{x}=\displaystyle \frac{{GM}\dot{M}}{(1+z)R}=4\pi {d}^{2}{f}_{x},\end{eqnarray} \tag{ 1 }$$

where M, R, ${(1+z)=(1-2{GM}/{c}^{2}R)}^{-1/2}$, and f_x are the neutron star mass, radius, surface redshift, and X-ray flux respectively, we estimate an average mass accretion rate, $\dot{M}$, for the September dwells with mHz oscillations of $\dot{M}=0.18{\dot{M}}_{\mathrm{Edd}}$, where we have assumed M = 1.4 M_⊙, R = 10 km, d = 5.7 kpc, and ${\dot{M}}_{\mathrm{Edd}}=1.7\times {10}^{-8}{M}_{\odot }$ yr⁻¹. This estimate also assumes that the accretion luminosity is radiated isotropically.

Previous studies have shown that mHz oscillations can be associated with bursting activity. For example, Altamirano et al. (2008) found episodes of mHz oscillations in 4U 1636−536 with decreasing frequency, and when the frequency dropped below about 7.5 mHz, the oscillations faded and bursting resumed; however, such behavior was not associated with all X-ray bursts (see their Figure 3). Indeed, the data for 4U 1636−536 suggest that intervals of mHz oscillations can come and go without an obvious trigger, as long as the source is in the appropriate accretion state, indicated by its position on the banana branch or in the transition to the banana branch (see Altamirano et al. 2008, Figure 1). For GS 1826 we find no evidence for mHz oscillations after the X-ray burst on June 22, which at face value is consistent with the behavior seen in 4U 1636−536. However, we emphasize that the present data are very sparse around the time of the X-ray burst. The only thing we can say with certainty is that we did not detect mHz oscillations in any of the NICER dwells preceding the X-ray burst or immediately after it.

Our results suggest the following long-term behavior in GS 1826. Prior to 2014 the source was almost exclusively in the hard state, with typical accretion rates lower than indicated by the present NICER observations. For example, comparing long-term variations in the persistent X-ray flux—a reasonable proxy for mass accretion rate—we find from Galloway et al. (2004) that the 2.5–25 keV persistent (absorbed) flux from RXTE observations was 2.2 × 10⁻⁹ erg cm⁻² s⁻¹ on 2002 July 29. Based on the NICER spectral modeling for the 2017 September epochs described above we estimate a flux in the same band of 4.5 × 10⁻⁹ erg cm⁻² s⁻¹, a factor of two or so higher. In making this estimate, we note that we have extrapolated the best-fit NICER model (Table 1) to 25 keV, that is, outside of NICER's nominal bandpass. However, the (absorbed) flux measured by NICER in the 2.5–9 keV band is 3.6 × 10⁻⁹ erg cm⁻² s⁻¹, which is already significantly larger than the 2.5–25 keV RXTE flux. Since the 2.5–9 keV NICER flux measurement is a firm lower bound to the 2.5–25 keV flux, this strongly suggests that the flux measured by NICER in 2017 September in the same band is significantly higher than the 2002 July flux measured with RXTE. This provides further support to the notion that the apparent long-term increase in X-ray flux likely represents an increasing mass accretion rate.

While in the hard state, the source was extensively observed with RXTE, but no mHz QPOs were ever observed. As the source continued to brighten and soften, the accretion rate became large enough to enter the accretion regime where mHz oscillations are observed in other LMXBs (see, for example, Altamirano et al. 2008, Figure 1). This overall behavior appears consistent with the observations of the mHz QPOs in 4U 1636−536, 4U 1608−52, and Aql X-1 and suggests we are seeing the same phenomenon in GS 1826. More observations of GS 1826 in its soft, higher accretion rate state will be needed to determine if it shows the full range of mHz oscillation phenomenology evident in other sources, including a closer connection with X-ray bursts.

There have been some reports of QPOs not linked with marginally stable burning, but with similar frequencies. In some of those cases the accretor is known or strongly suspected to be a black hole and not a neutron star. Examples are the QPOs identified in the black hole candidates H1743−322 (Altamirano & Strohmayer 2012), and LMC X-1 (Alam et al. 2014). In the case of H1743−322, the QPOs are seen in a hard spectral state, whereas the mHz QPOs reported here occur in the soft state of GS 1826. For LMC X-1 the observed QPO frequencies are higher, at ≈27 mHz, than any previously reported mHz oscillation associated with marginally stable burning. Based on this it appears unlikely to us that the mHz oscillations observed in GS 1826 are related to these QPOs.

Several accreting millisecond X-ray pulsars have also occasionally shown episodes of QPOs with mHz frequencies during some outbursts, an example being the 8 mHz "flaring" observed in IGR J00291+5934 (Ferrigno et al. 2017). In those QPOs, however, the modulation amplitude was >10% (rms) for photon energies less than 2 keV, which is substantially larger than is found for the mHz oscillations linked to marginally stable burning. Moreover, the energy spectrum softened at the peak of these flares, which is opposite to the behavior indicated in GS 1826. This argues rather strongly that the mHz oscillations reported here from GS 1826 are not related to these mHz "flaring" QPOs.

As discussed in Section 2.1 the oscillations reported here from GS 1826 have a rather high quality factor, that is, they appear to be rather coherent over the intervals they are observed. We emphasize, however, that these intervals (≈2 ks) are relatively short compared to the hours-long timescales over which mHz QPO frequencies have been observed to drift in several sources (Altamirano et al. 2008; Lyu et al. 2015). Additionally, current modeling of marginally stable burning can produce quite regular trains of mHz pulses, with the limiting factor here likely being temporal variations of the mass accretion rate (Heger et al. 2007b; Keek et al. 2014). Based on these arguments we think that the relatively high coherence seen in the mHz oscillations reported here is not inconsistent with their connection to marginally stable burning.

Figure 9 shows the energy dependence of the fractional (sinusoidal) modulation amplitude of the mHz oscillations in GS 1826. From about 1–3 keV an increasing trend is rather clear, with an apparent flattening above 3 keV. This appears to be at least qualitatively similar to the behavior seen in IGR J17480−2446—which shows an increasing trend in rms amplitude from about 2.5–10 keV based on RXTE data (Linares et al. 2012)—but opposite to that reported for 4U 1608−52 and 4U 1636−536 by Revnivtsev et al. (2001), also using RXTE. However, some caution is warranted here due to the different bandpasses of NICER and RXTE, as well as the extent to which uncertainties in backgrounds could influence derived amplitudes. A closer comparison of the energy dependence of the mHz oscillations reported here with NICER, and those detected with XMM-Newton should prove useful, as these instruments have more closely matched bandpasses.

We extracted spectra as a function of mHz pulsation phase, and find that they are consistent with the modulation being produced by a variation in the temperature of a thermal (blackbody) component at an effectively constant emitting area, though we emphasize again that this is not a unique spectral interpretation (see discussion in Section 3 above). Assuming a spherical emitting source and isotropic emission, the fitted blackbody normalization implies a radius of R ≈ 14 km. Given uncertainties in the source distance, the exact nature of the continuum describing the persistent emission, the extent to which disk reprocessing is present, as well as our neglect of color corrections, this value should only be taken as roughly indicative of the surface area associated with the mHz oscillations, and definitely not a precise measurement. The evidence that the oscillating flux can be associated with a thermal (blackbody) component with an area approximately consistent with that for a neutron star provides support for the interpretation of the oscillations as due to marginally stable nuclear burning; however, more data will be needed to establish this definitively.

We note that Stiele et al. (2016) reached different conclusions from a study of the phase-resolved spectra of mHz oscillations from 4U 1636−536. They argued that the flux modulation resulted from a variation in the emitting area at approximately constant temperature. While this is at odds with current theory, which finds that the temperature oscillates (Heger et al. 2007b), we emphasize that current theoretical calculations of marginally stable burning are all one-dimensional, so they essentially do not explore the possibility of variations in the burning area.

Here, we point out that there are several physical effects that suggest that lateral variations in the burning should be included in future modeling. For example, Heger et al. (2007b) showed that the oscillations associated with marginal stability occur in a very narrow range of mass accretion rates, $\dot{m}$, for a given surface gravity. They found the oscillations were present in a range of only 1% of the critical accretion rate. Note that it is the local accretion rate that is relevant for the nuclear burning. For a fast spinning neutron star the variations in effective surface gravity from equator to pole can result in a latitudinal change in the local accretion rate. Cooper & Narayan (2007) used this effect to explore the extent to which ignition of thermonuclear bursts could preferentially occur at different latitudes. However, this process would appear to be relevant for marginally stable burning as well, as a greater than 1% variation in the local accretion rate with latitude could influence the range of stellar latitudes for which marginally stable oscillations could occur. Indeed, this appears to be a possible mechanism for limiting the marginally stable burning to a belt in latitude (Lyu et al. 2016). As a rough guide, a star rotating at 582 Hz (close to the inferred spin frequency for 4U 1636−536), and with a mass of 1.4 M_⊙ and equatorial radius of 10 km, would have an ≈11% increase in effective surface gravity from equator to pole, using the results of AlGendy & Morsink (2014). This is certainly large enough to have a significant effect on the local accretion rate at the level required to influence marginally stable burning. Interestingly, for a spin rate of 11 Hz, appropriate for IGR J17480−2446 in Terzan 5, the variations in surface gravity are <0.004% equator to pole, and thus latitudinal variations would likely be negligible for such a slow spinner. This could perhaps account for some of the observed differences in the mHz oscillations of IGR J17480−2446 and fast spinners like 4U 1608−52, 4U 1636−536, and Aql X-1. As pointed out by Heger et al. (2007b), the critical accretion rate for marginally stable oscillations is itself a function of surface gravity, so the occurrence of oscillations, and their properties, will in principle depend on both effects; however, this discussion highlights that such latitudinal dependencies of the burning are likely, at least for fast spinners, and should be considered in future modeling of marginally stable oscillations.

Heger et al. (2007b) argued that the oscillation frequency associated with marginally stable burning should be sensitive to the neutron star surface gravity, which is proportional to M/R², and the mass fraction of hydrogen in the accreted matter (see their Figure 9). While their results on this score were based on one-zone models, the good overall agreement between predictions of their one-zone models and multi-zone hydrodynamic calculations with the KEPLER code (Heger et al. 2007b) suggests that these dependencies are robust. This is particularly promising for GS 1826 because of the extensive theoretical modeling already done on the "clocked" bursts (Heger et al. 2007a; Meisel 2018), which suggests a near solar composition for the accreted fuel. Interestingly, the bursts observed so far at the higher accretion rates associated with the soft spectral state have shorter durations, and several have reached higher peak fluxes, as is evidenced by the occurrence of PRE. These are both good indications of depletion of hydrogen in the fuel compared to the hard state, "clocked" bursts. This suggests that additional theoretical modeling of the soft state bursts, coupled with further comparisons to the "clocked" bursts, could place tighter constraints on the hydrogen fraction present in the fuel during episodes of mHz oscillations. It might then be possible to provide more robust limits on the neutron star surface gravity from measurements of the mHz oscillations.

This work was supported by NASA through the NICER mission and the Astrophysics Explorers Program. This research also made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. D.A. acknowledges support from the Royal Society. S.G. acknowledges the support of the Centre National d'Etudes Spatiales (CNES). We thank the anonymous referee for valuable efforts.

Facilities: NICER - , ADS - , HEASARC. -