# On the Frequency Evolution of X-ray Brightness Oscillations During Thermonuclear X-ray Bursts: Evidence for Coherent Oscillations

## 1 Introduction

Millisecond oscillations in the X-ray brightness during thermonuclear bursts, “burst oscillations”, have been observed from six low mass X-ray binaries (LMXB) with the Rossi X-ray Timing Explorer (RXTE) (see Strohmayer, Swank & Zhang et al. 1998 for a review). The presence of large amplitudes near burst onset combined with spectral evidence for localized thermonuclear burning suggests that these oscillations are caused by rotational modulation of thermonuclear inhomogeneities (see Strohmayer, Zhang & Swank 1997). The asymptotic pulsation frequency in the cooling tails of bursts from 4U 1728-34 are stable over year timescales, also supporting a coherent mechanism such as rotational modulation (Strohmayer et al. 1998a).

An intriguing aspect of these oscillations is the frequency evolution evident during many bursts. The frequency is observed to increase in the cooling tail, reaching a plateau or asymptotic limit (see Strohmayer et al. 1998a). However, Strohmayer (1999) has recently discovered an episode of spin down in the cooling tail of a burst from 4U 1636-53. Evidence of frequency change has been seen in five of the six burst oscillation sources and appears to be commonly associated with the physical process responsible for the pulsations. Strohmayer et. al (1997) have argued this evolution results from angular momentum conservation of the thermonuclear shell. The thermonuclear flash expands the shell, increasing its rotational moment of inertia and slowing its spin rate. Near burst onset the shell is thickest and thus the observed frequency lowest. The shell then spins back up as it recouples to the bulk of the neutron star as it cools. This scenario is viable as long as the shell decouples from the bulk of the neutron star during the thermonuclear flash and then comes back into co-rotation with it over the $`10`$ s of the burst fall-off. Calculations indicate that the $`10`$ m thick pre-burst shell expands to $`30`$ m during the flash (see Joss 1978; Bildsten 1995), which gives a frequency shift due to angular momentum conservation of $`2\nu _{spin}(20\mathrm{m}/R)`$, where $`\nu _{spin}`$ and $`R`$ are the stellar spin frequency and radius, respectively. For the several hundred Hz spin frequencies inferred from burst oscillations this gives a shift of $`2`$ Hz, similar to that observed.

In bursts where frequency drift is evident the drift broadens the peak in the power spectrum and produces quality values $`Q\nu _0/\mathrm{\Delta }\nu _{FWHM}300`$. In some bursts a relatively short train of pulses is observed during which there is no strong evidence for a varying frequency. A burst such as this from KS 1743-26 with 524 Hz oscillations yielded the highest coherence of $`Q900`$ yet reported in a burst oscillation (see Smith, Morgan & Bradt 1997).

In this Letter we investigate the time dependence of the frequency observed in bursts from 4U 1728-34 and 4U 1702-429. We show that in the cooling tails of bursts the pulse trains are effectively coherent. We show that with accurate modeling of the drift quality factors as high as $`Q4,000`$ are achieved in some bursts. We investigate the functional form of the frequency drift and show that a simple exponential “chirp” model works remarkably well. We use this model to search for significant power at the harmonics and first subharmonic of the strongest oscillation frequency in each source. Such searches are important in establishing whether the strongest oscillation frequency is the stellar spin frequency or its first harmonic, as appears now to be the case for 4U 1636-53 (see Miller 1999). The detection of harmonic signals or limits on them is also important in obtaining constraints on the stellar compactness (see Miller & Lamb 1998, Strohmayer et al. 1998b). We note that Zhang et al. (1998) have previously reported on a model for the frequency evolution during a burst from Aql X-1.

## 2 Modelling the Frequency Drift

To investigate the frequency evolution in burst data we use the $`Z_n^2`$ statistic (see Buccheri et al. 1983)

$$Z_n^2=2/N\underset{k=1}{\overset{n}{}}\left(\underset{j=1}{\overset{N}{}}\mathrm{cos}(k\varphi _j)\right)^2+\left(\underset{j=1}{\overset{N}{}}\mathrm{sin}(k\varphi _j)\right)^2.$$
(1)   
Here $`N`$ is the total number of photons in the time series, $`\varphi _j`$ are the phases of each photon derived from a frequency model, $`\nu (t)`$, vis. $`\varphi _j=2\pi _0^{t_j}\nu (t^{})𝑑t^{}`$, and $`n`$ is the total number of harmonics added together. For the burst oscillations, which are highly sinusoidal, we will henceforth restrict ourselves to $`n=1`$. This statistic is particularly suited to event mode data, since no binning is introduced. $`Z_1^2`$ has the same statistical properties as the well known Leahy normalized power spectrum, which for a Poisson process is distributed as $`\chi ^2`$ with 2 degrees of freedom. All of the bursts discussed here were observed with the Proportional Counter Array (PCA) onboard RXTE and sampled with 125 $`\mu `$s (1/8192 s) resolution.

For $`\nu (t)`$ we have investigated a number of functional forms, including; $`\nu (t)=\nu _0`$ (a constant frequency), $`\nu (t)=\nu _0(1+d_\nu t)`$ (a linearly increasing frequency), and $`\nu (t)=\nu _0(1\delta _\nu \mathrm{exp}(t/\tau ))`$ (an exponential “chirp”). For a given data set and frequency model we vary the parameters so as to maximize the $`Z_1^2`$ statistic. We then compare the maximum values from different models to judge which is superior in a statistical sense. Our aim is to both constrain the functional form of the frequency evolution and to determine whether the pulse train during all or a portion of a burst is coherent, or not. We judge the coherence of a given model by computing the quality factor $`Q\nu _0/\mathrm{\Delta }\nu _{fwhm}`$ from the width of the peak in a plot of $`Z_1^2`$ vs the frequency parameter $`\nu _0`$. We also compare the peak width to that expected for a coherent pulsation in data of the same length. A pulsation in a time series of finite extent produces a broadened peak in a power spectrum. The well known window function, $`W(\nu )=|\mathrm{sin}(\pi \nu T)/\pi \nu |^2`$, gives a width of $`\mathrm{\Delta }1/T`$, where $`T`$ is the length of the data. We also confirm that for a successful frequency model the integrated power under the $`Z_1^2`$ peak is consistent with that calculated assuming no frequency evolution.

### 2.1 Linear and Exponential Frequency Drift

To begin we demonstrate how a linear increase in frequency yields a significant improvement in the $`Z_1^2`$ statistic compared with a constant frequency model. We use the burst from 4U 1702-429 observed on July 26, 1997 at 14:04:19 UT, which we refer to as burst A (see Figure 4 in Markwardt, Strohmayer & Swank 1999). We used a 5.25 s interval during this burst to investigate the frequency evolution. In figure 1a we show results from our calculations of $`Z_1^2`$ for the constant frequency model (top panel) and the model with a linearly increasing frequency (bottom panel). In both cases the ordinate corresponds to the frequency parameter $`\nu _0`$ defined in the models. For the linear frequency model we found that $`Z_1^2`$ was maximized with $`d_\nu =1.264\times 10^3`$ s<sup>-1</sup>. Including the linear drift increased $`Z_1^2`$ from 88.48 to 271.4, a dramatic improvement of $`183`$ obtained with only 1 additional degree of freedom. The resulting $`Z_1^2`$ peak is also substantially narrower (see figure 1a), leaving no doubt that the pulsation frequency is increasing during this time interval.

The frequency evolution during bursts can also be explored with dynamic power spectra. Several such spectra have been presented elsewhere (see Strohmayer et al. 1998a, Strohmayer, Swank, & Zhang 1998, etc.). A striking behavior is that the pulsation frequency reaches an asymptotic limit in many bursts. Motivated by this behavior we investigated a simple exponential “chirp” model with a limiting frequency, $`\nu (t)=\nu _0(1\delta _\nu \mathrm{exp}(t/\tau ))`$. This model has three parameters, the limiting frequency $`\nu _0`$, the fractional frequency change, or “bite”, $`\delta _\nu `$, and the relaxation timescale, $`\tau `$. We fit this model to burst A and find a maximum $`Z_1^2`$ of 342.9, an increase of 71.5 in $`Z_1^2`$ over the linear frequency model. This is also a dramatically significant improvement in $`Z_1^2`$. We fit the peak in $`Z_1^2`$ vs. $`\nu _0`$ obtained with the chirp model to a gaussian in order to determine its width. Figure 1b shows the resulting fit. The peak is well described by a gaussian with a width $`\mathrm{\Delta }\nu _{fwhm}=0.201`$ Hz, which gives $`Q=\nu _0/\mathrm{\Delta }\nu _{fwhm}=1,641`$. We can compare this with the width caused by windowing, which for a 5.25 s interval gives a width (FWHM) of $`0.17`$ Hz.

We used the chirp model to investigate a sample of bursts from 4U 1702-429 and 4U 1728-34. We do not present here a systematic description of all observed bursts, rather, we demonstrate the main results with several illustrative examples. A burst from 4U 1702-429 observed on July 30, 1997 at 12:11:58 UT (burst B) revealed a $`12`$ s interval during which oscillations were detected. Our results using the chirp model for this burst are summarized in figure 2. Panel (a) shows a contour plot of the time evolution of the $`Z_1^2`$ statistic through the burst. It was computed by calculating $`Z_1^2`$ on a grid of constant frequency values using 2 s intervals with a new interval starting every 0.25 s, that is, assuming no frequency evolution. The burst countrate profile (solid histogram) and best fitting exponential chirp model (heavy solid line) are overlaid. The extent of the model curve defines the time interval used to fit the chirp model. The best model tracks the dynamic $`Z_1^2`$ contours remarkably well. Panel (b) compares $`Z_1^2`$ vs. $`\nu _0`$ for the constant frequency (dashed curve) and chirp models (solid curve). We again fit a gaussian to the peak calculated with the chirp model and find a width $`\mathrm{\Delta }\nu _{fwhm}=0.086`$ Hz, which yields $`Q=\nu _0/\mathrm{\Delta }\nu _{fwhm}=3,848`$ for this burst. This compares with a width of $`0.071`$ for a windowed pulsation of duration 12.5 s.

We carried out similar analyses to investigate the frequency evolution in bursts from 4U 1728-34. We again found that the chirp model provides a remarkably useful description of the frequency drift. Table 1 summarizes our results using the chirp model for several bursts from both 4U 1702-429 and 4U 1728-34.

We find the peaks obtained with the chirp model are only modestly broader than those expected for a coherent pulsation of the same length. Some of this additional width is likely due to the fact that pulsations are not present during the entirety of each interval examined. It is also likely that the chirp model is not the exact functional form of the frequency evolution, this is suggested by the broader wings of the $`Z_1^2`$ peaks computed for several bursts, however, the success of such a simple model argues strongly that the pulsations during the cooling tails of these bursts are phase coherent.

## 3 Harmonics and Subharmonics

Pulsations from a rotating hotspot can be used to place constraints on neutron star compactnesses (see Strohmayer et al. 1998b; Miller & Lamb 1998; and Miller 1999). The pulsation amplitude is constrained by the strength of gravitational light deflection. An observed amplitude places an upper limit on the compactness, $`GM/c^2R`$, because too compact stars cannot achieve the observed modulation amplitude. Further, an upper limit on the harmonic content places a lower limit on the compactness, since less compact stars produce more harmonic content, and at some limit the harmonics should become detectable.

In some models for the kHz QPO observed in the accretion driven X-ray flux from neutron star LMXB, the QPO frequency separation is closely related to the stellar spin frequency inferred from burst oscillations (see Miller, Lamb & Psaltis 1998; Strohmayer et al. 1996). In two sources, burst oscillations are seen with frequencies close to twice the kHz QPO frequency separation (Wijnands et al. 1997; Wijnands & van der Klis 1997; Mendez, van der Klis & van Paradijs ). Miller (1999) has reported evidence for a significant 290 Hz subharmonic of the strong 580 Hz pulsation seen in 4U 1636-539 (Zhang et al. 1996), suggesting that the strongest signal observed during bursts may actually be the first harmonic of the spin frequency, not the spin frequency itself. Based on these new results and the evidence for a beat frequency interpretation it is important to search for the subharmonic of the strongest signal detected during bursts.

We have shown that frequency drift during bursts can greatly smear out the signal power. We have also shown that simple models can recover a coherent peak. By modelling the drift we can make a much more sensitive search for harmonics. Moreover, we can coherently add signals from different bursts by first modelling their frequency evolution and then computing a total $`Z_1^2`$ by phase aligning each burst. We note that this procedure will also coherently add together power at any higher harmonics of the known signal. However, there will be a $`\pi `$ phase ambiguity of any signal at the first subharmonic (see Miller 1999).

We have added coherently the 330 Hz signals in all five bursts from 4U 1702-429 seen during our 1997, July observations (see Markwardt, Strohmayer & Swank 1999). We fit the chirp model to oscillations in each burst and then computed a total $`Z_1^2`$ by phase aligning them. Figure 3 shows the results of this analysis. The top panel shows the total $`Z_1^2`$ power at 330 Hz obtained by adding the bursts coherently. The peak value is $`1,400`$ and demonstrates that we have succesfully added the bursts coherently. The highest power for any burst individually was $`487`$. The two lower panels show the power at the first and second harmonics of the 330 Hz signal. We find no evidence for a significant signal at these or higher harmonics. To search for a signal at the 165 Hz sub-harmonic we computed a total $`Z_1^2`$ for each of 16 different combinations of the phases from each of the five bursts. Since there is a 2-fold ambiguity when coherently adding a subharmonic signal from two bursts, with a total of 5 we have $`2^4=16`$ possible combinations. We found no significant power at the subharmonic. We performed a similar analysis using 4 bursts from 4U 1728-34 which showed strong oscillations in their cooling tails, again we found no significant harmonic or subharmonic signals. The 90 % confidence upper limits on the signal power, $`Z_1^2`$, at the first harmonic in bursts from 4U 1702-429 and 4U 1728-34 are 5.8 and 1.8, respectively. These correspond to lower limits on the ratio of power at the fundamental to power at the first harmonic, $`h`$, of 242 and 556, respectively.

