Computing Remnant Properties from Eccentric SXS Simulations
This notebook loads eccentric 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, passing the reference eccentricity viae_ref.
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
ecc = md.get("reference_eccentricity")
try:
ecc = float(ecc)
except (TypeError, ValueError):
ecc = None
info = dict(
q=float(q),
ecc=ecc,
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 (low ecc)": "SXS:BBH:3971",
"non-spinning (high ecc)": "SXS:BBH:3972",
"precessing": "SXS:BBH:3712",
}
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}, 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 SXS:BBH:3971 (non-spinning (low ecc))...
Loading SXS simulations using latest tag 'v3.0.0', published at 2025-05-14T18:17:30Z.
q = 1.000, ecc = 0.0407
chi1 = [ 0. -0. -0.], chi2 = [-0. -0. -0.]
Time range: [460.1, 5363.5] M
Loading SXS:BBH:3972 (non-spinning (high ecc))...
q = 1.000, ecc = 0.4032
chi1 = [-0. -0. 0.], chi2 = [-0. 0. -0.]
Time range: [1139.9, 11242.7] M
Loading SXS:BBH:3712 (precessing)...
q = 1.252, ecc = 0.0729
chi1 = [-0.016 -0.331 -0.036], chi2 = [0.188 0.291 0.732]
Time range: [487.0, 3320.5] M
Loading SXS simulations using latest tag 'v3.0.0', published at 2025-05-14T18:17:30Z.
q = 1.000, ecc = 0.0407
chi1 = [ 0. -0. -0.], chi2 = [-0. -0. -0.]
Time range: [460.1, 5363.5] M
Loading SXS:BBH:3972 (non-spinning (high ecc))...
q = 1.000, ecc = 0.4032
chi1 = [-0. -0. 0.], chi2 = [-0. 0. -0.]
Time range: [1139.9, 11242.7] M
Loading SXS:BBH:3712 (precessing)...
q = 1.252, ecc = 0.0729
chi1 = [-0.016 -0.331 -0.036], chi2 = [0.188 0.291 0.732]
Time range: [487.0, 3320.5] 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, low eccentricity
[4]:
t, h_dict, info = data["non-spinning (low ecc)"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['non-spinning (low ecc)']}: q = {info['q']:.3f}, ecc = {info['ecc']:.4f}")
print(f" 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:3971: q = 1.000, ecc = 0.0407
E_initial = 0.007166, L_initial = 1.118046
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03949449 M
Peak luminosity : 0.00097531
Remnant mass : 0.95333993 M
Remnant spin (dimensionless) : 0.70637563
Remnant spin vector (x,y,z) : (0.00000000, -0.00000000, 0.70637563)
Remnant kick velocity : 0.00000000 c
Remnant kick velocity : 0.00 km/s
Remnant kick vector (x,y,z) : (-0.00000000, 0.00000000, 0.00000000) c
Remnant displacement (x,y,z) : (-0.00000009, 0.00000048, 0.00000051) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.953340 0.951621
Remnant |chi| 0.706376 0.686648
Remnant spin (z) 0.706376 0.686648
Kick velocity [c] 0.000000 0.000000
Spin vector [ 4.0653e-11 -1.4600e-10 7.0638e-01]
[-6.8953e-11 1.6464e-10 6.8665e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/76828473.py:6: 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-spinning, high eccentricity
[5]:
t, h_dict, info = data["non-spinning (high ecc)"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['non-spinning (high ecc)']}: q = {info['q']:.3f}, ecc = {info['ecc']:.4f}")
print(f" 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:3972: q = 1.000, ecc = 0.4032
E_initial = 0.004994, L_initial = 1.198603
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.000
Initial mass : 1.00000000 M
Total energy radiated : 0.04258557 M
Peak luminosity : 0.00102816
Remnant mass : 0.95242072 M
Remnant spin (dimensionless) : 0.70022783
Remnant spin vector (x,y,z) : (0.00000000, -0.00000000, 0.70022783)
Remnant kick velocity : 0.00000001 c
Remnant kick velocity : 0.00 km/s
Remnant kick vector (x,y,z) : (0.00000001, -0.00000001, 0.00000000) c
Remnant displacement (x,y,z) : (0.00000325, -0.00000187, 0.00000004) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.952421 0.951787
Remnant |chi| 0.700228 0.688662
Remnant spin (z) 0.700228 0.688662
Kick velocity [c] 0.000000 0.000000
Spin vector [ 2.4232e-11 -5.0568e-11 7.0023e-01]
[-2.5855e-11 4.3417e-11 6.8866e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/1384759221.py:6: 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, eccentric
[6]:
t, h_dict, info = data["precessing"]
E_init, L_init = nr_initial_conditions(info)
print(f"{sims['precessing']}: q = {info['q']:.3f}, ecc = {info['ecc']:.4f}")
print(f" 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:3712: q = 1.252, ecc = 0.0729
E_initial = 0.008075, L_initial = 1.030219
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.252
Initial mass : 1.00000000 M
Total energy radiated : 0.04747640 M
Peak luminosity : 0.00114379
Remnant mass : 0.94444837 M
Remnant spin (dimensionless) : 0.77238009
Remnant spin vector (x,y,z) : (0.02686384, -0.03576397, 0.77238009)
Remnant kick velocity : 0.00075243 c
Remnant kick velocity : 225.57 km/s
Remnant kick vector (x,y,z) : (0.00056890, -0.00026982, -0.00041195) c
Remnant displacement (x,y,z) : (0.06677721, -0.03921815, -0.06489974) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.944448 0.944022
Remnant |chi| 0.773674 0.754679
Remnant spin (z) 0.772380 0.753761
Kick velocity [c] 0.000752 0.000821
Spin vector [ 0.0269 -0.0358 0.7724]
[ 0.0269 -0.0257 0.7538] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/2075057483.py:6: 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. For eccentric orbits, the reference eccentricity is passed via e_ref so that the PN expressions (SEOBNRv5EHM-based, valid for eccentric orbits and aligned spins to 3PN) account for the eccentricity.
3a. Non-spinning, low eccentricity
[7]:
t, h_dict, info = data["non-spinning (low ecc)"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]),
e_ref=info["ecc"])
print(f"{sims['non-spinning (low ecc)']}: 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)
SXS:BBH:3971: e_ref = 0.0407
PN E_initial = 0.007415, PN L_initial = 1.102606
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03949449 M
Peak luminosity : 0.00097531
Remnant mass : 0.95309092 M
Remnant spin (dimensionless) : 0.68974664
Remnant spin vector (x,y,z) : (0.00000000, -0.00000000, 0.68974664)
Remnant kick velocity : 0.00000000 c
Remnant kick velocity : 0.00 km/s
Remnant kick vector (x,y,z) : (-0.00000000, 0.00000000, 0.00000000) c
Remnant displacement (x,y,z) : (-0.00000009, 0.00000048, 0.00000051) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.953091 0.951621
Remnant |chi| 0.689747 0.686648
Remnant spin (z) 0.689747 0.686648
Kick velocity [c] 0.000000 0.000000
Spin vector [ 4.0674e-11 -1.4607e-10 6.8975e-01]
[-6.8953e-11 1.6464e-10 6.8665e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/2539476334.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"],
3b. Non-spinning, high eccentricity
[8]:
t, h_dict, info = data["non-spinning (high ecc)"]
calc = GWRemnantCalculator(time=t, h_dict=h_dict, q=info["q"],
chi1=list(info["chi1"]), chi2=list(info["chi2"]),
e_ref=info["ecc"])
print(f"{sims['non-spinning (high ecc)']}: 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)
SXS:BBH:3972: e_ref = 0.4032
PN E_initial = 0.005595, PN L_initial = 1.163696
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.000
Initial mass : 1.00000000 M
Total energy radiated : 0.04258557 M
Peak luminosity : 0.00102816
Remnant mass : 0.95181990 M
Remnant spin (dimensionless) : 0.66258257
Remnant spin vector (x,y,z) : (0.00000000, -0.00000000, 0.66258257)
Remnant kick velocity : 0.00000001 c
Remnant kick velocity : 0.00 km/s
Remnant kick vector (x,y,z) : (0.00000001, -0.00000001, 0.00000000) c
Remnant displacement (x,y,z) : (0.00000325, -0.00000187, 0.00000004) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.951820 0.951787
Remnant |chi| 0.662583 0.688662
Remnant spin (z) 0.662583 0.688662
Kick velocity [c] 0.000000 0.000000
Spin vector [ 2.4262e-11 -5.0632e-11 6.6258e-01]
[-2.5855e-11 4.3417e-11 6.8866e-01] (SXS)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/66046558.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, eccentric
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"]),
e_ref=info["ecc"])
print(f"{sims['precessing']}: 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)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23618/806739006.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:3712: e_ref = 0.0729
PN E_initial = 0.008480, PN L_initial = 1.012780
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 1.252
Initial mass : 1.00000000 M
Total energy radiated : 0.04747640 M
Peak luminosity : 0.00114379
Remnant mass : 0.94404335 M
Remnant spin (dimensionless) : 0.75347566
Remnant spin vector (x,y,z) : (0.02688690, -0.03579467, 0.75347566)
Remnant kick velocity : 0.00075275 c
Remnant kick velocity : 225.67 km/s
Remnant kick vector (x,y,z) : (0.00056915, -0.00026994, -0.00041212) c
Remnant displacement (x,y,z) : (0.06680590, -0.03923495, -0.06492752) M
==================================================
Property gw_remnant SXS metadata
---------------------------------------------------------
Remnant mass [M] 0.944043 0.944022
Remnant |chi| 0.754804 0.754679
Remnant spin (z) 0.753476 0.753761
Kick velocity [c] 0.000753 0.000821
Spin vector [ 0.0269 -0.0358 0.7535]
[ 0.0269 -0.0257 0.7538] (SXS)
[ ]: