KerrGeoPy

KerrGeoPy is a python implementation of the KerrGeodesics Mathematica library. It is intended for use in computing orbital trajectories for extreme-mass-ratio inspirals (EMRIs). It implements the analytical solutions for plunging orbits from Dyson and van de Meent, as well as solutions for stable orbits from Fujita and Hikida. The library also provides a set of methods for computing constants of motion and orbital frequencies. See the documentation for more information.

Installation

Install using Anaconda

conda install -c conda-forge kerrgeopy

or using pip

pip install kerrgeopy

Note

This library uses functions introduced in scipy 1.8, so it may also be necessary to update scipy by running pip install scipy -U, although in most cases this should be done automatically by pip. Certain plotting and animation functions also make use of features introduced in matplotlib 3.7 and rely on ffmpeg, which can be easily installed using homebrew or anaconda.

Contributing

For contribution guidelines, see CONTRIBUTING.

Stable Bound Orbits

KerrGeoPy computes orbits in Boyer-Lindquist coordinates $(t,r,\theta,\phi)$. Let $M$ to represent the mass of the primary body and let $J$ represent its angular momentum. Working in geometrized units where $G=c=1$, stable bound orbits are parametrized using the following variables:

$a$ - spin of the primary body
$p$ - orbital semilatus rectum
$e$ - orbital eccentricity
$x$ - cosine of the orbital inclination

$$ a = \frac{J}{M^2}, \quad\quad p = \frac{2r_{\text{min}}r_{\text{max}}}{M(r_{\text{min}}+r_{\text{max}})}, \quad\quad e = \frac{r_{\text{max}}-r_{\text{min}}}{r_{\text{max}}+r_{\text{min}}}, \quad\quad x = \cos{\theta_{\text{inc}}} $$

Note that $a$ and $x$ are restricted to values between -1 and 1, while $e$ is restricted to values between 0 and 1. Retrograde orbits are represented using a negative value for $a$ or for $x$. Polar orbits, marginally bound orbits, and orbits around an extreme Kerr black hole are not supported.

First, construct a StableOrbit using the four parameters described above.

import kerrgeopy as kg
from math import cos, pi

orbit = kg.StableOrbit(0.999,3,0.4,cos(pi/6))

Plot the orbit from $\lambda = 0$ to $\lambda = 10$ using the plot() method

fig, ax = orbit.plot(0,10)

Next, compute the time, radial, polar and azimuthal components of the trajectory as a function of Mino time using the trajectory() method. By default, the time and radial components of the trajectory are given in geometrized units and are normalized using $M$ so that they are dimensionless.

t, r, theta, phi = orbit.trajectory()

import numpy as np
import matplotlib.pyplot as plt

time = np.linspace(0,20,200)

plt.figure(figsize=(20,4))

plt.subplot(1,4,1)
plt.plot(time, t(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$t(\lambda)$")

plt.subplot(1,4,2)
plt.plot(time, r(time))
plt.xlabel("$\lambda$")
plt.ylabel("$r(\lambda)$")

plt.subplot(1,4,3)
plt.plot(time, theta(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\theta(\lambda)$")

plt.subplot(1,4,4)
plt.plot(time, phi(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\phi(\lambda)$")

Orbital Properties

Use the constants_of_motion() method to compute the dimensionless energy, angular momentum and Carter constant. By default, constants of motion are given in geometrized units where $G=c=1$ and are scale-invariant, meaning that they are normalized according to the masses of the two bodies as follows:

$$ \mathcal{E} = \frac{E}{\mu}, \quad \mathcal{L} = \frac{L}{\mu M}, \quad \mathcal{Q} = \frac{Q}{\mu^2 M^2} $$

Here, $M$ is the mass of the primary body and $\mu$ is the mass of the secondary body.

Frequencies of motion can be computed in Mino time using the mino_frequencies() method and in Boyer-Lindquist time using the fundamental_frequencies() method. As with constants of motion, the frequencies returned by both methods are given in geometrized units and are normalized by $M$ so that they are dimensionless.

from IPython.display import display, Math

E, L, Q = orbit.constants_of_motion()

upsilon_r, upsilon_theta, upsilon_phi, gamma = orbit.mino_frequencies()

omega_r, omega_theta, omega_phi = orbit.fundamental_frequencies()


display(Math(fr"a = {orbit.a} \quad p = {orbit.p} \quad e = {orbit.e} \quad x = {orbit.x}"))

display(Math(fr"\mathcal{{E}} = {E:.3f} \quad \mathcal{{L}} = {L:.3f} \quad \mathcal{{Q}} = {Q:.3f}"))

display(Math(fr"""\Upsilon_r = {upsilon_r:.3f} \quad 
             \Upsilon_\theta = {upsilon_theta:.3f} \quad 
             \Upsilon_\phi = {upsilon_phi:.3f} \quad 
             \Gamma = {gamma:.3f}"""))

display(Math(fr"""\Omega_r = {omega_r:.3f} \quad
            \Omega_\theta = {omega_theta:.3f} \quad
            \Omega_\phi = {omega_phi:.3f}"""))

$\displaystyle a = 0.999 \quad p = 3 \quad e = 0.4 \quad x = 0.8660254037844387$

$\displaystyle \mathcal{E} = 0.877 \quad \mathcal{L} = 1.903 \quad \mathcal{Q} = 1.265$

$\displaystyle \Upsilon_r = 1.145 \quad \Upsilon_\theta = 2.243 \quad \Upsilon_\phi = 3.118 \quad \Gamma = 20.531$

$\displaystyle \Omega_r = 0.056 \quad \Omega_\theta = 0.109 \quad \Omega_\phi = 0.152$

Plunging Orbits

Plunging orbits are parametrized using the spin parameter and the three constants of motion.

$a$ - spin of the primary body
$\mathcal{E}$ - Energy
$\mathcal{L}$ - $z$-component of angular momentum
$\mathcal{Q}$ - Carter constant

It is assumed that all orbital parameters are given in geometrized units where $G=c=1$ and are normalized according to the masses of the two bodies as follows:

$$ a = \frac{J}{M^2}, \quad \mathcal{E} = \frac{E}{\mu}, \quad \mathcal{L} = \frac{L}{\mu M}, \quad \mathcal{Q} = \frac{Q}{\mu^2 M^2} $$

Construct a PlungingOrbit by passing in these four parameters.

orbit = kg.PlungingOrbit(0.9, 0.94, 0.1, 12)

As with stable orbits, the components of the trajectory can be computed using the trajectory() method

t, r, theta, phi = orbit.trajectory()

import numpy as np
import matplotlib.pyplot as plt

time = np.linspace(0,20,200)

plt.figure(figsize=(20,4))

plt.subplot(1,4,1)
plt.plot(time, t(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$t(\lambda)$")

plt.subplot(1,4,2)
plt.plot(time, r(time))
plt.xlabel("$\lambda$")
plt.ylabel("$r(\lambda)$")

plt.subplot(1,4,3)
plt.plot(time, theta(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\theta(\lambda)$")

plt.subplot(1,4,4)
plt.plot(time, phi(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\phi(\lambda)$")

Alternative Parametrizations

Use the from_constants() class method to construct a StableOrbit from the spin parameter and constants of motion $(a,E,L,Q)$

orbit = kg.StableOrbit.from_constants(0.9, 0.95, 1.6, 8)

Use the Orbit class to construct an orbit from the spin parameter $a$, initial position $(t_0,r_0,\theta_0,\phi_0)$ and initial four-velocity $(u^t_0,u^r_0,u^{\theta}_0,u^{\phi}_0)$

stable_orbit = kg.StableOrbit(0.999,3,0.4,cos(pi/6))

x0 = stable_orbit.initial_position
u0 = stable_orbit.initial_velocity

orbit = kg.Orbit(0.999,x0,u0)

t, r, theta, phi = orbit.trajectory()

time = np.linspace(0,20,200)

plt.figure(figsize=(20,4))

plt.subplot(1,4,1)
plt.plot(time, t(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$t(\lambda)$")

plt.subplot(1,4,2)
plt.plot(time, r(time))
plt.xlabel("$\lambda$")
plt.ylabel("$r(\lambda)$")

plt.subplot(1,4,3)
plt.plot(time, theta(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\theta(\lambda)$")

plt.subplot(1,4,4)
plt.plot(time, phi(time))
plt.xlabel("$\lambda$")
plt.ylabel(r"$\phi(\lambda)$")

Citation

If you use this software, please cite our article in the Journal of Open Source Software.

@article{kerrgeopy, 
doi = {10.21105/joss.06587}, 
url = {https://doi.org/10.21105/joss.06587}, 
year = {2024}, 
publisher = {The Open Journal}, 
volume = {9}, 
number = {98}, 
pages = {6587}, 
author = {Seyong Park and Zachary Nasipak}, 
title = {KerrGeoPy: A Python Package for Computing Timelike Geodesics in Kerr Spacetime}, 
journal = {Journal of Open Source Software} 
} 

Authors

• Seyong Park
• Zach Nasipak