Im Folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse Probekörper $m << M$ Masse schwarzes Loch) im Raum bewegt und sagt voraus, dass diese Bewegung sich stets entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht. Zunächst wird das Python-Modul "EinsteinPy" eingebunden, welches auf dem Modul SymPy basiert und symbolische Berechnungen in der Allgemeinen Relativitätstheorie relativ einfach möglich macht.
Definition der Koordinaten und der kovarianten Raumzeit-Metrik der Schwarzschildmetrik: $$ g_{\mu\nu}=\left( \begin{array}{ccc} 1-\frac{2\,M}{r} & 0 & 0 & 0\\ 0& -\frac{1}{1-\frac{2\,M}{r}}& 0&0 \\ 0& 0& -r^2& 0\\ 0& 0& 0& -r^2 \hbox{sin}^2(\theta)\end{array} \right)$$
from sympy import *
init_printing()
from einsteinpy.symbolic import *
t, r, theta, phi, M = symbols('t, r, theta, phi, M')
Metric = diag((1-2*M/r), -1/(1-2*M/r), -r**2, -r**2*sin(theta)**2).tolist()
g = MetricTensor(Metric, [t, r, theta, phi])
g.tensor()
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)
chr_llu = chr.change_config('llu')
print(chr_llu.config)
chr_llu.tensor()[1,1,1]
llu
Der Riemann Tensor: $$ R^{\mu}{}\!_{\nu \rho \sigma} = \frac{\partial \Gamma^{\mu}_{\nu \sigma}}{\partial x^{\rho}} - \frac{\partial \Gamma^{\mu}_{\nu \rho}}{\partial x^{\sigma}} + \Gamma^{\alpha}_{\nu \sigma}\Gamma^{\mu}_{\rho \alpha} - \Gamma^{\alpha}_{\nu \rho}\Gamma^{\mu}_{\sigma \alpha} $$
Hier speziell $$ R^{0}{}\!_{2 0 2} = R^{t}{}\!_{\theta t \theta} $$
rm = RiemannCurvatureTensor.from_christoffels(chr)
print(rm.config)
rm[0,2,0,2]
ulll
Der Ricci Tensor:
$$ R_{\mu \nu} = R^{\sigma}{}\!_{\mu \sigma \nu} = \frac{\partial \Gamma^{\sigma}_{\mu \nu}}{\partial x^{\sigma}} - \frac{\partial \Gamma^{\sigma}_{\mu \sigma}}{\partial x^{\nu}} + \Gamma^{\rho}_{\mu \nu}\Gamma^{\sigma}_{\rho \sigma} - \Gamma^{\sigma}_{\mu \rho}\Gamma^{\rho}_{\nu \sigma} $$
ric = RicciTensor.from_metric(g)
print(ric.config)
ric.tensor()
ll
Geodätische Bewegung eines Probekörpers¶
Die Geodätengleichung lässt sich durch folgendes Variationsprinzip herleiten
\begin{equation} \int_{x^\mu(A)}^{x^\mu(B)} ds= \int_{x^\mu(A)}^{x^\mu(B)} \sqrt{g_{\mu\nu} dx^\mu dx^\nu}=\int_A^B \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}} d\lambda \rightarrow \hbox{Extremal} \end{equation}
, wobei man, bei festgehaltenen Anfangs- und Endpunkten ($x^\mu(A)$ und $x^\mu(B)$), nach den extremalsten Kurven $x^\mu(\lambda)$ in der Raumzeit sucht. Den Parameter $\lambda$ bezeichnet man als affiner Parameter, welcher z.B. die Eigenzeit $\tau$ des Körpers sein kann. Die Geodätengleichung ergibt sich dann mittels der Euler-Lagrange Gleichungen, wobei man als Lagrangedichte den Ausdruck $L = \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}}$, bzw. alternativ $L = g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}$ verwendet.
\begin{equation} \frac{d}{d\lambda} \left( \frac{\partial L}{\partial \frac{d x^\alpha}{d \lambda}} \right) - \frac{\partial L}{\partial x^\alpha} = 0 \quad \xrightarrow{\cdot \frac{1}{2}g^{\mu \alpha }} \quad \frac{d^2 x^\mu}{d\lambda^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} ~=~ 0 \end{equation}
$\Gamma^\mu_{\nu \rho}$ sind die Christoffel Symbole zweiter Art und $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$ des Probekörpers). Die Geodätengleichung ist ein System gekoppelter Differenzialgleichungen, das wir uns im Folgenden näher betrachten werden.
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 Differenzialgleichungen 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
Radiales Einfallen eines Probekörpers in ein schwarzes Loch¶
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$) und beschränken uns zusätzlich auf ein radial in das schwarze Loch einfallenden Probekörper ($\frac{d \phi}{d\tau}=0, \frac{d^2 \phi}{d\tau^2}=0$). Unter diesen Annahmen sind nur die ersten beiden Geodätengleichungen relevant und diese lauten dann:
GeodEq0 = GeodEq0.subs({(theta,pi/2),(diff(x2,tau),0),(diff(x3,tau),0)})
GeodEq0
GeodEq1 = GeodEq1.subs({(theta,pi/2),(diff(x2,tau),0),(diff(x3,tau),0)})
GeodEq1
bzw. ...
GeodEq0_a = Eq(diff(x0,tau,tau),solve(GeodEq0,diff(x0,tau,tau))[0])
GeodEq0_a
GeodEq1_a = Eq(diff(x1,tau,tau),solve(GeodEq1,diff(x1,tau,tau))[0])
GeodEq1_a
Es ist manchmal hilfreich auch eine Latex-Ausgabe der DGL zu erzeugen:
print(latex(GeodEq0_a))
print(latex(GeodEq1_a))
\frac{d^{2}}{d \tau^{2}} t{\left(\tau \right)} = \frac{2 M \frac{d}{d \tau} r{\left(\tau \right)} \frac{d}{d \tau} t{\left(\tau \right)}}{r \left(2 M - r\right)}
\frac{d^{2}}{d \tau^{2}} r{\left(\tau \right)} = \frac{M \left(- 4 M^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + 4 M r \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + r^{2} \left(\frac{d}{d \tau} r{\left(\tau \right)}\right)^{2} - r^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2}\right)}{r^{3} \left(- 2 M + r\right)}
Um dieses System von zwei Differenzialgleichungen zweiter Ordnung mit Python numerisch lösen zu können, müssen wir es zunächst in ein System von vier Differenzialgleichungen erster Ordnung umschreiben. Details zum numerischen Lösung von Differenzialgleichungen finden Sie in der Vorlesung Einführung in die Programmierung für Studierende der Physik in den Unterpunkten Differenzialgleichungen: Numerische Lösung von Anfangswertproblemen und Systeme von gekoppelten Differenzialgleichungen und Differenzialgleichungen zweiter Ordnung.
Wir definieren dazu formal $$ t(\tau) = y_1(\tau) \, , \quad r(\tau) = y_2(\tau) \, , \quad \frac{d t(\tau)}{d\tau} = y_3(\tau) \, , \quad \frac{d r(\tau)}{d\tau} = y_4(\tau) $$
und schreiben unser DGL-System wie folgt um:
$$ \begin{eqnarray} \frac{d^2 t(\tau)}{d\tau^2} &=& \frac{d y_3(\tau)}{d\tau} = f_3(\tau,y_1,y_2,y_3,y_4) = \frac{2 M \frac{d}{d \tau} r{\left(\tau \right)} \frac{d}{d \tau} t{\left(\tau \right)}}{r \left(2 M - r\right)} \\ \frac{d^2 r(\tau)}{d\tau^2} &=& \frac{d y_4(\tau)}{d\tau} = f_4(\tau,y_1,y_2,y_3,y_4) = \\ &=& \frac{M \left(- 4 M^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + 4 M r \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2} + r^{2} \left(\frac{d}{d \tau} r{\left(\tau \right)}\right)^{2} - r^{2} \left(\frac{d}{d \tau} t{\left(\tau \right)}\right)^{2}\right)}{r^{3} \left(- 2 M + r\right)} \end{eqnarray} $$
y1, y2, y3, y4 = symbols('y_1, y_2, y_3, y_4')
SysDGL_0 = Eq(diff(x0,tau,tau),GeodEq0_a.rhs.subs({(r,y2),(diff(x0,tau),y3),(diff(x1,tau),y4)}))
SysDGL_0
SysDGL_1 = Eq(diff(x1,tau,tau),GeodEq1_a.rhs.subs({(r,y2),(diff(x0,tau),y3),(diff(x1,tau),y4)}))
SysDGL_1
Diese beiden DGLs werden wir später zusätzlich mittels der Programmiersprache C++ lösen. Wir exportieren deshalb schon jetzt das zugrundeliegende System von DGLs als C++ Code:
from sympy.codegen.ast import Assignment
wert = symbols('wert')
print(ccode(Assignment(wert,SysDGL_0.rhs)))
print(ccode(Assignment(wert,SysDGL_1.rhs)))
wert = 2*M*y_3*y_4/(y_2*(2*M - y_2)); wert = M*(-4*pow(M, 2)*pow(y_3, 2) + 4*M*y_2*pow(y_3, 2) - pow(y_2, 2)*pow(y_3, 2) + pow(y_2, 2)*pow(y_4, 2))/(pow(y_2, 3)*(-2*M + y_2));
Es ist manchmal hilfreich auch eine Latex-Ausgabe der DGL zu erzeugen:
print(latex(SysDGL_0))
print(latex(SysDGL_1))
\frac{d^{2}}{d \tau^{2}} t{\left(\tau \right)} = \frac{2 M y_{3} y_{4}}{y_{2} \left(2 M - y_{2}\right)}
\frac{d^{2}}{d \tau^{2}} r{\left(\tau \right)} = \frac{M \left(- 4 M^{2} y_{3}^{2} + 4 M y_{2} y_{3}^{2} - y_{2}^{2} y_{3}^{2} + y_{2}^{2} y_{4}^{2}\right)}{y_{2}^{3} \left(- 2 M + y_{2}\right)}
Nun lösen wir die DGLs mittels Python und setzen im Folgenden $M=1$:
f_3 = lambdify((y1,y2,y3,y4), SysDGL_0.rhs.subs(M,1))
f_4 = lambdify((y1,y2,y3,y4), SysDGL_1.rhs.subs(M,1))
Wir berechnen nun die numerische Lösung des oben beschriebenen Systems bestehend aus zwei gekoppelten DGLs.
Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. Zusätzlich berechnen wir die numerische Lösung auch mittels der älteren Methode "integrate.odeint()", die ebenfalls in dem Python-Modul "scipy" definiert ist.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy import integrate
Wir definieren das System der gekoppelten vier Differenzialgleichungen:
def DGLsys(time, vy):
y1,y2,y3,y4 = vy
dy1 = y3
dy2 = y4
dy3 = f_3(y1,y2,y3,y4)
dy4 = f_4(y1,y2,y3,y4)
return np.array([dy1,dy2,dy3,dy4])
Zusätzlich auch für die ältere Methode "integrate.odeint()"
def DGLsys_odeint(vy, time):
y1,y2,y3,y4 = vy
dy1 = y3
dy2 = y4
dy3 = f_3(y1,y2,y3,y4)
dy4 = f_4(y1,y2,y3,y4)
return np.array([dy1,dy2,dy3,dy4])
Zum Vergleich lösen wir auch die Bewegungsgleichung nach Newton:
$$\frac{d^2r}{dt^2}=-\frac{1}{r^2}$$
def DGLsysNewton(time, vy):
r, y2 = vy
dr = y2
dy2 = -1/r**2
return np.array([dr,dy2])
Zusätzlich auch für die ältere Methode "integrate.odeint()"
def DGLsysNewton_odeint(vy, time):
r, y2 = vy
dr = y2
dy2 = -1/r**2
return np.array([dr,dy2])
In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel die Anfangswerte eines Probekörpers, der sich zur Zeit $t=\tau=0$ bei einem Radius von $r_0=10=5 R_S$ befindet und eine Anfangsgeschwindigkeit in radialer Richtung von $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ hat. 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). Wir beschreiben die Bewegung aus der Sichtweise eines im Unendlichen ruhenden Beobachters (Schwarzschild Koordinaten).
N = 10000
tval = np.linspace(0, 33.7, N+1)
tvalNewton = np.linspace(0, 35.12, N+1)
t0 = 0.0
r0 = 10.0
dr0 = 0.0
dt0 = 1/np.sqrt(1-2/r0)
initialval = np.array([t0,r0,dt0,dr0])
initialvalNewton = np.array([r0,dr0])
Loes = integrate.solve_ivp(DGLsys, [tval[0], tval[-1]], initialval, t_eval=tval)
LoesNewton = integrate.solve_ivp(DGLsysNewton, [tvalNewton[0], tvalNewton[-1]], initialvalNewton, t_eval=tvalNewton)
Loes_odeint = integrate.odeint(DGLsys_odeint, initialval, tval)
LoesNewton_odeint = integrate.odeint(DGLsysNewton_odeint, initialvalNewton, tvalNewton)
rcParams.update({
'figure.figsize' : [10,8],
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
})
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,57)
plt.plot(Loes.y[1],Loes.y[0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
plt.plot(LoesNewton.y[0],tval,c="red", linewidth=1.5, linestyle='-',label="Newton");
plt.legend(loc='upper right', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 5, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Zusätzlich auch für die ältere Methode "integrate.odeint()"
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,57)
plt.plot(Loes_odeint[:, 1],Loes_odeint[:, 0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
plt.plot(LoesNewton_odeint[:, 0],tval,c="red", linewidth=1.5, linestyle='-',label="Newton");
plt.legend(loc='upper right', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 5, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Der in das schwarze Loch fallende Probekörpers bleibt somit, für einen äußeren Beobachter, anscheinend am Ereignishorizont bei $R_S=2M$ stehen (siehe die blaue Kurve im oberen Raumzeitdiagramm). Seine Bewegung in radialer Richtung wird nahe dem Ereignishorizont immer langsamer und das Bild, welches ein äußerer Beobachter wahrnimmt, friert am Ereignishorizont ein. Man kann zeigen, dass das wahrgenommene Bild gleichzeitig unendlich rotverschoben wird, sodass es nahe dem Ereignishorizont, langsam den wahrnehmbaren, sichtbaren Bereich verlässt.
Eine gute und zugleich lehrreiche Verbildlichung dieser und weiterer seltsamen Eigenschaften von schwarzen Löchern findet sich in dem Essay Black Holes and the German Reichstag (Schwarze Löcher und der Deutsche Reichstag) welcher 2005 im Magazin Physics World erschien (siehe Black Holes and the German Reichstag).
Die untere Abbildung veranschaulicht zusätzlich zur Lösung des radialen Falls des Probekörpers einige radial nach Außen gesendete Lichtstrahlen (gepunktete rote Kurven).
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,80)
plt.plot(Loes.y[1],Loes.y[0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
positions = [8, 6, 4, 3, 2.3, 2.1, 2.05, 2.01]
for pos in positions:
rval = np.linspace(pos, 10, 1000)
it = np.where(Loes.y[1] < pos)
PhotonOut = 2*np.log(np.abs(rval-2))+rval-(2*np.log(np.abs(pos-2))+pos)+Loes.y[0][it[0][0]]
plt.plot(rval, PhotonOut, c="red", linewidth=1.5, linestyle=":")
plt.legend(loc='lower center', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 15, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
Zusätzlich auch für die ältere Methode "integrate.odeint()"
plt.cla()
plt.ylabel(r"$\rm Koordinaten \, Zeit \,\,t$")
plt.xlabel(r"$\rm Radius \,\,r$")
plt.xlim(0,10.1)
plt.ylim(0,80)
plt.plot(Loes_odeint[:, 1],Loes_odeint[:, 0],c="blue", linewidth=1.5, linestyle='-',label="Allgemeine Relativitätstheorie");
positions = [8, 6, 4, 3, 2.3, 2.1, 2.05, 2.01]
for pos in positions:
rval = np.linspace(pos, 10, 1000)
it = np.where(Loes_odeint[:, 1] < pos)
PhotonOut = 2*np.log(np.abs(rval-2))+rval-(2*np.log(np.abs(pos-2))+pos)+Loes_odeint[it[0][0],0]
plt.plot(rval, PhotonOut, c="red", linewidth=1.5, linestyle=":")
plt.legend(loc='lower center', fontsize=18)
plt.axvline(x=2, linewidth=1.5, color="black", linestyle=':')
plt.text(2, 15, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",
bbox=dict(boxstyle="round",ec="black",fc="white"));
