This notebook illustrates the use of functions in the sxs
module related to gravitational-wave memory, as described in Adding Gravitational Memory to Waveform Catalogs using BMS Balance Laws by Mitman, et al. The documentation for these functions can be found here.
import warnings
import numpy as np
import matplotlib.pyplot as plt
import sxs
Load the data directly from the files:
with warnings.catch_warnings():
warnings.simplefilter("ignore") # To quiet warnings about spin weight not being set
h = sxs.load("rhOverM_Extrapolated_N5_CoM.h5")
psi2 = sxs.load("r3Psi2OverM_Extrapolated_N5_CoM.h5")
h._metadata['spin_weight'] = -2
psi2._metadata['spin_weight'] = 0
peak_time = h.max_norm_time()
i1, i2 = 0, h.n_times
ia, ib = h.index_closest_to(peak_time-220), h.index_closest_to(peak_time+170)
Add electric component of null memory to the waveform, following Eq. (17b):
h_ð = sxs.waveforms.memory.add_memory(h[i1:i2])
Compute the various components of the BMS strain:
J_m = sxs.waveforms.memory.J_m(h[i1:i2], psi2[i1:i2])
J_ð = sxs.waveforms.memory.J_E(h[i1:i2])
J_NĖ = sxs.waveforms.memory.J_Nhat(h[i1:i2], psi2[i1:i2])
J_ð = sxs.waveforms.memory.J_J(h[i1:i2])
J = J_m + J_ð + J_NĖ + J_ð
Plot the results, as in Figs. 1 and 6 of the paper:
for ell, part in [(2, np.real), (3, np.imag)]:
fig, axes = plt.subplots(2, sharex=True, figsize=(10, 12))
axes[0].plot(h.t[ia:ib]-peak_time, part(h[ia:ib, h.index(ell, 0)].ndarray), color='black')
axes[0].plot(h_ð.t[ia:ib]-peak_time, part(h_ð[ia:ib, h_ð.index(ell, 0)].ndarray), ls='dashed', color='red')
axes[0].plot(J.t[ia:ib]-peak_time, part(J[ia:ib, J.index(ell, 0)].ndarray), ls='dashdot', color='green')
axes[0].set_title(rf'{part.__name__} part of $({ell}, 0)$')
axes[0].set_ylabel(r'$R/M\ h$')
axes[0].legend([r'$h$', r'$h + J_{ð}$', r'$J$'])
axes[0].set_xlim(-220, 170)
axes[1].plot(J_m.t[ia:ib]-peak_time, part(J_m[ia:ib, J_m.index(ell, 0)].ndarray), ls='solid', color='black')
axes[1].plot(J_ð.t[ia:ib]-peak_time, part(J_ð[ia:ib, J_ð.index(ell, 0)].ndarray), ls='dotted', color='blue')
axes[1].plot(J_NĖ.t[ia:ib]-peak_time, part(J_NĖ[ia:ib, J_NĖ.index(ell, 0)].ndarray), ls='dashed', color='red')
axes[1].plot(J_ð.t[ia:ib]-peak_time, part(J_ð[ia:ib, J_ð.index(ell, 0)].ndarray), ls='dashdot', color='green')
axes[1].legend(['$J_m$', '$J_ð$', '$J_NĖ$', '$J_ð$'])
axes[1].set_xlabel(r'$(u - u_{\mathrm{peak}}) / M$')