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.

In [1]:
import warnings
import numpy as np
import matplotlib.pyplot as plt
import sxs

Load the data directly from the files:

In [2]:
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):

In [3]:
h_𝓔 = sxs.waveforms.memory.add_memory(h[i1:i2])

Compute the various components of the BMS strain:

In [4]:
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:

In [5]:
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$')
In [ ]: