Computing Remnant Properties from Maya Catalog Waveforms
This notebook loads numerical relativity waveforms from the Maya (Georgia Tech) catalog, computes remnant properties with gw_remnant, and compares against the NR remnant properties from the simulation.
Unlike SXS simulations, Maya waveforms do not provide ADM energy/angular momentum metadata, so we rely on the built-in PN expressions for the initial energy and angular momentum.
Requires pip install gw_remnant[nr].
Contact: Tousif Islam [tousifislam24@gmail.com]
1. Setup
[1]:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
warnings.filterwarnings("ignore", category=RuntimeWarning)
import os
import numpy as np
from mayawaves.utils.catalogutils import Catalog, Parameter as P
from mayawaves.coalescence import Coalescence
from gw_remnant.gw_remnant_calculator import GWRemnantCalculator
lal.MSUN_SI != Msun
[2]:
SIMS = {
"non-spinning": "GT0446",
"aligned-spin": "GT0588",
"precessing": "GT0460",
"eccentric": "MAYA0951",
}
DLDIR = "/tmp/maya"
os.makedirs(DLDIR, exist_ok=True)
cat = Catalog()
have = {f.replace(".h5", "") for f in os.listdir(DLDIR) if f.endswith(".h5")}
todo = [sid for sid in SIMS.values() if sid not in have]
if todo:
print(f"Downloading {todo}...")
cat.download_waveforms(todo, DLDIR, lvcnr_format=False)
print("All simulations ready.")
All simulations ready.
[3]:
R_EXT = 200.0
LMAX = 8
def load_maya_simulation(name):
"""Load a Maya simulation and return (time, h_dict, info)."""
co = Coalescence(os.path.join(DLDIR, f"{name}.h5"))
# Extract waveform modes at the largest available extraction radius
h_dict = {}
for (l, m) in co.included_modes:
if l <= LMAX:
t, hp, hc = co.strain_for_mode(l, m, extraction_radius=R_EXT)
h_dict[(l, m)] = hp - 1j * hc
# Remove junk radiation: skip first 5% of the pre-merger waveform
ipk = np.argmax(np.abs(h_dict[(2, 2)]))
junk_idx = max(int(0.05 * ipk), 1)
t = t[junk_idx:]
h_dict = {k: v[junk_idx:] for k, v in h_dict.items()}
# Get catalog parameters for eccentricity
p = cat.get_parameters_for_simulation(name)
ecc = float(p[P.ECCENTRICITY])
info = dict(
q=co.mass_ratio,
chi1=list(co.primary_compact_object.initial_dimensionless_spin),
chi2=list(co.secondary_compact_object.initial_dimensionless_spin),
ecc=ecc,
rem_mass=co.final_compact_object.final_horizon_mass,
rem_chi=np.array(co.final_compact_object.final_dimensionless_spin),
)
co.close()
return t, h_dict, info
def compare_with_nr(calc, info):
"""Print a side-by-side comparison of gw_remnant vs Maya NR remnant properties."""
nr_chi = info["rem_chi"]
print(f"{'Property':<25} {'gw_remnant':>15} {'Maya NR':>15}")
print("-" * 57)
print(f"{'Remnant mass [M]':<25} {calc.remnant_mass:>15.6f} {info['rem_mass']:>15.6f}")
print(f"{'Remnant |chi|':<25} {np.linalg.norm(calc.remnant_spin_vector):>15.6f} {np.linalg.norm(nr_chi):>15.6f}")
print(f"{'Remnant spin (z)':<25} {calc.remnant_spin:>15.6f} {nr_chi[2]:>15.6f}")
print(f"{'Kick velocity [c]':<25} {calc.remnant_kick:>15.6f}")
print(f"{'Spin vector':<25} {np.array2string(calc.remnant_spin_vector, precision=4)}")
print(f"{'':25} {np.array2string(nr_chi, precision=4)} (Maya NR)")
[4]:
data = {}
for label, name in SIMS.items():
print(f"Loading {name} ({label})...")
t, h_dict, info = load_maya_simulation(name)
data[label] = (t, h_dict, info)
print(f" q = {info['q']:.3f}, ecc = {info['ecc']:.4f}")
print(f" chi1 = {np.round(info['chi1'], 3)}, chi2 = {np.round(info['chi2'], 3)}")
print(f" Time range: [{t[0]:.1f}, {t[-1]:.1f}] M")
print()
Loading GT0446 (non-spinning)...
q = 2.000, ecc = 0.0016
chi1 = [-0. -0. 0.], chi2 = [-0. -0. 0.]
Time range: [82.7, 1999.0] M
Loading GT0588 (aligned-spin)...
q = 2.000, ecc = 0.0074
chi1 = [ 0. -0. 0.601], chi2 = [ 0. -0. 0.601]
Time range: [106.6, 2499.0] M
Loading GT0460 (precessing)...
q = 2.000, ecc = 0.0080
chi1 = [0.425 0. 0.425], chi2 = [ 0. -0. 0.601]
Time range: [75.7, 1698.9] M
Loading MAYA0951 (eccentric)...
q = 2.000, ecc = 0.0387
chi1 = [-0. -0. 0.], chi2 = [-0. -0. 0.]
Time range: [69.7, 1604.5] M
q = 2.000, ecc = 0.0016
chi1 = [-0. -0. 0.], chi2 = [-0. -0. 0.]
Time range: [82.7, 1999.0] M
Loading GT0588 (aligned-spin)...
q = 2.000, ecc = 0.0074
chi1 = [ 0. -0. 0.601], chi2 = [ 0. -0. 0.601]
Time range: [106.6, 2499.0] M
Loading GT0460 (precessing)...
q = 2.000, ecc = 0.0080
chi1 = [0.425 0. 0.425], chi2 = [ 0. -0. 0.601]
Time range: [75.7, 1698.9] M
Loading MAYA0951 (eccentric)...
q = 2.000, ecc = 0.0387
chi1 = [-0. -0. 0.], chi2 = [-0. -0. 0.]
Time range: [69.7, 1604.5] M
2. Non-spinning quasi-circular (GT0446, q = 2)
[5]:
t, h_dict, info = data["non-spinning"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=info["chi1"], chi2=info["chi2"])
print(f"{SIMS['non-spinning']}: PN E_initial = {calc.E_initial:.6f}, PN L_initial = {calc.L_initial:.6f}")
calc.print_remnants()
compare_with_nr(calc, info)
GT0446: PN E_initial = 0.007998, PN L_initial = 0.905186
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03036453 M
Peak luminosity : 0.00077459
Remnant mass : 0.96163725 M
Remnant spin (dimensionless) : 0.64705216
Remnant spin vector (x,y,z) : (-0.00000000, -0.00000000, 0.64705216)
Remnant kick velocity : 0.00049345 c
Remnant kick velocity : 147.93 km/s
Remnant kick vector (x,y,z) : (0.00042115, 0.00025714, 0.00000000) c
Remnant displacement (x,y,z) : (0.13480465, 0.08521788, 0.00000000) M
==================================================
Property gw_remnant Maya NR
---------------------------------------------------------
Remnant mass [M] 0.961637 0.961202
Remnant |chi| 0.647052 0.623272
Remnant spin (z) 0.647052 0.623272
Kick velocity [c] 0.000493
Spin vector [-2.0578e-13 -1.3328e-15 6.4705e-01]
[-8.6293e-10 1.8724e-09 6.2327e-01] (Maya NR)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23621/2038426531.py:3: 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=t, h_dict=h_dict, q=info["q"],
3. Aligned-spin (GT0588, q = 2, chi1z = 0.6, chi2z = 0.6)
[6]:
t, h_dict, info = data["aligned-spin"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=info["chi1"], chi2=info["chi2"])
print(f"{SIMS['aligned-spin']}: PN E_initial = {calc.E_initial:.6f}, PN L_initial = {calc.L_initial:.6f}")
calc.print_remnants()
compare_with_nr(calc, info)
GT0588: PN E_initial = 0.007827, PN L_initial = 0.894864
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.05143106 M
Peak luminosity : 0.00111540
Remnant mass : 0.94074223 M
Remnant spin (dimensionless) : 0.88902800
Remnant spin vector (x,y,z) : (0.00000000, -0.00000000, 0.88902800)
Remnant kick velocity : 0.00009365 c
Remnant kick velocity : 28.08 km/s
Remnant kick vector (x,y,z) : (-0.00007839, 0.00005124, 0.00000000) c
Remnant displacement (x,y,z) : (-0.06646915, -0.00213525, 0.00000000) M
==================================================
Property gw_remnant Maya NR
---------------------------------------------------------
Remnant mass [M] 0.940742 0.939592
Remnant |chi| 0.889028 0.838695
Remnant spin (z) 0.889028 0.838695
Kick velocity [c] 0.000094
Spin vector [ 9.3928e-12 -1.6034e-11 8.8903e-01]
[-6.0832e-09 2.2952e-08 8.3870e-01] (Maya NR)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23621/3298260311.py:3: 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=t, h_dict=h_dict, q=info["q"],
4. Precessing (GT0460, q = 2)
Note: The PN expressions assume the initial orbital angular momentum is aligned with the z-axis. For precessing binaries this assumption breaks down, so the remnant spin is slightly off even though the remnant mass and kick velocity remain accurate.
[7]:
t, h_dict, info = data["precessing"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=info["chi1"], chi2=info["chi2"])
print(f"{SIMS['precessing']}: PN E_initial = {calc.E_initial:.6f}, PN L_initial = {calc.L_initial:.6f}")
calc.print_remnants()
compare_with_nr(calc, info)
GT0460: PN E_initial = 0.009531, PN L_initial = 0.826227
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.04126588 M
Peak luminosity : 0.00089644
Remnant mass : 0.94920289 M
Remnant spin (dimensionless) : 0.79947742
Remnant spin vector (x,y,z) : (0.27472730, 0.00414207, 0.79947742)
Remnant kick velocity : 0.00269611 c
Remnant kick velocity : 808.27 km/s
Remnant kick vector (x,y,z) : (-0.00040168, -0.00010377, 0.00266400) c
Remnant displacement (x,y,z) : (-0.09322832, -0.02204109, 0.44337114) M
==================================================
Property gw_remnant Maya NR
---------------------------------------------------------
Remnant mass [M] 0.949203 0.945084
Remnant |chi| 0.845374 0.804466
Remnant spin (z) 0.799477 0.791574
Kick velocity [c] 0.002696
Spin vector [0.2747 0.0041 0.7995]
[ 0.1434 -0.0039 0.7916] (Maya NR)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23621/461026616.py:3: 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=t, h_dict=h_dict, q=info["q"],
5. Non-spinning eccentric (MAYA0951, q = 2, e = 0.04)
[8]:
t, h_dict, info = data["eccentric"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=info["chi1"], chi2=info["chi2"],
e_ref=info["ecc"])
print(f"{SIMS['eccentric']}: e_ref = {info['ecc']:.4f}")
print(f" PN E_initial = {calc.E_initial:.6f}, PN L_initial = {calc.L_initial:.6f}")
calc.print_remnants()
compare_with_nr(calc, info)
MAYA0951: e_ref = 0.0387
PN E_initial = 0.009052, PN L_initial = 0.861188
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03001894 M
Peak luminosity : 0.00077694
Remnant mass : 0.96092890 M
Remnant spin (dimensionless) : 0.61477410
Remnant spin vector (x,y,z) : (-0.00000000, -0.00000000, 0.61477410)
Remnant kick velocity : 0.00050057 c
Remnant kick velocity : 150.07 km/s
Remnant kick vector (x,y,z) : (0.00048055, 0.00014012, 0.00000000) c
Remnant displacement (x,y,z) : (0.09038768, 0.02593011, 0.00000000) M
==================================================
Property gw_remnant Maya NR
---------------------------------------------------------
Remnant mass [M] 0.960929 0.961318
Remnant |chi| 0.614774 0.623432
Remnant spin (z) 0.614774 0.623432
Kick velocity [c] 0.000501
Spin vector [-3.9591e-13 -2.7957e-15 6.1477e-01]
[-8.7386e-10 1.9332e-09 6.2343e-01] (Maya NR)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23621/1632622103.py:3: 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=t, h_dict=h_dict, q=info["q"],
[ ]: