In diesem Jupyter Notebook betrachteten wir, wie sich eine, nur gravitativ wechselwirkende Ansammlung von Teilchen (Staub), aufgrund der gegenseitigen gravitativen Anziehungskraft zusammenzieht, allmählich dichter wird um dann schließlich zu einem schwarzen Loch zu kollabieren. Die mathematisch analytische Betrachtung dieses Problems wurde erstmals im Jahre 1939 von J.R. Oppenheimer und H. Snyder in der Zeitschrift Physical Review publiziert (siehe Oppenheimer, J. Robert, and Hartland Snyder. "On continued gravitational contraction." Physical Review 56.5 (1939): 455.) und wird in vielen Lehrbüchern der Allgemeinen Relativitätstheorie unter Verwendung unterschiedlicher Formulierungen behandelt:
- Allgemeine Relativitätstheorie von Torsten Fließbach (Kapitel 47 und 48)
- Gravitation and Cosmology von Steven Weinberg (Seite 342 - 350)
- Gravitation von Charles W. Misner, Kip S. Thorne und John Archibald Wheeler (Kapitel 32)
- Relativistic hydrodynamics von Luciano Rezzolla und Olindo Zanotti (Kapitel 12.4)
Die zugrundeliegende allgemein-relativistische Gleichung, die es bei diesem Problem zu lösen gilt, ist die Einstein Gleichung (hier in kontra-kovarianter Darstellung)
$$ G^\mu{}\!_\nu = R^\mu{}\!_\nu - \frac{1}{2}g^\mu{}\!_\nu R = \frac{8\pi G}{c^4} T^\mu{}\!_\nu \quad,$$wobei hier der kosmologische Term $\Lambda$ vernachlässigt wurde, da es sich bei dem kollabierenden System um einen räumlich sehr kleinen Bereich der gesamten Raumzeit handelt. Zusätzlich folgt aus der Einsteingleichung die kovariante Erhaltung des Energie-Impulses. Diese sogenannten hydrodynamischen Gleichungen sind durch die folgenden vier Gleichungen definiert:
$$ \nabla\!_\mu G^\mu{}\!_\nu = 0 \quad \rightarrow \quad \nabla\!_\mu T^\mu{}\!_\nu = 0 \quad , $$wobei die kovariante Ableitung eines Tensors zweiter Stufe wie folgt definiert ist:
$$ \nabla\!_\alpha T^\mu{}\!_\nu = \partial_\alpha T^\mu{}\!_\nu + \Gamma^\mu_{\alpha \rho} T^\rho{}\!_\nu - \Gamma^\rho_{\alpha \nu} T^\mu{}\!_\rho \quad . $$Genau genommen, betrachten wir im Folgenden ein sphärisch symmetrisches System aus Staubmaterie und unterteilen die Raumzeit in einen inneren Bereich, der mit Staubmaterie gefüllt ($T^\mu{}\!_\nu \neq 0$) ist und einen äußeren Bereich ($T^\mu{}\!_\nu \equiv 0 $), der aufgrund des Birkoff Theorems durch die Schwarzschildmetrik gegeben ist. Zusätzlich nehmen wir an, dass die Materie im inneren Bereich der kollabierenden Kugel homogen verteilt ist.
Da sich zwangsläufig während des Kollapses der Ereignishorizont des entstehenden schwarzen Loches formt, ist es ratsam Koordinaten zu verwenden, die diese Koordinatensingularität am Horizont vermeiden. Eine Möglichkeit sind die "Gaußschen Koordinaten", bei denen die Zeitkoordinate die Eigenzeit $\tau$ ist. Für den inneren Bereich der Raumzeit verwenden wir die Friedmann-Lemaître-Robertson-Walker-Metrik (FLRW-Metrik),
$$ g_{\mu\nu}=\left( \begin{array}{ccc} c^2 & 0 & 0 & 0\\ 0& - \frac{R(\tau)^2}{1-k\, r^2} & 0&0 \\ 0& 0& -R(\tau)^2 r^2& 0\\ 0& 0& 0& -R(\tau)^2 r^2 \hbox{sin}^2(\theta)\\ \end{array} \right) \quad \forall \quad r \leq R(\tau) \quad , $$welche ursprünglich die homogene, isotrope zeitliche Entwicklung der Galaxiencluster des Universums beschrieben hatte. Die Funktion $R(\tau)$ bezeichnet den Radius der kollabierenden homogenen Kugel zur Eigenzeit $\tau$ und der Parameter $k$ kennzeichnet den Krümmungsparameter der Raumzeit, der in unserem Fall einen positiven Wert annimmt ($k > 0$). Die griechischen, raumzeitlichen Indices $\mu, \nu, \alpha, \beta$ laufen von 0..3 und entsprechen den folgenden sphärischen Koordinaten: $x^\mu = \left( x^0, x^1, x^2, x^3 \right) = \left( c \, \tau, r, \theta, \phi \right)$. $\tau$ hingegen stellt kein raumzeitlicher Index dar, sondern beschreibt die Eigenzeit des frei fallenden Staubteilchens.
Die Materie setzen wir als eine ideale Flüssigkeit an, die homogen verteilt ist und nur gravitativ wechselwirkt (Druck $p = 0$). Der Energie-Impuls Tensor besitzt im Innenraum somit die folgende einfache Struktur
$$ T^\mu{}\!_\nu=\left( \begin{array}{ccc} c^2 \, \varrho(\tau) & 0 & 0 & 0\\ 0& 0& 0&0 \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \end{array} \right) \quad , $$wobei die Funktion $\varrho(\tau)$ die Energiedichte der Materie beschreibt. Im Außenraum verschwindet der Energie-Impuls Tensor identisch.
Die raumzeitliche Struktur im Inneren der Materie erhält man mittels der Einstein Gleichung, die im betrachteten Fall ein System von vier gekoppelten Differentialgleichungen darstellt
$$ \begin{eqnarray} G^\tau{}\!_\tau &=& R^\tau\!_\tau - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\tau{}\!_\tau = \frac{8\pi G}{c^2} \varrho(\tau) \\ G^r\!_r &=& R^r\!_r - \frac{1}{2} R = \frac{8\pi G}{c^4} T^r\!_r = 0\\ G^\theta{}\!_\theta &=& R^\theta\!_\theta - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\theta{}\!_\theta = 0\\ G^\phi{}\!_\phi &=& R^\phi{}\!_\phi - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\phi{}\!_\phi = 0 \quad. \end{eqnarray} $$Wir berechnen uns im Folgenden die zugrundeliegenden Differentialgleichungen des kollabierenden Systems.
from sympy import *
init_printing()
from einsteinpy.symbolic import *
Wir definieren die FLRW-Metrik
$$ g_{\mu\nu}=\left( \begin{array}{ccc} c^2 & 0 & 0 & 0\\ 0& - \frac{R(\tau)^2}{1-k\, r^2} & 0&0 \\ 0& 0& -R(\tau)^2 r^2& 0\\ 0& 0& 0& -R(\tau)^2 r^2 \hbox{sin}^2(\theta)\\ \end{array} \right) \quad , $$wobei die Funktion $R(\tau)$ als eine "sympy"-Funktion definiert wurde.
tau, r, theta, phi, Pi, k, c, G = symbols('tau, r, theta, phi, pi, k, c, G')
R = Function('R')(tau)
Metric = diag(c**2, -R**2/(1-k*r**2), -R**2*r**2, -R**2*r**2*sin(theta)**2).tolist()
g = MetricTensor(Metric, [tau, 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)
print(chr.config)
simplify(chr.tensor()[1,1,1])
Die einzelnen Komponenten der Christoffel Symbole kann man anhand der folgenden Matrixform ablesen.
simplify(chr.tensor())
Die Christoffel Symbole in kovarianter Form (erster Art): $$ \Gamma_{\mu \nu \sigma} = \Gamma^{\rho}_{\mu \nu} g_{\rho \sigma} = \frac{1}{2} \left( g_{\nu \sigma| \mu} + g_{\mu \sigma| \nu} - g_{\mu \nu| \sigma} \right)$$
Hier speziell $$ \Gamma_{1 1 1} = \Gamma_{r r r}$$
chr_lll = chr.change_config('lll')
print(chr_lll.config)
chr_lll.tensor()[1,1,1]
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]
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.simplify()
Der Ricci Skalar ergibt sich aus der Kontraktion des Ricci Tensors: $R = R^{\mu}{}_{\!\mu} = g^{\mu \nu}R_{\nu \mu}$
ric_S = RicciScalar.from_riccitensor(ric)
ric_S.simplify()
Der Einstein Tensor:$$ G_{\mu \nu} = R_{\mu \nu} - \frac{1}{2}g_{\mu \nu}R $$
einst = EinsteinTensor.from_metric(g)
print(einst.config)
einst.simplify()
Die Materie setzen wir als eine ideale Flüssigkeit an, die homogen verteilt ist und nur gravitativ wechselwirkt (Druck $p = 0$). Der Energie-Impuls Tensor besitzt somit die folgende einfache Struktur
$$ T^\mu{}\!_\nu=\left( \begin{array}{ccc} c^2 \, \varrho(\tau) & 0 & 0 & 0\\ 0& 0& 0&0 \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \end{array} \right) \quad , $$wobei die Funktion $\varrho(\tau)$ die Energiedichte der Materie beschreibt.
rho = Function('varrho')(tau)
T_munu = diag(c**2 * rho, 0, 0, 0).tolist()
T = Tensor(T_munu,config="ul")
T.tensor()
Die Einstein Gleichung $$ G^\mu{}\!_\nu = R^\mu{}\!_\nu - \frac{1}{2}g^\mu{}\!_\nu R = \frac{8\pi G}{c^4} T^\mu{}\!_\nu$$
stellt demnach (in dem betrachteten Fall) ein System von vier gekoppelten Differentialgleichungen zweiter Ordnung dar
$$ \begin{eqnarray} G^\tau{}\!_\tau &=& R^\tau\!_\tau - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\tau{}\!_\tau = \frac{8\pi G}{c^2} \varrho(\tau) \\ G^r\!_r &=& R^r\!_r - \frac{1}{2} R = \frac{8\pi G}{c^4} T^r\!_r = 0\\ G^\theta{}\!_\theta &=& R^\theta\!_\theta - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\theta{}\!_\theta = 0\\ G^\phi{}\!_\phi &=& R^\phi{}\!_\phi - \frac{1}{2} R = \frac{8\pi G}{c^4} T^\phi{}\!_\phi = 0 \quad. \end{eqnarray} $$einst_ul = einst.change_config('ul')
EinstGl_0=Eq(einst_ul[0,0],8*Pi*G/c**4 * T[0,0])
EinstGl_0
EinstGl_1 = Eq(einst_ul[1,1],8*Pi*T[1,1])
EinstGl_1
EinstGl_2 = Eq(einst_ul[2,2],8*Pi*T[2,2])
EinstGl_2
EinstGl_3 = Eq(einst_ul[3,3],8*Pi*T[3,3])
EinstGl_3
Die letzten drei Gleichungen sind identisch ($G^r\!_r = G^\theta{}\!_\theta = G^\phi{}\!_\phi$).
Die kovariante Erhaltung des Energie-Impulses wird nun näher betrachtet. Diese sogenannten hydrodynamischen Gleichungen sind durch die folgenden vier Gleichungen definiert:
$$ \nabla\!_\mu G^\mu{}\!_\nu = 0 \quad \rightarrow \quad \nabla\!_\mu T^\mu{}\!_\nu = 0 \quad , $$wobei die kovariante Ableitung eines Tensors zweiter Stufe wie folgt definiert ist:
$$ \nabla\!_\alpha T^\mu{}\!_\nu = \partial_\alpha T^\mu{}\!_\nu + \Gamma^\mu_{\alpha \rho} T^\rho{}\!_\nu - \Gamma^\rho_{\alpha \nu} T^\mu{}\!_\rho \quad . $$und im Folgenden mittels Python analytisch berechnet wird.
Für die partielle Ableitung erhält man:
dT=GenericVector([diff(T[0,0],tau), diff(T[1,1],r), diff(T[2,2],theta), diff(T[3,3],phi)],[tau, r, theta, phi], config='l',parent_metric=Metric)
dT.tensor()
Der zweite Term $+ \Gamma^\mu_{\mu \rho} T^\rho{}\!_\nu$ berechnet sich zu:
tensorcontraction(tensorproduct(chr.tensor(),T.tensor()),(0, 1),(2, 3))
Der dritte Term $- \, \Gamma^\rho_{\mu \nu} T^\mu{}\!_\rho$ ist identisch Null:
-tensorcontraction(tensorproduct(chr.tensor(),T.tensor()),(0, 4),(1, 3))
Für die kovariante Ableitung erhält man schließlich:
$$ \nabla\!_\mu T^\mu{}\!_\nu = \partial_\mu T^\mu{}\!_\nu + \Gamma^\mu_{\mu \rho} T^\rho{}\!_\nu - \Gamma^\rho_{\mu \nu} T^\mu{}\!_\rho \quad , $$DT = dT.tensor() + tensorcontraction(tensorproduct(chr.tensor(),T.tensor()),(0, 1),(2, 3)) - tensorcontraction(tensorproduct(chr.tensor(),T.tensor()),(0, 4),(1, 3))
DT
wobei lediglich die Zeitkomponente des Vektors $\nabla\!_\mu T^\mu{}\!_0 = \nabla\!_\mu T^\mu{}\!_\tau$ von null verschieden ist und man folgende allgemeinrelativistische Gleichung der Massen/Energieerhaltung erhält:
Eq_E = Eq(simplify(DT[0]),0)
Eq_E
Diese einfache Differentialgleichung lässt sich in $\frac{d}{d\tau} \left( \varrho(\tau) \cdot R(\tau)^3 \right) = 0$ umformen
Eq_E1 = Eq(simplify(diff((rho*R**3),tau)),0)
Eq_E1
und somit gilt $\varrho(\tau) \cdot R(\tau)^3 = \rm const \quad \forall \, \tau \in \mathbb{R}$.
Diese Gleichung stellt die Massenerhaltung der kollabierenden Kugel dar und wir definieren die Konstante $\rm M_0 := \frac{4 \, \pi}{3} \cdot \varrho(\tau) \cdot R(\tau)^3 $ welche die gravitative Masse der Kugel darstellt. Wir substituieren nun $\varrho(\tau) = \frac{3 , \rm M_0}{4 \pi \, R(\tau)^3}$ in der $\tau\tau$-Komponente der Einsteingleichung.
M0 = symbols('M_0')
EinstGl_0a = simplify(EinstGl_0.subs(rho,3*M0/(4 * Pi * R**3)))
EinstGl_0a
Wir lösen die obere Gleichung nach $\left( \frac{d R(\tau)}{d\tau} \right)^2$ auf
EinstGl_0b=Eq(diff(R,tau)**2,solve(EinstGl_0a,diff(R,tau)**2)[0])
EinstGl_0b
und multiplizieren beide Seiten mit $\frac{M_0}{2}$.
EinstGl_0c = Eq(EinstGl_0b.lhs*M0/2,expand(EinstGl_0b.rhs*M0/2))
EinstGl_0c
Die obere Gleichung repräsentiert die Energieerhaltung des kollabierenden Systems:
$$\underbrace{\frac{M_0}{2} \left( \frac{d R(\tau)}{d\tau} \right)^2}_{\hbox{kinetische Energie}} - \underbrace{\frac{G M_0}{R(\tau)}}_{\hbox{potentielle Energie}} = - \frac{M_0 c^2 k}{2} = \hbox{const} \quad,$$wobei der erste Term die kinetische Energie der Massenverteilung darstellt und der zweite Term als potentielle Energie interpretiert werden kann.
Wir nehmen im Folgenden an, dass die kollabierende Kugel zur Zeit $\tau=0$ in Ruhe ist ( $\left. \frac{d R(\tau)}{d\tau} \right|_{\tau=0} = 0$ ) und einen Radius von $\rm R_0 := R(\tau=0)$ besitzt. Substituiert man diese Anfangswerte in den oberen Energiesatz, so erhält man den folgenden Ausdruck für den Krümmungsparameter $k$.
R0 = symbols('R_0')
Eq_k = Eq(k,solve(EinstGl_0c.subs(diff(R,tau),0).subs(R,R0),k)[0])
Eq_k
Der Energiesatz vereinfacht sich somit zur folgender Dfferentialgleichung erster Ordnung.
EinstGl_0d = Eq(diff(R,tau)**2,solve(EinstGl_0c.subs(k,Eq_k.rhs),diff(R,tau)**2)[0])
EinstGl_0d
Die Lösung dieser Differentialgleichung
$$ \left( \frac{d R(\tau)}{d\tau} \right)^2 = c^2 \,k \frac{R_0-R(\tau)}{R(\tau)} $$ist eine Zykloide mit der Parameterdarstellung (siehe Fließbach, Kapitel 47)
$$ \begin{eqnarray} \tau = \tau(\psi) &=& \frac{R_0}{2 c \sqrt{k}} \left( \psi + {\rm sin}(\psi) \right) \\ R = R(\psi) &=& \frac{R_0}{2} \left( 1 + {\rm cos}(\psi) \right) \quad , \end{eqnarray} $$wobei $\psi \in [0,\pi]$ den Parameter darstellt, der die Zykloide parametrisiert und $\psi = \pi $ den Endzustand des Kollapses repräsentiert, bei dem sich die gesamte Kugel auf einen singulären Punkt im Zentrum zusammengezogen hat.
Wir stellen uns im Folgenden die Lösung grafisch dar, wobei wir $c=1$ und $G=1$ setzen und die gravitative Masse der Kugel auf $M_0=5$ setzen. Der Anfangsradius $R_0$ sei dreimal so groß wie der theoretische Schwarzschildradius $R_S=2 M_0$ ($R_0 = 3 R_S = 6 M_0 = 30$); daraus folgt $\rightarrow k=\frac{1}{3}$.
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
params = {
'figure.figsize' : [12,8],
# 'text.usetex' : True,
'axes.titlesize' : 16,
'axes.labelsize' : 18,
'xtick.labelsize' : 16 ,
'ytick.labelsize' : 16
}
matplotlib.rcParams.update(params)
set_M0 = 5
set_R0 = 3*2*set_M0
set_c = 1
set_G = 1
set_k = float(Eq_k.rhs.subs(G,set_G).subs(c,set_c).subs(M0,set_M0).subs(R0,set_R0))
anz_pts=10000
psi_val = np.linspace(0, np.pi, anz_pts)
psi_val1 = np.linspace(np.pi, 1.5*np.pi, anz_pts)
x = set_R0/(2*set_c*np.sqrt(set_k))*(psi_val + np.sin(psi_val))
y = set_R0/2*(1 + np.cos(psi_val))
x1 = set_R0/(2*set_c*np.sqrt(set_k))*(psi_val1 + np.sin(psi_val1))
y1 = set_R0/2*(1 + np.cos(psi_val1))
plt.ylabel(r"$\rm R(\tau) $")
plt.xlabel(r"$\rm Eigenzeit \,\,\tau$")
plt.plot(x,y,c="blue", linewidth=1.5, linestyle='-');
plt.plot(x1,y1,c="black", linewidth=1.5, linestyle=':');
Die Entstehung der punktförmigen Singularität im Zentrum wird nach einer Eigenzeit
$$ T = \tau(\pi) = \frac{R_0}{2 c \sqrt{k}} \left( \pi + {\rm sin}(\pi) \right) = \frac{R_0 \, \pi}{2 c \sqrt{k}} = \frac{\pi}{2} \sqrt{\frac{R_0^3}{2 \, G \, M_0}} $$erreicht. In unserem Fall berechnet sich der Wert zu $T=81.621$, wobei die in der oberen Abbildung dargestellte gepunktete schwarze Kurve die unphysikalische Erweiterung der gefundenen Lösung mit $\psi > \pi$ darstellt.
T = float(set_R0/(2*set_c*np.sqrt(set_k))*(np.pi + np.sin(np.pi)))
T
Das entstehende Schwarze Loch besitzt zusätzlich eine Koordinatensingularität beim Schwarzschildradius $R_S=2 \, M_0$. Die Zeit $\tau_{S}$, die die kollabierende Kugel benötigt, um auf $R_S = R(\tau_{S})$ zu schrumpfen, berechnet sich zu
$$ R(\psi_S) = \frac{R_0}{2} \left( 1 + {\rm cos}(\psi_S) \right) = 2 \, M_0= R_S \quad \Rightarrow \quad \psi_{S} = \operatorname{acos}{\left(\frac{4 M_{0} - R_{0}}{R_{0}} \right)} $$und somit
$$ \tau_{S} = \tau_{S}(\psi_S) = \frac{R_{0} \left(\sqrt{1 - \frac{\left(4 M_{0} - R_{0}\right)^{2}}{R_{0}^{2}}} + \operatorname{acos}{\left(\frac{4 M_{0} - R_{0}}{R_{0}} \right)}\right)}{2 c \sqrt{k}} = \frac{\sqrt{1 - \frac{\left(4 M_{0} - R_{0}\right)^{2}}{R_{0}^{2}}} + \operatorname{acos}{\left(\frac{4 M_{0} - R_{0}}{R_{0}} \right)}}{\sqrt{\frac{8 G M_{0}}{R_{0}^3}}} $$In unserem Fall berechnet sich der Wert zu $\tau_{S}=74.135$,
psiS = symbols('psi_S')
tauS = symbols('tau_S')
Eq_S = Eq(R0/2*(1 + cos(psiS)) , 2*M0)
L_Eq_S = solve(Eq_S,psiS)[1]
Eq(psiS,L_Eq_S)
L_tauS = Eq(tauS,R0/(2*c*sqrt(k))*(L_Eq_S + sin(L_Eq_S)))
L_tauS
Eq(tauS,float(L_tauS.rhs.subs(R0,set_R0).subs(M0,set_M0).subs(k,set_k).subs(G,set_G).subs(c,set_c)))
Zweiter mathematischer Zugang zum Lösen der Dynamik des Kollapses¶
Wir wollen im Folgenden eine weitere Herleitung zeigen, die die $rr$-Komponente der Einsteingleichung verwendet. Die $rr$-Komponente der Einsteingleichung lautete:
EinstGl_1
Wir lösen hierzu die erste Einsteingleichung ($\tau\tau$-Komponente der Einsteingleichung) nach $\left( \frac{d R(\tau)}{d\tau} \right)$ auf
dR2 = solve(EinstGl_0a,(R.diff(tau))**2)[0]
Eq((R.diff(tau))**2,dR2)
und setzen den Ausdruck in die $rr$-Komponente der Einsteingleichung ein:
EinstGl_1a = EinstGl_1.subs((R.diff(tau))**2,dR2)
EinstGl_1a
und formen die Gleichung nach $\frac{d^2 R(\tau)}{d\tau^2}$ um:
DGL = Eq(R.diff(tau,tau),solve(EinstGl_1a,R.diff(tau,tau))[0])
DGL
Diese Differentialgleichung zweiter Ordnung beschreibt die zeitliche Entwicklung des Radius der kollabierenden Kugel. Man erkennt, dass die Differentialgleichung unabhängig vom Krümmungsparameter $k$ ist und dem klassischen Newtonschen Gravitationsgesetz entspricht. Im Folgenden werden wir diese DGL numerisch lösen, wobei wir die gravitative Masse der kollabierenden Kugel auf $\rm M_0 = 5$ festlegen und $G=1$ setzen.
from scipy import integrate
y1, y2 = symbols('y1, y2')
DGL_1 = DGL.subs(R,y1).subs(M0,set_M0).subs(G,set_G)
DGL_1
DGL_L=lambdify([tau,(y1,y2)], DGL_1.rhs)
Wir schreiben die Differentialgleichung zweiter Ordnung in ein System von zwei Differentialgeichungen erster Ordnung um:
def DGL_GravColl(y, tau):
y1, y2 = y
dy1 = y2
dy2 = DGL_L(tau,(y1,y2))
return np.array([dy1, dy2])
Als Anfangswerte setzen wir wieder den Radius der Kugel auf $\rm R_0 := R(\tau=0) = 100$ und nehmen an, dass sich die Kugel zum Anfangszeitpunkt in Ruhe befindet ( $\left. \frac{d R(\tau)}{d\tau} \right|_{\tau=0} = 0$ ).
pts = np.linspace(0, T-0.0003, anz_pts)
set_dR0=0.0
Randbedingungen = np.array([set_R0, set_dR0])
Loes_DGL = integrate.odeint(DGL_GravColl, Randbedingungen, pts)
Wir stellen uns die Lösung wieder grafisch dar, wobei wir zusätzlich auch $\frac{d R(\tau)}{d\tau}$ und $\varrho(\tau) = \frac{3 , \rm M_0}{4 \pi \, R(\tau)^3}$ visualisieren.
import matplotlib.gridspec as gridspec
params = {
'figure.figsize' : [16,5],
# 'text.usetex' : True,
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
}
matplotlib.rcParams.update(params)
fig = plt.figure(figsize=(16,5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.35)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax1.set_ylabel(r"$\rm R(\tau) $")
ax1.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax2.set_ylabel(r"$\rm dR/d\tau $")
ax2.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax3.set_ylabel(r"$\rm log(\varrho(\tau)) $")
ax3.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax1.plot(pts, Loes_DGL[:, 0], linewidth=1.5, c="blue")
ax2.plot(pts, Loes_DGL[:, 1], linewidth=1.5, c="blue")
ax3.plot(pts, np.log10(3*set_M0/(4*np.pi*Loes_DGL[:, 0]**3)), linewidth=1.5, c="blue")
ax1.set_ylim(0, set_R0+2)
ax2.set_ylim(-2.5, 0.1);
#ax3.set_ylim(-3, 4);
Aufgrund des Birkhoff-Theorems muss die Innenraum-Lösung der Metrik der kollabierenden Kugel in die äußere Schwarzschildmetrik am Sternrand stetig ineinander übergehen. Die kovariante Form der Raumzeit-Metrik der Schwarzschildmetrik lautet (hier mit $G=c=1$)
$$ g_{\mu\nu}=\left( \begin{array}{ccc} 1-\frac{2\,M}{\tilde{r}} & 0 & 0 & 0\\ 0& -\frac{1}{1-\frac{2\,M}{\tilde{r}}}& 0&0 \\ 0& 0& -\tilde{r}^2& 0\\ 0& 0& 0& -\tilde{r}^2 \hbox{sin}^2(\tilde{\theta})\end{array} \right) \quad \forall \quad \tilde{r} \geq \tilde{R}(t) \quad ,$$wobei $\tilde{R}(t)$ den Radius der Kugel in Schwarzschildkoordinaten repräsentiert.
Da diese Metrik nicht in Gaußscher Normalform ist und die Koordinaten $\tilde{x}^\mu = \left( c \, t, \tilde{r}, \tilde{\theta}, \tilde{\phi} \right)$ eine andere Bedeutung besitzen (Sichtweise eines im Unendlichen ruhenden Beobachters), muss man entweder die Innenraumlösung in diese Koordinaten überführen, oder die Schwarzschildmetrik in Gaußscher Normalform bringen. Beide Zugänge sind mathematisch anspruchsvoll und zeitintensiv und können in dieser Vorlesung nicht behandelt werden.
Um dennoch einen Eindruck zu erhalten, wie die kollabierende Kugel aus der Sichtweise eines im Unendlichen ruhenden Beobachters erscheint, nehmen wir an, dass auf der Oberfläche der kollabierenden Kugel Licht radial nach Außen emittiert wird und der im Unendlichen ruhende Beobachter dieses Licht beobachtet.
Zunächst zeigt ein Vergleich der Schwarzschildmetrik mit der FLRW-Metrik, dass $\tilde{\theta} = \theta$ und $\tilde{\phi} = \phi$ gilt. Zusätzlich gilt, dass der Flächeninhalt $\cal{A}$ der Oberfläche der Kugel von der verwendeten Metrik unabhängig sein muss. In der FLRW-Metrik gilt ${\cal{A}} = 4 \pi \, \left( R(\tau) \right)^2$ und in der Schwarzschildmetrik ${\cal \tilde{A}} = 4 \pi \, \left( \tilde{R}(t) \right)^2$, sodass $\tilde{R}(t) = R(\tau)$ gilt. Ferner ist $M$ die gravitative Masse, die das asymptotische Gravitationsfeld bestimmt, die wir im Folgenden mit dem gefundenen Ausdruck $\rm M_0 = \frac{4 \, \pi}{3} \cdot \varrho(\tau) \cdot R(\tau)^3 $ gleichsetzen ($M=M_0$).
Wir betrachten im Folgenden ein Staubteilchen, dass sich äußeren Rand der Kugel befindet und beschreiben seine geodätische Bewegung aus Sichtweise eines im Unendlichen ruhende Beobachters. Zur Zeit $t=\tau=0$ soll sich das Teilchen auf der äquatorialen Ebene der Kugel ($\tilde{\theta} = \theta = \frac{\pi}{2}, {\rm cos}(\theta)=0, {\rm sin}(\theta)=1$) befinden und wir berechnen seine geodätische Bewegung mittels der Geodätengleichung
$$ \frac{d^2 \tilde{x}^\mu}{d\lambda^2} + \Gamma^\mu_{\nu \rho} \frac{d \tilde{x}^\nu}{d\lambda} \frac{d \tilde{x}^\rho}{d\lambda} ~=~ 0 $$in vorgegebener Schwarzschild-Raumzeit, wobei $\lambda$ ein affiner Parameter ist, welcher die Bewegung des Staubteilchens parametrisiert, und $\Gamma^\mu_{\nu \rho}$ die Christoffel Symbole zweiter Art darstellen.
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 \tilde{x}^\nu}{d\lambda} \frac{d \tilde{x}^\rho}{d\lambda} \\ && \frac{d^2 \tilde{r}}{d\lambda^2} = - \Gamma^1_{\nu \rho} \frac{d \tilde{x}^\nu}{d\lambda} \frac{d \tilde{x}^\rho}{d\lambda}\\ && \frac{d^2 \tilde{\theta}}{d\lambda^2} = - \Gamma^2_{\nu \rho} \frac{d \tilde{x}^\nu}{d\lambda} \frac{d \tilde{x}^\rho}{d\lambda}\\ && \frac{d^2 \tilde{\phi}}{d\lambda^2} = - \Gamma^3_{\nu \rho} \frac{d \tilde{x}^\nu}{d\lambda} \frac{d \tilde{x}^\rho}{d\lambda} \quad . \end{eqnarray} $$Wir implementieren im Folgenden die Geodätengleichung in Python, wobei wir der Einfachheit halber die $\tilde{ }$ über den Koordinaten der Schwarzschildmetrik weglassen.
t, r, theta, phi, M = symbols('t, r, theta, phi, M')
Metric = diag(c**2*(1-2*G*M/(c**2*r)), -1/(1-2*G*M/(c**2*r)), -r**2, -r**2*sin(theta)**2).tolist()
g = MetricTensor(Metric, [t, r, theta, phi])
chr = ChristoffelSymbols.from_metric(g)
T_lambda = symbols('lambda')
x0 = Function('t')(T_lambda)
x1 = Function('r')(T_lambda)
x2 = Function('theta')(T_lambda)
x3 = Function('phi')(T_lambda)
xmu=GenericVector([x0, x1, x2, x3],[t, r, theta, phi], config='u',parent_metric=Metric)
Geod1=diff(xmu.tensor(),T_lambda,T_lambda)
Geod2T=tensorproduct(chr.tensor(),diff(xmu.tensor(),T_lambda),diff(xmu.tensor(),T_lambda))
Geod2=tensorcontraction(Geod2T, (1, 3),(2,4))
Geod=Geod1+Geod2
Es folgen die vier Gleichungen des Systems der Geodätengleichung:
GeodEq0=Eq(Geod[0],0)
GeodEq0
GeodEq1=Eq(Geod[1],0)
GeodEq1
GeodEq2=Eq(Geod[2],0)
GeodEq2
GeodEq3=Eq(Geod[3],0)
GeodEq3
Da die Bewegung des Teilchens sich radial und auf der äquatorialen Ebene abspielen wird, gilt $\frac{d \phi}{d\lambda}=0, \frac{d^2 \phi}{d\lambda^2}=0$ und $\frac{d \theta}{d\lambda}=0, \frac{d^2 \theta}{d\lambda^2}=0$. Dies hat zur Folge, dass die letzen beiden Gleichungen des DGL-Systems identisch verschwinden.
Um dieses System von zwei Differentialgleichungen zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir es zunächst in ein System von vier Differentialgleichungen erster Ordnung umschreiben. Wir definieren dazu formal: $y_1=\frac{dt}{d\lambda}$, $y_2=\frac{dr}{d\lambda}$, $y_3=\frac{dy_1}{d\lambda}=\frac{d^2t}{d\lambda^2}$ und $y_4=\frac{dy_2}{d\lambda}=\frac{d^2r}{d\lambda^2}$.
y1, y2, y3, y4 = symbols('y_1, y_2, y_3, y_4')
F0=solve(GeodEq0,diff(x0,T_lambda,T_lambda))[0]
F1=solve(GeodEq1.subs({(theta,pi/2),(diff(x2,T_lambda),0),(diff(x3,T_lambda),0)}),diff(x1,T_lambda,T_lambda))[0]
y3=lambdify((r,y1,y2), F0.subs({(diff(x0,T_lambda),y1),(diff(x1,T_lambda),y2),(M,set_M0)}).subs(G,set_G).subs(c,set_c))
y4=lambdify((r,y1,y2), F1.subs({(diff(x0,T_lambda),y1),(diff(x1,T_lambda),y2),(M,set_M0)}).subs(G,set_G).subs(c,set_c))
def DGLsys(vy, time):
t, r, y1, y2 = vy
dt = y1
dr = y2
dy1 = y3(r,y1,y2)
dy2 = y4(r,y1,y2)
return np.array([dt,dr,dy1,dy2])
Wir setzen die Parameterwerte auf die gleichen Werte, die wir bei der Kollapsrechnung festlegten ($M_0=5$, $c=1$ und $G=1$). Zur Zeit $t=\lambda=0$ befindet sich das Staubteilchen bei einem Radius von $r(0) = R_0 =30$ und Winkel $\phi_0=0$. Seine Anfangsgeschwindigkeit in radialer Richtung sei $\left. \frac{dr}{d\lambda} \right|_{\lambda=0}=0$. Die weitere nötige Anfangsbedingung ( $\left. \frac{dt}{d\lambda} \right|_{\tau=0}$ ) erhalten wir mittels der Bedingung $\left( \frac{ds}{d\lambda} \right)^2= g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda} = 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 Teilchens). Durch die Bedingung $\left( \frac{ds}{d\lambda} \right)^2 = 1$ wird gleichzeitig die Interpretation des affiner Parameters $\lambda$ festgelegt, der dann die Eigenzeit $\tau$ des fallenden Staubteilchens darstellt. Die Gleichung für das infinitesimale Weglängenelement $\left( \frac{ds}{d\lambda} \right)^2 = 1$ lautet (beachte $\frac{d\phi}{d\lambda} = \frac{d\theta}{d\lambda} =0$).
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),(dphi,0)}),1)
Eq_ds
Wir lösen diese Gleichung nach $dt$ auf und erhalten zwei Lösungen:
solve(Eq_ds,dt)
Zur Zeit $t=\lambda=0$ gilt $\left. \frac{dr}{d\lambda} \right|_{\lambda=0}=0$, sodass sich hier die Lösungen noch vereinfachen.
Anf_dt = solve(Eq_ds.subs(dr,0),dt)
Anf_dt
Einsetzen der Parameterwerte liefert für die positive, zweite Lösung:
set_Anf_dt = Anf_dt[1].subs(G,set_G).subs(c,set_c).subs(M,set_M0).subs(r,set_R0)
set_Anf_dt
Bei der nummerischen Lösung ist jetzt zu beachten, dass durch die äußere Sichtweise auf das Problem, die Schwarzschild-Singularität am Ereignishorizont ein Problem liefern wird. Wir können somit das System der Differentialgleichungen lediglich im Bereich $\lambda \in [0,\tau_{S}[ \, = \, [0,74.135[$ lösen.
anz_pts=100001
lambda_val = np.linspace(0, 74.13, anz_pts)
t0=0.0
r0=set_R0
dr0=0.0
dt0=float(set_Anf_dt)
initialval = np.array([t0,r0,dt0,dr0])
Loes = integrate.odeint(DGLsys, initialval, lambda_val)
params = {
'figure.figsize' : [12,8],
# 'text.usetex' : True,
'axes.titlesize' : 16,
'axes.labelsize' : 18,
'xtick.labelsize' : 16 ,
'ytick.labelsize' : 16
}
matplotlib.rcParams.update(params)
plt.cla()
x_max=160
plt.xlabel(r"$\rm Zeit $")
plt.ylabel(r"$\rm Radius \,\,R$")
plt.xlim(0,x_max)
plt.ylim(0,set_R0+1)
plt.plot(Loes[:, 0],Loes[:, 1],c="blue", linewidth=1.5, linestyle='-',label="Schwarzschild Zeit")
plt.plot(pts,Loes_DGL[:, 0],c="red", linewidth=1.5, linestyle='-',label="Eigenzeit")
plt.plot([0,x_max],[2*set_M0,2*set_M0],c="black", linewidth=1.5, linestyle=':')
plt.legend(loc='upper right', fontsize=18)
plt.text(20, 10, r"$R_S=2M$", size=20, rotation=0.,ha="center", va="center",bbox=dict(boxstyle="round",ec="black",fc="white"));
Das in das schwarze Loch fallende Staubteilchen 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. In der oberen Abbildung erkennt man deutlich, dass der Standpunkt des Beobachters eine bedeutende Rolle spielt. Während, aus der Sichtweise des kollabierende Körpers (bzw. aus der Sichtweise des Staubkorns auf seiner Oberfläche) der Ereignishorizont übertreten werden kann, friert hingegen die kollabierende Kugel, aus dem Betrachtungswinkel eines im Unendlichen ruhenden Beobachters (Schwarzschild Zeit) ein und das Staubkorn übertritt nie die Grenze bei $2 \, M$.
Die Entstehung des Ereignishorizontes kann leider in dieser Vorlesung nicht im Detail besprochen werden. Es sei nur erwähnt, dass dieser sich, ab einem gewissen Zeitpunkt im Zentrum des Sterns bildet, er sich dann immer weiter ausdehnt, bis schließlich die gesamte Kugel hinter dem Horizont ist. Wichtig bei der mathematischen Herleitung sind die sogenannten trapped surfaces und im Speziellen der outermost trapped surface (apparent horizon), welcher lokal für jede raumartige Hypersphäre definiert werden kann. Roger Penrose konnte dies im Jahre 1965 zeigen:
Penrose, Roger. "Gravitational collapse and space-time singularities." Physical Review Letters 14.3 (1965): 57. siehe auch Penrose, Roger. "“Golden Oldie”: gravitational collapse: the role of general relativity." General Relativity and Gravitation 34 (2002): 1141-1165.Enden möchte ich diese Vorlesung mit einem Zitat von Theodor W. Adorno, dem bedeutenden Philosophen der "Frankfurter Schule", das in dem besprochenen Kontext (Bedeutung des Beobachtungsstandpunktes in der Allgemeinen Relativitätstheorie) irgendwie passend ist:
Philosophie, wie sie im Angesicht der Verzweiflung einzig noch zu verantworten ist, wäre der Versuch, alle Dinge so zu betrachten, wie sie vom Standpunkt der Erlösung aus sich darstellten.