Das Mehrkörperproblem in der Astronomie/ Hierarchische Algorithmen/ Quasikontinuierliche Massenverteilung (für Fortgeschrittene)

Konstruktion

Bei allen bisher vorgestellten Lösungsverfahren wurde ein Mehrkörpersystem als eine Ansammlung diskreter Massenpunkte betrachtet, wobei eine möglichst genaue Bestimmung der Bahnen der einzelnen Mitglieder zumeist im Vordergrund stand. Wie jedoch schon angedeutet, kommt es bei aus vielen Objekten bestehenden Ensembles wie Sternhaufen oder Galaxien aber meist gar nicht auf eine präzise Kenntnis ohnehin kaum prognostizierbarer individueller Orbits an, sondern nur auf eine möglichst korrekte Wiedergabe der räumlichen Verteilung der Körper und ihrer typischen Geschwindigkeiten. Zudem ist das Modell diskreter Massenpunkte für Systeme, die von dunkler Materie dominiert scheinen, wegen deren diffusen Verteilung oft nicht angemessen.

Als Alternative zu diskreten Körpern bietet es sich an, mit einer kontinuierlichen Massenverteilung zu arbeiten. Strenggenommen ist dies jedoch gar nicht möglich, da zu diesem Zweck die an einer bestimmten Stelle vorhandene Masse mit beliebig genauer räumlicher Auflösung ermittelt werden müsste. Tatsächlich kann eine Materieverteilung auch ohne Verwendung diskreter Massenpunkte lediglich mit einer gewissen Körnigkeit angegeben werden, so dass man von einer quasikontinuierlichen Verteilung spricht.

Um eine solche Verteilung darzustellen, teilt man das zu untersuchende System in kleine Würfel einer festen Kantenlänge $a$ ein. Für jeden dieser Kuben wird die darin enthaltene Masse und die mittlere Dichte mittels Division durch das Volumen $a^{3}$ bestimmt. Dieser Wert wird der Position der Würfelmitte zugeordnet, nicht dem Schwerpunkt der im Oktanten vorhandenen Einzelmassen. Auf diese Weise gewinnt man Dichten, die zwar weiterhin nur für diskrete Orte gelten. Diese weisen jedoch untereinander jeweils den gleichen Abstand voneinander auf.

Konstruktion einer quasikontinuierlichen Massenverteilung für ein System diskreter Massenpunkte

Im Gegensatz zum Barnes-Hut-Algorithmus werden keine Würfel unterschiedlicher Kantenlänge entsprechend einer Raumhierarchie betrachtet. Im Folgenden wird aber gezeigt, dass eine aus $N$ diskreten Körpern hervorgehende quasikontinuierliche Dichteverteilung ebenfalls eine Kräfteberechnung mit einem Aufwand proportional zu $N\log N$ erlaubt. Mathematisch ist dies im Vergleich zu den bislang vorgestellten Algorithmen allerdings viel anspruchsvoller, so dass sich dieses Unterkapitel vor allem an fortgeschrittene Leser richtet.

Potential einer kontinuierlichen Dichteverteilung

Während für ein System diskreter Massenpunkte die dort herrschenden Kräfte unmittelbar angegeben werden können, muss man für eine kontinuierliche Materieverteilung den Umweg über die potentielle Energie $E$ bzw. das Potential $\Phi$ nehmen, um zu einem Verfahren des Typs $N\log N$ zu gelangen. Unter dem Potential versteht man die pro Masse vorhandene potentielle Energie. Eine Punktmasse $m$ im Abstand $r$ von einer Punktmasse $M$ weist die potentielle Energie $E=-GmM/r$ auf, das Potential von $M$ lautet dementsprechend $\Phi =-GM/r$.

Um einen allgemeingültigen Zusammenhang zwischen Dichte und Potential zu gewinnen, muss man die von einer kleinen Masse $dM$ ausgehende Gravitationsbeschleunigung ${\vec {da}}$ sowie den durch diese bewirkten Fluss $d\Psi$ durch die Oberfläche $A$ des von $dM$ eingenommenen Volumens $V$ betrachten. Für die Beschleunigung gilt, wobei ${\vec {e}}_{r}$ den Einheitsvektor zwischen $dM$ und einem beliebigen Probekörper darstellt:

${\vec {da}}=-{\frac {GdM}{r^{2}}}{\vec {e}}_{r}$

Dieser Beschleunigung entspricht ein durch eine Kugelschale mit Radius $r$ hindurch tretender Fluss:

$d\Psi =\oint _{A}{\vec {da}}\cdot {\vec {dA}}=-\oint _{A}{\frac {GdM}{r^{2}}}dA$

Die letzte Umformung nutzt aus, dass die beiden Vektoren ${\vec {da}}$ und ${\vec {dA}}$ in die gleiche Richtung weisen. Das Integral über $dA$ liefert die Oberfläche der Kugelschale $4\pi r^{2}$ multipliziert mit dem Faktor $GdM/r^{2}$ und damit die Beziehung:

$d\Psi =-4\pi GdM$

$dM$ wiederum kann geschrieben werden als $\rho ({\vec {r}})dV$ und damit $d\Psi$ als $-4\pi G\rho ({\vec {r}})dV$. Daraus folgt schließlich für den eine beliebige Oberfläche eines beliebigen Volumens der Masse $M$ passierenden Fluss:

$\Psi =-4\pi G\int _{V}\rho ({\vec {r}})dV$

Nach dem Satz von Gauß kann das Flussintegral andererseits folgendermaßen umgeformt werden:

$\Psi =\oint _{A}{\vec {a}}\cdot {\vec {dA}}=\int _{V}\nabla \cdot {\vec {a}}dV$

Da über ein beliebiges Volumen integriert wird, müssen die für die Berechnung des Flusses benutzten Integranden gleich sein, so dass gilt:

$\nabla \cdot {\vec {a}}=-4\pi G\rho ({\vec {r}})$

Für ein konservatives Kraftfeld wie die Gravitation ist die von diesem ausgeübte Beschleunigung durch den Gradienten $-\nabla \Phi$ des Potentials gegeben. Dies liefert den gesuchten Zusammenhang zwischen Potential und Dichte, welcher als Poisson-Gleichung bezeichnet wird.

$\Delta \Phi =4\pi G\rho ({\vec {r}})$

In kartesischen Koordinaten lautet diese:

${\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial y^{2}}}+{\frac {\partial ^{2}\Phi }{\partial z^{2}}}=4\pi G\rho ({\vec {r}})$

Lösung der Poisson-Gleichung durch diskrete Fourier-Transformation

Eine Beziehung wie die Poisson-Gleichung, welche Ableitungen nach mehreren Variablen enthält, wird als partielle Differentialgleichung bezeichnet. Analog zur Beobachtung der zeitlichen Entwicklung eines Systems durch diskrete Zeitschritte wird eine solche Gleichung gelöst, indem die Differentiale durch diskrete Raumschritte ersetzt werden. Aus einer kontinuierlichen Materieverteilung wird somit eine durch ein regelmäßiges Würfelgitter beschriebene diskrete Dichteverteilung, wie sie bereits in der Einleitung dieses Unterkapitels beschrieben worden ist.

Mittels diskreter Raumschritte gewinnt man aus der Differentialgleichung für jeden Würfel des Gitters eine lineare Differenzengleichung, d. h. für $N$ Würfel ein System von $N$ derartigen Gleichungen. Im Prinzip kann ein solches durch klassische Verfahren wie das Gaußsche Eliminationsverfahren gelöst werden. Der Aufwand hierzu ist jedoch proportional zu $N^{2}$.

Wie nun erörtert wird, kann die Poisson-Gleichung auch auf Grundlage einer diskreten Fourier-Transformation gelöst werden. Für eine solche existieren sehr effiziente Algorithmen, welche als schnelle Fourier-Transformationen bezeichnet werden und wie gewünscht einen Aufwand proportional zu $N\log N$ aufweisen.

Eindimensional

Das Prinzip der Lösung mittels diskreter Fourier-Transformation soll zunächst für eine eindimensionale Dichteverteilung $\rho (x)$ erläutert werden. Aus den Würfeln werden dann Liniensegmente einer festen Länge $\Delta x$. Um den Laplace-Operator durch Differenzen wiederzugeben, muss man für eine gegebene Stützstelle $x$ zumindest deren nächste Nachbarn $x-\Delta x$ und $x+\Delta x$ betrachten. Diesen entsprechen Potentialdifferenzen $\Phi (x)-\Phi (x-\Delta x)$ und $\Phi (x+\Delta x)-\Phi (x)$ und damit 1.Ableitungen $(\Phi (x)-\Phi (x-\Delta x))/\Delta x$ und $(\Phi (x+\Delta x)-\Phi (x))/\Delta x$. Um die 2.Ableitung zu erhalten, muss man wiederum die Differenz der 1.Ableitungen bilden und diese durch $\Delta x$ dividieren. Dies liefert schließlich die Darstellung:

${\frac {\Phi (x+\Delta x)-2\Phi (x)+\Phi (x-\Delta x)}{(\Delta x)^{2}}}=4\pi G\rho (x)$

Da die Stützstellen untereinander den gleichen Abstand haben, genügt es für die folgende Diskussion, diese einfach mit einem laufenden Index $j$ durchzunummerieren. Mit $N$ Stützstellen, wobei im Folgenden $N$ gerade sein soll, lautet die Fouriertransformierte des Potentials:

${\tilde {\Phi }}_{n}={\frac {1}{N}}\sum _{j=-N/2}^{N/2-1}e^{-2\pi inj/N}\Phi _{j}$

Die Rücktransformation erfolgt durch die Vorschrift:

$\Phi _{j}=\sum _{n=-N/2}^{N/2-1}e^{2\pi inj/N}{\tilde {\Phi }}_{n}$

Die Beziehungen für die Dichte $\rho _{j}$ und ihre Fouriertransformierte ${\tilde {\rho }}_{n}$ sind völlig analog. Um die Verhältnisse an einer Stützstelle $j$ zu erfassen, werden noch die Potentiale $\Phi _{j+1}$ und $\Phi _{j-1}$ benötigt. Diese lassen sich sehr elegant mit Hilfe der Rücktransformation ausdrücken, da die Summe jeweils über $n$ läuft und ${\tilde {\Phi }}_{n}$ nicht von $j$ abhängt:

$\Phi _{j+1}=\sum _{n=-N/2}^{N/2-1}e^{2\pi in/N}e^{2\pi inj/N}{\tilde {\Phi }}_{n}$

$\Phi _{j-1}=\sum _{n=-N/2}^{N/2-1}e^{-2\pi in/N}e^{2\pi inj/N}{\tilde {\Phi }}_{n}$

Setzt man nun die Rücktransformationen in die Differenzengleichung ein, so steht auf beiden Seiten eine über $n$ laufende Summe. Damit die Gleichheit erfüllt wird, müssen jeweils die Summanden gleich sein. Transformiertes Potential und Dichte folgen damit der Beziehung:

${\frac {(e^{2\pi in/N}-2+e^{-2\pi in/N}){\tilde {\Phi }}_{n}}{(\Delta x)^{2}}}=4\pi G{\tilde {\rho }}_{n}$

Verwendet man noch den für die komplexe Exponentialfunktion gültigen Zusammenhang $e^{ix}=\sin x+i\cos x$ und löst nach dem Potential auf, so erhält man:

${\tilde {\Phi }}_{n}={\frac {(\Delta x)^{2}}{8\pi G}}{\frac {{\tilde {\rho }}_{n}}{\sin(2\pi n/N)-1}}$

Die Fouriertransformierte der Dichte liefert unmittelbar die Fouriertransformierte des Potentials, so dass kein Gleichungssystem gelöst werden muss.

Dreidimensional

Ausgehend vom eindimensionalen Fall ist das Vorgehen für das dreidimensionale Problem nun klar zu verstehen. Wie bereits erläutert, werden in der Praxis die diskreten Raumschritte $\Delta x$, $\Delta y$ und $\Delta z$ in Richtung der drei Koordinatenachsen gleich groß gewählt, das zu behandelnde System in Würfel konstanter Kantenlänge $a$ aufgeteilt. Die Differenzengleichung für solch einen Würfel lautet:

${\frac {\Phi (x+a,y,z)+\Phi (x-a,y,z)+\Phi (x,y+a,z)+\Phi (x,y-a,z)+\Phi (x,y,z+a)+\Phi (x,y,z-a)-6\Phi (x,y,z)}{a^{2}}}=4\pi G\rho (x,y,z)$

Die nun erforderliche dreidimensionale Fouriertransformation vereinfacht sich erheblich, wenn man davon ausgeht, dass in allen drei Koordinatenrichtungen die gleiche Anzahl $N$ von Stützstellen vorhanden ist, insgesamt also $N^{3}$ Würfel (man kann dies jederzeit sicherstellen, indem man gegebenenfalls die Verteilung mit Oktanten der Dichte 0 ergänzt). Um diese zu nummerieren, werden jetzt drei Indizes $j$, $k$ und $l$ benötigt. Die Transformation und Rücktransformation des Potentials (für die Dichte gelten gleichartige Formeln) lauten dann:

${\tilde {\Phi }}_{n,o,p}={\frac {1}{N}}\sum _{j=-N/2}^{N/2-1}\sum _{k=-N/2}^{N/2-1}\sum _{l=-N/2}^{N/2-1}e^{-2\pi i(nj+ok+pl)/N}\Phi _{j,k,l}$

$\Phi _{j,k,l}=\sum _{n=-N/2}^{N/2-1}\sum _{o=-N/2}^{N/2-1}\sum _{p=-N/2}^{N/2-1}e^{2\pi i(nj+ok+pl)/N}{\tilde {\Phi }}_{n,o,p}$

Für das Einsetzen in die Differenzengleichung sind zusätzlich die Größen $\Phi _{j+1,k,l}$, $\Phi _{j-1,k,l}$, $\Phi _{j,k+1,l}$, $\Phi _{j,k-1,l}$, $\Phi _{j,k,l+1}$ und $\Phi _{j,k,l-1}$ zu betrachten. Abermals gelingt dies mit Hilfe der Rücktransformation, da ${\tilde {\Phi }}_{n,o,p}$ nicht von $j$,$k$,$l$ abhängt. Das Resultat leuchtet unmittelbar ein:

${\frac {(e^{2\pi in/N}+e^{-2\pi in/N}+e^{2\pi io/N}+e^{-2\pi io/N}+e^{2\pi ip/N}+e^{-2\pi ip/N}-6){\tilde {\Phi }}_{n,o,p}}{a^{2}}}=4\pi G{\tilde {\rho }}_{n,o,p}$

${\tilde {\Phi }}_{n,o,p}={\frac {a^{2}}{8\pi G}}{\frac {{\tilde {\rho }}_{n,o,p}}{\sin(2\pi n/N)+\sin(2\pi o/N)+\sin(2\pi p/N)-3}}$

Beschleunigung durch eine kontinuierliche Dichteverteilung

Hat man das zu einer Dichteverteilung $\rho (x,y,z)$ dazugehörige Potential $\Phi (x,y,z)$ gefunden, so folgt daraus unmittelbar die Beschleunigung ${\vec {a}}(x,y,z)$. Wie schon bei der Herleitung der Poisson-Gleichung erwähnt, ist ${\vec {a}}$ durch den Gradienten von $\Phi$ gegeben;

${\vec {a}}=-\nabla \Phi$

In kartesischen Koordinaten lautet dieser:

${\begin{pmatrix}a_{x}\\a_{y}\\a_{z}\end{pmatrix}}=-{\begin{pmatrix}\partial \Phi /\partial x\\\partial \Phi /\partial y\\\partial \Phi /\partial z\end{pmatrix}}$

Ein weiteres Mal müssen die Ableitungen durch kleine Raumschritte ersetzt werden, wozu wieder für jede Stützstelle die nächsten Nachbarn herangezogen werden. Aus $\partial \Phi /\partial x$ wird so $(\Phi (x+a,y,z)-\Phi (x-a,y,z))/2a$, entsprechende Beziehungen gelten für die y- und z-Richtung:

${\begin{pmatrix}a_{x}\\a_{y}\\a_{z}\end{pmatrix}}=-{\frac {1}{2a}}{\begin{pmatrix}\Phi (x+a,y,z)-\Phi (x-a,y,z)\\\Phi (x,y+a,z)-\Phi (x,y-a,z)\\\Phi (x,y,z+a)-\Phi (x,y,z-a)\end{pmatrix}}$

Um die zeitliche Entwicklung des Systems zu verfolgen, muss man ein Stück weit wieder den Standpunkt diskreter Massenpunkte einnehmen. Jeder Körper wird jetzt als Probeteilchen betrachtet, das sich mit der aus der Poisson-Gleichung folgenden Beschleunigung innerhalb der quasikontinuierlichen Dichteverteilung bewegt. Mit jedem Zeitschritt nehmen die Massenpunkte neue Positionen ein, wodurch die daraus abgeleitete Dichteverteilung sich allmählich ändert. Die Aktualisierung der Positionen und Geschwindigkeiten kann mit den üblichen Methoden wie Hermite-Polynome, Leapfrog-Verfahren usw. erfolgen.

Die von der Poisson-Gleichung gelieferten Beschleunigungen sind nur an den Stützstellen gegeben, wohingegen die einzelnen Körper beliebige Orte einnehmen können. Um die dort herrschenden Beschleunigungen zu bestimmen, muss man zwischen den einem Massenpunkt nächstgelegenen Stützpunkten in jeder Achsen-Richtung interpolieren.

Für die Festlegung eines Zeitschritts $\Delta t$ bietet sich nun unmittelbar das Kriterium von Zemp und anderen (2007) an, das wie bereits geschildert $\Delta t$ mit der lokalen Dichte verknüpft:

$\Delta t=0.03{\frac {1}{\sqrt {G\rho }}}$

Die Einführung einer quasikontinuierlichen Dichteverteilung bringt automatisch eine erhebliche Glättung der Gravitation mit sich, so dass kein Plummerradius festgelegt werden muss, um unrealistische Beschleunigungsspitzen abzufangen. Andererseits können nur solche Strukturen korrekt erfasst werden, deren Größe zumindest derjenigen eines Würfels des Dichtegitters entspricht.