Allgemeine Relativitätstheorie mit dem Computer¶

General Theory of Relativity on the Computer¶

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main (Sommersemester 2021)¶

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske¶

Frankfurt am Main 04.04.2021¶

Erster Vorlesungsteil: Allgemeine Relativitätstheorie mit Python¶

Bewegung eines Probekörpers um ein rotierendes schwarzes Loch (Kerr Metrik)¶

Im Folgenden betrachten wir die Bewegung eines massiven Probekörpers um ein rotierendes schwarzes Loch und lösen die Geodätengleichung in vorgegebener Kerr-Raumzeit (in Boyer-Lindquist Koordinaten). Die kovarianten Kerr Raumzeit-Metrik eines rotierenden schwarzen Lochs der Masse M und Rotation a ($ a \in [-1,1]$ ist ein spezifischer Drehimpuls $a=J/M$ und wird als der sogenannte Kerr-Rotationsparameter bezeichnet) besitzt in Boyer-Lindquist Koordinaten folgendes Aussehen: $$ \begin{eqnarray} &g_{\mu\nu}=\left( \begin{array}{ccc} g_{tt}(r,\theta) & 0 & 0 & g_{t\phi}(r,\theta)\\ 0& g_{rr}(r,\theta)& 0&0 \\ 0& 0& g_{\theta\theta}(r,\theta)& 0\\ g_{\phi t}(r,\theta)& 0& 0& g_{\phi\phi}(r,\theta)\\ \end{array} \right)& \\ & g_{tt}(r,\theta)=\left( \frac{1-2\,M\,r}{\rho^2} \right)\,\,, \, g_{t\phi}(r,\theta)=\frac{2aMr\hbox{sin}^2(\theta)}{\rho^2} &\\ &g_{rr}(r,\theta)=-\frac{\rho^2}{\Delta}\,\,, \quad g_{\theta\theta}(r,\theta)=-\rho^2\,\,, \, & \\ &g_{\phi\phi}(r,\theta)=-\left( \frac{r^2+a^2+2 M r a^2 \hbox{sin}^2(\theta)}{\rho^2} \right)\hbox{sin}^2(\theta)\,\,,&\\ &\rho^2=r^2+a^2 \hbox{cos}^2(\theta) \,\,, \quad \Delta=r^2-2Mr+a^2& \end{eqnarray} $$

Wir definieren zunächst die kovarianten Raumzeit-Metrik eines rotierenden schwarzen Lochs der Masse M und Rotation a in Boyer-Lindquist Koordinaten (a ist der spezifische Drehimpuls ($a=J/M$) und wird als der sogenannte Kerr-Rotationsparameter bezeichnet):

from sympy import *
init_printing()
from einsteinpy.symbolic import *

t, r, theta, phi, M, a = symbols('t, r, theta, phi, M, a')
rho2=r**2+(a*cos(theta))**2
Delta=r**2-2*M*r+a**2
Metric = Matrix([[(1-2*M*r/rho2), 0, 0, (2*a*M*r*(sin(theta))**2)/rho2], [0, -rho2/Delta, 0, 0], [0, 0, -rho2, 0], [(2*a*M*r*(sin(theta))**2)/rho2, 0, 0, -(r**2+a**2+(2*M*r*a**2*(sin(theta))**2)/rho2)*(sin(theta))**2]]).tolist()
g = MetricTensor(Metric, [t, r, theta, phi])

g.tensor()

Rotierende schwarze Löcher besitzen, im Gegensatz zu nicht-rotierenden schwarzen Löchern, eine Ringsingularität und eine kompliziertere Struktur der Ereignishorizonte.

Ereignishorizonte ereignen sich formal an den Raumzeitpunkten, bei denen die radiale Komponente der zugrunde liegenden Metrik singulär wird ($g_{rr}\rightarrow\infty$ (bzw. $g^{rr}=0$).

EqHorizonte=Eq(g.inv()[1,1],0)
EqHorizonte

Für $a \neq 0$ erhält man zwei Lösungen, die man gewöhnlich mit den Symbolen $r_{+}$ und $r_{-}$ bezeichnet

$$ r_{+}= M + \sqrt{M^2 - a^2}\,, \quad r_{-}= M - \sqrt{M^2 - a^2}\quad, $$

und die für $a = 0$ in den Schwarzschild-Limes des Ereignishorizontes eines nicht-rotierenden schwarzen Lochs übergehen ($a=0 \rightarrow r_{+}=r_{-}=2\,M$).

LoesEqHorizonte=solve(EqHorizonte,r)
LoesEqHorizonte

Zusätzlich existieren die Flächen der stationären Grenze (stationary limit surfaces) und die der unendlichen Rotverschiebung; sie sind durch $g_{tt}=0$ bestimmt und werden gewöhnlich mit den Symbolen $r_{S^+}$ und $r_{S^-}$ bezeichnet: $$ r_{S^+}= M + \sqrt{M^2 - a^2{{\rm cos}^2(\theta)}}\,, \,\, r_{S^-}= M - \sqrt{M^2 - a^2{{\rm cos}^2(\theta)}} $$

EqUnRot=Eq(g.tensor()[0,0],0)
EqUnRot

LoesEqUnRot=solve(EqUnRot,r)
LoesEqUnRot

Die Horizontstruktur und die Flächen der stationären Grenze (unendliche Rotverschiebung) eines rotierenden schwarzen Lochs veranschaulichen wir in den folgenden Bildern (Kerr-Rotationparameter $a \in [0.5, \,0.95, \,0.99]$).

import numpy as np
import matplotlib.pyplot as plt 
import matplotlib
import matplotlib.gridspec as gridspec
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
params = {
    'figure.figsize'    : [16,5],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)

PlotHplus=lambdify((a), LoesEqHorizonte[1].subs({(M,1)}))
PlotHminus=lambdify((a), LoesEqHorizonte[0].subs({(M,1)}))
PlotUnRotplus=lambdify((a,theta), LoesEqUnRot[1].subs({(M,1)}))
PlotUnRotminus=lambdify((a,theta), LoesEqUnRot[0].subs({(M,1)}))

ntheta, nphi = np.linspace(0, 1.5 * np.pi, 400), np.linspace(0, np.pi, 400)
THETA, PHI = np.meshgrid(ntheta, nphi)

xylim=2.0
fig = plt.figure(figsize=(16,4.5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.25)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1],projection='3d')
ax3 = plt.subplot(gs[2],projection='3d')

seta=0.7
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax1.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax1.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax1.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax1.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

seta=0.95
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax2.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax2.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax2.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax2.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

seta=0.999
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax3.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax3.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax3.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax3.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

ax1.set_title(r"$\rm a=0.7$")
ax1.set_ylabel(r"$\rm x$")
ax1.set_xlabel(r"$\rm y$")
ax1.set_zlabel(r"$\rm z$")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(-xylim, xylim)
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.xaxis.pane.set_edgecolor('white')
ax1.yaxis.pane.set_edgecolor('white')
ax1.zaxis.pane.set_edgecolor('white')
ax1.view_init(azim=-55, elev=20)

ax2.set_title(r"$\rm a=0.95$")
ax2.set_ylabel(r"$\rm x$")
ax2.set_xlabel(r"$\rm y$")
ax2.set_zlabel(r"$\rm z$")
ax2.set_xlim(-xylim, xylim)
ax2.set_ylim(-xylim, xylim)
ax2.set_zlim(-xylim, xylim)
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.xaxis.pane.set_edgecolor('white')
ax2.yaxis.pane.set_edgecolor('white')
ax2.zaxis.pane.set_edgecolor('white')
ax2.view_init(azim=-55, elev=20)

ax3.set_title(r"$\rm a=0.99$")
ax3.set_ylabel(r"$\rm x$")
ax3.set_xlabel(r"$\rm y$")
ax3.set_zlabel(r"$\rm z$")
ax3.set_xlim(-xylim, xylim)
ax3.set_ylim(-xylim, xylim)
ax3.set_zlim(-xylim, xylim)
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.xaxis.pane.set_edgecolor('white')
ax3.yaxis.pane.set_edgecolor('white')
ax3.zaxis.pane.set_edgecolor('white')
ax3.view_init(azim=-55, elev=20)
#ax1.grid(False)

Rotierende schwarze Löcher besitzen eine kompliziertere Struktur ihrer Ereignishorizonte und die Flächen unendlicher Rotverschiebung und die Grenzflächen stationärer Bewegungen sind im Allgemeinen nicht identisch mit den Horizonten. Die oberen Abbildungen zeigen dies für drei unterschiedlich schnell rotierende schwarze Löcher (links $a=0.7$, Mitte $a=0.95$ und rechts $a=0.99$). Bei nichtverschwindender Rotation existieren zwei Ereignishorizonte (siehe obere Abbildung, r₊: graue Fläche und r_-: grüne Fläche) und ebenso existieren zwei Flächen unendlicher Rotverschiebung (r_{s₊</sub>: äußere rote Fläche und r_{s_-</sub>: blaue Fläche). Der Bereich zwischen der äußeren roten Fläche unendlicher Rotverschiebung und dem äußeren grauen Ereignishorizont bei r₊ nennt man Ergosphäre; hier können Testteilchen nicht mehr entgegengesetzt der Rotationsrichtung des schwarzen Lochs rotieren, stationäre Bewegungen sind somit nicht mehr möglich. Die unteren Abbildungen zeigen die gleichen Bilder bei einem anderen Viewpoint.}}

xylim=2.0
fig = plt.figure(figsize=(16,4.5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.25)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1],projection='3d')
ax3 = plt.subplot(gs[2],projection='3d')

seta=0.7
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax1.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax1.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax1.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax1.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

seta=0.95
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax2.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax2.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax2.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax2.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

seta=0.999
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)

R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)

R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)

R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

ax3.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
ax3.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
ax3.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
ax3.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)

ax1.set_title(r"$\rm a=0.7$")
ax1.set_ylabel(r"$\rm x$")
ax1.set_xlabel(r"$\rm y$")
ax1.set_zlabel(r"$\rm z$")
ax1.set_xlim(-xylim, xylim)
ax1.set_ylim(-xylim, xylim)
ax1.set_zlim(-xylim, xylim)
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.xaxis.pane.set_edgecolor('white')
ax1.yaxis.pane.set_edgecolor('white')
ax1.zaxis.pane.set_edgecolor('white')
ax1.view_init(azim=90, elev=0)

ax2.set_title(r"$\rm a=0.95$")
ax2.set_ylabel(r"$\rm x$")
ax2.set_xlabel(r"$\rm y$")
ax2.set_zlabel(r"$\rm z$")
ax2.set_xlim(-xylim, xylim)
ax2.set_ylim(-xylim, xylim)
ax2.set_zlim(-xylim, xylim)
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.xaxis.pane.set_edgecolor('white')
ax2.yaxis.pane.set_edgecolor('white')
ax2.zaxis.pane.set_edgecolor('white')
ax2.view_init(azim=90, elev=0)

ax3.set_title(r"$\rm a=0.99$")
ax3.set_ylabel(r"$\rm x$")
ax3.set_xlabel(r"$\rm y$")
ax3.set_zlabel(r"$\rm z$")
ax3.set_xlim(-xylim, xylim)
ax3.set_ylim(-xylim, xylim)
ax3.set_zlim(-xylim, xylim)
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.xaxis.pane.set_edgecolor('white')
ax3.yaxis.pane.set_edgecolor('white')
ax3.zaxis.pane.set_edgecolor('white')
ax3.view_init(azim=90, elev=0)
#ax1.grid(False)

Wir stellen diese Abbildungen in einer Animation dar:

import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from IPython.display import HTML

params = {
    'figure.figsize'    : [8,8],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)

anzframes=13
xylim=2.0
fig = plt.figure(figsize=(15,7))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1],projection='3d')
Lseta=[0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.98,0.99,0.999,0.9999]

def init():    
    seta=Lseta[0]
    settitle='a='+str(round(seta, 3))
    R = PlotHplus(seta)
    HplusX = R * np.sin(THETA) * np.cos(PHI)
    HplusY = R * np.sin(THETA) * np.sin(PHI) 
    HplusZ = R * np.cos(THETA)
    R = PlotHminus(seta)
    HminusX = R * np.sin(THETA) * np.cos(PHI)
    HminusY = R * np.sin(THETA) * np.sin(PHI) 
    HminusZ = R * np.cos(THETA)
    R = PlotUnRotplus(seta,THETA)
    UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
    UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
    UnRotplusZ = R * np.cos(THETA)
    R = PlotUnRotminus(seta,THETA)
    UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
    UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
    UnRotminusZ = R * np.cos(THETA)
    ax1.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
    ax1.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
    ax1.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
    ax1.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)
    ax2.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
    ax2.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
    ax2.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
    ax2.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)
    ax1.set_title(settitle)
    ax1.set_ylabel(r"$\rm x$")
    ax1.set_xlabel(r"$\rm y$")
    ax1.set_zlabel(r"$\rm z$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(-xylim, xylim)
    ax1.set_xticks([-2,-1,0,1,2])
    ax1.set_yticks([-2,-1,0,1,2])
    ax1.set_zticks([-2,-1,0,1,2])
    ax1.xaxis.pane.fill = False
    ax1.yaxis.pane.fill = False
    ax1.zaxis.pane.fill = False
    ax1.xaxis.pane.set_edgecolor('white')
    ax1.yaxis.pane.set_edgecolor('white')
    ax1.zaxis.pane.set_edgecolor('white')
    ax2.set_title(settitle)
    ax2.set_ylabel(r"$\rm x$")
    ax2.set_xlabel(r"$\rm y$")
    ax2.set_zlabel(r"$\rm z$")
    ax2.set_xlim(-xylim, xylim)
    ax2.set_ylim(-xylim, xylim)
    ax2.set_zlim(-xylim, xylim)
    ax2.set_xticks([-2,-1,0,1,2])
    ax2.set_yticks([-2,-1,0,1,2])
    ax2.set_zticks([-2,-1,0,1,2])
    ax2.xaxis.pane.fill = False
    ax2.yaxis.pane.fill = False
    ax2.zaxis.pane.fill = False
    ax2.xaxis.pane.set_edgecolor('white')
    ax2.yaxis.pane.set_edgecolor('white')
    ax2.zaxis.pane.set_edgecolor('white')
    ax1.view_init(azim=-55, elev=20)
    ax2.view_init(azim=90, elev=0)
    return fig,

def animate(i):
    ax1.cla()
    ax2.cla()
    seta=Lseta[i]
    R = PlotHplus(seta)
    HplusX = R * np.sin(THETA) * np.cos(PHI)
    HplusY = R * np.sin(THETA) * np.sin(PHI) 
    HplusZ = R * np.cos(THETA)
    R = PlotHminus(seta)
    HminusX = R * np.sin(THETA) * np.cos(PHI)
    HminusY = R * np.sin(THETA) * np.sin(PHI) 
    HminusZ = R * np.cos(THETA)
    R = PlotUnRotplus(seta,THETA)
    UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
    UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
    UnRotplusZ = R * np.cos(THETA)
    R = PlotUnRotminus(seta,THETA)
    UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
    UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
    UnRotminusZ = R * np.cos(THETA)
    ax1.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
    ax1.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
    ax1.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
    ax1.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)
    ax2.plot_surface(HminusX, HminusY, HminusZ, color="green", linewidth=0, alpha=0.5)
    ax2.plot_surface(HplusX, HplusY, HplusZ, color="grey", linewidth=0, alpha=0.5)
    ax2.plot_wireframe(UnRotplusX, UnRotplusY, UnRotplusZ, rstride=10, cstride=10, linewidth=0.5,antialiased=True,color="red", alpha=0.5)
    ax2.plot_wireframe(UnRotminusX, UnRotminusY, UnRotminusZ, rstride=10, cstride=10, linewidth=1,antialiased=True,color="blue", alpha=0.5)
    settitle='a='+str(round(seta, 3))
    ax1.set_title(settitle)
    ax1.set_ylabel(r"$\rm x$")
    ax1.set_xlabel(r"$\rm y$")
    ax1.set_zlabel(r"$\rm z$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(-xylim, xylim)
    ax1.set_xticks([-2,-1,0,1,2])
    ax1.set_yticks([-2,-1,0,1,2])
    ax1.set_zticks([-2,-1,0,1,2])
    ax1.xaxis.pane.fill = False
    ax1.yaxis.pane.fill = False
    ax1.zaxis.pane.fill = False
    ax1.xaxis.pane.set_edgecolor('white')
    ax1.yaxis.pane.set_edgecolor('white')
    ax1.zaxis.pane.set_edgecolor('white')
    ax2.set_title(settitle)
    ax2.set_ylabel(r"$\rm x$")
    ax2.set_xlabel(r"$\rm y$")
    ax2.set_zlabel(r"$\rm z$")
    ax2.set_xlim(-xylim, xylim)
    ax2.set_ylim(-xylim, xylim)
    ax2.set_zlim(-xylim, xylim)
    ax2.set_xticks([-2,-1,0,1,2])
    ax2.set_yticks([-2,-1,0,1,2])
    ax2.set_zticks([-2,-1,0,1,2])
    ax2.xaxis.pane.fill = False
    ax2.yaxis.pane.fill = False
    ax2.zaxis.pane.fill = False
    ax2.xaxis.pane.set_edgecolor('white')
    ax2.yaxis.pane.set_edgecolor('white')
    ax2.zaxis.pane.set_edgecolor('white')
    ax1.view_init(azim=-55, elev=20)
    ax2.view_init(azim=90, elev=0)
    return fig,

ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=600)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Die Visualisierung kann man auch wieder mittels der "Plotly Python Open Source Graphing Library" erstellen ( siehe https://plotly.com/python/ ).

import plotly.graph_objects as go

seta=0.98
R = PlotHplus(seta)
HplusX = R * np.sin(THETA) * np.cos(PHI)
HplusY = R * np.sin(THETA) * np.sin(PHI) 
HplusZ = R * np.cos(THETA)
R = PlotHminus(seta)
HminusX = R * np.sin(THETA) * np.cos(PHI)
HminusY = R * np.sin(THETA) * np.sin(PHI) 
HminusZ = R * np.cos(THETA)
R = PlotUnRotplus(seta,THETA)
UnRotplusX = R * np.sin(THETA) * np.cos(PHI)
UnRotplusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotplusZ = R * np.cos(THETA)
R = PlotUnRotminus(seta,THETA)
UnRotminusX = R * np.sin(THETA) * np.cos(PHI)
UnRotminusY = R * np.sin(THETA) * np.sin(PHI) 
UnRotminusZ = R * np.cos(THETA)

PHminus = go.Surface(x=HminusX, y=HminusY, z=HminusZ, colorscale="greens", showscale=False)
PHplus = go.Surface(x=HplusX, y=HplusY, z=HplusZ, colorscale="greys", showscale=False)
PUnRotplus = go.Surface(x=UnRotplusX, y=UnRotplusY, z=UnRotplusZ, colorscale="reds", showscale=False)
PUnRotminus = go.Surface(x=UnRotminusX, y=UnRotminusY, z=UnRotminusZ, colorscale="blues", showscale=False)

layout = go.Layout(autosize=True,
                  width=700, height=700,
                  margin=dict(l=15, r=50, b=65, t=90),
                  scene_aspectmode='cube',
                  scene = dict(
                    xaxis_showbackground=False,
                    yaxis_showbackground=False,
                    zaxis_showbackground=False,
                    xaxis_gridcolor="lightgrey",
                    yaxis_gridcolor="lightgrey",
                    zaxis_gridcolor="lightgrey",
                    xaxis_range=[-2,2],
                    yaxis_range=[-2,2],
                    zaxis_range=[-2,2],
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='z')
            )

data=[PHminus,PHplus,PUnRotplus,PUnRotminus]
fig=go.Figure(data=data, layout=layout)

camera = dict(eye=dict(x=1, y=-1.8, z=1.8))
fig.update_layout(scene_camera=camera)
fig.show()

Ein rotierendes schwarzes Loch zieht die Raumzeit mit sich mit. Die Rotationsfrequenz mit der die raumzeitlichen Struktur mitgeführt wird, nennt man ''Frame dragging'' Frequenz $\omega$; sie beschreibt wie viel sich die $\phi$-Koordinate pro Koordinatenzeit $t$ verändert: $$ \begin{equation} \omega(r,\theta)=\frac{d\phi}{dt} = \frac{\frac{d\phi}{d\tau}}{\frac{dt}{d\tau}} = \frac{u^\phi}{u^t} = \frac{g^{t\phi}}{g^{tt}} \end{equation} $$

FrameDrag=g.inv()[0,3]/g.inv()[0,0]
FrameDrag

Die folgende Abbildung zeigt die Frame Dragging Frequenz in der äquatorialen Ebene ($\theta=\pi/2$) für $M=1$, $a=0.95$ (green) und $a=0.2$ (rot).

FrameDragP=lambdify((a,r), FrameDrag.subs({(M,1),(theta,pi/2)}))
FrameDrag.subs({(M,1),(theta,pi/2)})

params = {
    'figure.figsize'    : [8,8],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)

rval = np.linspace(0.001, 10, 10001)
plt.cla()
plt.ylabel(r"$\rm Frame \, Dragging \,Frequenz \,\, \omega$")
plt.xlabel(r"$\rm Radius \,\, r $")
plt.plot(rval,FrameDragP(0.2,rval),c="red", linewidth=1.5, linestyle='-');
plt.plot(rval,FrameDragP(0.95,rval),c="green", linewidth=1.5, linestyle='-');
plt.xlim(1, 6)
plt.ylim(0, 0.6);

Dieser Mitführungs-Effekt der Raumzeit (Frame-Dragging-Effekt, bzw. Lense-Thirring-Effekt) ist eine generelle Eigenschaft von rotierenden Körpern in der allgemeinen Relativitätstheorie, die in der Newtonschen Gravitationstheorie nicht auftritt. Dieses Frame-Dragging der raumzeitlichen Struktur wurde bereits im Jahr 1918 von dem Mathematiker Josef Lense und dem Physiker Hans Thirring vorhergesagt.

Die obere Abbildung zeigt, dass der Frame-Dragging-Effekt in der Nähe eines schnell rotierenden schwarzen Loches sehr groß werden kann. Der Effekt ist jedoch auch bei langsam rotierenden Körpern existent und konnte im Jahre 2004 mittels der NASA Satelliten-Mission Gravity Probe B am Beispiel des durch die Erdrotation verursachten Mitführungs-Effekt der Raumzeit experimentell nachgewiesen werden (siehe auch Gravity Probe B: Testing Einsteins Universe).

Die Ergebnisse des Gravity Probe B Experiments wurden im Jahre 2011 veröffentlicht:

Everitt, CW Francis, et al. "Gravity probe B: final results of a space experiment to test general relativity." Physical Review Letters 106.22 (2011): 221101 (see also Clifford Will's Physics Viewpoint)

Geodätische Bewegung eines Probekörpers in vorgegebener Kerr-Raumzeit¶

Im Folgenden wird die Geodätengleichung in vorgegebener Kerr-Raumzeit (Boyer-Lindquist Koordinaten) betrachtet: $$ \frac{d^2 x^\mu}{d\tau^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau} ~=~ 0 $$ Die Geodätengleichung stellt ein System gekoppelter nichtlinearer Differentialgleichungen dar

$$ \begin{eqnarray} && \frac{d^2 t}{d\lambda^2} = - \Gamma^0_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \\ && \frac{d^2 r}{d\lambda^2} = - \Gamma^1_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \theta}{d\lambda^2} = - \Gamma^2_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \phi}{d\lambda^2} = - \Gamma^3_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \quad , \end{eqnarray} $$

wobei $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$), t, r, $\theta$ und $\phi$ die Koordinaten und $\Gamma^\mu_{\nu \rho}$ die Christoffel Symbole zweiter Art darstellen.

Die Chistoffel Symbole (zweiter Art): $$ \Gamma^{\sigma}_{\mu \nu} = \frac{1}{2}g^{\sigma \alpha} \left( g_{\nu \alpha| \mu} + g_{\mu \alpha| \nu} - g_{\mu \nu| \alpha} \right)$$

Hier speziell $$ \Gamma^{1}_{1 1} = \Gamma^{r}_{r r}$$

chr = ChristoffelSymbols.from_metric(g)
print(chr.config)
chr.tensor()[1,1,1]

ull

Wir definieren zunächst den ersten Term der Geodätengleichung $\frac{d^2 x^\mu}{d\tau^2}$:

tau = symbols('tau')
x0 = Function('t')(tau)
x1 = Function('r')(tau)
x2 = Function('\\theta')(tau)
x3 = Function('\\phi')(tau)

xmu=GenericVector([x0, x1, x2, x3],[t, r, theta, phi], config='u',parent_metric=Metric)
xmu.tensor()

Geod1=diff(xmu.tensor(),tau,tau)
Geod1

und nun den zweiten Term $\Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau}$:

Geod2T=tensorproduct(chr.tensor(),diff(xmu.tensor(),tau),diff(xmu.tensor(),tau))
Geod2=tensorcontraction(Geod2T, (1, 3),(2,4))
Geod2

Die Geodätengleichung ist somit ein System von nichtlinearen gekoppelten Differentialgleichungen zweiter Ordnung. Die entsprechenden vier Gleichungen lauten:

Geod=Geod1+Geod2
GeodEq0=Eq(Geod[0],0)
GeodEq0

GeodEq1=Eq(Geod[1],0)
GeodEq1

GeodEq2=Eq(Geod[2],0)
GeodEq2

GeodEq3=Eq(Geod[3],0)
GeodEq3

Beispiel eines radial in ein rotierendes schwarzes Loch einfallenden Probekörpers¶

Wir lassen im Folgenden nur ebene Bewegungen zu ($\theta=\frac{\pi}{2}, {\rm cos}(\theta)=0, {\rm sin}(\theta)=1, \frac{d \theta}{d\tau}=0, \frac{d^2 \theta}{d\tau^2}=0$). Unter diesen Annahmen sind nur drei der vier Geodätengleichungen relevant und diese lauten dann:

GeodEq0=GeodEq0.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq0

GeodEq1=GeodEq1.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq1

