This Jupyter notebook requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw. To install it, simply run sage -pip install kerrgeodesic_gw
version()
'SageMath version 9.4, Release Date: 2021-08-22'
We are interested in the gravitational wave emission from a particle of mass $\mu$ on the prograde ISCO of a Kerr black hole of mass $M$ and angular momentum parameter $a = 0.90 M$. The radius $r_0$ (in terms of Boyer-Lindquist coordinates) of the prograde ISCO is obtained as follows:
from kerrgeodesic_gw import KerrBH
a = 0.90
r0 = KerrBH(a).isco_radius()
r0
2.32088304176189
The function h_plus_particle
gives the rescaled gravitational wave in the $+$ mode, $(r/\mu)\, h_+$, at the retarded time $u = t - r_*$:
from kerrgeodesic_gw import h_plus_particle
theta = pi/4
phi = 0
u = 0
h_plus_particle(a, r0, u, theta, phi)
0.5582496798737479
while the function h_cross_particle
gives $(r/\mu)\, h_\times$:
from kerrgeodesic_gw import h_cross_particle
h_cross_particle(a, r0, u, theta, phi)
-0.05710956283550757
Plot of the waveform on a time interval of twice the orbital period:
from kerrgeodesic_gw import plot_h_particle
umin = 0
umax = 2/KerrBH(a).orbital_frequency(r0)
plot_h_particle(a, r0, theta, phi, umin, umax, plot_points=400)
The corresponding spectrum:
from kerrgeodesic_gw import plot_spectrum_particle
plot_spectrum_particle(a, r0, theta, legend_label=r'$h_+$') \
+ plot_spectrum_particle(a, r0, theta, mode='x', color='red',
legend_label=r'$h_\times$', offset=0.08)
Same thing for $\theta=\frac{\pi}{2}$ (edge-on view):
theta = pi/2
plot_h_particle(a, r0, theta, phi, umin, umax, plot_points=400)
plot_spectrum_particle(a, r0, theta, legend_label=r'$h_+$') \
+ plot_spectrum_particle(a, r0, theta, mode='x', color='red',
legend_label=r'$h_\times$', offset=0.08)
Same thing for $\theta=0$ (face-on view):
theta = 0
plot_h_particle(a, r0, theta, phi, umin, umax, plot_points=400)
plot_spectrum_particle(a, r0, theta, legend_label=r'$h_+$') \
+ plot_spectrum_particle(a, r0, theta, mode='x', color='red',
legend_label=r'$h_\times$', offset=0.08)
The factor $\mu/r$, with $\mu = 1 M_\odot$ and $r$ equal to the distance to Sgr A*:
from kerrgeodesic_gw import astro_data
mu_ov_r = astro_data.solar_mass_m / astro_data.SgrA_distance_m
mu_ov_r
5.893577116499728e-18
The mass of Sgr A* in seconds:
BH_time_scale = astro_data.SgrA_mass_s
BH_time_scale
20.19522512635793
The strain spectral sensitivity of LISA detector:
from kerrgeodesic_gw import lisa_detector
hn = lisa_detector.strain_sensitivity
plot(hn, (1e-5, 1), plot_points=2000, ymin=1e-20, ymax=1e-14,
axes_labels=[r"$f\ [\mathrm{Hz}]$",
r"$S(f)^{1/2} \ \left[\mathrm{Hz}^{-1/2}\right]$"],
scale='loglog', gridlines='minor', frame=True, axes=False)
The effective power spectral density of LISA noise:
psd = lisa_detector.power_spectral_density
The signal-to-noise ratio for one day of observations is:
from kerrgeodesic_gw import signal_to_noise_particle
theta = pi/4
t_obs = 24*3600 # 1 day in seconds
signal_to_noise_particle(a, r0, theta, psd, t_obs,
BH_time_scale, scale=mu_ov_r)
103718.23722663395
Plot of the SNR as a function of the orbital radius $r_0$:
rmin = KerrBH(a).isco_radius()
rmax = 50
fsnr = lambda r0: signal_to_noise_particle(a, r0, theta, psd, t_obs,
BH_time_scale, scale=mu_ov_r)
graph = plot(fsnr, (rmin, rmax), plot_points=50,
thickness=1.5, scale='semilogy',
axes_labels=[r"$r_0/M$",
r"SNR$\times \left(\frac{1\, M_\odot}{\mu}\right) \left(\frac{1\, {\rm d}}{T}\right) ^{1/2}$"],
gridlines='minor', frame=True, axes=False)
graph