[1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%load_ext snakeviz

import glob
import matplotlib.pyplot as plt
from astropy.table import Table
import numpy as np
from hendrics.io import load_events, save_pds, load_pds
from stingray.gti import create_gti_from_condition, gti_border_bins, time_intervals_from_gtis, cross_two_gtis
from stingray.utils import show_progress
from stingray.fourier import avg_cs_from_events, avg_pds_from_events, poisson_level, get_average_ctrate
from stingray import AveragedPowerspectrum, AveragedCrossspectrum, EventList

Spectral timing exploration

In this page, we will run a quicklook timing/spectral timing analysis of a NICER observation of an accreting black hole. We will not give a scientific interpretation, just pure exploration, showing off Stingray’s capabilities.

Load events and plot light curve

Let us take a look at the light curve

[2]:
def load_and_cleanup_events(fname):
    """Load data and apply GTIs"""
    events = EventList.read(fname)
    lc = events.to_lc(dt=1)
    lc.apply_gtis()
    plt.figure()
    plt.plot(lc.time, lc.counts)

    new_gti = create_gti_from_condition(lc.time, lc.counts > 0, safe_interval=1)
    lc.gti = new_gti
    lc.apply_gtis()
    plt.plot(lc.time, lc.counts)
    plt.title("Light curve")
    plt.xlabel(f"Time (s from {events.mjdref})")
    plt.ylabel(f"Counts/bin")
    events.gti = new_gti
    return events

fname = "data.hdf5"

# Load events. NB. for standard FITS files from HEASARC, it would be
# events = EventList.read(fname, "hea")
events = EventList.read(fname)
events.fname = fname

# Create light curve and apply GTIs
lc = events.to_lc(dt=1)
lc.apply_gtis()
plt.figure()
plt.plot(lc.time, lc.counts)
plt.title("Light curve")
plt.xlabel(f"Time (s from {events.mjdref})")
plt.ylabel(f"Counts/bin")

[2]:
Text(0, 0.5, 'Counts/bin')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_2_1.png

The light curve has zeros where it should not have them. Let us define new GTIs avoiding these zeros.

[3]:
new_gti = create_gti_from_condition(lc.time, lc.counts > 0, safe_interval=1)
lc.gti = new_gti
lc.apply_gtis()

[4]:
plt.figure()
plt.plot(lc.time, lc.counts)
plt.title("Light curve")
plt.xlabel(f"Time (s from {events.mjdref})")
plt.ylabel(f"Counts/bin")

[4]:
Text(0, 0.5, 'Counts/bin')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_5_1.png

Much better. Hereafter, we will use these new GTIs to process the data.

[5]:
events.gti = new_gti

Calculate periodogram and cross spectrum

Let us now take a look at the periodogram and the cross spectrum

[6]:
# Calculate the periodogram in fractional rms normalization
segment_size=50
dt=0.001
norm="frac"
pds = AveragedPowerspectrum.from_events(events, segment_size=segment_size, dt=dt, norm=norm)

# Calculate the Poisson noise level
ctrate = get_average_ctrate(events.time, events.gti, segment_size)
noise = poisson_level(norm, meanrate=ctrate)

# Rebin the periodogam
pds_reb = pds.rebin_log(0.02)
285it [00:00, 305.45it/s]
[7]:
plt.figure()

plt.plot(pds.freq, pds.power, drawstyle="steps-mid", color="grey", alpha=0.5)
plt.plot(pds_reb.freq, pds_reb.power, drawstyle="steps-mid", color="k")
plt.axhline(noise)
plt.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2 Hz^{-1}}$")

plt.figure()
plt.plot(pds.freq, (pds.power - noise) * pds.freq, drawstyle="steps-mid", color="grey", alpha=0.5)
plt.plot(pds_reb.freq, (pds_reb.power - noise) * pds_reb.freq, drawstyle="steps-mid", color="k")
plt.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2}$")

[7]:
Text(0, 0.5, '$\\mathrm{(rms / mean)^2}$')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_10_1.png
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_10_2.png

And the cross spectrum between the bands 0.3–5 keV and 5–12 keV

[8]:
ref_band = [5, 12]
sub_band = [0.3, 5]
events_ref = events.filter_energy_range(ref_band)
events_sub = events.filter_energy_range(sub_band)

cs = AveragedCrossspectrum.from_events(events_sub, events_ref, segment_size=segment_size, dt=dt, norm=norm)
cs_reb = cs.rebin_log(0.02)

285it [00:01, 202.98it/s]
[9]:
plt.figure()
# plt.plot(cs.freq, cs.power * cs.freq, drawstyle="steps-mid", color="grey", alpha=0.5)
plt.plot(cs_reb.freq, cs_reb.power * cs_reb.freq, drawstyle="steps-mid", color="k")
plt.errorbar(cs_reb.freq, cs_reb.power * cs_reb.freq, yerr=cs_reb.power_err * cs_reb.freq, fmt="none", color="r")
plt.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2}$")

[9]:
Text(0, 0.5, '$\\mathrm{(rms / mean)^2}$')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_13_1.png

With the cross spectrum we can explore the time lags versus frequency

[10]:
# Use shorter segments, rebin a little more heavily
cs = AveragedCrossspectrum.from_events(events_sub, events_ref, segment_size=4, dt=0.01, norm=norm)
cs_reb = cs.rebin_log(0.05)

lag, lag_e = cs_reb.time_lag()
plt.figure()
plt.errorbar(cs_reb.freq, lag, yerr=lag_e, fmt="o")
plt.xlabel("Frequency (Hz)")
plt.ylabel(f"Time lag ({sub_band[0]:g}-{sub_band[1]:g} keV vs {ref_band[0]:g}-{ref_band[1]:g} keV, in seconds)")
plt.axhline(0, ls="--")
plt.loglog()
plt.ylim([1e-4, None]);
# plt.xlim([None, 80])
# plt.legend();
3644it [00:00, 3733.07it/s]
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_15_1.png

Another interesting thing to measure is the coherence at different frequencies

[11]:
coh, coh_e = cs_reb.coherence()
plt.figure()
plt.errorbar(cs_reb.freq, coh, yerr=coh_e, fmt="o")
plt.xlabel("Frequency (Hz)")
plt.ylabel(f"Coherence ({sub_band[0]:g}-{sub_band[1]:g} keV vs {ref_band[0]:g}-{ref_band[1]:g} keV)")
plt.axhline(0, ls="--")
plt.loglog()
# plt.ylim([1e-4, None]);
# plt.xlim([None, 80])
# plt.legend();
[11]:
[]
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_17_1.png

Spectral timing

Now let us explore the spectral timing properties of this observation, with no physical interpretation, just for the sake of data exploration.

[12]:
from stingray.varenergyspectrum import CountSpectrum, CovarianceSpectrum, RmsSpectrum, LagSpectrum

Lags and coherence

Let us start with the lag spectrum with respect to energy, in different frequency bands. This might be confusing for people coming from other wavelengths, so let us specify that

  • “frequency” refers to the frequency of the variability.

  • “energy” refers to the photon energy.

The photons at 0.3-12 keV are modulated by oscillations and other stochastic noise up to ~100 Hz (see section above). As an example, we will now analyze the spectral timing properties using the variability up to 1 Hz and between 4 and 10 Hz.

[13]:
energy_spec = np.geomspace(0.3, 10, 31)
segment_size = 4
bin_time = 0.03
freq_interval = [0, 1]

# If not specified, the reference energy band is the whole band.

lagspec_0_1 = LagSpectrum(events, freq_interval=freq_interval,
                          segment_size=segment_size, bin_time=bin_time,
                          energy_spec=energy_spec)
energies = np.mean(lagspec_0_1.energy_intervals, axis=1)


100%|██████████| 30/30 [00:34<00:00,  1.15s/it]
[14]:
plt.figure()
plt.errorbar(energies, lagspec_0_1.spectrum, yerr=lagspec_0_1.spectrum_error, fmt='o', label="0-1 Hz")
plt.xlabel("Energy (keV)")
plt.ylabel("Time Lag (s)")
plt.semilogx()
[14]:
[]
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_24_1.png
[15]:
lagspec_4_10 = LagSpectrum(events, freq_interval=[4, 10],
                           segment_size=segment_size, bin_time=bin_time,
                           energy_spec=energy_spec)
energies = np.mean(lagspec_4_10.energy_intervals, axis=1)
energies_err = np.diff(lagspec_4_10.energy_intervals, axis=1).flatten() / 2


100%|██████████| 30/30 [00:35<00:00,  1.18s/it]
[16]:
plt.figure()
plt.errorbar(energies, lagspec_4_10.spectrum, xerr=energies_err, yerr=lagspec_4_10.spectrum_error, fmt='o', label="4-10 Hz")
plt.errorbar(energies, lagspec_0_1.spectrum, xerr=energies_err, yerr=lagspec_0_1.spectrum_error, fmt='o', label="0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Time lag (s)")
[16]:
Text(0, 0.5, 'Time lag (s)')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_26_1.png
[17]:
plt.figure()
plt.errorbar(energies, lagspec_4_10.spectrum * (2 * np.pi * 7),
             xerr=energies_err, yerr=lagspec_4_10.spectrum_error * (2 * np.pi * 7), fmt='o', label="4-10 Hz")
plt.errorbar(energies, lagspec_0_1.spectrum * (2 * np.pi * 0.5),
             xerr=energies_err, yerr=lagspec_0_1.spectrum_error * (2 * np.pi * 0.5), fmt='o', label="0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Phase lag (s)")
[17]:
Text(0, 0.5, 'Phase lag (s)')
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_27_1.png

Interesting: the low-frequency variability has much longer time lags than the high-frequency variability, but the phase lags are on the same order of magnitude.

Covariance and RMS spectrum

[18]:
covspec_0_1 = CovarianceSpectrum(events, freq_interval=[0, 1],
                                 segment_size=segment_size, bin_time=bin_time,
                                 energy_spec=energy_spec, norm="abs")
covspec_4_10 = CovarianceSpectrum(events, freq_interval=[4, 10],
                                  segment_size=segment_size, bin_time=bin_time,
                                  energy_spec=energy_spec, norm="abs")

100%|██████████| 30/30 [00:37<00:00,  1.26s/it]
100%|██████████| 30/30 [00:45<00:00,  1.51s/it]
[19]:
plt.figure()
plt.errorbar(energies, covspec_4_10.spectrum,
             xerr=energies_err, yerr=covspec_4_10.spectrum_error, fmt='o', label="4-10 Hz")
plt.errorbar(energies, covspec_0_1.spectrum,
             xerr=energies_err, yerr=covspec_0_1.spectrum_error, fmt='o', label="0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Absolute Covariance (counts / s)");
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_31_0.png

This covariance, plotted this way, mostly tracks the number of counts in each energy bin. To get an unfolded covariance, we need to use the response of the instrument. Another way is to plot the fractional covariance, normalizing by the number of counts in each bin.

To do this, we calculate the Count Spectrum and divide by it.

[20]:
countsp = CountSpectrum(events, energy_spec=energy_spec)
30it [00:05,  5.15it/s]
[21]:
plt.figure()
plt.errorbar(energies, covspec_4_10.spectrum / countsp.spectrum,
             xerr=energies_err, yerr=covspec_4_10.spectrum_error / countsp.spectrum, fmt='o', label="4-10 Hz")
plt.errorbar(energies, covspec_0_1.spectrum / countsp.spectrum,
             xerr=energies_err, yerr=covspec_0_1.spectrum_error / countsp.spectrum, fmt='o', label="0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Normalized Covariance (1 / s)");
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_35_0.png

Alternatively, we can calculate the Covariance Spectrum in fractional rms normalization

[22]:
covspec_0_1 = CovarianceSpectrum(events, freq_interval=[0, 1],
                                 segment_size=segment_size, bin_time=bin_time,
                                 energy_spec=energy_spec, norm="frac")
covspec_4_10 = CovarianceSpectrum(events, freq_interval=[4, 10],
                                  segment_size=segment_size, bin_time=bin_time,
                                  energy_spec=energy_spec, norm="frac")

100%|██████████| 30/30 [00:38<00:00,  1.28s/it]
100%|██████████| 30/30 [00:33<00:00,  1.13s/it]
[23]:
plt.figure()
plt.errorbar(energies, covspec_4_10.spectrum,
             xerr=energies_err, yerr=covspec_4_10.spectrum_error, fmt='o', label="4-10 Hz")
plt.errorbar(energies, covspec_0_1.spectrum,
             xerr=energies_err, yerr=covspec_0_1.spectrum_error, fmt='o', label="0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Fractional Covariance");
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_38_0.png

This should largely be equivalent to the RMS spectrum

[24]:
rmsspec_0_1 = RmsSpectrum(events, freq_interval=[0, 1],
                          segment_size=segment_size, bin_time=bin_time,
                          energy_spec=energy_spec, norm="frac")
rmsspec_4_10 = RmsSpectrum(events, freq_interval=[4, 10],
                           segment_size=segment_size, bin_time=bin_time,
                           energy_spec=energy_spec, norm="frac")

100%|██████████| 30/30 [00:10<00:00,  2.76it/s]
100%|██████████| 30/30 [00:11<00:00,  2.69it/s]
[25]:
plt.figure()
plt.errorbar(energies, covspec_4_10.spectrum,
             xerr=energies_err, yerr=covspec_4_10.spectrum_error, fmt='o', label="Cov. 4-10 Hz", alpha=0.5)
plt.errorbar(energies, covspec_0_1.spectrum,
             xerr=energies_err, yerr=covspec_0_1.spectrum_error, fmt='o', label="Cov. 0-1 Hz", alpha=0.5)
plt.errorbar(energies, rmsspec_4_10.spectrum,
             xerr=energies_err, yerr=rmsspec_4_10.spectrum_error, fmt='o', label="RMS 4-10 Hz")
plt.errorbar(energies, rmsspec_0_1.spectrum,
             xerr=energies_err, yerr=rmsspec_0_1.spectrum_error, fmt='o', label="RMS 0-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Fractional RMS");
../../_images/notebooks_Spectral_Timing_Spectral_Timing_Exploration_41_0.png

QED, except that the error bars in some points look underestimated. It is always recommended to test error bars with simulations, in any case, as analytic formulas are based on a series of assumptions (in particular, on the coherence) that might not be correct in real life.

[ ]: