Library for computing stable and plunging geodesics in Kerr spacetime
Project description
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 the latest version of kerrgeopy using pip
pip install kerrgeopy
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)$")
Authors
- Seyong Park
- Zach Nasipak
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for kerrgeopy-0.9.1-py3-none-any.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 7a49da84a71b11ba251566a74692955523e8faaf224aa7de1b164141ec1d568d |
|
MD5 | a60bacb6bc4b01fbb98f4608639f0d5a |
|
BLAKE2b-256 | 6d9ba5948c1a2d87083cade3c859f9daea1f5bfe555d2cd6e74511131669d438 |