import numpy as np
import matplotlib.pyplot as plt

Application: X-ray Burst Timing

The Rossi X-ray Timing Explorer (RXTE) has very high time resolution and has been used to observe a large number of X-ray sources. We’ll look at the data from the low-mass X-ray binary 4U 1728-34—this is an X-ray burst system.

In Strohmayer et al. 1996 it was shown that the neutron star spin rate can be seen in the Fourier transform of the lightcurve of the burst. Here we repeat the analysis.

We thank Tod Strohmayer for sharing the data from that paper. : 4u1728_burstdata.txt

We’ll read in the data—although it is stored in a multidimensional fashion in the file, it is actually a time-series, so we’ll flatten the data into a 1-d array.

data = np.loadtxt("4u1728_burstdata.txt").flatten()

N = len(data)
N

262144

We can start out by plotting all of the data.

fig, ax = plt.subplots()
ax.plot(data)

[<matplotlib.lines.Line2D at 0x7fdd2c602aa0>]

We’ll bin the data into fewer samples to build up the signal to noise. We can do this by reshaping it into a 2-d array and then summing over the columns.

binned_data = data.reshape(N // 256, 256).sum(axis=1)

fig, ax = plt.subplots()
ax.plot(binned_data)

[<matplotlib.lines.Line2D at 0x7fdd2c6f79a0>]

Now we see a much clearer X-ray lightcurve.

Next we want to take the FFT and look to see if there are any interesting frequencies with a lot of power

c_k = np.fft.rfft(data)

We want to get the physical frequencies corresponding to the Fourier coefficients

T = 32.0
kfreq = np.fft.rfftfreq(N) * N / T

Now we can plot the power spectrum. We normalize by \(2/N\)

fig, ax = plt.subplots()
ax.plot(kfreq, np.abs(c_k)**2 * 2 / N)

ax.set_xscale("log")
ax.set_yscale("log")

Here we see a signal around 300 Hz. We’ll bin the Fourier data to increase the signal to noise

c_k_binned = np.abs(c_k[1:]).reshape(int(len(c_k)/8), 8).mean(axis=1)
kfreq_binned = kfreq[1:].reshape(int(len(kfreq)/8), 8).mean(axis=1)

fig, ax = plt.subplots()
ax.plot(kfreq_binned, c_k_binned**2 * 2 / N)

ax.set_xlim(345, 380)
ax.set_ylim(0, 20)

(0.0, 20.0)

Now we see the strong signal at around 363 Hz—this is the same frequency identified as the neutron star rotation rate in the original paper.