Gravitational radiation from circular equatorial orbits around a Kerr black hole

This Jupyter notebook requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw. To install it, simply run sage -pip install kerrgeodesic_gw

In [1]:
version()
Out[1]:
'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:

In [2]:
from kerrgeodesic_gw import KerrBH
a = 0.90
r0 = KerrBH(a).isco_radius()
r0
Out[2]:
2.32088304176189

Waveforms

The function h_plus_particle gives the rescaled gravitational wave in the $+$ mode, $(r/\mu)\, h_+$, at the retarded time $u = t - r_*$:

In [3]:
from kerrgeodesic_gw import h_plus_particle
theta = pi/4
phi = 0
u = 0
h_plus_particle(a, r0, u, theta, phi)
Out[3]:
0.5582496798737479

while the function h_cross_particle gives $(r/\mu)\, h_\times$:

In [4]:
from kerrgeodesic_gw import h_cross_particle
h_cross_particle(a, r0, u, theta, phi)
Out[4]:
-0.05710956283550757

Plot of the waveform on a time interval of twice the orbital period:

In [5]:
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)
Out[5]:

The corresponding spectrum:

In [6]:
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)
Out[6]:

Same thing for $\theta=\frac{\pi}{2}$ (edge-on view):

In [7]:
theta = pi/2
plot_h_particle(a, r0, theta, phi, umin, umax, plot_points=400)
Out[7]:
In [8]:
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)
Out[8]:

Same thing for $\theta=0$ (face-on view):

In [9]:
theta = 0
plot_h_particle(a, r0, theta, phi, umin, umax, plot_points=400)
Out[9]:
In [10]:
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)
Out[10]:

Signal-to-noise ratio in LISA detector for emission from Sgr A*

The factor $\mu/r$, with $\mu = 1 M_\odot$ and $r$ equal to the distance to Sgr A*:

In [11]:
from kerrgeodesic_gw import astro_data
mu_ov_r = astro_data.solar_mass_m / astro_data.SgrA_distance_m
mu_ov_r
Out[11]:
5.893577116499728e-18

The mass of Sgr A* in seconds:

In [12]:
BH_time_scale = astro_data.SgrA_mass_s
BH_time_scale
Out[12]:
20.19522512635793

The strain spectral sensitivity of LISA detector:

In [13]:
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)
Out[13]:

The effective power spectral density of LISA noise:

In [14]:
psd = lisa_detector.power_spectral_density

The signal-to-noise ratio for one day of observations is:

In [15]:
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)
Out[15]:
103718.23722663395

Plot of the SNR as a function of the orbital radius $r_0$:

In [16]:
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
Out[16]: