Computing Remnant Properties from pySEOBNR Waveforms
This notebook generates waveforms using the pySEOBNR effective-one-body models (SEOBNRv5HM, SEOBNRv5PHM, SEOBNRv5EHM), computes remnant properties with gw_remnant, and compares against NR surrogate remnant fits.
Requires pip install gw_remnant[eob].
Contact: Tousif Islam [tousifislam24@gmail.com]
1. Setup
[1]:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import numpy as np
import lal
import surfinBH
import gwtools
from pyseobnr import GenerateWaveform
from gw_remnant.gw_remnant_calculator import GWRemnantCalculator
lal.MSUN_SI != Msun
[2]:
fit_3dq8 = surfinBH.LoadFits('NRSur3dq8Remnant')
fit_7dq4 = surfinBH.LoadFits('NRSur7dq4Remnant')
Loaded NRSur3dq8Remnant fit.
Loaded NRSur7dq4Remnant fit.
[3]:
M_TOTAL = 70.0 # total mass in Msun (arbitrary, cancels in geometric units)
DISTANCE = 100.0 # distance in Mpc (arbitrary, cancels in geometric units)
def generate_seobnr(approximant, q, chi1, chi2, e_ref=0.0):
"""Generate a pySEOBNR waveform and return in geometric units (M=1)."""
m1 = M_TOTAL * q / (1 + q)
m2 = M_TOTAL / (1 + q)
params = {
"mass1": m1, "mass2": m2,
"spin1x": chi1[0], "spin1y": chi1[1], "spin1z": chi1[2],
"spin2x": chi2[0], "spin2y": chi2[1], "spin2z": chi2[2],
"approximant": approximant,
"f22_start": 20.0,
"distance": DISTANCE,
"deltaT": 1.0 / 4096.0,
}
if e_ref > 0:
params["eccentricity"] = e_ref
gen = GenerateWaveform(params)
t_si, hlm_si = gen.generate_td_modes()
# Convert SI -> geometric units (M=1)
Mtot_s = M_TOTAL * lal.MTSUN_SI
D_m = DISTANCE * 1e6 * lal.PC_SI
Mtot_m = M_TOTAL * lal.MRSUN_SI
t = t_si / Mtot_s
h_dict = {(l, m): -h * D_m / Mtot_m for (l, m), h in hlm_si.items()}
print(f"{approximant} time grid: [{t[0]:.1f}, {t[-1]:.1f}] M")
return t, h_dict
def compare_remnants(calc, fit, q, chi1, chi2, omega0=None, fit_name="fit"):
"""Print a side-by-side comparison of gw_remnant vs a surfinBH remnant fit."""
if omega0 is not None:
mf, chif, vf, *_ = fit.all(q, chi1, chi2, omega0=omega0)
else:
mf, chif, vf, *_ = fit.all(q, chi1, chi2)
print(f"{'Property':<25} {'gw_remnant':>15} {fit_name:>20}")
print("-" * 62)
print(f"{'Remnant mass [M]':<25} {calc.remnant_mass:>15.6f} {mf:>20.6f}")
print(f"{'Remnant spin (z)':<25} {calc.remnant_spin:>15.6f} {chif[2]:>20.6f}")
print(f"{'Remnant |chi|':<25} {np.linalg.norm(calc.remnant_spin_vector):>15.6f} {np.linalg.norm(chif):>20.6f}")
print(f"{'Kick velocity [c]':<25} {calc.remnant_kick:>15.6f} {np.linalg.norm(vf):>20.6f}")
print(f"{'Spin vector':<25} {np.array2string(calc.remnant_spin_vector, precision=4)}")
print(f"{'':25} {np.array2string(chif, precision=4)} ({fit_name})")
2. Non-spinning (SEOBNRv5HM, q = 3)
[4]:
q = 3.0
chi1 = [0, 0, 0]
chi2 = [0, 0, 0]
times, h = generate_seobnr('SEOBNRv5HM', q, chi1, chi2)
SEOBNRv5HM time grid: [-2329.5, 347.8] M
[5]:
calc_ns = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
calc_ns.print_remnants()
compare_remnants(calc_ns, fit_3dq8, q, chi1, chi2, fit_name='NRSur3dq8Remnant')
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 3.000
Initial mass : 1.00000000 M
Total energy radiated : 0.02140288 M
Peak luminosity : 0.00050797
Remnant mass : 0.97180206 M
Remnant spin (dimensionless) : 0.53837510
Remnant spin vector (x,y,z) : (-0.00000000, 0.00000000, 0.53837510)
Remnant kick velocity : 0.00050084 c
Remnant kick velocity : 150.15 km/s
Remnant kick vector (x,y,z) : (0.00046378, 0.00018906, 0.00000000) c
Remnant displacement (x,y,z) : (0.16253014, 0.05912293, 0.00000000) M
==================================================
Property gw_remnant NRSur3dq8Remnant
--------------------------------------------------------------
Remnant mass [M] 0.971802 0.971233
Remnant spin (z) 0.538375 0.540600
Remnant |chi| 0.538375 0.540600
Kick velocity [c] 0.000501 0.000513
Spin vector [-7.4205e-20 6.1728e-21 5.3838e-01]
[0. 0. 0.5406] (NRSur3dq8Remnant)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23625/3118498222.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_ns = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
3. Aligned-spin (SEOBNRv5HM, q = 2, chi1z = 0.5, chi2z = -0.3)
[6]:
q = 2.0
chi1 = [0, 0, 0.5]
chi2 = [0, 0, -0.3]
times, h = generate_seobnr('SEOBNRv5HM', q, chi1, chi2)
SEOBNRv5HM time grid: [-2257.0, 372.2] M
[7]:
calc_al = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
calc_al.print_remnants()
compare_remnants(calc_al, fit_3dq8, q, chi1, chi2, fit_name='NRSur3dq8Remnant')
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03913652 M
Peak luminosity : 0.00092268
Remnant mass : 0.95265535 M
Remnant spin (dimensionless) : 0.75606931
Remnant spin vector (x,y,z) : (-0.00000000, -0.00000000, 0.75606931)
Remnant kick velocity : 0.00019383 c
Remnant kick velocity : 58.11 km/s
Remnant kick vector (x,y,z) : (0.00019118, 0.00003196, -0.00000000) c
Remnant displacement (x,y,z) : (0.07562675, 0.00941414, -0.00000000) M
==================================================
Property gw_remnant NRSur3dq8Remnant
--------------------------------------------------------------
Remnant mass [M] 0.952655 0.951674
Remnant spin (z) 0.756069 0.762036
Remnant |chi| 0.756069 0.762036
Kick velocity [c] 0.000194 0.000277
Spin vector [-1.8888e-20 -5.4327e-21 7.5607e-01]
[0. 0. 0.762] (NRSur3dq8Remnant)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23625/604263903.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_al = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
4. Precessing (SEOBNRv5PHM, 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.
[8]:
q = 2.0
chi1 = [0.4, 0.0, 0.3]
chi2 = [0.0, 0.3, -0.2]
times, h = generate_seobnr('SEOBNRv5PHM', q, chi1, chi2)
SEOBNRv5PHM time grid: [-2163.8, 365.5] M
[9]:
calc_pr = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
calc_pr.print_remnants()
omega0 = abs(np.gradient(0.5 * gwtools.phase(h[(2, 2)])) / np.gradient(times))[0]
compare_remnants(calc_pr, fit_7dq4, q, chi1, chi2, omega0=omega0,
fit_name='NRSur7dq4Remnant')
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03532244 M
Peak luminosity : 0.00085427
Remnant mass : 0.95652174 M
Remnant spin (dimensionless) : 0.70962809
Remnant spin vector (x,y,z) : (0.14346747, 0.00813686, 0.70962809)
Remnant kick velocity : 0.00031231 c
Remnant kick velocity : 93.63 km/s
Remnant kick vector (x,y,z) : (-0.00024042, -0.00007480, -0.00018478) c
Remnant displacement (x,y,z) : (-0.08045489, -0.02178497, -0.06562347) M
==================================================
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23625/469808673.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_pr = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2)
Loaded NRSur7dq4 model
Property gw_remnant NRSur7dq4Remnant
--------------------------------------------------------------
Remnant mass [M] 0.956522 0.954669
Remnant spin (z) 0.709628 0.709370
Remnant |chi| 0.724031 0.723153
Kick velocity [c] 0.000312 0.002488
Spin vector [0.1435 0.0081 0.7096]
[0.1382 0.0254 0.7094] (NRSur7dq4Remnant)
Loaded NRSur7dq4 model
Property gw_remnant NRSur7dq4Remnant
--------------------------------------------------------------
Remnant mass [M] 0.956522 0.954669
Remnant spin (z) 0.709628 0.709370
Remnant |chi| 0.724031 0.723153
Kick velocity [c] 0.000312 0.002488
Spin vector [0.1435 0.0081 0.7096]
[0.1382 0.0254 0.7094] (NRSur7dq4Remnant)
5. Eccentric (SEOBNRv5EHM, q = 2, e = 0.1)
[10]:
q = 2.0
chi1 = [0, 0, 0]
chi2 = [0, 0, 0]
e_ref = 0.1
times, h = generate_seobnr('SEOBNRv5EHM', q, chi1, chi2, e_ref=e_ref)
SEOBNRv5EHM time grid: [-1958.1, 352.4] M
[11]:
calc_ecc = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2,
e_ref=e_ref)
calc_ecc.print_remnants()
compare_remnants(calc_ecc, fit_3dq8, q, chi1, chi2, fit_name='NRSur3dq8Remnant')
==================================================
Remnant Properties Summary
==================================================
Mass ratio : 2.000
Initial mass : 1.00000000 M
Total energy radiated : 0.03039380 M
Peak luminosity : 0.00075077
Remnant mass : 0.96167657 M
Remnant spin (dimensionless) : 0.62364773
Remnant spin vector (x,y,z) : (-0.00000000, -0.00000000, 0.62364773)
Remnant kick velocity : 0.00047126 c
Remnant kick velocity : 141.28 km/s
Remnant kick vector (x,y,z) : (-0.00036052, 0.00030349, -0.00000000) c
Remnant displacement (x,y,z) : (-0.10889326, 0.10845557, -0.00000000) M
==================================================
Property gw_remnant NRSur3dq8Remnant
--------------------------------------------------------------
Remnant mass [M] 0.961677 0.961195
Remnant spin (z) 0.623648 0.623424
Remnant |chi| 0.623648 0.623424
Kick velocity [c] 0.000471 0.000445
Spin vector [-1.0200e-20 -4.5642e-22 6.2365e-01]
[0. 0. 0.6234] (NRSur3dq8Remnant)
/var/folders/9k/xlxjfz4d2cz1063spty94p0r0000gq/T/ipykernel_23625/495517983.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_ecc = GWRemnantCalculator(time=times, h_dict=h, q=q, chi1=chi1, chi2=chi2,
[ ]: