kerrgeodesic_gw
¶This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
It requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw (version $\geq$ 0.3). To install the latter, simply run
sage -pip install kerrgeodesic_gw
in a terminal.
version()
'SageMath version 9.2, Release Date: 2020-10-24'
First, we set up the notebook to use LaTeX-formatted display:
%display latex
and we ask for CPU demanding computations to be performed in parallel on 8 processes:
Parallelism().set(nproc=8)
A Kerr black bole is entirely defined by two parameters $(m, a)$, where $m$ is the black hole mass and $a$ is the black hole angular momentum divided by $m$.
In this notebook, we shall set $m=1$ and we denote the angular momentum parameter $a$ by the symbolic variable a
, using a0
for a specific numerical value:
a = var('a')
a0 = 0.998
The spacetime object is created as an instance of the class KerrBH
:
from kerrgeodesic_gw import KerrBH
M = KerrBH(a)
print(M)
Kerr spacetime M
The Boyer-Lindquist coordinate $r$ of the event horizon:
rH = M.event_horizon_radius()
rH
rH0 = rH.subs({a: a0})
rH0
The method boyer_lindquist_coordinates()
returns the chart of Boyer-Lindquist coordinates BL
and allows the user to instanciate the Python variables (t, r, th, ph)
to the coordinates $(t,r,\theta,\phi)$:
BL.<t, r, th, ph> = M.boyer_lindquist_coordinates()
BL
The metric tensor is naturally returned by the method metric()
:
g = M.metric()
g.display()
We set $\mu=1$ and pick some values of $E$, $L$ and $Q$, with $E<1$ to ensure that we are dealing with a bound geodesic:
mu = 1
E = 0.9
#L = 1.9
L = 2.
Q = 1.3
r = var('r')
R(r) = ((E^2 - 1)*r^4 + 2*r^3 + (a0^2*(E^2 - 1) - Q - L^2)*r^2
+ 2*(Q + (L - a0*E)^2)*r - a0^2*Q)
graph = plot(-R(r), (r, -1.5, 7.5), thickness=1.5,
axes_labels=[r'$r/m$', r'$-\mathcal{R}(r)/m^4$'])
graph
rp = find_root(R(r), 1.5, 4)
rp
ra = find_root(R(r), 5, 10)
ra
graph += line([(rp, 0), (ra, 0)], thickness=2.5, color='red') \
+ text(r'$r_{\rm p}$', (1.05*rp, -2.5), color='red', fontsize=16) \
+ text(r'$r_{\rm a}$', (1.02*ra, -2.5), color='red', fontsize=16)
graph += line([(rH0, -20), (rH0, 30)], thickness=2, color='black')
graph
zoom = plot(-R(r), (r, -0.1, 2.5), thickness=1.5,
axes_labels=[r'$r/m$', r'$-\mathcal{R}(r)/m^4$'])
zoom += line([(rp, 0), (2.5, 0)], thickness=2.5, color='red') \
+ text(r'$r_{\rm p}$', (1.01*rp, 0.2), color='red', fontsize=12)
zoom += line([(rH0, -0.8), (rH0, 1.8)], thickness=2, color='black')
zoom
graph = graph.inset(zoom, (0.7, 0.7, 0.3, 0.3), fontsize=8)
graph
graph.save("gek_R_potential.pdf")
Let us choose the initial point $P$ for the geodesic:
P = M.point((0, (rp + ra)/2, pi/2, 0), name='P')
print(P)
Point P on the Kerr spacetime M
A geodesic is constructed by providing the range $[\lambda_{\rm min},\lambda_{\rm max}]$ of the affine parameter $\lambda$, the initial point and either
We shall also specify some numerical value for the Kerr spin parameter $a$.
Here, we choose $\lambda\in[0, 600\, m]$, the option (ii) and $a=0.998 \,m$, where $m$ in the black hole mass::
lmax = 600 # lambda_max
Li = M.geodesic([0, lmax], P, mu=mu, E=E, L=L, Q=Q, a_num=0.998,
name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector:
The curve was correctly set. Parameters appearing in the differential system defining the curve are [a].
print(Li)
Geodesic Li of the Kerr spacetime M
The numerical integration of the geodesic equation is performed via integrate()
, by providing the integration step $\delta\lambda$ in units of $m$:
Li.integrate(step=0.005)
We can then plot the geodesic:
Li.plot(plot_points=2000, thickness=1.5)
Actually, many options can be passed to the method plot()
. For instance to a get a 3D spacetime diagram:
Li.plot(coordinates='txy', thickness=2)
or to get the trace of the geodesic in the $(x,y)$ plane:
graph = Li.plot(coordinates='xy', plot_points=2000, thickness=1.5,
axes_labels=[r'$x/m$', r'$y/m$'])
graph
graph.save("gek_timelike_xy.pdf")
or else to get the trace in the $(x,z)$ plane:
graph = Li.plot(coordinates='xz', plot_points=2000, thickness=1.5,
axes_labels=[r'$x/m$', r'$z/m$'])
graph
graph.save("gek_timelike_xz.pdf")
th = var('th', latex_name=r'\theta')
V(th) = cos(th)^2 * (a0^2*(1 - E^2) + L^2/sin(th)^2)
V(th)
graph = plot(V(th), (th, 0.7, pi-0.7), axes_labels=[r'$\theta$', r'$V(\theta)$'])
graph += line([(0, Q), (pi, Q)], color='red')
graph
th_min = find_root(V(th) - Q, 0.5, pi/2)
th_min
Analytic formula for $\theta_{\rm min}$:
A2 = a0^2*(mu^2 - E^2)
sthm = (A2 - L^2 - Q + sqrt((L^2 + Q - A2)^2 + 4*L^2*A2))/(2*A2)
th_min0 = arcsin(sqrt(sthm))
th_min0
th_min - th_min0
th_min_deg = n(th_min/pi*180)
th_min_deg
th_inc_deg = 90 - th_min_deg
th_inc_deg
p_K = 2*ra*rp/(ra + rp)
p_K
e_K = (ra - rp)/(ra + rp)
e_K
Let us check that the values of $\mu$, $E$, $L$ and $Q$ evaluated at $\lambda=300 \, m$ are equal to those at $\lambda=0$ up to the numerical accuracy of the integration scheme:
Li.check_integrals_of_motion(lmax)
quantity | value | initial value | diff. | relative diff. |
Decreasing the integration step leads to smaller errors:
Li.integrate(step=0.001)
Li.check_integrals_of_motion(lmax)
quantity | value | initial value | diff. | relative diff. |
We choose a ingoing null geodesic in the equatorial plane with $L = -6 E < 0$, starting at the point of Boyer-Lindquist coordinates $(t,r,\theta,\phi) = (0, 12, \pi/2, 0)$:
lmax = 13.07
Li = M.geodesic([0, lmax], M((0,12,pi/2,0)), mu=0, E=1, L=-6, Q=0,
r_increase=False, a_num=a0,
name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector:
The curve was correctly set. Parameters appearing in the differential system defining the curve are [a].
Li.integrate(step=0.00002)
graph = Li.plot(coordinates='xy', color='green', prange=[0, lmax], plot_points=40000,
thickness=1.5, display_tangent=True, scale=0.0001,
color_tangent='green', plot_points_tangent=4,
axes_labels=[r'$x/m$', r'$y/m$'])
graph
graph.save('gek_winding_null.pdf')
Li.check_integrals_of_motion(0.9999*lmax)
quantity | value | initial value | diff. | relative diff. |
- |
lmax = 26.41
Li = M.geodesic([0, lmax], M((0,15,pi/2,0)), mu=1, E=1, L=0, Q=0,
r_increase=False, a_num=a0,
name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector:
The curve was correctly set. Parameters appearing in the differential system defining the curve are [a].
Li.integrate(step=0.0002)
graph = Li.plot(coordinates='xy', color='red', prange=[0, lmax], plot_points=40000,
thickness=1.5, display_tangent=True, scale=0.0001,
color_tangent='red', plot_points_tangent=6,
axes_labels=[r'$x/m$', r'$y/m$'])
show(graph, ymin=-2, ymax=3)
graph.save('gek_frame_dragging.pdf')
Li.check_integrals_of_motion(0.9999*lmax)
quantity | value | initial value | diff. | relative diff. |
- |