Computing Remnant Properties with gw_remnant
Using a pre-computed waveform
This notebook demonstrates how to compute remnant black hole properties from a gravitational waveform stored in a .npy file, compare with the NRSur3dq8Remnant surrogate fit, and visualize the results.
Contact: Tousif Islam [tousifislam24@gmail.com]
1. Load gw_remnant
[1]:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import numpy as np
import gw_remnant
from gw_remnant.gw_remnant_calculator import GWRemnantCalculator
lal.MSUN_SI != Msun
[2]:
help(GWRemnantCalculator)
Help on class GWRemnantCalculator in module gw_remnant.gw_remnant_calculator:
class GWRemnantCalculator(gw_remnant.gw_utils.gw_plotter.GWPlotter, gw_remnant.remnant_calculators.peak_luminosity_calculator.PeakLuminosityCalculator, gw_remnant.remnant_calculators.remnant_spin_calculator.AngularMomentumCalculator, gw_remnant.remnant_calculators.trajectory_calculator.TrajectoryCalculator, gw_remnant.remnant_calculators.kick_velocity_calculator.LinearMomentumCalculator, gw_remnant.remnant_calculators.remnant_mass_calculator.RemnantMassCalculator, gw_remnant.remnant_calculators.initial_energy_momenta.InitialEnergyMomenta)
| GWRemnantCalculator(time: 'np.ndarray', h_dict: 'dict[tuple[int, int], np.ndarray]', q: 'float', chi1: 'np.ndarray | list[float] | None' = None, chi2: 'np.ndarray | list[float] | None' = None, e_ref: 'float | None' = None, E_initial: 'float | None' = None, L_initial: 'float | None' = None, M_initial: 'float' = 1, use_filter: 'bool' = False) -> 'None'
|
| Calculator for remnant properties of binary black hole mergers.
|
| This class computes remnant properties and time evolution of physical quantities
| from gravitational waveform data. It combines functionality from multiple
| calculator classes to provide comprehensive analysis of binary black hole mergers.
|
| Computed quantities include:
| - Energy flux and radiated energy
| - Linear momentum and kick velocity evolution
| - Angular momentum and spin evolution
| - Final remnant mass, spin, and kick velocity
| - Peak luminosity
|
| All calculations are performed in geometric units where G=c=1, with masses
| in units of total mass M and time in units of M.
|
| Args:
| time (np.ndarray): Array of time values in geometric units (M)
| h_dict (dict): Dictionary of complex waveform modes with (l,m) tuple keys,
| e.g., {(2,2): h_22(t), (3,3): h_33(t), ...}
| q (float): Mass ratio q = m1/m2, where m1 >= m2
| chi1 (list or np.ndarray): Spin vector [sx, sy, sz] for primary black
| hole at the start of the waveform, in dimensionless units. Default is None
| chi2 (list or np.ndarray): Spin vector [sx, sy, sz] for secondary black
| hole at the start of the waveform, in dimensionless units. Default is None
| e_ref (float): Eccentricity at the reference time. User must provide
| accurate value; code does not validate. Default is None
| E_initial (float): Initial energy of the binary in units of total mass M.
| If None, computed using PN expressions. Set to 0 to inspect energy changes
| relative to reference. Set to known value (e.g., from NR simulation) for
| absolute tracking. Default is None
| L_initial (float): Initial angular momentum of the binary in units of M^2.
| If None, computed using PN expressions. Set to 0 to inspect angular momentum
| changes relative to reference. Set to known value (e.g., from NR simulation)
| for absolute tracking. Default is None
| M_initial (float): Initial total mass in units of M. Default is 1
| use_filter (bool): Whether to apply filtering to computed quantities.
| Default is False
|
| Attributes:
| h_dot (dict): Time derivative of waveform modes
| E_dot (np.ndarray): Energy flux as a function of time
| E_rad (float): Total radiated energy
| Eoft (np.ndarray): Cumulative radiated energy as a function of time
| Moft (np.ndarray): Binary mass as a function of time
| remnant_mass (float): Final remnant black hole mass
| P_dot (np.ndarray): Time derivative of linear momentum vector [3 x N_times]
| Poft (np.ndarray): Linear momentum vector as a function of time [3 x N_times]
| voft (np.ndarray): Center of mass velocity vector as a function of time [3 x N_times]
| kickoft (np.ndarray): Kick velocity magnitude as a function of time
| remnant_kick (float): Final kick velocity magnitude
| remnant_kick_vector (np.ndarray): Final kick velocity vector (3,) in units of c
| J_dot (np.ndarray): Time derivative of angular momentum vector [3 x N_times]
| Joft (np.ndarray): Angular momentum vector as a function of time [3 x N_times]
| spinoft (np.ndarray): Dimensionless spin z-component as a function of time
| remnant_spin (float): Final remnant dimensionless spin z-component
| spin_vector_oft (np.ndarray): Dimensionless spin vector [3 x N_times] (x, y, z)
| remnant_spin_vector (np.ndarray): Final remnant dimensionless spin vector (3,)
| L_peak (float): Peak luminosity
| peak_kick (float): Peak kick velocity
| xoft (np.ndarray): Center-of-mass displacement vector [N_times x 3] in units of M
| remnant_displacement (np.ndarray): Final center-of-mass displacement vector (3,)
|
| Inherits From:
| GWPlotter: Plotting utilities for visualizing results
| PeakLuminosityCalculator: Peak luminosity calculations
| AngularMomentumCalculator: Angular momentum evolution
| TrajectoryCalculator: Center-of-mass trajectory
| LinearMomentumCalculator: Linear momentum and kick velocity
| RemnantMassCalculator: Mass and energy calculations
| InitialEnergyMomenta: Initial condition calculations
|
| Example:
| >>> import numpy as np
| >>> time = np.arange(-1000, 100, 0.1)
| >>> h_dict = {(2,2): h_22_data, (3,3): h_33_data}
| >>> calc = GWRemnantCalculator(time, h_dict, q=2.0,
| ... chi1=[0, 0, 0.5])
| >>> calc.print_remnants()
| >>> calc.plot_mass_energy()
|
| Method resolution order:
| GWRemnantCalculator
| gw_remnant.gw_utils.gw_plotter.GWPlotter
| gw_remnant.remnant_calculators.peak_luminosity_calculator.PeakLuminosityCalculator
| gw_remnant.remnant_calculators.remnant_spin_calculator.AngularMomentumCalculator
| gw_remnant.remnant_calculators.trajectory_calculator.TrajectoryCalculator
| gw_remnant.remnant_calculators.kick_velocity_calculator.LinearMomentumCalculator
| gw_remnant.remnant_calculators.remnant_mass_calculator.RemnantMassCalculator
| gw_remnant.remnant_calculators.initial_energy_momenta.InitialEnergyMomenta
| builtins.object
|
| Methods defined here:
|
| __init__(self, time: 'np.ndarray', h_dict: 'dict[tuple[int, int], np.ndarray]', q: 'float', chi1: 'np.ndarray | list[float] | None' = None, chi2: 'np.ndarray | list[float] | None' = None, e_ref: 'float | None' = None, E_initial: 'float | None' = None, L_initial: 'float | None' = None, M_initial: 'float' = 1, use_filter: 'bool' = False) -> 'None'
| Initialize self. See help(type(self)) for accurate signature.
|
| get_remnant_properties(self) -> 'dict[str, float]'
| Return remnant properties as a dictionary.
|
| Returns:
| dict: Dictionary containing:
| - 'mass_ratio': Mass ratio q = m1/m2
| - 'M_initial': Initial total mass
| - 'E_rad': Total radiated energy
| - 'L_peak': Peak luminosity
| - 'remnant_mass': Final remnant mass
| - 'remnant_spin': Final dimensionless spin (z-component)
| - 'remnant_spin_x': Final dimensionless spin x-component
| - 'remnant_spin_y': Final dimensionless spin y-component
| - 'remnant_spin_z': Final dimensionless spin z-component
| - 'remnant_kick': Final kick velocity in units of c
| - 'remnant_kick_kmps': Final kick velocity in km/s
| - 'remnant_kick_x/y/z': Final kick velocity vector components in units of c
| - 'peak_kick': Peak kick velocity in units of c
| - 'remnant_displacement_x/y/z': Final center-of-mass displacement
| components in units of M
|
| print_remnants(self) -> 'None'
| Print summary of remnant properties.
|
| Displays key remnant quantities including mass ratio, initial mass,
| total radiated energy, peak luminosity, and final remnant mass, spin,
| and kick velocity. All quantities are printed in geometric units.
|
| ----------------------------------------------------------------------
| Methods inherited from gw_remnant.gw_utils.gw_plotter.GWPlotter:
|
| plot_angular_momentum(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot angular momentum components as a function of time.
|
| Creates a 3-row, 2-column figure showing the x, y, and z components of the
| angular momentum vector. Each row is split into pre-merger (t <= -500M, 60% width)
| and post-merger (t > -500M, 40% width) regions for better visualization of
| different phases of the evolution.
|
| The angular momentum components are in units of total mass squared M^2.
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| plot_kick_velocity(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot the magnitude of kick velocity as a function of time.
|
| Creates a single-row, 2-column figure showing the time evolution of the recoil
| (kick) velocity magnitude imparted to the remnant black hole due to asymmetric
| gravitational wave emission. The plot is split into pre-merger (t <= -500M, 60% width)
| and post-merger (t > -500M, 40% width) regions.
|
| Kick velocity is given in units of the speed of light c.
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| plot_linear_momentum(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot linear momentum components as a function of time.
|
| Creates a 3-row, 2-column figure showing the x, y, and z components of the
| linear momentum vector. Each row is split into pre-merger (t <= -500M, 60% width)
| and post-merger (t > -500M, 40% width) regions for better visualization of
| different phases of the evolution.
|
| The momentum components are in units of total mass M.
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| plot_mass_energy(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot mass and energy evolution as a function of time.
|
| Creates a 4-row, 2-column figure showing:
| - Top row: Real part of (2,2) mode waveform derivative
| - Second row: Energy flux (time derivative of energy)
| - Third row: Total energy as a function of time
| - Bottom row: Mass as a function of time
|
| Each row is split into pre-merger (t <= -500M, 60% width) and post-merger
| (t > -500M, 40% width) regions for better visualization.
|
| All quantities are plotted in geometric units where G=c=1.
| Time is in units of total mass M.
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| plot_spin_vector(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot dimensionless spin vector components as a function of time.
|
| Creates a 3-row, 2-column figure showing the x, y, and z components of the
| dimensionless spin vector. Each row is split into pre-merger (t <= -500M, 60% width)
| and post-merger (t > -500M, 40% width) regions for better visualization of
| different phases of the evolution.
|
| The spin components are dimensionless (chi = S/M^2).
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| plot_trajectory(self, save_path: 'str | None' = None) -> 'matplotlib.figure.Figure'
| Plot the center-of-mass trajectory of the binary.
|
| Creates a two-panel figure:
| - Left: the center-of-mass path projected onto the orbital (x-y) plane,
| x(t) = \int v(t) dt, with the start and final (remnant) points marked.
| - Right: the displacement magnitude |x(t)| as a function of time, split
| into pre-merger (t <= -500M) and post-merger (t > -500M) regions.
|
| Positions are in units of total mass M (geometric units, G=c=1).
| See arXiv:1802.04276 and arXiv:1510.03386.
|
| Args:
| save_path (str or None): If provided, save figure to this path.
|
| Returns:
| matplotlib.figure.Figure: The figure object.
|
| ----------------------------------------------------------------------
| Data descriptors inherited from gw_remnant.remnant_calculators.initial_energy_momenta.InitialEnergyMomenta:
|
| __dict__
| dictionary for instance variables (if defined)
|
| __weakref__
| list of weak references to the object (if defined)
2. Load waveform data
We load a pre-computed waveform from q8_NR.npy. The file contains a dictionary with a time array t and complex waveform modes keyed by (l, m) tuples.
[3]:
data = np.load("q8_NR.npy", allow_pickle=True)
q = 8
wf = {}
for key in data.item():
if key == 't':
time = data.item()[key]
else:
wf[key] = data.item()[key]
print(f"Time range: [{time[0]:.1f}, {time[-1]:.1f}] M")
print(f"Available modes: {list(wf.keys())}")
Time range: [-5000.0, 49.9] M
Available modes: [(2, 2), (2, -2), (2, 1), (2, -1), (3, 1), (3, -1), (3, 2), (3, -2), (3, 3), (3, -3), (4, 2), (4, -2), (4, 3), (4, -3), (4, 4), (4, -4)]
3. Compute remnants using gw_remnant
[4]:
calc = GWRemnantCalculator(time=time, h_dict=wf, q=q)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23586/3950355643.py:1: UserWarning: Tips: If you are using NR waveforms, ensure that they do not contain junk radiation, as that is known to corrupt the remnant property estimation.
calc = GWRemnantCalculator(time=time, h_dict=wf, q=q)
[5]:
calc.print_remnants()
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 8.000
Initial mass : 1.00000000 M
Total energy radiated : 0.00702148 M
Peak luminosity : 0.00012866
Remnant mass : 0.98956868 M
Remnant spin (dimensionless) : 0.30699791
Remnant spin vector (x,y,z) : (0.00000000, 0.00000000, 0.30699791)
Remnant kick velocity : 0.00025734 c
Remnant kick velocity : 77.15 km/s
Remnant kick vector (x,y,z) : (0.00005844, -0.00025062, 0.00000000) c
Remnant displacement (x,y,z) : (-0.00363199, -0.01214107, 0.00000000) M
==================================================
[6]:
props = calc.get_remnant_properties()
props
[6]:
{'mass_ratio': 8,
'M_initial': 1,
'E_rad': np.float64(0.0070214795341573665),
'L_peak': np.float64(0.00012865542051061172),
'remnant_mass': np.float64(0.9895686824890562),
'remnant_spin': np.float64(0.3069979106066966),
'remnant_spin_x': np.float64(0.0),
'remnant_spin_y': np.float64(0.0),
'remnant_spin_z': np.float64(0.3069979106066966),
'remnant_kick': np.float64(0.0002573401394491285),
'remnant_kick_kmps': np.float64(77.148632947517),
'remnant_kick_x': np.float64(5.8442257495459254e-05),
'remnant_kick_y': np.float64(-0.00025061614056267673),
'remnant_kick_z': np.float64(0.0),
'peak_kick': np.float64(0.00031984873237122857),
'remnant_displacement_x': np.float64(-0.0036319857377641878),
'remnant_displacement_y': np.float64(-0.012141067720693373),
'remnant_displacement_z': np.float64(0.0)}
4. Compare remnant properties with NRSur3dq8Remnant
Requires surfinBH (install via pip install gw_remnant[surrogates]).
[7]:
import surfinBH
from gw_remnant.gw_utils import waveform_generator
fit = surfinBH.LoadFits('NRSur3dq8Remnant')
rem = waveform_generator.compute_nrsur3dq8_remnant(fit, q=q,
chi1=[0, 0, 0],
chi2=[0, 0, 0])
Loaded NRSur3dq8Remnant fit.
============================================================
Remnant predictions from NRSur3dq8Remnant
============================================================
Final mass : 0.989306 ± 0.000133 M
Final spin |chi| : 0.306742 ± 0.000039
Final spin vector : (0.000000, 0.000000, 0.306742)
Kick |v| : 0.000276 ± 0.000025 c
Kick vector : (2.735904e-04, -3.549509e-05, 0.000000e+00) c
============================================================
[8]:
print(f"{'Property':<25} {'gw_remnant':>15} {'NRSur3dq8Remnant':>20}")
print("-" * 62)
print(f"{'Remnant mass [M]':<25} {calc.remnant_mass:>15.6f} {rem['final_mass']:>20.6f}")
print(f"{'Remnant spin':<25} {calc.remnant_spin:>15.6f} {rem['final_spin']:>20.6f}")
print(f"{'Kick velocity [c]':<25} {calc.remnant_kick:>15.6f} {rem['kick_velocity']:>20.6f}")
Property gw_remnant NRSur3dq8Remnant
--------------------------------------------------------------
Remnant mass [M] 0.989569 0.989306
Remnant spin 0.306998 0.306742
Kick velocity [c] 0.000257 0.000276
5. Plot remnant quantities
[9]:
calc.plot_mass_energy();
[10]:
calc.plot_linear_momentum();
[11]:
calc.plot_angular_momentum();
[12]:
calc.plot_kick_velocity();
[13]:
calc.plot_spin_vector();
[14]:
calc.plot_trajectory();
[ ]: