gwEccEvNS : a simple model for eccentricity evolution in time domain for non-spinning binaries

Presented by Islam and Venumadhav (https://arxiv.org/abs/2502.02739)

import sys
# install gwModels in editable mode from the parent directory
!{sys.executable} -m pip install -e ../ --no-deps --quiet

import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import gwModels
gwModels.utils.set_rcparams()

lal.MSUN_SI != Msun

Example evolution

# load eccentric NR data from SXS:BBH:1355 and circular NR data from SXS:BBH:0180
wfs = np.load('../gwModels/data/SXSBBH1355.npy', allow_pickle=True)[()]

t_ecc = wfs['t_ecc']
hecc_dict = wfs['hecc_dict']
t_cir = wfs['t_cir']
hcir_dict = wfs['hcir_dict']

# compute eccentricity from waveform data
obj = gwModels.ecc_measures.ComputeEccentricity(t_ecc = t_ecc,
                                                h_ecc_dict = {'h_l2m2': hecc_dict['h_l2m2']},
                                                t_cir = t_cir,
                                                h_cir_dict = {'h_l2m2': hcir_dict['h_l2m2']},
                                                q=1,
                                                distance_btw_peaks=2000, # choice for the distance between peak
                                                fit_funcs_orders=['3PN_m1over8', '3PN_m1over8'], # default choices for the fit orders using PN+pseudo PN expressions
                                                ecc_prefactor=2/3) # default choice for the Newtonian limit

# calculate eccentricity from gwEccEvNS model using q
et_model = gwModels.dynamics.gwEccEvNS_model(t=obj.time_xi,
                                             q=1,
                                             e0=obj.ecc_xi[0])

... gwModels eccentricity at t_ref=-2638.56 : 0.07996

plt.figure(figsize=(6,4))
plt.plot(obj.time_xi, obj.ecc_xi, c='C0', lw=3, label='from gwModels')
plt.plot(obj.time_xi, et_model, c='k', ls='--', label='$\\rm gwEccEvNS$')
plt.xlabel('$t/M$')
plt.ylabel('$\\rm Eccentricity$')
plt.ylim(0.00,0.08)
plt.xlim(xmin=-2500,xmax=50)
plt.legend(frameon=False)
plt.tight_layout()
plt.show()