Intro to Waveforms

This tutorial demonstrates how to use PyART’s waveform class

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import seaborn as sns

# logging configuration, default INFO level
from PyART.logging_config import setup_logging
setup_logging()

# Set up Seaborn aesthetics
sns.set_context('talk') 
sns.set_theme(font_scale=1.2)
sns.set_style('ticks')

# Update matplotlib rcParams
rcParams.update(
    {
        'text.usetex': False,
        'font.family': 'stixgeneral',
        'mathtext.fontset': 'stix',
        'axes.grid': True,
        'grid.linestyle': ':',
        'grid.color': '#bbbbbb',
        'axes.linewidth': 1,
    }
)

We now create an empty Waveform object and inspect its attributes. This class is inherited by all waveform catalogs and models, meaning that they will all have these attributes and methods.
For more information, see the Waveform class documentation.

from PyART.waveform import Waveform

# Create an empty Waveform object
waveform = Waveform()

attributes = [attr for attr in dir(waveform) if not attr.startswith('_') and not callable(getattr(waveform, attr))]
methods = [method for method in dir(waveform) if not method.startswith('_') and callable(getattr(waveform, method))]

print("Attributes of the Waveform object:")
for attr in attributes:
    print(f"- {attr}")

print("\nMethods of the Waveform object:")
for method in methods:
    print(f"- {method}")

Attributes of the Waveform object:
- domain
- dothlm
- dyn
- f
- hc
- hlm
- hp
- kind
- psi4lm
- t
- t_psi4
- u
- units

Methods of the Waveform object:
- compute_dothlm
- compute_hphc
- compute_psi4lm
- cut
- dynamics_from_hlm
- ej_from_hlm
- extract_merger_ringdown_qts
- find_max
- integrate_data
- interpolate_hlm
- phase_shift
- plot
- plot_modes
- to_SI
- to_frequency
- to_geom
- to_time

We now load a waveform from the SXS catalog and plot its modes \(h_{\ell m}\).
They are stored as a dictionary of the form: hlm[(l,m)]. Each dictionary entry has the amplitude: A and phase p of the mode, as well as the real and imag part of \(h_{\ell,m}\).

Note that, by default, the modes are not interpolated to a uniform time grid. Therefore, before plotting, we:

• Interpolate to a uniform time grid with waveform.interpolate_hlm

• Shift the time array so that the (2,2) mode peaks at \(t=0M\), using the find_max method.

from PyART.catalogs import sxs

# Utility function to check if an array is uniformly spaced
def is_uniform(arr, rtol=1e-10, atol=1e-12):
    delta = np.diff(arr)
    return np.allclose(delta, delta[0], rtol=rtol, atol=atol)

sxs_waveform = sxs.Waveform_SXS(ID='0180',
                                download=True,
                                ignore_deprecation=True,
                                downloads=["hlm", "metadata", "horizons"],
                                load=["hlm", "metadata", "horizons"]
                                )


# verify that the time array is uniform
print("Is the time array uniform?", is_uniform(sxs_waveform.u))

# interpolate to uniform time gri, dT = 1M
t_interp, hlm_interp = sxs_waveform.interpolate_hlm(1)

# overwrite time array and modes in the waveform object
sxs_waveform._t = t_interp
sxs_waveform._u = t_interp
sxs_waveform._hlm = hlm_interp

# verify that the time array is uniform
print("Is the time array uniform?", is_uniform(sxs_waveform.u))

t_max, _, _, _, mrg_idx = sxs_waveform.find_max(return_idx=True, mode=(2,2))

fig, ax = plt.subplots()
modes = [(2, 2), (2, 1), (3, 3), (4, 4), (2, -2), (2, -1), (3, -3), (4, -4)]
for mode in modes:
    ax.semilogy(sxs_waveform.u-t_max, sxs_waveform.hlm[mode]['A'], label=f'{mode}')

# legend, out of box
ax.legend(ncol=1, bbox_to_anchor=(1.05, 1), loc='upper left')
ax.axvline(0, color='k', linestyle='--', alpha=0.5)
ax.set_xlabel('t/M')
ax.set_ylabel(r'$h_{\ell m}/\nu$')
plt.show()

2026-07-29 05:43:06 The path ../dat/SXS/SXS_BBH_0180 does not exist, contains no 'Lev*'directory, or does not contain the requested order.

2026-07-29 05:43:06 Downloading the simulation from the SXS catalog.

2026-07-29 05:43:07 Setting the download (cache) directory to ../dat/SXS/

2026-07-29 05:43:07 Loading SXS simulations using latest tag 'v3.0.0', published at 2025-05-14T18:17:30Z.

/opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
2026-07-29 05:43:07 Downloading to /home/runner/work/PyART/PyART/docs/source/dat/SXS/simulations_v3.0.0.temp.json:

  0%|          | 0/5589954 [00:00<?, ?it/s]

25.3MB [00:00, 265MB/s]                    

31.1MB [00:00, 267MB/s]

2026-07-29 05:43:19 Loaded SXS simulation SXS:BBH:0180.

2026-07-29 05:43:20 Downloading to /home/runner/work/PyART/PyART/docs/source/dat/SXS/SXS:BBH:0180v2.0/Lev4:Strain_N2.h5:

  0%|          | 0/1390909 [00:00<?, ?it/s]

  5%|▍         | 64.0k/1.33M [00:00<00:04, 300kB/s]

  9%|▉         | 128k/1.33M [00:00<00:03, 322kB/s] 

 24%|██▎       | 320k/1.33M [00:00<00:01, 629kB/s]

 42%|████▏     | 576k/1.33M [00:00<00:00, 908kB/s]

 57%|█████▋    | 768k/1.33M [00:00<00:00, 947kB/s]

 71%|███████   | 960k/1.33M [00:01<00:00, 969kB/s]

 90%|████████▉ | 1.19M/1.33M [00:01<00:00, 1.09MB/s]

100%|██████████| 1.33M/1.33M [00:01<00:00, 940kB/s] 

2026-07-29 05:43:31 Saved hlm data.

2026-07-29 05:43:31 Downloading to /home/runner/work/PyART/PyART/docs/source/dat/SXS/SXS:BBH:0180v2.0/Lev4:Horizons.h5:

  0%|          | 0/3484192 [00:00<?, ?it/s]

  2%|▏         | 64.0k/3.32M [00:00<00:09, 346kB/s]

  4%|▍         | 128k/3.32M [00:00<00:07, 420kB/s] 

  9%|▉         | 320k/3.32M [00:00<00:04, 734kB/s]

 17%|█▋        | 576k/3.32M [00:00<00:02, 1.00MB/s]

 23%|██▎       | 768k/3.32M [00:00<00:02, 1.22MB/s]

 32%|███▏      | 1.06M/3.32M [00:00<00:01, 1.42MB/s]

 38%|███▊      | 1.25M/3.32M [00:01<00:01, 1.28MB/s]

 45%|████▌     | 1.50M/3.32M [00:01<00:01, 1.31MB/s]

 51%|█████     | 1.69M/3.32M [00:01<00:01, 1.44MB/s]

 56%|█████▋    | 1.88M/3.32M [00:01<00:01, 1.28MB/s]

 62%|██████▏   | 2.06M/3.32M [00:01<00:01, 1.19MB/s]

 68%|██████▊   | 2.25M/3.32M [00:02<00:00, 1.14MB/s]

 73%|███████▎  | 2.44M/3.32M [00:02<00:00, 1.11MB/s]

 83%|████████▎ | 2.75M/3.32M [00:02<00:00, 1.30MB/s]

 90%|█████████ | 3.00M/3.32M [00:02<00:00, 1.32MB/s]

 96%|█████████▌| 3.19M/3.32M [00:02<00:00, 1.23MB/s]

100%|██████████| 3.32M/3.32M [00:02<00:00, 1.19MB/s]

2026-07-29 05:43:34 Saved horizons data.

2026-07-29 05:43:34 Saved metadata.

Is the time array uniform? False
Is the time array uniform? True

We now use the modes to compute the energy and angular momentum radiated as a function of time. To do so we use the ej_from_hlm method of the Waveform class.

Note that this method requires knowledge of the initial ADM mass and angular momentum of the system, which are stored as attributes of the Waveform object (and automatically loaded, for most catalogs).

M_adm_0 = sxs_waveform.metadata['E0byM']
J_adm_0 = sxs_waveform.metadata['Jz0']
m1 = sxs_waveform.metadata['m1']
m2 = sxs_waveform.metadata['m2']

eb, e, jorb = sxs_waveform.ej_from_hlm(M_adm_0, J_adm_0, m1, m2, modes=modes)

# identify merger
e_mrg = eb[mrg_idx]
j_mrg = jorb[mrg_idx]

fig, ax = plt.subplots()
ax.plot(jorb, eb)
ax.set_xlabel(r'$J/M^2$')
ax.set_ylabel(r'$E_b/M$')
ax.scatter(j_mrg, e_mrg, color='k', label='Merger')
ax.set_title('Binding energy vs Angular Momentum')
plt.show()

We will now look at the dynamics. This is not available for all catalogs/models. In the case of SXS waveforms, we do load Horizon data when we download the waveform and can access it.

The horizon data is stored in waveform.dyn. The time is stored once, in dyn['t'], and the vectors are stored without it, with shape (N, 3):

key	meaning
t	time
m1, m2	Christodoulou masses
chi1, chi2	dimensionless spin vectors \(\vec{\chi}_i = \vec{S}_i/m_i^2\)
S1_mag, S2_mag	dimensionful spin magnitudes \(|\vec{S}_i|\)
x1, x2	coordinate positions in the inertial frame

If the binary merges, the common horizon is loaded too. It only forms at merger, so it lives on its own, shorter time array t_remnant, alongside m_remnant, chi_remnant, S_remnant_mag and x_remnant.

# plot the mass / spin evolution
fig, axs = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
axs[0].semilogy(sxs_waveform.dyn['t'], abs(sxs_waveform.dyn['m1']-m1), label=r'$m_1$')
axs[0].semilogy(sxs_waveform.dyn['t'], abs(sxs_waveform.dyn['m2']-m2), label=r'$m_2$', linestyle='--')
axs[1].semilogy(sxs_waveform.dyn['t'], sxs_waveform.dyn['S1_mag'], label=r'$S_1$')
axs[1].semilogy(sxs_waveform.dyn['t'], sxs_waveform.dyn['S2_mag'], label=r'$S_2$', linestyle='--')
axs[0].set_ylabel(r'$\Delta m_i$')
axs[1].set_ylabel(r'$|S_i|$')
axs[1].set_xlabel('t/M')
axs[0].legend()
axs[1].legend()
plt.show()

# the dimensionless spins are vectors: chi_i = S_i / m_i^2
fig, ax = plt.subplots(figsize=(12, 4))
for i, comp in enumerate(['x', 'y', 'z']):
    ax.plot(sxs_waveform.dyn['t'], sxs_waveform.dyn['chi1'][:, i], label=rf'$\chi_1^{comp}$')
ax.set_xlabel('t/M')
ax.set_ylabel(r'$\chi_1$')
ax.set_title('Dimensionless spin components of BH 1')
ax.legend()
plt.show()

# the common horizon forms at merger, so it has its own (shorter) time array
if 'm_remnant' in sxs_waveform.dyn:
    chi_rem = np.linalg.norm(sxs_waveform.dyn['chi_remnant'][-1])
    print(f"remnant mass: {sxs_waveform.dyn['m_remnant'][-1]:.6f} "
          f"(metadata: {sxs_waveform.metadata['Mf']:.6f})")
    print(f"remnant spin: {chi_rem:.6f} "
          f"(metadata: {sxs_waveform.metadata['af']:.6f})")

# plot the trajectories
xA = sxs_waveform.dyn['x1']
xB = sxs_waveform.dyn['x2']
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(xA[:, 0], xA[:, 1], label='BH 1')
ax.plot(xB[:, 0], xB[:, 1], label='BH 2')
ax.set_xlabel('x/M')
ax.set_ylabel('y/M')
ax.set_title('Black Hole Trajectories')
ax.legend()
plt.show()

remnant mass: 0.951615 (metadata: 0.951615)
remnant spin: 0.686432 (metadata: 0.686430)

The waveform class itself also has a method to plot modes and/or polarizations, to get quick visualizations:

help(sxs_waveform.plot)

Help on method plot in module PyART.waveform:

plot(quantity, show=False, ax=None, labels_on=True, **kwargs) method of PyART.catalogs.sxs.Waveform_SXS instance
    Plot the specified quantity for rapid visualization.
    
    Parameters
    ----------
    quantity: str
        quantity to plot, e.g. 'hp', 'hc', 'hlm', 'dothlm', 'psi4lm', 'dyn'
    
    show: bool
        if True, show the plot immediately, otherwise return the axes object for further customization.
    
    ax: matplotlib.axes.Axes
        If provided, the plot will be drawn on this axes object. If not provided, a new figure and axes will be created.
    
    labels_on: bool or str
        If True, set both x and y labels. If 'x', set only x label. If 'y', set only y label. If False, do not set any labels.
    
    kwargs: dict
        additional keyword arguments to pass to the plotting function.
        The following special keyword arguments are also recognized:
        - mode: tuple
            For 'hlm', 'dothlm', or 'psi4lm', specify the (l,m) mode to plot. Default is (2,2).
        - dyn_quantities: list
            For 'dyn', specify the dynamics quantities to plot. Can be one or two quantities. If two are provided,
            the first will be plotted against the second. If one is provided, it will be plotted against time.
            Default is to plot r vs t
    
    Returns
    -------
    If show is False, returns the matplotlib.axes.Axes object containing the plot. Otherwise, nothing is returned and the plot is displayed.

# we compute the polarizations & r for default plots
sxs_waveform.compute_hphc(phi=0., i=np.pi/3)
sxs_waveform.dyn['r'] = np.sqrt((xA[:, 0]-xB[:, 0])**2 + (xA[:, 1]-xB[:, 1])**2 + (xA[:, 2]-xB[:, 2])**2)

# plot pol, modes and r(t)
fig, ax = plt.subplots(4, 1, figsize=(6, 10), sharex=True)
for i, quantity in enumerate(['hp', 'hc', 'hlm', 'dyn']):
    sxs_waveform.plot(quantity, show=False, ax=ax[i])