This Jupyter notebook illustrates some applications of SageMath functionalities in general relativity, specifically in deriving and displaying the horizons and ergosurfaces of the Kerr spacetime. Most of the involved tools are part of the SageManifolds project.
Rogerio T. Cavalcanti
It requires the SageMath version at least equal to 9.2.
version()
'SageMath version 9.4, Release Date: 2021-08-22'
%display latex
Parallelism().set(nproc=8)
a = var('a', domain='positive')
M.<t, r, th, ph> = manifolds.Kerr(m=1, a=a, coordinates='BL')
BL = M.default_chart()
g = M.metric()
g.display_comp()
singular_ring = {r:0, th:pi/2}
horizons = solve(1/g[1,1].expr()==0,r,solution_dict=True)
horizons
inner_horizon, outer_horizon = horizons
K = M.vector_field(1,0,0,0, name='K')
K.display()
Checking if $K$ is a Killing vector field $(\mathcal{L}_{_K}g = 0)$
g.lie_derivative(K) == 0
g(K,K).display()
ergosurfaces = solve(g(K,K).expr(),r,solution_dict=True)
ergosurfaces
inner_ergo, outer_ergo = ergosurfaces
List of surfaces
surfaces_param = [outer_ergo,outer_horizon,inner_horizon,inner_ergo,singular_ring]
In rational polinomial coordinates all components of the Kerr metric are rational polynomials, which in principle make it easyer to handle. We are going to use such coordinates for checking the Kretschmann scalar over the horizon and ergosurfaces of the spacetime.
RP.<t, r, ch, ph> = M.chart(r't:(-oo,+oo) r:(0,+oo) ch:(-1,1):\chi ph:(-pi,pi):periodic:\phi')
Transition map from Boyer–Lindquist coordinates to rational polynomial coordinates and its inverse.
BL_to_RP = BL.transition_map(RP, [t, r, cos(th), ph])
BL_to_RP.display()
BL_to_RP.set_inverse(t, r, acos(ch), ph)
BL_to_RP.inverse().display()
g.display_comp(RP.frame(),RP)
Setting the default chart and frame
M.set_default_chart(RP)
M.set_default_frame(RP.frame())
Riemann tensor, Ricci tensor and Kretschmann scalar
%time Riem = g.riemann()
CPU times: user 5.26 s, sys: 634 ms, total: 5.89 s Wall time: 36.6 s
%time Ric = g.ricci()
CPU times: user 4.43 s, sys: 65.2 ms, total: 4.49 s Wall time: 2.6 s
Ric == 0
%time R_up = Riem.up(g)
CPU times: user 1.54 s, sys: 236 ms, total: 1.78 s Wall time: 21.2 s
%time R_down = Riem.down(g)
CPU times: user 404 ms, sys: 92.4 ms, total: 497 ms Wall time: 3.52 s
%time Kretschmann_scalar = R_up['^{abcd}']*R_down['_{abcd}']
CPU times: user 49.8 s, sys: 891 ms, total: 50.7 s Wall time: 31.7 s
Kretschmann_scalar.display()
Getting and factoring the symbolic expression in the default chart
K_scalar = Kretschmann_scalar.expr().factor()
K_scalar
Singular Ring $(r=0,\chi=0)$
K_scalar.subs(r=0)
K_scalar.subs(ch=0)
Outer ergosurface, outer horizon and inner horizon at $\chi = 0$
for k in ['outer_ergo','outer_horizon','inner_horizon']:
print(k)
display(K_scalar.subs(eval(k)).subs({cos(th):ch}).subs(ch=0))
outer_ergo
outer_horizon
inner_horizon
Inner ergosurface for $\chi \neq 0$ (the inner ergosurface coincides with the singular ring at $\chi = 0$)
K_inner_ergo = K_scalar.subs(inner_ergo).subs({cos(th):ch}).canonicalize_radical()
K_inner_ergo
Series expansion up to $O(\chi^5)$
K_inner_ergo.series(ch,5)
Setting the default chart and frame back to Boyer–Lindquist
M.set_default_chart(BL)
M.set_default_frame(BL.frame())
The Kerr original coordinates will be used as an intermediate step for introducing the Kerr-Schild coordinates.
Kr.<u, r, th, vph> = M.chart(r'u:(-oo,+oo) r:(0,+oo) th:(0,pi):\theta vph:(-pi,pi):periodic:\varphi')
f(r) = r/(a^2+r^2-2*r)
assume(a<1)
F(r) = integral(f(r),r)
Kr_to_BL = Kr.transition_map(BL, [u-2*F(r), r, th, vph-a*F(r)])
Kr_to_BL.display()
Kr_to_BL.inverse().display()
Showing the change of frame from BL to Kerr
M.change_of_frame(BL.frame(),Kr.frame())[:]
KS.<u,x,y,z> = M.chart()
Change of coordinates from Kerr to Kerr-Schild
Kr_to_KS = Kr.transition_map(KS, [u, (r*cos(vph) - a*sin(vph))*sin(th),
(r*sin(vph) + a*cos(vph))*sin(th),
r*cos(th)])
Kr_to_KS.display()
Parametrization of the surfaces in Kerr-Schild coordinates
surfaces_KS = [vector([s.subs(param) for s in Kr_to_KS(u, r, th, ph)[1:]]) for param in surfaces_param]
# plotting data
surfs_data = {
'outer_ergo': {'name': 'Outer ergosurface',
'color': 'gray',
'z_label': 1,
'param_index': 0,
'phi1': 7*pi/5,
'plot_points_factor': 1},
'outer_orizon': {'name': 'Outer horizon',
'color': colormaps.Set1(1)[:3],
'z_label': .5,
'param_index': 1,
'phi1': 7*pi/5,
'plot_points_factor': 1},
'inner_horizon': {'name': 'Inner horizon',
'color': colormaps.Set1(4)[:3],
'z_label': 0,
'param_index': 2,
'phi1': 6*pi/5,
'plot_points_factor': .7},
'inner_ergo': {'name': 'Inner ergosurface',
'color': colormaps.Set1(3)[:3],
'z_label': -0.5,
'param_index': 3,
'phi1': 2*pi,
'plot_points_factor': .7}}
Python function for plotting the surfaces
def kerr_surfaces(surf, a=.99, print_labels=True, plot_points=30, mesh=True, **kwargs):
if a >1 or a < 0:
print("choose 'a' between 0 and 1")
return None
# Labels
if print_labels:
Ker_BH = text3d('Kerr black hole', (-2,-5,2.5), fontsize='200%', fontfamily='serif', fontweight='bold')
sep_line = text3d(r'___', (-2,-5,1.7), fontsize='160%', fontfamily='serif', fontweight='bold')
a_label = text3d('a = ' + str(a.n(digits=5)), (-2,-5,2), fontsize='170%')
s_ring_label = text3d('Singular ring', (-2,-5,-1), color='red', fontsize='170%', fontfamily='serif')
labels = sum([text3d(S['name'],
(-2,-5,S['z_label']),
color=S['color'],
fontsize='170%',
fontfamily='serif') for S in surfs_data.values()])
labels += Ker_BH + a_label + s_ring_label + sep_line
else: labels = Graphics()
# Surfaces
s_ring = parametric_plot3d(surf[4].subs(a=a),(ph,0,2*pi), color='red', thickness=4)
plots = sum([parametric_plot3d(surf[S['param_index']].subs(a=a),
(th,0,pi),
(ph,0,S['phi1']),
color=S['color'],
mesh=mesh,
plot_points=S['plot_points_factor']*plot_points,
frame=False,
**kwargs) for S in surfs_data.values()])
plots += s_ring
return plots+labels
kerr_surfaces(surfaces_KS, 0.999, viewpoint=[[-0.6557,-0.5284,-0.5394],112.41])
Dark theme
kerr_surfaces(surfaces_KS, .9, viewpoint=[[-0.6557,-0.5284,-0.5394],112.41], theme='dark')
We can also see the surfaces immersed in the Euclidian space $\mathbb{E}^3$.
E.<x,y,z> = EuclideanSpace(3)
spherical.<r, th, ph> = E.spherical_coordinates()
Differential map from Kerr coordinates to the Euclidean space
Kr_to_E = M.diff_map(E, {(Kr, spherical): [r,th,ph]},
name='Kr_to_E',
latex_name=r'\Phi_{_{\text{Kerr} \to \mathbb{E}^3}}')
Kr_to_E.display()
Coordinates in Euclidean space
coordinates = Kr_to_E(M.point((u,r,th,vph), chart=Kr)).coordinates()
Surfaces in Euclidean space
surfaces_E = [vector([s.subs(param) for s in coordinates]) for param in surfaces_param]
kerr_surfaces(surfaces_E, .96, viewpoint=[[-0.8499,-0.3478,-0.396],91.88])
We now create an animation by varying the parameter $a$.
theme = 'light' #'light' or 'dark'
frames1 = [kerr_surfaces(surfaces_KS,k, theme=theme,
viewpoint=[[-0.6557,-0.5284,-0.5394],112.41]) for k in srange(0,.5,.1)]
frames2 = [kerr_surfaces(surfaces_KS,k, theme=theme,
viewpoint=[[-0.6557,-0.5284,-0.5394],112.41]) for k in srange(.5,.95,.075)]
frames3 = [kerr_surfaces(surfaces_KS,k, theme=theme,
viewpoint=[[-0.6557,-0.5284,-0.5394],112.41]) for k in srange(.95,.9999,.005)]
frames4 = [kerr_surfaces(surfaces_KS,k, theme=theme,
viewpoint=[[-0.6557,-0.5284,-0.5394],112.41]) for k in srange(.9999,1,.000045)]
frames = frames1+frames2+frames3+frames4
animate(frames).interactive()