Calculate LISA SNR and merger time due to GW of J0256+5934

J0526+5934_in_LISA.py

import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt

# LEGWORK is a python package that can calculate LISA SNRs - https://legwork.readthedocs.io/
# Citation info at: https://legwork.readthedocs.io/en/latest/cite.html
import legwork as lw

# take 100,000 random samples based on uncertainties
N = 100000

# parameters and uncertainties from Rebassa-Mansergas+2024
m_1 = np.random.normal(0.257, 0.049, N) * u.Msun
m_2 = np.random.normal(0.71, 0.09, N) * u.Msun
p_orb = np.random.normal(1232, 0.66, N) * u.s
ecc = np.zeros(N)
inclination = np.random.normal(65, 7, N) * u.deg

# parallax and position from GAIA DR3 (for obj 282679289838317184)
parallax = np.random.normal(1.18, 0.0910, N) * u.mas
dist = (1 / parallax.to(u.arcsec).value) * u.pc
position = SkyCoord(ra=np.repeat(81.5434157692, N) * u.deg, dec=np.repeat(59.5792549926, N) * u.deg)

# setup source with repeated values (and varied polarisations)
sources = lw.source.Source(m_1=m_1, m_2=m_2, f_orb=1/p_orb, ecc=ecc, dist=dist,
                           position=position, inclination=inclination)

# get SNR and merger time
sources.get_snr()

# plot SNR and add mean and standard deviation
fig, ax = sources.plot_source_variables("snr", show=False)
ax.set_title("J0526+5934 in LISA", fontsize=24)
ax.axvline(np.mean(sources.snr), color="lightgrey", linestyle="--")
ax.axvspan(np.mean(sources.snr) - np.std(sources.snr), np.mean(sources.snr) + np.std(sources.snr),
           color="lightgrey", alpha=0.3)

# annotate plot with mean and standard deviation of SNR in top right corner
ax.text(0.95, 0.95, f"SNR = {np.mean(sources.snr):.2f} $\\pm$ {np.std(sources.snr):.2f}",
        va="top", ha="right", transform=ax.transAxes, fontsize=16)

# annotate merger time in the same way
t_merge = sources.get_merger_time().to(u.Myr)
ax.text(0.95, 0.85, f"$t_{{\\rm merge, GW}}$ = {np.mean(t_merge.value):.2f} $\\pm$ {np.std(t_merge):.2f}",
        va="top", ha="right", transform=ax.transAxes, fontsize=16)

plt.show()

print(f"J0526+5934 SNR in LISA is {np.mean(sources.snr):.2f} +/- {np.std(sources.snr):.2f}")
print(f"J0526+5934 merger time due to GW is {np.mean(t_merge):.2f} +/- {np.std(t_merge):.2f}")

# >>> J0526+5934 SNR in LISA is 13.79 +/- 3.57
# >>> J0526+5934 merger time due to GW is 3.21 Myr +/- 0.72 Myr