# A Search for Aperiodic Millisecond Variability in Cygnus X-1

## 1 Introduction

Identifying and understanding short timescale variability in cosmic sources has repeatedly led to a better understanding of their fundamental nature and of the important physical processes present. For example, during the last decade, studies of fast time variability in low-mass X-ray binaries concentrated on quasi-periodic oscillations (QPOs) and various noise components frequencies up to 2 kHz. The study of QPOs and associated noise components led to a qualitatively more complete understanding of the accretion processes and the various omnipresent instabilities (van der Klis (1997)). Similar breakthroughs in our understanding of black hole candidates have yet to occur.

The dynamical and radiative timescales in the inner disk of accreting black hole candidates are predicted to be in the millisecond range. The thin disk models of Wallinder, Kato & Abramowicz (1992), scaled to stellar mass black holes, have local thermal and acoustic timescales $`<`$ 1 ms, and quasiperiodic variability in the X-ray emission is predicted at these timescales. Bao & $`Ø`$stgaard (1995) have numerically modeled orbiting shots in a geometrically thin accretion disk around a black hole including all relativistic effects. The shots, or “hot spots”, were simulated to radiate photons isotropically in their proper rest frames. For various shot distributions and different inclination angles, they find that the power density spectrum (PDS) exhibits a “cutoff” at the Keplerian frequency corresponding to the inner edge of the accretion disk. This cutoff is present for both optically thick and thin disks. The accretion model of Nowak & Wagoner (1995) also predicts a sharp cutoff in the PDS falling as $`f^5`$ for frequencies greater than the Keplerian frequency of the inner edge of the disk. The $`f^5`$ dependence arises from the three-dimensional hydrodynamic turbulent flow interior to the edge of the disk.

Advection-dominated disk models for black holes (Chakrabarti & Titarchuck (1996); Narayan (1996)) provide a self-consistent explanation for the energy spectra of hard and soft states of black hole candidate binary systems. In these types of models a standing shock can develop in the accretion flow at about 10 — 30 Schwarzschild radii (R<sub>Sch</sub>). The location of the shock defines an effective inner edge for both the disk and the halo components that can lead to abrupt changes in the PDS. Depending on the mass and angular momentum of the black hole these effects are predicted to be in the 3 — 100 Hz range. There are also interesting theoretical predictions of QPOs at a few hundred Hz arising from the special character of BH accretion (Perez et al.  (1997)).

Millisecond variability in Cyg X-1 has been reported twice. Rothschild et al. (1974) reported millisecond bursts in an observation of Cyg X-1 obtained with a rocket experiment. These bursts appeared as excess counts over that expected from Poisson statistics assuming that the Poisson expectation remains constant. However, the leakage of variability at lower frequencies ($`10`$ Hz) into the higher frequencies of interest ($`1000`$ Hz) invalidates this assumption (Press & Schechter (1974); Weisskopf & Sotherland (1978)). Indeed, when the pre-1978 literature is carefully reviewed, the analysis of Cyg X-1 timing spectra from a number of experiments show no conclusive evidence for millisecond variability (Weisskopf & Sotherland (1978)). More recent results show no model-independent evidence for millisecond variability (Lochner, Swank & Szymokowiak (1989)), except in the context of the shot model (Lochner, Swank & Szymokowiak ; Negoro et al.  (1995)). Lochner et al. (1991) used the phase portrait idea to determine parameters of a shot noise model. Using data from HEAO A-2 and EXOSAT, they find evidence for characteristic shot durations lasting from milliseconds to a few seconds. However, this analysis is model-dependent.

Meekins et al. (1984) (hereafter M84 ) worked to untangle leakage effects from slower time scales and developed a $`\chi ^2`$ method to claim detected variability at timescales of 0.3 – 3.0 millisecond in a HEAO A-1 observation of Cyg X-1 with 8 $`\mu `$s resolution. This result has been prominently quoted as the only strong evidence for variability in Cyg X-1 at millisecond time scales; for example, see van der Klis (1995) and Liang (1998). In this paper we present a reanalysis of HEAO A-1 data and an analysis of Rossi X-ray Timing Explorer (RXTE) data that contradicts the apparently clear observation in the PDS of millisecond power from Cyg X-1 of M84. From our analyses, we conclude that the observed millisecond power in the M84 PDS is due to either the known HEAO A-1 reset problem (Wood et al.  (1984)) or a previously unknown instrumental effect.

## 2 Analyses

### 2.1 Observations

We have analyzed archival observations of Cyg X-1 and the supernova remnant Cas A from the High Energy Astrophysics Observatory A-1 experiment (HEAO A-1). We have also analyzed new observations of Cyg X-1 made by the Rossi X-ray Timing Explorer (RXTE). The HEAO A-1 observations were made on 1978 May 7 while Cyg X-1 was in the hard (low) state. Cas A was observed on 1978 August 2. Data from this presumably Poisson source were used to model the response of the detector and to search for instrumental effects. Both sets of A-1 data were recorded using the high bit rate mode described below. The RXTE observations occurred on 1996 June 8, 1996 June 17, 1996 June 27 and 1996 July 12 as part of our approved RXTE Guest Investigator program. Similar observations have been previously published (Cui et al.  (1997); Belloni et al.  (1996)). Cyg X-1 was in its soft (high) state at the time of the RXTE observations. Table 1 shows the observation dates and times for all of the data used in our analyses.

We have used two techniques to analyze the HEAO A-1 data: the relative integral power method of M84 and the standard Fast Fourier transform (FFT) power spectrum method. The RXTE observations were analyzed using the FFT power spectrum method only. For the HEAO A-1 analyses, we have derived new methods to correct the data for both dead time and instrumental effects. For the RXTE analysis we use the standard RXTE dead time correction (Zhang et al. (1995)).

### 2.2 Analysis of the HEAO A-1 Data of Cas A and Cyg X-1

The HEAO A-1 data were recorded in the high bit rate (HBR) mode and consisted of a series of zeros and ones. A zero indicated no photons in the previous 8 $`\mu `$s and a one indicated that at least one photon was detected in the 8 $`\mu `$s interval. There was no energy information in this mode. The energy range covered by these observations is about 1 to 25 keV.

The data for Cas A were Fourier transformed to search for deviations from the expected Poisson source spectrum. For a Poisson source measured by a detector with no dead time, the expected spectrum is flat with a value 2, using the Leahy normalization (Leahy et al.  (1983)). Introducing dead time into the system slightly reduces the normalization value, but the shape remains relatively flat in the region in which we are interested. The observed Fourier power for Cas A-1 was not flat, but instead showed a broad “knee” in the spectrum, as shown in Figure 1. The distribution of differences in photon arrival times ($`\mathrm{\Delta }t_\gamma `$) showed a kink in the expected offset exponential distribution. This effect was modeled under the assumption that it was a previously uncorrected instrumental effect. The effect may have been unique to the HBR data or general to the HEAO A-1 data, but it would have been difficult to observe in the well-studied 5 and 320 ms binned data modes. In the HBR PDS it was apparent only at frequencies above about 100 Hz, which is the Nyquist frequency for the 5 ms data. It was not possible to determine the times between individual events for the 5 and 320 ms modes since the data were binned; therefore the kink in the $`\mathrm{\Delta }t_\gamma `$ distribution is not likely to be observable. We searched the 5 ms data for this effect and did not find any indications of it. Eadie et al. (1971) note that the hyperexponential function is applicable in situations where there is a mixture of exponential processes. We find that an offset hyperexponential function is a good representation of the HEAO A-1 HBR $`\mathrm{\Delta }t_\gamma `$ distribution:

$$f_H(t)=U(t\tau )(p_1\rho _1e^{\rho _1(t\tau )}+(1p_1)\rho _2e^{\rho _2(t\tau )})$$
(1)   
where $`U(t\tau )`$ is the Heavyside step function, $`\tau `$ is the dead time, $`\rho _1`$ and $`\rho _2`$ are the count rates for two Poisson processes and $`p_1`$ is the probability of generating an $`\mathrm{\Delta }t_\gamma `$ from the first Poisson process. The model was fitted to the Cas A data and the results are shown in Figure 2.

### 2.3 Power Spectrum Analysis of HEAO A-1 Data

We analyzed the Cyg X-1 data using a Fourier transform and observed a spectrum similarly shaped to the Cas A PDS. We again interpreted this as a manifestation of either the instrument reset problem or of a previously unreported instrumental effect. We determined the effective Poisson noise floor in the presence of the instrumental effect for a non-Poisson source, Cyg X-1. Using the following procedure we fit Equation 1 to the Cyg X-1 $`\mathrm{\Delta }t_\gamma `$ distribution. Random $`\mathrm{\Delta }t_\gamma `$ were drawn from Equation 1 as defined by the above fit parameters and accumulated to generate absolute times. This Monte Carlo time series was Fourier transformed in the same manner as the data were. The $`\chi ^2`$ of the Monte Carlo PDS to the data PDS was calculated for frequencies above 100 Hz. This procedure was repeated for a grid of parameter values whose origin was defined by the initial fit to the $`\mathrm{\Delta }t_\gamma `$ distribution. The resulting best fit PDS defined the effective Poisson noise floor. The resultant noise corrected PDS for the HEAO data is shown in Fig. 3. Fig. 4 shows the region of the PDS above 10 Hz to better examine the power at high frequencies. No statistically significant power above the noise floor is observed above 25 Hz. This is consistent with our assumption that power above 100 Hz is attributable to Poisson noise.

### 2.4 Relative integral power analysis

M84 derived a statistic, that they called the relative integral power, to quantify aperiodic variability. The details of their approach are described in Section III of their paper.

The new statistic, $`P_{rel}`$, of M84 is defined as the total discrete Fourier transform power of the mean subtracted time series divided by the square of the total number of x-ray counts, $`N^2`$, in the time series of length $`T`$ divided into $`m`$ equal length bins.

$$P_{rel}\frac{1}{N^2}\underset{j=\frac{m}{2}1}{\overset{\frac{m}{2}}{}}(|a_j|^2a_0^2)=\frac{\chi ^2}{N}$$
(2)   
where the $`a_j`$ are the standard Fourier coefficients

$$a_j=\underset{k=0}{\overset{m1}{}}x_ke^{2\pi ijk/m},$$
(3)   
$`x_k`$ is the number of events in the $`k`$th time bin, and $`a_j`$ is the Fourier coefficient at frequency $`f_j(=j/T)`$.

The distribution of power variability over all possible frequencies forms the Fourier power spectrum. With the Leahy normalization, this power spectrum is given by Leahy et al.(1983);

$$P_j=\frac{2|a_j|^2}{N},j=1,\mathrm{},m/21$$
(4)   
where $`N`$ is the total number of x-ray counts observed in the time interval 0$``$T.

The M84 analysis did not include corrections for dead time or instrumental effects. We have derived an approximation to the relative integral power that allows for simple corrections to equation 1 of M84 for these effects. We define $`\chi ^2`$ as,

$$\chi ^2=\underset{j=1}{\overset{\frac{m}{2}1}{}}P_j+\frac{1}{2}P_{\frac{m}{2}}$$
(5)   
where $`m`$ is the number of time bins in the data segment under consideration. As in M84, the 9 minutes of HEAO A-1 data is divided into $`L`$ contiguous data segments of width $`\mathrm{\Delta }t_{seg}`$, each containing $`m`$ (= 10) time bins, and various quantities were calculated. This is repeated with $`T\mathrm{\Delta }t_{seg}`$ = 0.3 ms, 1 ms, 3 ms, 10 ms, 30 ms and 100 ms.

The average of the $`\chi ^2`$ in Equation 5 over the entire ensemble of 10 bin data segments for a given $`\mathrm{\Delta }t_{seg}`$ can be approximated by

$$\chi ^2\frac{m1}{2}P$$
(6)   
where $`P`$ is the average Leahy-normalized power over the entire ensemble of 10 bin data segments ($`L`$ 9 minutes/$`\mathrm{\Delta }t_{seg}`$) and the set of frequencies $`f_j=1/\mathrm{\Delta }t_{seg},2/\mathrm{\Delta }t_{seg},\mathrm{}m/2\mathrm{\Delta }t_{seg}`$. Using equation 6 for the average $`\chi ^2`$ in equation 1 of M84 yields

$$P_{rel}\frac{\frac{m1}{2}[PP_{noise}]}{(N1)}$$
(7)   
for each $`\mathrm{\Delta }t_{seg}`$. The expected noise is simply proportional to the Poisson floor.

The M84 analysis of the HEAO A-1 HBR observation of Cyg X-1 found an excess of variability at time scales above 1 ms and a sharp cutoff at about 1 ms . We have reanalyzed the same observation using their method without dead time and instrumental corrections and have found excellent consistency with their results. There are some minor discrepancies which can be attributed to differences in the bin offsets and the use of a different digitization of the original analog tape. The comparison of our analysis to that of M84 is shown in Figure 5.

By ignoring instrumental effects, M84 chose a value of 2 for the noise floor, where 2 is the value at all frequencies of the Fourier transform of a Poisson source in the Leahy normalization. Using Equation 1 as the probability distribution for the noise floor, we have calculated the expected noise floor in Equation 7. We binned the data in a different manner than the original M84 analysis. The uncorrected results with the new binning are shown in Figure 6. The shape and normalization are in good agreement with that of the original M84 work, shown in Figure 5. Figure 7 shows the results of correcting for standard Poisson dead time. Note that the normalization is increased and that the peak is broader, enhancing the effect observed by M84. The results of our analysis, which has been corrected for dead time and instrumental effects, of the HEAO A-1 observation of Cyg X-1 are shown in Figure 8. We see no evidence for the previously reported rise in the relative integral power, once the corrections for the previously uncorrected instrumental effects or the manifestation of the known reset problem are applied.

### 2.5 Power spectrum Analysis of RXTE Data

