Computing Remnant Properties from SXS Numerical Relativity Simulations
This notebook loads quasi-circular SXS simulations, computes remnant properties with gw_remnant, and compares against the remnant properties in the SXS simulation metadata.
Two approaches are demonstrated:
Section 2: Supply
E_initialandL_initialfrom the NR metadata (exact initial conditions).Section 3: Let
gw_remnantcomputeE_initialandL_initialfrom its built-in PN expressions.
Requires pip install gw_remnant[nr].
Contact: Tousif Islam [tousifislam24@gmail.com]
1. Setup
[1]:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import numpy as np
import sxs
from gw_remnant.gw_remnant_calculator import GWRemnantCalculator
lal.MSUN_SI != Msun
[2]:
LMAX = 8
def load_sxs_simulation(name):
"""Load an SXS simulation and return (time, h_dict, metadata_info)."""
sim = sxs.load(name)
md = sim.metadata
h = sim.h
# discard junk radiation before reference time
reference_time = md.get("reference_time", h.t[0])
reference_index = h.index_closest_to(reference_time)
h = h[reference_index:]
t = np.array(h.t, float)
h_dict = {}
for l in range(2, LMAX + 1):
for m in range(-l, l + 1):
h_dict[(l, m)] = np.array(h.data[:, h.index(l, m)], complex)
m1 = md.get("reference_mass1") or md.get("initial_mass1")
m2 = md.get("reference_mass2") or md.get("initial_mass2")
q = m1 / m2 if (m1 and m2) else md.get("reference_mass_ratio")
if q is not None and q < 1:
q = 1.0 / q
info = dict(
q=float(q),
M_ADM=float(md["initial_ADM_energy"]),
J_ADM=np.array(md["initial_ADM_angular_momentum"], float),
chi1=np.array(md["initial_dimensionless_spin1"], float),
chi2=np.array(md["initial_dimensionless_spin2"], float),
rem_mass=float(md["remnant_mass"]),
rem_chi=np.array(md["remnant_dimensionless_spin"], float),
kick=np.array(md.get("remnant_velocity", [0, 0, 0]), float),
)
return t, h_dict, info
def nr_initial_conditions(info):
"""Compute E_initial and L_initial from NR metadata."""
q = info["q"]
X1, X2 = q / (1 + q), 1 / (1 + q)
S1z = info["chi1"][2] * X1**2
S2z = info["chi2"][2] * X2**2
E_initial = 1.0 - info["M_ADM"]
L_initial = info["J_ADM"][2] - S1z - S2z
return E_initial, L_initial
def compare_with_nr(calc, info):
"""Print a side-by-side comparison of gw_remnant vs SXS metadata."""
nr_chi = info["rem_chi"]
nr_kick = np.linalg.norm(info["kick"])
print(f"{'Property':<25} {'gw_remnant':>15} {'SXS metadata':>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} {nr_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)} (SXS)")
[3]:
sims = {
"non-spinning": "SXS:BBH:2265",
"non-precessing": "SXS:BBH:2014",
"precessing": "SXS:BBH:2000",
}
data = {}
for label, name in sims.items():
print(f"Loading {name} ({label})...")
t, h_dict, info = load_sxs_simulation(name)
data[label] = (t, h_dict, info)
print(f" q = {info['q']:.3f}, 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 SXS:BBH:2265 (non-spinning)...
Loading SXS simulations using latest tag 'v3.0.0', published at 2025-05-14T18:17:30Z.
q = 3.000, chi1 = [ 0. -0. -0.], chi2 = [ 0. -0. -0.]
Time range: [1077.2, 30739.8] M
Loading SXS:BBH:2014 (non-precessing)...
q = 4.000, chi1 = [ 0. -0. 0.8], chi2 = [-0. 0. 0.4]
Time range: [820.8, 5231.0] M
Loading SXS:BBH:2000 (precessing)...
q = 3.999, chi1 = [-0.163 0.782 0.05 ], chi2 = [ 0.782 -0.107 -0.132]
Time range: [1218.8, 5097.4] M
Loading SXS simulations using latest tag 'v3.0.0', published at 2025-05-14T18:17:30Z.
q = 3.000, chi1 = [ 0. -0. -0.], chi2 = [ 0. -0. -0.]
Time range: [1077.2, 30739.8] M
Loading SXS:BBH:2014 (non-precessing)...
q = 4.000, chi1 = [ 0. -0. 0.8], chi2 = [-0. 0. 0.4]
Time range: [820.8, 5231.0] M
Loading SXS:BBH:2000 (precessing)...
q = 3.999, chi1 = [-0.163 0.782 0.05 ], chi2 = [ 0.782 -0.107 -0.132]
Time range: [1218.8, 5097.4] M
2. Using E_initial and L_initial from NR metadata
Here we supply the initial binding energy and orbital angular momentum directly from the SXS metadata:
E_initial = 1 - M_ADML_initial = J_ADM_z - S1_z - S2_z
This gives the most accurate remnant mass since it uses the exact ADM energy.
2a. Non-spinning quasi-circular
[4]:
t, h_dict, info = data["non-spinning"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['non-spinning']}: q = {info['q']:.3f}, E_initial = {E_init:.6f}, L_initial = {L_init:.6f}")
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]),
E_initial=E_init, L_initial=L_init)
calc.print_remnants()
compare_with_nr(calc, info)
SXS:BBH:2265: q = 3.000, E_initial = 0.003840, L_initial = 0.972388
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 3.000
Initial mass : 1.00000000 M
Total energy radiated : 0.02482047 M
Peak luminosity : 0.00052717
Remnant mass : 0.97133946 M
Remnant spin (dimensionless) : 0.54414793
Remnant spin vector (x,y,z) : (0.00000004, 0.00000016, 0.54414793)
Remnant kick velocity : 0.00056293 c
Remnant kick velocity : 168.76 km/s
Remnant kick vector (x,y,z) : (0.00031969, -0.00046335, 0.00000006) c
Remnant displacement (x,y,z) : (0.04215534, -0.08285644, 0.00002058) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.971339 0.971102
Remnant |chi| 0.544148 0.540609
Remnant spin (z) 0.544148 0.540609
Kick velocity [c] 0.000563 0.000582
Spin vector [3.8320e-08 1.5622e-07 5.4415e-01]
[ 5.3327e-08 -7.7351e-08 5.4061e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/953974870.py:5: 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"],
2b. Non-precessing quasi-circular
[5]:
t, h_dict, info = data["non-precessing"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['non-precessing']}: q = {info['q']:.3f}, E_initial = {E_init:.6f}, L_initial = {L_init:.6f}")
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]),
E_initial=E_init, L_initial=L_init)
calc.print_remnants()
compare_with_nr(calc, info)
SXS:BBH:2014: q = 4.000, E_initial = 0.005468, L_initial = 0.646792
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 4.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03939614 M
Peak luminosity : 0.00073047
Remnant mass : 0.95513612 M
Remnant spin (dimensionless) : 0.89469251
Remnant spin vector (x,y,z) : (-0.00000004, -0.00000033, 0.89469251)
Remnant kick velocity : 0.00006424 c
Remnant kick velocity : 19.26 km/s
Remnant kick vector (x,y,z) : (0.00003287, -0.00005519, 0.00000004) c
Remnant displacement (x,y,z) : (-0.00345658, -0.01176818, 0.00000834) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.955136 0.954579
Remnant |chi| 0.894693 0.881767
Remnant spin (z) 0.894693 0.881767
Kick velocity [c] 0.000064 0.000040
Spin vector [-3.9407e-08 -3.2589e-07 8.9469e-01]
[1.6820e-07 6.9549e-07 8.8177e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/1814872700.py:5: 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"],
2c. Precessing quasi-circular
[6]:
t, h_dict, info = data["precessing"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['precessing']}: q = {info['q']:.3f}, E_initial = {E_init:.6f}, L_initial = {L_init:.6f}")
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]),
E_initial=E_init, L_initial=L_init)
calc.print_remnants()
compare_with_nr(calc, info)
SXS:BBH:2000: q = 3.999, E_initial = 0.004981, L_initial = 0.687246
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 3.999
Initial mass : 1.00000000 M
Total energy radiated : 0.01918366 M
Peak luminosity : 0.00041586
Remnant mass : 0.97583504 M
Remnant spin (dimensionless) : 0.57633725
Remnant spin vector (x,y,z) : (-0.08356523, 0.42361271, 0.57633725)
Remnant kick velocity : 0.00090260 c
Remnant kick velocity : 270.59 km/s
Remnant kick vector (x,y,z) : (-0.00013259, 0.00010672, -0.00088641) c
Remnant displacement (x,y,z) : (-0.02164086, -0.00085242, -0.12443661) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.975835 0.975265
Remnant |chi| 0.720136 0.691529
Remnant spin (z) 0.576337 0.559751
Kick velocity [c] 0.000903 0.000970
Spin vector [-0.0836 0.4236 0.5763]
[-0.0868 0.3967 0.5598] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/4040783995.py:5: 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. Using native PN E_initial and L_initial
Here we let gw_remnant compute the initial energy and angular momentum from its built-in post-Newtonian expressions (SEOBNRv5EHM-based, valid for aligned spins to 3PN + non-spinning circular to 5PN). No E_initial or L_initial is passed.
3a. Non-spinning quasi-circular
[7]:
t, h_dict, info = data["non-spinning"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(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)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/2826678238.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"],
SXS:BBH:2265: PN E_initial = 0.003887, PN L_initial = 0.967556
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 3.000
Initial mass : 1.00000000 M
Total energy radiated : 0.02482047 M
Peak luminosity : 0.00052717
Remnant mass : 0.97129257 M
Remnant spin (dimensionless) : 0.53907838
Remnant spin vector (x,y,z) : (0.00000004, 0.00000016, 0.53907838)
Remnant kick velocity : 0.00056296 c
Remnant kick velocity : 168.77 km/s
Remnant kick vector (x,y,z) : (0.00031970, -0.00046337, 0.00000006) c
Remnant displacement (x,y,z) : (0.04215739, -0.08286044, 0.00002058) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.971293 0.971102
Remnant |chi| 0.539078 0.540609
Remnant spin (z) 0.539078 0.540609
Kick velocity [c] 0.000563 0.000582
Spin vector [3.8323e-08 1.5624e-07 5.3908e-01]
[ 5.3327e-08 -7.7351e-08 5.4061e-01] (SXS)
3b. Non-precessing quasi-circular
[8]:
t, h_dict, info = data["non-precessing"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]))
print(f"{sims['non-precessing']}: PN E_initial = {calc.E_initial:.6f}, PN L_initial = {calc.L_initial:.6f}")
calc.print_remnants()
compare_with_nr(calc, info)
SXS:BBH:2014: PN E_initial = 0.005803, PN L_initial = 0.632742
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 4.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03939614 M
Peak luminosity : 0.00073047
Remnant mass : 0.95480117 M
Remnant spin (dimensionless) : 0.87990845
Remnant spin vector (x,y,z) : (-0.00000004, -0.00000033, 0.87990845)
Remnant kick velocity : 0.00006426 c
Remnant kick velocity : 19.26 km/s
Remnant kick vector (x,y,z) : (0.00003288, -0.00005521, 0.00000004) c
Remnant displacement (x,y,z) : (-0.00345770, -0.01177226, 0.00000834) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.954801 0.954579
Remnant |chi| 0.879908 0.881767
Remnant spin (z) 0.879908 0.881767
Kick velocity [c] 0.000064 0.000040
Spin vector [-3.9435e-08 -3.2612e-07 8.7991e-01]
[1.6820e-07 6.9549e-07 8.8177e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/3571571940.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"],
3c. Precessing quasi-circular
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. When directly supplying NR initial conditions (Section 2), this issue does not arise.
[9]:
t, h_dict, info = data["precessing"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(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)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23613/2529988132.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"],
SXS:BBH:2000: PN E_initial = 0.005430, PN L_initial = 0.667319
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 3.999
Initial mass : 1.00000000 M
Total energy radiated : 0.01918366 M
Peak luminosity : 0.00041586
Remnant mass : 0.97538668 M
Remnant spin (dimensionless) : 0.55592193
Remnant spin vector (x,y,z) : (-0.08364207, 0.42400225, 0.55592193)
Remnant kick velocity : 0.00090302 c
Remnant kick velocity : 270.72 km/s
Remnant kick vector (x,y,z) : (-0.00013265, 0.00010677, -0.00088682) c
Remnant displacement (x,y,z) : (-0.02165076, -0.00085276, -0.12449380) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.975387 0.975265
Remnant |chi| 0.704147 0.691529
Remnant spin (z) 0.555922 0.559751
Kick velocity [c] 0.000903 0.000970
Spin vector [-0.0836 0.424 0.5559]
[-0.0868 0.3967 0.5598] (SXS)
[ ]:
[ ]: