Gravitational wave emission by a single orbiting particle

This is a TEST notebook for the functions h_plus_particle and h_cross_particle. It does not make use of the function plot_h_particle from the kerrgeodesic_gw package; instead the plots are performed by loops calling h_plus_particle and h_cross_particle at each step. This is much less efficient than using plot_h_particle.

In [1]:
version()
Out[1]:
'SageMath version 9.4, Release Date: 2021-08-22'
In [2]:
%display latex
In [3]:
# For benchmarking
import time
system_time0 = time.time()
CPU_time0 = time.perf_counter()
In [4]:
from kerrgeodesic_gw import (h_plus_particle, h_cross_particle)

Directory to store the figure files (created if it does not exist already):

In [5]:
import os
figdir = "figures_gw_single_particle/"
if not os.path.exists(figdir):
    os.makedirs(figdir)
In [6]:
def ordital_period(a, r0):
    return RDF(2*pi*(r0^(3/2) + a))

Waveform from near the ISCO for $a=0$

Plot for comparison with Fig. 1 of E. Poisson, PRD 47, 1497 (1993):

In [7]:
a = 0.
r0 = 6.2 # near ISCO value used by Poisson (1993)
theta, phi = 4*pi/9, 0
norm_Poisson93 = -0.065^(2/3) # the minus sign arises from a different metric signature
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5, 
                               algorithm_Zinf='1.5PN') / norm_Poisson93
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5,
                                algorithm_Zinf='1.5PN') / norm_Poisson93
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$t/M$', r'$r h / \mu/(M\Omega)^{2/3}$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0,\quad r_0=6.2M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig1_Poisson93.pdf")
graph
Out[7]:

Waveform from the ISCO for $a=0$

In [8]:
a = 0.
r0 = 6.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$t/M$', r'$r h / \mu$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0,\quad r_0=6M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5, 
                               algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5, 
                                algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5,
              legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), color='red', linestyle=':', thickness=1.5, 
              legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r6.pdf")
graph
Out[8]:
In [9]:
show(graph, figsize=7)

Waveform from $r_0=12M$ for $a=0$

In [10]:
a = 0.
r0 = 12.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$t/M$', r'$r h / \mu$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0,\quad r_0=12M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5, 
                               algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5, 
                                algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5, 
              legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), linestyle=':', color='red', thickness=1.5, 
              legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r12.pdf")
graph
Out[10]:

Waveform from $r_0=50M$ for $a=0$

In [11]:
a = 0.
r0 = 50.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$t/M$', r'$r h / \mu$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0,\quad r_0=50M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5, 
                               algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5, 
                                algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5,
              legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), linestyle=':', color='red',thickness=1.5,
              legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r50.pdf")
graph
Out[11]:

Waveform from the ISCO for $a=0.5 M$

In [12]:
a = 0.5
r0 = 4.234
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=7)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=7)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$t/M$', r'$r h / \mu$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0.5M,\quad r_0=4.234M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "waveform_a0.5_r4.234.pdf")
graph
Out[12]:

Waveform from the ISCO for $a=0.9 M$

In [13]:
a = 0.9
r0 = 2.321
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', plot_points=200,
             axes_labels=[r'$t/M$', r'$r h / \mu$'], 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0.9M,\quad r_0=2.321M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$',
              plot_points=200)
graph.save(figdir + "waveform_a0.9_r2.321.pdf")
graph
Out[13]:

Comparison with Detweiler (1978)

Comparison with Fig. 8 of S.L. Detweiler, ApJ 225, 687 (1978) (given that Detweiler is using $-u = r_*-t$ on the $x$-axis, we use $-t$ as argument of h_plus_particle)

In [14]:
a = 0.5
r0 = 4.45
theta, phi = pi/2, 0
hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=7)  # NB: t --> -t
hc = lambda t: h_cross_particle(a, r0, -t, theta, phi, l_max=7)
graph = plot(hp, (0, 125), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$(r_*-t)/M$', r'$r h / \mu$'], ymin=-1, ymax=1, 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0.5M,\quad r_0=4.45M,\quad \theta=\frac{\pi}{2}$")
graph += plot(hc, (0, 125), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig8_Detweiler78_lmax7.pdf")
graph
Out[14]:

Comparison with Fig. 9 of S.L. Detweiler, ApJ 225, 687 (1978):

In [15]:
a = 0.9
r0 = 2.4
theta, phi = pi/2, 0
hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=7)  # NB: t --> -t
hc = lambda t: h_cross_particle(a, r0, -t, theta, phi, l_max=7)
graph = plot(hp, (0, 60), thickness=1.5, legend_label=r'$h_+$', 
             axes_labels=[r'$(r_*-t)/M$', r'$r h / \mu$'], ymin=-1, ymax=1, 
             gridlines=True, frame=True, axes=False, 
             title=r"$a=0.9M,\quad r_0=2.4M,\quad \theta=\frac{\pi}{2}$")
graph += plot(hc, (0, 60), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig8_Detweiler78_lmax_7.pdf")
graph
Out[15]:
In [16]:
a = 0.9
r0 = 2.4
theta, phi = pi/2, 0
graph = Graphics()
for l_max in [5, 7, 10]:
    hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=l_max)  # NB: t --> -t
    graph += plot(hp, (0, 60), thickness=1.5, color=hue((l_max-5)/8), 
                  legend_label=r"$\ell_{{\rm max}} = {}$".format(l_max),
                  axes_labels=[r'$(r_*-t)/M$', r'$r h_+ / \mu$'], ymin=-1., ymax=1., 
                  gridlines=True, frame=True, axes=False, 
                  title=r"$a=0.9M,\quad r_0=2.4M,\quad \theta=\frac{\pi}{2}$")
graph.save(figdir + "compar_fig9_Detweiler78.pdf")
graph
Out[16]:
In [17]:
print("System time elapsed since the start of this notebook: {} s".format(
                                                 time.time() - system_time0))
print("CPU time elapsed since the start of this notebook:    {} s".format(
                                                 time.perf_counter() - CPU_time0))
System time elapsed since the start of this notebook: 53.653218507766724 s
CPU time elapsed since the start of this notebook:    53.653356110999994 s