## 4 Discussion

In this work we have concentrated on the pulsations in the cooling tails of bursts. Bursts also show pulsations during the rising phase (Strohmayer, Zhang, & Swank 1997). We have not yet been able to show that the pulsations which begin near burst onset can be phase connected to those in the cooling tail with a simple model. To fully address this interesting question will require more sophisticated modelling than we have employed here. We will address this question in future work.

With the chirp model we find magnitudes of the frequency shift, $`\nu _0\delta _\nu `$ of $`23`$ Hz. These values are consistent with simple estimates based on angular momentum conservation using theoretical values for the pre- and post-burst thickness of bursting shells (Bildsten 1995). For the frequency relaxation timescale, $`\tau `$, we find a range of values from 1.7 to 4 s. Interestingly, different bursts from the same source can show markedly different decay timescales. For example, the two bursts from 4U 1728-34 summarized in table 1 show similar values for $`\nu _0`$ and $`\delta _\nu `$, but have decay timescales $`\tau `$ which differ by almost a factor of two. Of these two bursts, burst C had both a substantially greater peak flux and fluence. This seems consistent with the idea that the frequency increase is due to hydrostatic settling of the shell as it radiates away its thermal energy, however, study of more bursts is required to firmly establish such a connection.

If the angular momentum conservation argument is correct, it implies the existence of a shear layer in the neutron star atmosphere. In the chirp model the total amount of phase shearing is simply, $`\varphi _{shear}=\nu _0\delta _\nu \tau (1e^{T/\tau })`$, where $`T`$ is the length of the data interval. For the bursts examined here we find $`\varphi _{shear}48`$, so that the shell “slips” this many revolutions over the underlying neutron star during the duration of pulsations. The dynamics of this shear layer are no doubt complex. Given the physical conditions in the shell; the shear flows are characterized by a large Reynolds number, it is likely that dissipation of the shear velocity and recoupling will be dominated by turbulent momentum transport. Magnetic fields may also play a role as well. Shear layers can be unstable to Kelvin-Helmholtz instability, however, Bildsten (1998) has suggested that the shear may be stabilized by either thermal buoyancy or the mean molecular weight contrast. We urge new theoretical investigations to explore the mechanisms of recoupling to determine if such a shear layer can survive long enough to explain the persistence of pulsations for $`10`$ s, as well as the observed relaxation timescale.

For rotational modulation of a hotspot, the ratio, $`h`$, of signal power at the fundamental to that at the first harmonic is a function of the stellar compactness (see Miller & Lamb 1998), so that measurement of $`h`$ can be used to constrain the compactness. More compact stars have less harmonic content in their pulses and therefore larger $`h`$. Since the pulsations in the cooling tails of bursts are likely caused by a broad brightness anisotropy on the neutron star surface, and not a point spot, it will require more realistic modelling of such an emission geometry to use the limits on harmonic content derived here to constrain the stellar compactness. We will perform such modelling in future work.

## 5 Figure Captions

Figure 1a: $`Z_1^2`$ vs frequency parameter $`\nu _0`$ computed from a 5.25 s interval in burst A. The top panel was calculated with no frequency modulation, $`d_\nu =0`$, while the bottom panel was computed with a linear frequency increase of magnitude $`d_\nu =1.264\times 10^3`$ s<sup>-1</sup>.

Figure 1b: $`Z_1^2`$ peak and gaussian fit to burst A using the best fitting parameters of the chirp model. The solid line shows the gaussian fit to the peak. The derived peak centroid and width (FWHM) give $`Q\nu _0/\delta _\nu =1,641`$.

Figure 2a: Dynamic $`Z_1^2`$ spectrum for burst B. The contours show loci of constant $`Z_1^2`$ and were computed using 2 s intervals with a new interval starting every 0.25 s. The calculation was done on a grid of constant frequency points with no frequency modulation. The PCA countrate profile is shown (solid histogram) as well as the best fitting chirp model (heavy solid line). The interval used to fit the chirp model is denoted by the extent of the model curve.

Figure 2b: $`Z_1^2`$ vs frequency parameter $`\nu _0`$ computed from the time interval during burst B marked in figure 2a. The dashed curve shows $`Z_1^2`$ computed with no frequency modulation, ie. $`\delta _\nu =0`$, while the solid curve was computed with the best fitting chirp model.

Figure 3: Results of the search for harmonic signals by coherently adding 330 Hz signals in five bursts from 4U 1702-429. The top panel shows the total 330 Hz signal computed by adding all five bursts coherently. The two lower panels show $`Z_1^2`$ in the vicinty of the 1st and 2nd harmonics of the 330 Hz signal. There is no significant power at either harmonic.