In this tutorial, we will run a quicklook spectrotemporal analysis of a NICER observation of one epoch of the 2018 outburst of the accreting black hole MAXI 1820+070, largely reproducing the results from, e.g., Wang et al. 2021, De Marco et al. 2021. We will not give a scientific interpretation, just pure exploration.
We will use the Stingray software package, at the version specified in the installation process.
Let us first install the correct software version. From the shell,
$ pip install stingray pyfftw
The source code is available in the official Github repository
[1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
import copy
import glob
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy.modeling import models
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
from stingray.modeling.parameterestimation import PSDLogLikelihood
params = {
'font.size': 7,
'xtick.major.size': 0,
'xtick.minor.size': 0,
'xtick.major.width': 0,
'xtick.minor.width': 0,
'ytick.major.size': 0,
'ytick.minor.size': 0,
'ytick.major.width': 0,
'ytick.minor.width': 0,
'figure.figsize': (6, 4),
"axes.grid" : True,
"grid.color": "grey",
"grid.linewidth": 0.3,
"grid.linestyle": ":",
"axes.grid.axis": "y",
"axes.grid.which": "both",
"axes.axisbelow": False,
'axes.labelsize': 8,
'xtick.labelsize': 8,
'ytick.labelsize': 8,
'legend.fontsize': 8,
'legend.title_fontsize': 8,
'figure.dpi': 300, # the left side of the subplots of the figure
'figure.subplot.left': 0.195, # the left side of the subplots of the figure
'figure.subplot.right': 0.97, # the right side of the subplots of the figure
'figure.subplot.bottom': 0.145, # the bottom of the subplots of the figure
'figure.subplot.top': 0.97, # the top of the subplots of the figure
'figure.subplot.wspace': 0.2, # the amount of width reserved for space between subplots,
# expressed as a fraction of the average axis width
'figure.subplot.hspace': 0.2, # the amount of height reserved for space between subplots,
# expressed as a fraction of the average axis height
}
mpl.rcParams.update(params)
Load events and plot light curve¶
Let us take a look at the light curve. We load the NICER event list into a stingray.EventList object, and create a stingray.Lightcurve from it.
[2]:
fname = "ni1200120106_0mpu7_cl_bary.evt.gz"
events = EventList.read(fname, "hea")
events.fname = fname
/home/pupperemeritus/anaconda3/lib/python3.9/site-packages/stingray/io.py:235: UserWarning: Column energy not found
warnings.warn('Column ' + a + ' not found')
[3]:
# Create light curve and apply GTIs
lc_raw = events.to_lc(dt=1)
lc_raw.apply_gtis()
plt.figure()
plt.plot(lc_raw.time, lc_raw.counts, color="k")
plt.title("Light curve")
plt.xlabel(f"Time (s from {events.mjdref})")
plt.ylabel(f"Counts/bin")
[3]:
Text(0, 0.5, 'Counts/bin')
The light curve seems reasonably clean, with no need for further cleaning. Otherwise, we would have to filter out, e.g. flares or intervals with zero counts, doing something along the lines of:
new_gti = create_gti_from_condition(lc_raw.time, lc_raw.counts > 0, safe_interval=1)
lc = copy.deepcopy(lc_raw)
lc.gti = new_gti
lc.apply_gtis()
plt.figure()
plt.plot(lc_raw.time, lc_raw.counts, color="grey", alpha=0.5, label="Raw")
plt.plot(lc.time, lc.counts, color="k", label="Cleaned")
plt.title("Light curve")
plt.xlabel(f"Time (s from {events.mjdref})")
plt.ylabel(f"Counts/bin")
plt.legend();
events.gti = new_gti
Calculate periodogram and cross spectrum¶
Let us now take a look at the periodogram and the cross spectrum. The periodogram will be obtained with Bartlett’s method: splitting the light curve into equal-length segments, calculating the periodogram in each, and then averaging them into the final periodogram.
We will use the fractional rms normalization (sometimes referred to as the Belloni, or Miyamoto, normalization, from the papers Belloni & Hasinger 1990, Miyamoto et al. 1992). The background contribution is negligible and will be ignored.
Note: since the fractional rms normalization uses the mean count rate, the final result changes slightly if the normalization is applied in the single periodograms from each light curve segment, with the count rate of each chunk, or on the averaged periodogram, using the average count rate of the full light curve. We choose the second option (note the use_common_mean=True).
We will first plot the periodogram as is, in units of \((\mathrm{rms/mean)^2\,Hz^{-1}}\).
Then, from the periodogram, we will subtract the theoretical Poisson noise level of \(2/\mu\), where \(\mu\) is the mean count rate in the observation, and we will multiply the powers by the frequency, to have the periodogram in units of \((\mathrm{rms/mean)^2}\)
In both cases, we will rebin the periodogram geometrically, averaging more bins at larger frequencies, in order to lower the noise level.
[4]:
# Calculate the periodogram in fractional rms normalization.
# Length in seconds of each light curve segment
segment_size=50
# Sampling time of the light curve: 1ms, this will give a Nyquist
# frequency of 0.5 / dt = 500 Hz.
dt=0.001
# Fractional rms normalization
norm="frac"
pds = AveragedPowerspectrum.from_events(
events, segment_size=segment_size, dt=dt,
norm=norm, use_common_mean=True)
# Calculate the mean count rate
ctrate = get_average_ctrate(events.time, events.gti, segment_size)
# Calculate the Poisson noise level
noise = poisson_level(norm, meanrate=ctrate)
# Rebin the periodogam
pds_reb = pds.rebin_log(0.02)
65it [00:00, 65.69it/s]
[5]:
plt.figure()
plt.plot(pds.freq, pds.power, drawstyle="steps-mid", color="grey", alpha=0.5, label="PDS")
plt.plot(pds_reb.freq, pds_reb.power, drawstyle="steps-mid", color="k", label="Rebinned PDS")
plt.axhline(noise, ls=":", label="Poisson noise level")
plt.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2 Hz^{-1}}$");
plt.legend()
plt.figure()
plt.plot(pds.freq, (pds.power - noise) * pds.freq, drawstyle="steps-mid", color="grey", alpha=0.5, label="PDS")
plt.plot(pds_reb.freq, (pds_reb.power - noise) * pds_reb.freq, drawstyle="steps-mid", color="k", label="Rebinned PDS")
plt.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2}$");
plt.legend();
We will now do the same with the cross spectrum between the bands 0.3–5 keV and 5–12 keV.
In this case, there is no need to subtract the Poisson noise level, as it is zero in the cross spectrum, provided that the energy bands do not overlap.
[6]:
ref_band = [1.5, 3]
sub_band = [0.5, 1]
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)
65it [00:00, 112.32it/s]
/home/pupperemeritus/anaconda3/lib/python3.9/site-packages/stingray/fourier.py:720: RuntimeWarning: invalid value encountered in sqrt
dRe = dIm = dG = np.sqrt(power_over_2n * (seg_power - frac))
/home/pupperemeritus/anaconda3/lib/python3.9/site-packages/stingray/fourier.py:722: RuntimeWarning: invalid value encountered in sqrt
dphi = np.sqrt(power_over_2n * (seg_power / (Gsq - bsq) -
/home/pupperemeritus/anaconda3/lib/python3.9/site-packages/stingray/crossspectrum.py:2761: UserWarning: Some error bars in the Averaged Crossspectrum are invalid.Defaulting to sqrt(2 / M) in Leahy norm, rescaled to the appropriate norm.
warnings.warn(
[7]:
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.loglog()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2}$");
/home/pupperemeritus/.local/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1333: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
/home/pupperemeritus/.local/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1333: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
Periodogram modeling¶
This periodogram has a number of broad components, that can be approximated by Lorentzian curves. Let us try to model it.
[8]:
pds = AveragedPowerspectrum.from_events(events, segment_size=segment_size, dt=dt, norm="leahy")
pds_reb = pds.rebin_log(0.02)
65it [00:00, 72.65it/s]
We will model the periodogram using the maximum likelihood estimation from Barret & Vaughan 2012.
For periodograms averaged over \(L\) independent segments and \(M\) independent neighbouring frequencies,
where \(\theta\) are the model parameters, \(P_j\) are the periodogram values, \(S_j\) the model of the underlying signal, \(c(2ML)\) is a factor independent of \(P_j\) or \(S_j\), and thus unimportant to the parameter estimation problem considered here (it only scales the likelihood, but does not change its shape).
For non-uniformly binned periodograms, the factor \(ML\) should go inside the sum:
This is the formula that we will apply here.
Let us now create an initial model that more or less describes the periodogram
[9]:
fit_model = models.Lorentz1D(x_0=0.04, fwhm=0.15, amplitude=7000) + \
models.Lorentz1D(x_0=0.2, fwhm=3, amplitude=300)
plt.figure()
plt.plot(pds_reb.freq, (pds_reb.power - 2) * pds_reb.freq, drawstyle="steps-mid", color="k", label="Rebinned PDS")
plt.plot(pds.freq, fit_model(pds.freq) * pds.freq, color="r", label="Starting Model")
for mod in fit_model:
plt.plot(pds.freq, mod(pds.freq) * pds.freq, color="r", ls=":")
plt.semilogx()
plt.xlim([pds.freq[0], pds.freq[-1]])
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathrm{(rms / mean)^2}$");
plt.legend();
plt.ylim([0, None])
[9]:
(0.0, 488.21599547079995)
We will now add a constant at the Poisson noise level (2 in Leahy normalization) and fit using the Maximum Likelihood estimation in stingray
[10]:
from stingray.modeling import PSDParEst
fit_model = models.Const1D(amplitude=2) + fit_model
parest = PSDParEst(pds_reb, fitmethod="L-BFGS-B", max_post=False)
loglike = PSDLogLikelihood(
pds_reb.freq, pds_reb.power, fit_model, m=pds_reb.m)
res = parest.fit(loglike, fit_model.parameters)
fitmod = res.model
# The Poisson noise level was the first parameter.
poisson = fitmod.parameters[0]
print(res.p_opt)
[1.95227938e+00 6.97518942e+03 4.11961192e-02 1.42093997e-01
2.98070633e+02 4.06300000e-01 2.65743398e+00]
[11]:
plt.figure()
gs = plt.GridSpec(2, 1, hspace=0)
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1], sharex=ax0)
ax0.plot(pds_reb.freq, (pds_reb.power - poisson) * pds_reb.freq, drawstyle="steps-mid", color="k", label="Rebinned PDS")
ax0.plot(pds.freq, (fitmod(pds.freq) - poisson) * pds.freq, color="r", label="Best Model")
for mod in fitmod[1:]:
ax0.plot(pds.freq, mod(pds.freq) * pds.freq, color="r", ls=":")
ax0.set_xlabel("Frequency (Hz)")
ax0.set_ylabel(r"$\mathrm{(rms / mean)^2}$");
ax0.legend();
ax1.plot(pds_reb.freq, (pds_reb.power - poisson) * pds_reb.freq, drawstyle="steps-mid", color="k", label="Rebinned PDS")
ax1.plot(pds.freq, (fitmod(pds.freq) - poisson) * pds.freq, color="r", label="Best Model")
for mod in fitmod[1:]:
ax1.plot(pds.freq, mod(pds.freq) * pds.freq, color="r", ls=":")
ax1.set_xlabel("Frequency (Hz)")
ax1.set_ylabel(r"$\mathrm{(rms / mean)^2}$");
ax1.loglog()
ax1.set_ylim([1e-1, None]);
ax1.set_xlim([pds.freq[0], pds.freq[-1]]);
Lags and coherence¶
With the cross spectrum we can explore the time lags versus frequency
[12]:
# Use shorter segments, rebin a little more heavily
cs = AveragedCrossspectrum.from_events(events_sub, events_ref, segment_size=2, dt=0.01, norm=norm)
cs_reb = cs.rebin_log(0.4)
lag, lag_e = cs_reb.time_lag()
2627it [00:00, 2906.20it/s]
[13]:
plt.figure()
plt.errorbar(cs_reb.freq, lag, yerr=lag_e, fmt="o", color="k")
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.semilogx()
# plt.ylim([1e-4, None]);
# plt.xlim([None, 80])
# plt.legend();
[13]:
[]
Another interesting thing to measure is the coherence at different frequencies
[14]:
coh, coh_e = cs_reb.coherence()
plt.figure()
plt.errorbar(cs_reb.freq, coh, yerr=coh_e, fmt="o", color="k")
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();
[14]:
[]
Spectral timing¶
Now let us explore the spectral timing properties of this observation, with no physical interpretation, just for the sake of data exploration.
[15]:
from stingray.varenergyspectrum import CountSpectrum, CovarianceSpectrum, RmsSpectrum, LagSpectrum
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.
From Kara+2019, figure 3
[16]:
energy_spec = np.geomspace(0.5, 10, 41)
segment_size = 10
bin_time = 0.001
freq_interval = [3, 30]
ref_band=[0.5, 10]
# If not specified, the reference energy band is the whole band.
lagspec_3_30 = LagSpectrum(events, freq_interval=freq_interval,
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, ref_band=ref_band)
energies = lagspec_3_30.energy
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:57<00:00, 1.44s/it]
[17]:
plt.figure()
plt.errorbar(energies, lagspec_3_30.spectrum * 1e4, yerr=lagspec_3_30.spectrum_error * 1e4, fmt='o', label="3-30 Hz", color="k")
plt.xlabel("Energy (keV)")
plt.ylabel("Time Lag ($10^{-4}$ s)")
plt.xlim([0.5, 10])
plt.semilogx()
[17]:
[]
[18]:
lagspec_01_1 = LagSpectrum(events, freq_interval=[0.1, 1],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, ref_band=ref_band)
energies = lagspec_01_1.energy
energies_err = np.diff(lagspec_01_1.energy_intervals, axis=1).flatten() / 2
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:54<00:00, 1.37s/it]
[19]:
plt.figure()
plt.errorbar(energies, lagspec_01_1.spectrum, xerr=energies_err, yerr=lagspec_01_1.spectrum_error, fmt='o', label="0.1-1 Hz")
plt.errorbar(energies, lagspec_3_30.spectrum, xerr=energies_err, yerr=lagspec_3_30.spectrum_error, fmt='o', label="3-30 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Time lag (s)")
[19]:
Text(0, 0.5, 'Time lag (s)')
[20]:
freq_01_1 = (1 + 0.1) / 2 * 2 * np.pi
freq_3_30 = (3 + 30) / 2 * 2 * np.pi
plt.figure()
plt.errorbar(energies, lagspec_01_1.spectrum * freq_01_1 , xerr=energies_err, yerr=lagspec_01_1.spectrum_error * freq_01_1, fmt='o', label="0.1-1 Hz")
plt.errorbar(energies, lagspec_3_30.spectrum * freq_3_30, xerr=energies_err, yerr=lagspec_3_30.spectrum_error * freq_3_30, fmt='o', label="3-30 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Phase lag (rad)")
[20]:
Text(0, 0.5, 'Phase lag (rad)')
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¶
[21]:
covspec_3_30 = CovarianceSpectrum(events, freq_interval=[3, 30],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="abs", ref_band=ref_band)
covspec_01_1 = CovarianceSpectrum(events, freq_interval=[0.1, 1],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="abs", ref_band=ref_band)
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:55<00:00, 1.40s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:55<00:00, 1.40s/it]
[22]:
plt.figure()
plt.errorbar(energies, covspec_3_30.spectrum,
xerr=energies_err, yerr=covspec_3_30.spectrum_error, fmt='o', label="3-30 Hz")
plt.errorbar(energies, covspec_01_1.spectrum,
xerr=energies_err, yerr=covspec_01_1.spectrum_error, fmt='o', label="0.1-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Absolute Covariance (counts / s)");
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.
[23]:
countsp = CountSpectrum(events, energy_spec=energy_spec)
40it [00:08, 4.47it/s]
[24]:
plt.figure()
plt.errorbar(energies, covspec_3_30.spectrum / countsp.spectrum,
xerr=energies_err, yerr=covspec_3_30.spectrum_error / countsp.spectrum, fmt='o', label="3-30 Hz")
plt.errorbar(energies, covspec_01_1.spectrum / countsp.spectrum,
xerr=energies_err, yerr=covspec_01_1.spectrum_error / countsp.spectrum, fmt='o', label="0.1-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Normalized Covariance (1 / s)");
Alternatively, we can calculate the Covariance Spectrum in fractional rms normalization
[25]:
covspec_01_1 = CovarianceSpectrum(events, freq_interval=[0.1, 1],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="frac")
covspec_3_30 = CovarianceSpectrum(events, freq_interval=[3, 30],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="frac")
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [01:00<00:00, 1.50s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:59<00:00, 1.50s/it]
[26]:
plt.figure()
plt.errorbar(energies, covspec_01_1.spectrum,
xerr=energies_err, yerr=covspec_01_1.spectrum_error, fmt='o', label="0.1-1 Hz")
plt.errorbar(energies, covspec_3_30.spectrum,
xerr=energies_err, yerr=covspec_3_30.spectrum_error, fmt='o', label="3-30 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Fractional Covariance");
This should largely be equivalent to the RMS spectrum
[27]:
rmsspec_01_1 = RmsSpectrum(events, freq_interval=[0.1, 1],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="frac")
rmsspec_3_30 = RmsSpectrum(events, freq_interval=[3, 30],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="frac")
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:13<00:00, 3.03it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:13<00:00, 2.96it/s]
[28]:
plt.figure()
plt.errorbar(energies, covspec_3_30.spectrum,
xerr=energies_err, yerr=covspec_3_30.spectrum_error, fmt='o', label="Cov. 3-30 Hz", alpha=0.5)
plt.errorbar(energies, covspec_01_1.spectrum,
xerr=energies_err, yerr=covspec_01_1.spectrum_error, fmt='o', label="Cov. 0.1-1 Hz", alpha=0.5)
plt.errorbar(energies, rmsspec_3_30.spectrum,
xerr=energies_err, yerr=rmsspec_3_30.spectrum_error, fmt='o', label="RMS 3-30 Hz")
plt.errorbar(energies, rmsspec_01_1.spectrum,
xerr=energies_err, yerr=rmsspec_01_1.spectrum_error, fmt='o', label="RMS 0.1-1 Hz")
plt.legend()
plt.semilogx()
plt.xlabel("Energy (keV)")
plt.ylabel("Fractional RMS");
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.
[29]:
from stingray.varenergyspectrum import LagSpectrum
covspec_3_30 = CovarianceSpectrum(events, freq_interval=[3, 30],
segment_size=segment_size, bin_time=bin_time,
energy_spec=energy_spec, norm="frac")
100%|█████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:59<00:00, 1.49s/it]
[30]:
def variable_for_value(value):
for n,v in globals().items():
if id(v) == id(value):
return n
return None
for func in [lagspec_3_30, lagspec_01_1, covspec_01_1, covspec_3_30]:
name = variable_for_value(func)
func.write(name + ".csv", fmt="ascii")
[ ]: