# Gravitational radiation by a single particle¶

Gravitational radiation by a particle on a circular orbit around a Kerr black hole

The gravitational wave emitted by a particle of mass $$\mu$$ on a circular orbit in the equatorial plane of a Kerr black hole of mass $$M$$ and angular momentum parameter $$a$$ is given by the formula:

(1)$h_+ - i h_\times = \frac{2\mu}{r} \, \sum_{\ell=2}^{\infty} \sum_{{\scriptstyle m=-\ell\atop \scriptstyle m\not=0}}^\ell \frac{Z^\infty_{\ell m}(r_0)}{(m\omega_0)^2} \, _{-2}S^{am\omega_0}_{\ell m}(\theta,\phi) \, e^{-i m \phi_0} e^{- i m \omega_0 (t-r_*)}$

where

• $$h_+ = h_+(t,r,\theta,\phi)$$ and $$h_\times = h_\times(t,r,\theta,\phi)$$, $$(t,r,\theta,\phi)$$ being the Boyer-Lindquist coordinates of the observer

• $$r_*$$ is the tortoise coordinate corresponding to $$r$$

• $$r_0$$ is the Boyer-Lindquist radius of the particle’s orbit

• $$\phi_0$$ is some constant phase factor

• $$\omega_0$$ is the orbital angular velocity

• $$Z^\infty_{\ell m}(r_0)$$ is a solution of the radial component of the Teukolsky equation (cf. Zinf())

• $$_{-2}S^{am\omega_0}_{\ell m}(\theta,\phi)$$ is the spin-weighted spheroidal harmonic of weight $$-2$$ (cf. spin_weighted_spheroidal_harmonic())

According to Eq. (1), the Fourier-series expansion of the waveform $$(h_+,h_\times)$$ received at the location $$(t,r,\theta,\phi)$$ is

(2)$h_{+,\times}(t, r, \theta, \phi) = \sum_{m=1}^{+\infty} \left[ A_m^{+,\times}(r,\theta) \cos(m\psi) + B_m^{+,\times}(r,\theta)\sin(m\psi) \right] ,$

where

$\psi := \omega_0 (t - r_*) - \phi + \phi_0 ,$

$$\omega_0$$ being the orbital angular velocity of the particle and $$r_*$$ the tortoise coordinate corresponding to $$r$$.

Note that the dependence of the Fourier coefficients $$A_m^{+,\times}(r,\theta)$$ and $$B_m^{+,\times}(r,\theta)$$ with respect to $$r$$ is simply $$\mu/r$$, where $$\mu$$ is the particle’s mass, i.e. we may consider the following rescaled Fourier coefficients, which depend on $$\theta$$ only:

${\bar A}_m^{+,\times}(\theta) := \frac{r}{\mu} A_m^{+,\times}(r,\theta) \quad\mbox{and}\quad {\bar B}_m^{+,\times}(\theta) := \frac{r}{\mu} B_m^{+,\times}(r,\theta)$

According to Eqs. (1) and (2), we have

(3)${\bar A}_m^+(\theta) = \frac{2}{(m\omega_0)^2} \sum_{\ell=2}^{\infty} \mathrm{Re}\left( Z^\infty_{\ell m}(r_0) \right) \left[ (-1)^\ell\, {}_{-2}S^{-am\omega_0}_{\ell,- m}(\theta,0) + {}_{-2}S^{am\omega_0}_{\ell m}(\theta,0) \right]$
(4)${\bar B}_m^+(\theta) = \frac{2}{(m\omega_0)^2} \sum_{\ell=2}^{\infty} \mathrm{Im}\left( Z^\infty_{\ell m}(r_0) \right) \left[ (-1)^\ell\, {}_{-2}S^{-am\omega_0}_{\ell,- m}(\theta,0) + {}_{-2}S^{am\omega_0}_{\ell m}(\theta,0) \right]$
(5)${\bar A}_m^\times(\theta) = \frac{2}{(m\omega_0)^2} \sum_{\ell=2}^{\infty} \mathrm{Im}\left( Z^\infty_{\ell m}(r_0) \right) \left[ (-1)^\ell\, {}_{-2}S^{-am\omega_0}_{\ell,- m}(\theta,0) - {}_{-2}S^{am\omega_0}_{\ell m}(\theta,0) \right]$
(6)${\bar B}_m^\times(\theta) = \frac{2}{(m\omega_0)^2} \sum_{\ell=2}^{\infty} \mathrm{Re}\left( Z^\infty_{\ell m}(r_0) \right) \left[ (-1)^{\ell+1}\, {}_{-2}S^{-am\omega_0}_{\ell,- m}(\theta,0) + {}_{-2}S^{am\omega_0}_{\ell m}(\theta,0) \right]$

This module implements the following functions:

REFERENCES:

• S. A. Teukolsky, Astrophys. J. 185, 635 (1973)

• S. Detweiler, Astrophys. J. 225, 687 (1978)

• M. Shibata, Phys. Rev. D 50, 6297 (1994)

• D. Kennefick, Phys. Rev. D 58, 064012 (1998)

• S. A. Hughes, Phys. Rev. D 61, 084004 (2000) [doi:10.1103/PhysRevD.61.084004]

• E. Gourgoulhon, A. Le Tiec, F. Vincent, N. Warburton, Astron. Astrophys. 627, A92 (2019) [doi:10.1051/0004-6361/201935406]; Arxiv 1903.02049

kerrgeodesic_gw.gw_particle.decay_time(a, r_init, r_final, l_max=None, m_min=1, approximation=None, quad_epsrel=1e-06, quad_limit=500)

Return the time spent in the migration from a circular orbit of radius r_init to that of radius r_final, induced by gravitational radiation reaction.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r_init – Boyer-Lindquist radius of the initial orbit (in units of $$M$$)

• r_final – Boyer-Lindquist radius of the final orbit (in units of $$M$$)

• l_max – (default: None) upper bound in the summation over the harmonic degree $$\ell$$; if None, l_max is determined automatically from the available tabulated data

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• approximation – (default: None) string describing the computational method; allowed values are

• None: exact computation

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• 'quadrupole' (only for a=0): quadrupole approximation (0-post-Newtonian); see h_particle_quadrupole()

• quad_epsrel – (default: 1.e-6) relative error tolerance in the computation of the integral giving the decay time

• quad_limit – (default: 500) upper bound on the number of subintervals used in the adaptive algorithm to compute the integral (this corresponds to the argument limit of SciPy’s function quad)

OUTPUT:

• rescaled decay time $$T (\mu/M^2)$$, where $$M$$ is the BH mass and $$\mu$$ the mass of the orbiting particle.

EXAMPLES:

Time to migrate from $$r_0=10 M$$ to $$r_0=6 M$$ around a Schwarzschild black hole:

sage: from kerrgeodesic_gw import decay_time
sage: decay_time(0, 10, 6)  # tol 1.0e-13
90.80876605028857


Let us check that at large radius, there is a good agreement with the quadrupole formula:

sage: decay_time(0, 50, 6)  # tol 1.0e-13
121893.29664651724
sage: decay_time(0, 50, 6, approximation='quadrupole')  # tol 1.0e-13
122045.0

kerrgeodesic_gw.gw_particle.h_amplitude_particle_fourier(m, a, r0, theta, l_max=10, algorithm_Zinf='spline')

Return the amplitude Fourier mode of a given order $$m$$ of the rescaled gravitational wave emitted by a particle in circular orbit around a Kerr black hole.

The rescaled Fourier mode of order $$m$$ received at the location $$(t,r,\theta,\phi)$$ is

$\frac{r}{\mu} h_m^{+,\times} = {\bar A}_m^{+,\times} \cos(m\psi) + {\bar B}_m^{+,\times} \sin(m\psi)$

where $$\mu$$ is the particle mass and $$\psi := \omega_0 (t-r_*) - \phi$$, $$\omega_0$$ being the orbital frequency of the particle and $$r_*$$ the tortoise coordinate corresponding to $$r$$ and $${\bar A}_m^{+,\times}$$ and $${\bar B}_m^{+,\times}$$ are given by Eqs. (3)-(6) above.

The $$+$$ and $$\times$$ amplitudes of the Fourier mode $$m$$ are defined respectively by

$\frac{r}{\mu} |h_m^+| := \sqrt{({\bar A}_m^+)^2 + ({\bar B}_m^+)^2} \quad\mbox{and}\quad \frac{r}{\mu} |h_m^\times| := \sqrt{({\bar A}_m^\times)^2 + ({\bar B}_m^\times)^2}$

INPUT:

• m – positive integer defining the Fourier mode

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$ in Eqs. (3)-(6)

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

OUTPUT:

• tuple $$((r/\mu)|h_m^+|,\ (r/\mu)|h_m^\times|)$$ (cf. the above expression)

EXAMPLE:

For a Schwarzschild black hole ($$a=0$$):

sage: from kerrgeodesic_gw import h_amplitude_particle_fourier
sage: a = 0
sage: h_amplitude_particle_fourier(2, a, 6., pi/2)  # tol 1.0e-13
(0.27875846152963557, 1.5860176188287866e-16)
sage: h_amplitude_particle_fourier(2, a, 6., pi/4)  # tol 1.0e-13
(0.47180033963220214, 0.45008696580919527)
sage: h_amplitude_particle_fourier(2, a, 6., 0)  # tol 1.0e-13
(0.6724377101568336, 0.6724377101568336)
sage: h_amplitude_particle_fourier(2, a, 6., pi/4, l_max=5)  # tol 1.0e-13
(0.47179830286565255, 0.4500948389153302)
sage: h_amplitude_particle_fourier(2, a, 6., pi/4, l_max=5,  # tol 1.0e-13
....:                              algorithm_Zinf='1.5PN')
(0.5381495951380861, 0.5114366815383188)


For a rapidly rotating Kerr black hole ($$a=0.95 M$$):

sage: a = 0.95
sage: h_amplitude_particle_fourier(2, a, 6., pi/4)  # tol 1.0e-13
(0.39402068296301823, 0.37534143024659444)
sage: h_amplitude_particle_fourier(2, a, 2., pi/4)  # tol 1.0e-13
(0.7358730645589858, 0.7115113031184368)

kerrgeodesic_gw.gw_particle.h_cross_particle(a, r0, u, theta, phi, phi0=0, l_max=10, m_min=1, algorithm_Zinf='spline')

Return the rescaled $$h_\times$$-part of the gravitational radiation emitted by a particle in circular orbit around a Kerr black hole.

The computation is based on Eq. (2) above.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• u – retarded time coordinate of the observer (in units of $$M$$): $$u = t - r_*$$, where $$t$$ is the Boyer-Lindquist time coordinate and $$r_*$$ is the tortoise coordinate

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• phi – Boyer-Lindquist azimuthal coordinate $$\phi$$ of the observer

• phi0 – (default: 0) phase factor

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

OUTPUT:

• the rescaled waveform $$(r / \mu) h_\times$$, where $$\mu$$ is the particle’s mass and $$r$$ is the Boyer-Lindquist radial coordinate of the observer

EXAMPLES:

Let us consider the case $$a=0$$ (Schwarzschild black hole) and $$r_0=6 M$$ (emission from the ISCO). For $$\theta=\pi/2$$, we have $$h_\times=0$$:

sage: from kerrgeodesic_gw import h_cross_particle
sage: a = 0
sage: h_cross_particle(a, 6., 0., pi/2, 0.)  # tol 1.0e-13
1.0041370414185673e-16


while for $$\theta=\pi/4$$, we have:

sage: h_cross_particle(a, 6., 0., pi/4, 0.)  # tol 1.0e-13
0.275027796440582
sage: h_cross_particle(a, 6., 0., pi/4, 0., l_max=5)  # tol 1.0e-13
0.2706516303570341
sage: h_cross_particle(a, 6., 0., pi/4, 0., l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
0.2625307460899205


For an orbit of larger radius ($$r_0=12 M$$), the 1.5-post-Newtonian approximation is in better agreement with the exact computation:

sage: h_cross_particle(a, 12., 0., pi/4, 0.)
0.1050751824554463
sage: h_cross_particle(a, 12., 0., pi/4, 0., l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
0.10244926162224487


A plot of the waveform generated by a particle orbiting at the ISCO:

sage: hc = lambda u: h_cross_particle(a, 6., u, pi/4, 0.)
sage: plot(hc, (0, 200.), axes_labels=[r'$(t-r_*)/M$', r'$r h_\times/\mu$'],
....:      gridlines=True, frame=True, axes=False)
Graphics object consisting of 1 graphics primitive


Case $$a/M=0.95$$ (rapidly rotating Kerr black hole):

sage: a = 0.95
sage: h_cross_particle(a, 2., 0., pi/4, 0.)  # tol 1.0e-13
-0.2681353673743396


Assessing the importance of the mode $$m=1$$:

sage: h_cross_particle(a, 2., 0., pi/4, 0., m_min=2)  # tol 1.0e-13
-0.3010579420748449

kerrgeodesic_gw.gw_particle.h_cross_particle_fourier(m, a, r0, theta, l_max=10, algorithm_Zinf='spline')

Return the Fourier mode of a given order $$m$$ of the rescaled $$h_\times$$-part of the gravitational wave emitted by a particle in circular orbit around a Kerr black hole.

The rescaled Fourier mode of order $$m$$ received at the location $$(t,r,\theta,\phi)$$ is

$\frac{r}{\mu} h_m^\times = {\bar A}_m^\times \cos(m\psi) + {\bar B}_m^\times \sin(m\psi)$

where $$\mu$$ is the particle mass and $$\psi := \omega_0 (t-r_*) - \phi$$, $$\omega_0$$ being the orbital frequency of the particle and $$r_*$$ the tortoise coordinate corresponding to $$r$$ and $${\bar A}_m^\times$$ and $${\bar B}_m^\times$$ are given by Eqs. (5)-(6) above.

INPUT:

• m – positive integer defining the Fourier mode

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$ in Eqs. (5)-(6)

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

OUTPUT:

• tuple $$({\bar A}_m^\times, {\bar B}_m^\times)$$

EXAMPLES:

Let us consider the case $$a=0$$ first (Schwarzschild black hole), with $$m=2$$:

sage: from kerrgeodesic_gw import h_cross_particle_fourier
sage: a = 0


$$h_m^\times$$ is always zero in the direction $$\theta=\pi/2$$:

sage: h_cross_particle_fourier(2, a, 8., pi/2)  # tol 1.0e-13
(3.444996575846961e-17, 1.118234985040581e-16)


Let us then evaluate $$h_m^\times$$ in the direction $$\theta=\pi/4$$:

sage: h_cross_particle_fourier(2, a, 8., pi/4)  # tol 1.0e-13
(0.09841144532628172, 0.31201728756415015)
sage: h_cross_particle_fourier(2, a, 8., pi/4, l_max=5)  # tol 1.0e-13
(0.09841373523119075, 0.31202073689061305)
sage: h_cross_particle_fourier(2, a, 8., pi/4, l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
(0.11744823578781578, 0.34124272645755677)


Values of $$m$$ different from 2:

sage: h_cross_particle_fourier(3, a, 20., pi/4)  # tol 1.0e-13
(0.022251439699635174, -0.005354134279052387)
sage: h_cross_particle_fourier(3, a, 20., pi/4, l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
(0.017782177999686274, 0.0)
sage: h_cross_particle_fourier(1, a, 8., pi/4)  # tol 1.0e-13
(0.03362589948237155, -0.008465651545641889)
sage: h_cross_particle_fourier(0, a, 8., pi/4)
(0, 0)


The case $$a/M=0.95$$ (rapidly rotating Kerr black hole):

sage: a = 0.95
sage: h_cross_particle_fourier(2, a, 8., pi/4)  # tol 1.0e-13
(0.08843838892991202, 0.28159329265206867)
sage: h_cross_particle_fourier(8, a, 8., pi/4)  # tol 1.0e-13
(-0.00014588821622195107, -8.179557811364057e-06)
sage: h_cross_particle_fourier(2, a, 2., pi/4)  # tol 1.0e-13
(-0.6021994882885746, 0.3789513303450391)
sage: h_cross_particle_fourier(8, a, 2., pi/4)  # tol 1.0e-13
(0.01045760329050054, -0.004986913120370192)

kerrgeodesic_gw.gw_particle.h_particle_quadrupole(r0, u, theta, phi, phi0=0, mode='+')

Return the rescaled $$h_+$$ or $$h_\times$$ part of the gravitational radiation emitted by a particle in quasi-Newtonian circular orbit around a massive body, computed at the quadrupole approximation.

The computation is performed according to the following formulas:

$\begin{split}\begin{array}{l} \displaystyle h_+ = 2\, \frac{\mu}{r} \frac{M}{r_0} (1+\cos^2\theta) \cos\left[2\omega_0 (t - r) + 2(\phi_0-\phi)\right] \\ \displaystyle h_\times = 4\, \frac{\mu}{r} \frac{M}{r_0} \cos\theta \sin\left[2\omega_0 (t - r) + 2(\phi_0-\phi)\right] \end{array}\end{split}$

where $$M$$ is the mass of the central body, $$\mu\ll M$$ the mass of the orbiting particle, $$r_0$$ the orbital radius, $$\omega_0 = \sqrt{M/r_0^3}$$ the corresponding orbital angular velocity and $$(t, r, \theta,\phi)$$ the coordinates of the observer.

INPUT:

• r0 – radius of the orbit (in units of $$M$$, the mass of the central body)

• u – retarded time coordinate $$u = t - r$$ of the observer (in units of $$M$$)

• theta – colatitute $$\theta$$ of the observer

• phi – azimuthal coordinate $$\phi$$ of the observer

• phi0 – (default: 0) phase factor

• mode – (default: '+') string determining which GW polarization mode is considered; allowed values are '+' and 'x', for respectively $$h_+$$ and $$h_\times$$

OUTPUT:

• the rescaled waveform $$(r / \mu) h$$, where $$\mu$$ is the particle’s mass, $$r$$ is the radial coordinate of the observer and $$h$$ is either $$h_+$$ or $$h_\times$$ (depending on the value of mode)

EXAMPLES:

Values of $$h_+$$ for $$r_0 = 12 M$$:

sage: from kerrgeodesic_gw import (h_particle_quadrupole,
....:                              h_plus_particle, h_cross_particle)
sage: theta, phi = pi/3, 0
sage: r0 = 12.
sage: porb = n(2*pi*r0^1.5)  # orbital period
sage: u = porb/16
sage: h_particle_quadrupole(r0, u, theta, phi)  # tol 1.0e-13
0.1473139127471974
sage: h_plus_particle(0., r0, u, theta, phi)  # exact value, tol 1.0e-13
0.041745779012809306
sage: h_plus_particle(0., r0, u, theta, phi, l_max=5,
....:                 algorithm_Zinf='1.5PN')  # 1.5 PN approx, tol 1.0e-13
0.06396788402755646


Values of $$h_\times$$ for $$r_0 = 12 M$$:

sage: h_particle_quadrupole(r0, u, theta, phi, mode='x')  # tol 1.0e-13
0.1178511301977579
sage: h_cross_particle(0., r0, u, theta, phi)  # exact value, tol 1.0e-13
0.1351232731503482
sage: h_cross_particle(0., r0, u, theta, phi, l_max=5,
....:                  algorithm_Zinf='1.5PN')  # 1.5 PN approx, tol 1.0e-13
0.14924631043762673


Values of $$h_+$$ for $$r_0 = 50 M$$:

sage: r0 = 50.
sage: porb = n(2*pi*r0^1.5)  # orbital period
sage: u = porb/16
sage: h_particle_quadrupole(r0, u, theta, phi)  # tol 1.0e-13
0.03535533905932738
sage: h_plus_particle(0., r0, u, theta, phi, l_max=5)  # exact value, tol 1.0e-13
0.024850367000986223
sage: h_plus_particle(0., r0, u, theta, phi, l_max=5,
....:                 algorithm_Zinf='1.5PN')  # 1.5 PN approx, tol 1.0e-13
0.026007035506911764


The difference between the exact value and the one resulting from the quadrupole approximation looks large, but actually this results mostly from some dephasing, as one can see on a plot:

sage: hp_quad = lambda t: h_particle_quadrupole(r0, t, theta, phi)
sage: hc_quad = lambda t: h_particle_quadrupole(r0, t, theta, phi, mode='x')
sage: hp = lambda t: h_plus_particle(0., r0, t, theta, phi, l_max=5)
sage: hc = lambda t: h_cross_particle(0., r0, t, theta, phi, l_max=5)
sage: umax = 2*porb
sage: plot(hp_quad, (0, umax), legend_label=r'$h_+$ quadrupole',
....:      axes_labels=[r'$(t-r_*)/M$', r'$r h/\mu$'],
....:      title=r'$r_0=50 M$, $\theta=\pi/3$',
....:      gridlines=True, frame=True, axes=False) + \
....: plot(hp, (0, umax), color='red', legend_label=r'$h_+$ exact') + \
....: plot(hc_quad, (0, umax), linestyle='--',
....:      legend_label=r'$h_\times$ quadrupole') + \
....: plot(hc, (0, umax), color='red', linestyle='--',
....:      legend_label=r'$h_\times$ exact')
Graphics object consisting of 4 graphics primitives

kerrgeodesic_gw.gw_particle.h_particle_signal(a, r0, theta, phi, u_min, u_max, mode='+', nb_points=100, phi0=0, l_max=10, m_min=1, approximation=None, store=None)

Return a time sequence of the $$h_+$$ or the $$h_\times$$ part of the gravitational radiation from a particle in circular orbit around a Kerr black hole.

Note

It is more efficient to use this function than to perform a loop over h_plus_particle() or h_cross_particle(). Indeed, the Fourier modes, which involve the computation of spin-weighted spheroidal harmonics and of the functions $$Z^\infty_{\ell m}(r_0)$$, are evaluated once for all, prior to the loop on the retarded time $$u$$.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• phi – Boyer-Lindquist azimuthal coordinate $$\phi$$ of the observer

• u_min – lower bound of the retarded time coordinate of the observer (in units of the black hole mass $$M$$): $$u = t - r_*$$, where $$t$$ is the Boyer-Lindquist time coordinate and $$r_*$$ is the tortoise coordinate

• u_max – upper bound of the retarded time coordinate of the observer (in units of the black hole mass $$M$$)

• mode – (default: '+') string determining which GW polarization mode is considered; allowed values are '+' and 'x', for respectively $$h_+$$ and $$h_\times$$

• nb_points – (default: 100) number of points in the interval (u_min, u_max)

• phi0 – (default: 0) phase factor

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• approximation – (default: None) string describing the computational method; allowed values are

• None: exact computation

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• 'quadrupole' (only for a=0): quadrupole approximation (0-post-Newtonian); see h_particle_quadrupole()

• store – (default: None) string containing a file name for storing the time sequence; if None, no storage is attempted

OUTPUT:

• a list of nb_points pairs $$(u, r h/\mu)$$, where $$u$$ is the retarded time, $$h$$ is either $$h_+$$ or $$h_\times$$ depending on the parameter mode, $$\mu$$ is the particle mass and $$r$$ is the Boyer-Lindquist radial coordinate of the observer

EXAMPLES:

$$h_+$$ signal at $$\theta=\pi/2$$ from a particle at the ISCO of a Schwarzschild black hole ($$a=0$$, $$r_0=6M$$):

sage: from kerrgeodesic_gw import h_particle_signal
sage: h_particle_signal(0., 6., pi/2, 0., 0., 200., nb_points=9)  # tol 1.0e-13
[(0.000000000000000, 0.1536656546005028),
(25.0000000000000, -0.2725878162016865),
(50.0000000000000, 0.3525164756465054),
(75.0000000000000, 0.047367530900643974),
(100.000000000000, -0.06816472285771447),
(125.000000000000, -0.10904082076122341),
(150.000000000000, 0.11251491162759894),
(175.000000000000, 0.2819301792449237),
(200.000000000000, -0.24646401049292863)]


Storing the sequence in a file:

sage: h = h_particle_signal(0., 6., pi/2, 0., 0., 200.,
....:                       nb_points=9, store='h_plus.d')


The $$h_\times$$ signal, for $$\theta=\pi/4$$:

sage: h_particle_signal(0., 6., pi/4, 0., 0., 200., nb_points=9, mode='x')  # tol 1.0e-13
[(0.000000000000000, 0.275027796440582),
(25.0000000000000, -0.18713017721920192),
(50.0000000000000, 0.2133141583155321),
(75.0000000000000, -0.531073507307601),
(100.000000000000, 0.3968872953624949),
(125.000000000000, -0.4154274307718398),
(150.000000000000, 0.5790969355083798),
(175.000000000000, -0.24074783639714234),
(200.000000000000, 0.22869838143661578)]

kerrgeodesic_gw.gw_particle.h_plus_particle(a, r0, u, theta, phi, phi0=0, l_max=10, m_min=1, algorithm_Zinf='spline')

Return the rescaled $$h_+$$-part of the gravitational radiation emitted by a particle in circular orbit around a Kerr black hole.

The computation is based on Eq. (2) above.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• u – retarded time coordinate of the observer (in units of $$M$$): $$u = t - r_*$$, where $$t$$ is the Boyer-Lindquist time coordinate and $$r_*$$ is the tortoise coordinate

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• phi – Boyer-Lindquist azimuthal coordinate $$\phi$$ of the observer

• phi0 – (default: 0) phase factor

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

OUTPUT:

• the rescaled waveform $$(r / \mu) h_+$$, where $$\mu$$ is the particle’s mass and $$r$$ is the Boyer-Lindquist radial coordinate of the observer

EXAMPLES:

Let us consider the case $$a=0$$ (Schwarzschild black hole) and $$r_0=6 M$$ (emission from the ISCO):

sage: from kerrgeodesic_gw import h_plus_particle
sage: a = 0
sage: h_plus_particle(a, 6., 0., pi/2, 0.)  # tol 1.0e-13
0.1536656546005028
sage: h_plus_particle(a, 6., 0., pi/2, 0., l_max=5)  # tol 1.0e-13
0.157759938177291
sage: h_plus_particle(a, 6., 0., pi/2, 0., l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
0.22583887001798497


For an orbit of larger radius ($$r_0=12 M$$), the 1.5-post-Newtonian approximation is in better agreement with the exact computation:

sage: h_plus_particle(a, 12., 0., pi/2, 0.)  # tol 1.0e-13
0.11031251832047866
sage: h_plus_particle(a, 12., 0., pi/2, 0., l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
0.12935832450325302


A plot of the waveform generated by a particle orbiting at the ISCO:

sage: hp = lambda u: h_plus_particle(a, 6., u, pi/2, 0.)
sage: plot(hp, (0, 200.), axes_labels=[r'$(t-r_*)/M$', r'$r h_+/\mu$'],
....:      gridlines=True, frame=True, axes=False)
Graphics object consisting of 1 graphics primitive


Case $$a/M=0.95$$ (rapidly rotating Kerr black hole):

sage: a = 0.95
sage: h_plus_particle(a, 2., 0., pi/2, 0.)  # tol 1.0e-13
0.20326150400852214


Assessing the importance of the mode $$m=1$$:

sage: h_plus_particle(a, 2., 0., pi/2, 0., m_min=2)  # tol 1.0e-13
0.21845811047370495

kerrgeodesic_gw.gw_particle.h_plus_particle_fourier(m, a, r0, theta, l_max=10, algorithm_Zinf='spline')

Return the Fourier mode of a given order $$m$$ of the rescaled $$h_+$$-part of the gravitational wave emitted by a particle in circular orbit around a Kerr black hole.

The rescaled Fourier mode of order $$m$$ received at the location $$(t,r,\theta,\phi)$$ is

$\frac{r}{\mu} h_m^+ = {\bar A}_m^+ \cos(m\psi) + {\bar B}_m^+ \sin(m\psi)$

where $$\mu$$ is the particle mass and $$\psi := \omega_0 (t-r_*) - \phi$$, $$\omega_0$$ being the orbital frequency of the particle and $$r_*$$ the tortoise coordinate corresponding to $$r$$ and $${\bar A}_m^+$$ and $${\bar B}_m^+$$ are given by Eqs. (3)-(4) above.

INPUT:

• m – positive integer defining the Fourier mode

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$ in Eqs. (3)-(4)

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

OUTPUT:

• tuple $$({\bar A}_m^+, {\bar B}_m^+)$$

EXAMPLES:

Let us consider the case $$a=0$$ first (Schwarzschild black hole), with $$m=2$$:

sage: from kerrgeodesic_gw import h_plus_particle_fourier
sage: a = 0
sage: h_plus_particle_fourier(2, a, 8., pi/2)  # tol 1.0e-13
(0.2014580652208302, -0.06049343736886148)
sage: h_plus_particle_fourier(2, a, 8., pi/2, l_max=5)  # tol 1.0e-13
(0.20146097329552273, -0.060495372034569186)
sage: h_plus_particle_fourier(2, a, 8., pi/2, l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
(0.2204617125753912, -0.0830484439639611)


Values of $$m$$ different from 2:

sage: h_plus_particle_fourier(3, a, 20., pi/2)  # tol 1.0e-13
(-0.005101595598729037, -0.021302121442654077)
sage: h_plus_particle_fourier(3, a, 20., pi/2, l_max=5, algorithm_Zinf='1.5PN')  # tol 1.0e-13
(0.0, -0.016720919174427588)
sage: h_plus_particle_fourier(1, a, 8., pi/2)  # tol 1.0e-13
(-0.014348477201223874, -0.05844679575244101)
sage: h_plus_particle_fourier(0, a, 8., pi/2)
(0, 0)


The case $$a/M=0.95$$ (rapidly rotating Kerr black hole):

sage: a = 0.95
sage: h_plus_particle_fourier(2, a, 8., pi/2)  # tol 1.0e-13
(0.182748773431646, -0.05615306925896938)
sage: h_plus_particle_fourier(8, a, 8., pi/2)  # tol 1.0e-13
(-4.724709221209198e-05, 0.0006867183495228116)
sage: h_plus_particle_fourier(2, a, 2., pi/2)  # tol 1.0e-13
(0.1700402877617014, 0.33693580916655747)
sage: h_plus_particle_fourier(8, a, 2., pi/2)  # tol 1.0e-13
(-0.009367442995129153, -0.03555092085651877)

kerrgeodesic_gw.gw_particle.plot_h_particle(a, r0, theta, phi, u_min, u_max, plot_points=200, phi0=0, l_max=10, m_min=1, approximation=None, mode=('+', 'x'), color=None, linestyle=None, legend_label=('$h_+$', '$h_\\\\times$'), xlabel='$(t - r_*)/M$', ylabel=None, title=None, gridlines=True)

Plot the gravitational waveform emitted by a particle in circular orbit around a Kerr black hole.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the black hole mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• phi – Boyer-Lindquist azimuthal coordinate $$\phi$$ of the observer

• u_min – lower bound of the retarded time coordinate of the observer (in units of $$M$$): $$u = t - r_*$$, where $$t$$ is the Boyer-Lindquist time coordinate and $$r_*$$ is the tortoise coordinate

• u_max – upper bound of the retarded time coordinate of the observer (in units of $$M$$)

• plot_points – (default: 200) number of points involved in the sampling of the interval (u_min, u_max)

• phi0 – (default: 0) phase factor

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• approximation – (default: None) string describing the computational method; allowed values are

• None: exact computation

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• 'quadrupole' (only for a=0): quadrupole approximation (0-post-Newtonian); see h_particle_quadrupole()

• mode – (default: ('+', 'x')) string determining the plotted quantities: allowed values are '+' and 'x', for respectively $$h_+$$ and $$h_\times$$, as well as ('+', 'x') for plotting both polarization modes

• color – (default: None) a color (if mode = '+' or 'x') or a pair of colors (if mode = ('+', 'x')) for the plot(s); if None, the default colors are 'blue' for $$h_+$$ and 'red' for $$h_\times$$

• linestyle – (default: None) a line style (if mode = '+' or 'x') or a pair of line styles (if mode = ('+', 'x')) for the plot(s); if None, the default style is a solid line

• legend_label – (default: (r'$h_+$', r'$h_\times$')) labels for the plots of $$h_+$$ and $$h_\times$$; used only if mode is ('+', 'x')

• xlabel – (default: r'$(t - r_*)/M$') label of the $$x$$-axis

• ylabel – (default: None) label of the $$y$$-axis; if None, r'$(r/\mu)\, h_+$' is used for mode = '+', r'$(r/\mu)\, h_\times$' for mode = 'x' and r'$(r/\mu)\, h_{+,\times}$' for mode = ('+', 'x')

• title – (default: None) plot title; if None, the title is generated from a, r0 and theta (see the example below)

• gridlines – (default: True) indicates whether the gridlines are to be plotted

OUTPUT:

• a graphics object

EXAMPLES:

Plot of the gravitational waveform generated by a particle orbiting at the ISCO of a Kerr black hole with $$a=0.9 M$$:

sage: from kerrgeodesic_gw import plot_h_particle
sage: plot_h_particle(0.9, 2.321, pi/4, 0., 0., 70.)
Graphics object consisting of 2 graphics primitives


Plot of $$h_+$$ only, with some non-default options:

sage: plot_h_particle(0.9, 2.321, pi/4, 0., 0., 70., mode='+',
....:                 color='green', xlabel=r'$u/M$', gridlines=False,
....:                 title='GW from ISCO, $a=0.9M$')
Graphics object consisting of 1 graphics primitive

kerrgeodesic_gw.gw_particle.plot_spectrum_particle(a, r0, theta, mode='+', m_max=10, l_max=10, algorithm_Zinf='spline', color='blue', linestyle='-', thickness=2, legend_label=None, offset=0, xlabel=None, ylabel=None, title=None, gridlines=True)

Plot the Fourier spectrum of the gravitational radiation emitted by a particle in equatorial circular orbit around a Kerr black hole.

The Fourier spectrum is defined by the following sequence indexed by $$m$$:

$H^{+,\times}_m(\theta) := \frac{r}{\mu} \sqrt{(A_m^{+,\times}(\theta))^2 + (B_m^{+,\times}(\theta))^2}$

where the functions $$A_m^+$$, $$B_m^+$$, $$A_m^\times$$, $$B_m^\times$$ are defined by (2).

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• theta – Boyer-Lindquist colatitute $$\theta$$ of the observer

• mode – (default: '+') string determining which GW polarization mode is considered; allowed values are '+' and 'x', for respectively $$h_+$$ and $$h_\times$$

• m_max – (default: 10) maximal value of $$m$$

• l_max – (default: 10) upper bound in the summation over the harmonic degree $$\ell$$

• algorithm_Zinf – (default: 'spline') string describing the computational method for $$Z^\infty_{\ell m}(r_0)$$; allowed values are

• 'spline': spline interpolation of tabulated data

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• color – (default: 'blue') color of vertical lines

• linestyle – (default: '-') style of vertical lines

• legend_label – (default: None) legend label for this spectrum

• offset – (default: 0) horizontal offset for the position of the vertical lines

• xlabel – (default: None) label of the x-axis; if none is provided, the label is set to $$m$$

• ylabel – (default: None) label of the y-axis; if none is provided, the label is set to $$H_m^{+,\times}$$

• title – (default: None) plot title; if None, the title is generated from a, r0 and theta

• gridlines – (default: True) indicates whether the gridlines are to be plotted

OUTPUT:

• a graphics object

EXAMPLES:

Spectrum of gravitational radiation generated by a particle orbiting at the ISCO of a Schwarzschild black hole ($$a=0$$, $$r_0=6M$$):

sage: from kerrgeodesic_gw import plot_spectrum_particle
sage: plot_spectrum_particle(0, 6., pi/2)
Graphics object consisting of 10 graphics primitives

kerrgeodesic_gw.gw_particle.radiated_power_particle(a, r0, l_max=None, m_min=1, approximation=None)

Return the total (i.e. summed over all directions) power of the gravitational radiation emitted by a particle in circular orbit around a Kerr black hole.

The total radiated power (or luminosity) $$L$$ is computed according to the formula:

$L = \frac{\mu^2}{4\pi} \sum_{\ell=2}^{\infty} \sum_{{\scriptstyle m=-\ell\atop \scriptstyle m\not=0}}^\ell \frac{\left| Z^\infty_{\ell m}(r_0) \right| ^2}{(m\omega_0)^2}$

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• l_max – (default: None) upper bound in the summation over the harmonic degree $$\ell$$; if None, l_max is determined automatically from the available tabulated data

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• approximation – (default: None) string describing the computational method; allowed values are

• None: exact computation

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• 'quadrupole' (only for a=0): quadrupole approximation (0-post-Newtonian); see h_particle_quadrupole()

OUTPUT:

• rescaled radiated power $$L\, (M/\mu)^2$$, where $$M$$ is the BH mass and $$\mu$$ the mass of the orbiting particle.

EXAMPLES:

Power radiated by a particle at the ISCO of a Schwarzschild BH:

sage: from kerrgeodesic_gw import radiated_power_particle
0.000937262782177525


Power radiated by a particle at the ISCO of a rapidly rotating Kerr BH ($$a=0.98 M$$):

sage: radiated_power_particle(0.98, 1.61403)  # tol 1.0e-13
0.08629927053494096


Power computed according to various approximations:

sage: radiated_power_particle(0, 6., approximation='1.5PN')  # tol 1.0e-13
0.0010435864384751143
sage: radiated_power_particle(0, 6., approximation='quadrupole')  # tol 1.0e-13
0.000823045267489712


Let us check that at large radius, these approximations get closer to the actual result:

sage: radiated_power_particle(0, 50.)  # tol 1.0e-13
1.9624555021692417e-08
sage: radiated_power_particle(0, 50., approximation='1.5PN')  # tol 1.0e-13
1.958318200436188e-08
sage: radiated_power_particle(0, 50., approximation='quadrupole')  # tol 1.0e-13
2.0480000000000002e-08

kerrgeodesic_gw.gw_particle.secular_frequency_change(a, r0, l_max=None, m_min=1, approximation=None)

Return the gravitational-radiation induced change of the orbital frequency of a particle in circular orbit around a Kerr black hole.

INPUT:

• a – BH angular momentum parameter (in units of $$M$$, the BH mass)

• r0 – Boyer-Lindquist radius of the orbit (in units of $$M$$)

• l_max – (default: None) upper bound in the summation over the harmonic degree $$\ell$$; if None, l_max is determined automatically from the available tabulated data

• m_min – (default: 1) lower bound in the summation over the Fourier mode $$m$$

• approximation – (default: None) string describing the computational method; allowed values are

• None: exact computation

• '1.5PN' (only for a=0): 1.5-post-Newtonian expansion following E. Poisson, Phys. Rev. D 47, 1497 (1993) [doi:10.1103/PhysRevD.47.1497]

• 'quadrupole' (only for a=0): quadrupole approximation (0-post-Newtonian); see h_particle_quadrupole()

OUTPUT:

• rescaled fractional change in orbital frequency $$\dot{\omega}_0/\omega_0 \, (M^2/\mu)$$, where $$M$$ is the BH mass and $$\mu$$ the mass of the orbiting particle.

EXAMPLES:

Relative change in orbital frequency change at $$r_0 = 10 M$$ around a Schwarzschild black hole:

sage: from kerrgeodesic_gw import secular_frequency_change
sage: secular_frequency_change(0, 10)  # tol 1.0e-13
0.002701529506901975
sage: secular_frequency_change(0, 10, approximation='quadrupole')  # tol 1.0e-13
0.00192


At larger distance, the quadrupole approximation works better:

sage: secular_frequency_change(0, 50)  # tol 1.0e-13
3.048598175592674e-06
sage: secular_frequency_change(0, 50, approximation='quadrupole')  # tol 1.0e-13
3.072e-06


At the ISCO, $$\dot{\omega}_0$$ diverges:

sage: secular_frequency_change(0, 6)
+infinity


while the quadrupole approximation would have predict a finite value there:

sage: secular_frequency_change(0, 6, approximation='quadrupole')  # tol 1.0e-13
0.014814814814814817


Case of a Kerr black hole:

sage: secular_frequency_change(0.5, 6)  # tol 1.0e-13
0.020393023677356764