GeodEq3=GeodEq3.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq3

Um dieses System von drei Differentialgleichungen zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir es zunächst in ein System von sechs Differentialgleichungen erster Ordnung umschreiben. Wir definieren dazu formal: $y_1=\frac{dt}{d\tau}$, $y_2=\frac{dr}{d\tau}$, $y_3=\frac{d\phi}{d\tau}$, $y_4=\frac{dy_1}{d\tau}=\frac{d^2t}{d\tau^2}$, $y_5=\frac{dy_1}{d\tau}=\frac{dy_2}{d\tau}=\frac{d^2r}{d\tau^2}$ und $y_6=\frac{dy_3}{d\tau}=\frac{d^2\phi}{d\tau^2}$. Wir legen des weiteren die Masse des schwarzen Loches auf $M=1$ fest.

y1, y2, y3 = symbols('y_1, y_2, y_3')

F0=solve(GeodEq0,diff(x0,tau,tau))[0]
F1=solve(GeodEq1,diff(x1,tau,tau))[0]
F3=solve(GeodEq3,diff(x3,tau,tau))[0]

Eq(diff(x0,tau,tau),F0.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

Eq(diff(x1,tau,tau),F1.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

Eq(diff(x3,tau,tau),F3.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

Desweiteren setzen wir den Kerr-Rotationsparameter auf $a=0.95$.

seta=0.95

y4=lambdify((r,y1,y2,y3), F0.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1),(a,seta)}))

y5=lambdify((r,y1,y2,y3), F1.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1),(a,seta)}))

y6=lambdify((r,y1,y2,y3), F3.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1),(a,seta)}))

Wir definieren das System der gekoppelten sechs Differentialgleichungen:

def DGLsys(vy, time):
    t, r, phi, y1, y2, y3 = vy
    dt = y1
    dr = y2
    dphi = y3
    dy1 = y4(r,y1,y2,y3)
    dy2 = y5(r,y1,y2,y3)
    dy3 = y6(r,y1,y2,y3)
    return np.array([dt,dr,dphi,dy1,dy2,dy3])

import numpy as np
import matplotlib.pyplot as plt 
import matplotlib
from scipy import integrate

In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel einen radial in das rotierende schwarze Loch einfallenden Probekörper. Zur Zeit $t=\tau=0$ befindet sich der Probekörper bei einem Radius von $r_0=10$ und Winkel $\phi_0=0$. Seine Anfangsgeschwindigkeit in radialer Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ und seine Anfangsgeschwindigkeit in $\phi$Richtung sei auch $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$. Die weitere nötige Anfangsbedingung ( $\left. \frac{dt}{d\tau} \right|_{\tau=0}$ ) erhalten wir mittels der Bedingung $\left( \frac{ds}{d\tau} \right)^2= g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^\nu}{d\tau} = g_{\mu\nu} u^\mu u^\nu = 1$, die für massive Probekörper immer erfüllt sein muss ($u^\mu$ ist die Vierergeschwindigkeit des Probekörpers). Die Gleichung für das infinitesimale Weglängenelement $ds^2=1$ in der Ebene lautet:

dt, dr, dtheta, dphi = symbols('dt, dr, d\\theta, d\\phi')
dx=GenericVector([dt, dr, dtheta, dphi],[t, r, theta, phi], config='u',parent_metric=Metric)
dxdx=tensorproduct(dx.tensor(),dx.tensor())
dxdxT=Tensor(dxdx,config="uu")
ds2a=tensorproduct(g.tensor(),dxdxT.tensor())
ds2aT=Tensor(ds2a,config="lluu")
ds2=tensorcontraction(ds2aT.tensor(), (0, 2),(1, 3))
Eq_ds=Eq(ds2.subs({(theta,pi/2),(dtheta,0)}),1)
Eq_ds

Wir lösen diese Gleichung nach $dt$ auf:

LoesEq_ds=solve(Eq_ds,dt)
LoesEq_ds

Wir benutzen die zweite Lösung (positive Lösung) für unsere Anfangsbedingung und setzen den Kerr-Rotationsparameter $a=0.95 \quad \rightarrow$ $\left. \frac{dt}{d\tau} \right|_{\tau=0}=1.11803398874989$.

r0=10.0
t0=0.0
phi0=0.0
dr0=0.0
dphi0=0.0
dt0=float(LoesEq_ds[1].subs({(r,r0),(dr,dr0),(dphi,dphi0),(M,1),(a,seta)}))
dt0

numpoints=100001
tauval = np.linspace(0, 33.7975, numpoints)

r0=10.0
t0=0.0
phi0=0.0
dr0=0.0
dphi0=0.0
dt0=float(LoesEq_ds[1].subs({(r,r0),(dr,dr0),(dphi,dphi0),(M,1),(a,0.95)}))

initialval = np.array([t0,r0,phi0,dt0,dr0,dphi0])

Loes = integrate.odeint(DGLsys, initialval, tauval)

params = {
    'figure.figsize'    : [15,8],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)

endpoints=3
ax = plt.gca()
ax.cla() 
ax.set_ylabel(r"$\rm y \,[km]$")
ax.set_xlabel(r"$\rm x \,[km]$")
ax.plot(Loes[:numpoints-endpoints, 1]*np.cos(Loes[:numpoints-endpoints, 2]),Loes[:numpoints-endpoints, 1]*np.sin(Loes[:numpoints-endpoints, 2]),c="blue", linewidth=2.5, linestyle='-');
ax.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="red", s=60, marker='o');
PUnRotplus=plt.Circle((0, 0), PlotUnRotplus(0.95,np.pi/2), color='yellow',alpha=0.8)
ax.add_patch(PUnRotplus);
PHplus=plt.Circle((0, 0), PlotHplus(0.95), color='grey')
ax.add_patch(PHplus);
PHminus=plt.Circle((0, 0), PlotHminus(0.95), color='green')
ax.add_patch(PHminus);
plt.axis('equal');

Fällt ein Probekörper radial in ein rotierendes schwarzes Loch, so wird er aufgrund des Mitführungseffektes der rotierenden Raumzeit abgelenkt. Die obere Abbildung veranschaulicht diesen Effekt und stellt die Bewegung eines radial in ein rotierendes schwarzes Loch ($a=0.95$) einfallenden Probekörpers dar. Das rotierende schwarze Loch ist wie folgt visualisiert: Der Bereich der Ergosphäre des rotierenden schwarzen Loches (Bereich zwischen der äußeren Fläche unendlicher Rotverschiebung (r_{s₊</sub>) und dem äußeren Ereignishorizont (r₊)) ist gelb, der Bereich zwischen dem äußeren und inneren Ereignishorizont (Bereich zwischen r_- und r₊) ist grau und für den inneren grünen Bereich gilt r<r_-.}