We have analyzed four RXTE/PCA observations of Cyg X-1, see Table 1. The RXTE data were recorded with 4 $`\mu `$s time resolution. During all four observations the source was in the high (soft) state. Figure 9 shows the RXTE All Sky Monitor (ASM) light curve for Cyg X-1 around the time of our observations.

We binned the data into 50 microsecond bins. The light curves were divided into equal segments of 26 seconds. An FFT was performed on each data segment. The results were averaged over all segments and over equal logarithmic frequency intervals. We used the Leahy normalization for the PDS. The dead time corrected Poisson noise power was then subtracted from the PDS obtained to yield the remaining signal above the noise. To determine the Poisson noise floor, we calculated the Poisson power spectrum, correcting for nonparalyzable dead time using Equation 44 in Zhang et al. (1995) with a dead time of 10 microseconds (Zhang & Jahoda (1998)). Corrections were not made for the energy dependent dead time or for very large events. However, below 30 Hz, these corrections are not significant and can be ignored (Cui et al.  (1997)). Figure 10 shows the Poisson noise subtracted PDS for these observations. Figure 11 shows the 10–30 Hz region on a linear scale to show the behavior of the PDS as it approaches the limit imposed by the corrections.

Figure 12 shows the PDS’s for both RXTE and HEAO A-1 on the same axes.

## 3 Results

### 3.1 HEAO A-1

As reported above, we have discovered either an unknown instrumental effect or a manifestation of the known reset problem in the HEAO A-1 high bit rate data. This effect could not have been discovered using the binned 5 ms and 320 ms timing resolution modes of HEAO A-1 because their Nyquist frequencies are too low. The effect is observable only with the higher 8 $`\mu `$s time resolution of the HBR data or it is a mode-dependent problem that only occurs in the HBR data mode. Once this effect is taken into account, we observe no evidence for the rise or the sharp cutoff in the relative integral power for Cyg X-1 reported by M84.

We observe excess power with a 95% confidence level at frequencies below 25 Hz in the noise subtracted PDS from Cyg X-1 in its hard state. Above 30–40 Hz, the noise subtracted PDS is consistent with the null hypothesis. In the region where excess power is significant, we find that the spectral shape can be described by a power law spectrum with a break in the spectrum at 3 Hz. From 0.1 to 3 Hz the spectral index is $`1.20\pm 0.08`$ and above the 3 Hz the spectrum steepens to $`1.7\pm 0.2`$.

### 3.2 RXTE

Our results are consistent with previously published results (Cui et al.  (1997); Belloni et al.  (1996)). Below 30 Hz we find that the spectral shape can be described by a broken power law with the break occurring at about 10 Hz. Below 10 Hz the spectral index is 1.05$`\pm `$0.01 and between 10 and 30 Hz the spectral index steepens to 1.75$`\pm `$0.03.

The lack of corrections for very large events and energy dependent dead time in the standard RXTE corrections make it impossible to extend the search for excess power beyond about 30 Hz at this time. There is adequate data in the sample to extend the search to higher frequencies once these additional corrections are developed. Additional work on these corrections is necessary to exploit the full timing capabilities of RXTE to search for excess power using this method.

## 4 Conclusions

In light of our discovery of either a new instrumental effect or a correction for the known reset problem that accounts for the observed relative integral power of M84 there is no longer any evidence for model-independent aperiodic variability on millisecond timescales from Cyg X-1. This lack of observed variability does not rule out its existence. RXTE should have the capability to make measurements of excess power at millisecond time scales once the appropriate corrections are available.

Our results for the Cyg X-1 PDS’s are consistent with previous measurements in the hard and soft states. Previous measurements of Cyg X-1 often show a flat PDS below about 0.1 Hz, although the location of the maximum frequency varies from about 0.04 to 0.4 Hz. The minimum frequency studied here is only 0.1 Hz, and we see no indication of a flat PDS in our data. Our RXTE spectral shape determination is consistent with other RXTE observations (Cui et al.  (1997); Belloni et al.  (1996)) made around the same time. In both sets of observations the spectral shape is similar to that observed for most black-hole candidates in the same state (van der Klis (1995)).

In order to extend the range of the types of searches into the regime where they can start to impact the models discussed above, we must make several improvements to the data and the techniques. At present, the limiting factor the RXTE data is the lack of adequate background subtraction at high frequencies. This is being worked on (Zhang & Jahoda (1998)) but is not available. Cross checks on Poisson sources, as presented here using Cas A, are invaluable in searching for uncorrected dead time and instrument difficulties. To maximize the utility of the cross checks, the cross checking observations are best done in the same data taking modes and at around the same times as the the observations of the source under study. Unfortunately, the current available instrument data sets do not have these properties.

Work supported by Department of Energy contract DE-AC03-76SF00515, NASA RXTE Guest Investigator Grand, the Office of Naval Research, and Stanford University.