Das Mehrkörperproblem in der Astronomie/ Praktische Beispiele/ Planetensystem

Charakteristische Eigenschaften

Wie gerade bereits erwähnt besteht eine fundamentale Eigenschaft des Planetensystems darin, dass die Sonne fast ganz allein die auf einen Planeten einwirkende Schwerkraft bestimmt. So ist auf der Erde selbst die von Jupiter als massereichstem Planeten ausgehende Gravitation etwa 10000 Mal schwächer als diejenige der Sonne. Die Masse des Jupiters beträgt nur etwa 1/1000 Sonnenmasse, zudem ist selbst dessen minimale Entfernung etwa 4 Mal größer als unser Abstand zur Sonne. Die von den übrigen Planeten auf die Erde ausgeübten Kräfte sind noch weit geringer. Saturn weist weniger als 1/3 Jupitermasse auf und ist im Mittel etwa doppelt so weit wie jener von der Erde entfernt. Die Venus kommt unserem Planeten immer wieder recht nahe (bis auf etwa 30% der Entfernung zur Sonne), doch ihre Masse entspricht lediglich etwa 80% derjenigen der Erde selbst. Sowohl die von Saturn als auch Venus ausgehende Schwerkraft erreicht bei uns maximal etwa 1/100000 der Sonnengravitation - eine Größenordnung unter derjenigen des Jupiter.

Ein weiteres bedeutsames Charakteristikum des Sonnensystems sind die regelmäßigen Bewegungen der Planeten. Sie laufen auf nahezu perfekten Ellipsen um die Sonne und stehen relativ zu ihr in (fast) konstanten Zeitabständen immer wieder (fast) an den gleichen Positionen.

Im Gegensatz zu jedem anderen Mehrkörperensemble im Kosmos lassen sich für das Sonnensystem die für eine Simulation benötigten Startbedingungen mit hoher Präzision festlegen. Die Bahnen der Planeten sind für kurze Zeiträume sehr genau bekannt, ebenso ihre Massen. Astronomische Jahreskalender listen für jeden Planeten dessen Position am Himmel im Verlauf des Jahres auf, doch für eine Simulation sind stattdessen Orte und Geschwindigkeiten in kartesischen Koordinaten erforderlich. Ein solcher Satz einfach zu handhabender Anfangsbedingungen findet sich z.B. bei Ito und Tanikawa (2002) , welcher in folgender Tabelle zusammengestellt ist. Positionen sind in Astronomischen Einheiten, Geschwindigkeiten in km/s und Massen in millionstel Sonnenmassen angegeben. Die Sonne steht im Ursprung des Koordinatensystems. Pluto wurde erst einige Jahre nach der Veröffentlichung der Arbeit Itos und Tanikawas als Zwergplanet neu klassifiziert und ist daher in der Tabelle ebenfalls aufgenommen.

Planet	Masse	X	Y	Z	Vx	Vy	Vz
Merkur	0.166014	-0.126918	-0.447751	-0.031728	37.157	-10.930	-3.427
Venus	2.447838	-0.718938	-0.036831	0.021862	1.591	-35.126	-0.816
Erde	3.003490	-0.182443	0.966249	0.003394	-29.753	-5.628	-0.828
Mars	0.322715	1.391195	-0.005707	0.001958	1.003	26.295	0.771
Jupiter	954.791938	3.986089	2.960631	0.027793	-7.972	11.112	0.016
Saturn	285.885981	6.377297	6.606955	-0.145999	-7.471	6.694	0.039
Uranus	43.662517	14.509203	-13.654637	0.026676	4.611	4.653	0.117
Neptun	51.513890	16.941496	-24.901946	0.360724	4.452	3.102	-0.023
Pluto	0.007407	-9.869829	-28.059537	5.356149	5.289	-2.624	-1.117

Bahnelemente

Definition

Die Regelmäßigkeit der Planetenbahnen erlaubt es, diese durch einige wenige Größen zu beschreiben, anstatt mit einer großen Datenmenge von Positionen und Geschwindigkeiten in kartesischen Koordinaten hantieren zu müssen. Sie sind in folgender Zeichnung dargestellt.

Die Form einer Bahn wird durch die große Halbachse $a$ und die Exzentrizität $e$ gekennzeichnet. $a$ gibt die Ausdehnung des Orbits an und entspricht dem mittleren Abstand zur Sonne. $e$ charakterisiert die eigentliche Gestalt. Für $e$ = 0 ist die Bahn ein Kreis, für Werte zwischen 0 und 1 eine mit zunehmender Exzentrizität immer langgestreckter erscheinende Ellipse. Der minimale Abstand zur Sonne auf einer solchen beträgt $(1-e)a$, der maximale $(1+e)a$. Beläuft sich die Exzentrizität genau auf 1, ist die Bahn eine Parabel. Für $e$ größer als 1 liegt eine Hyperbel vor.

Die möglichen Bahnformen, die ein um eine stark dominierende Zentralmasse umlaufender Körper aufweisen kann, stehen mit dessen Energie in engem Zusammenhang. Bei Kreisen und Ellipsen ist seine kinetische Energie geringer als seine potentielle Energie, so dass er an die Zentralmasse gebunden ist. Im Fall einer Parabel halten die beiden Energien exakt einander die Waage, bei einer Hyperbel überwiegt die kinetische Energie. Die umlaufende Masse ist dann nicht an das Zentralgestirn gebunden. Sie passiert dieses ein einziges Mal, um danach wieder auf Nimmerwiedersehen zu verschwinden. Solche Orbits werden tatsächlich gelegentlich beobachtet, wenn ein Komet aus der Oortschen Wolke ins Innere des Sonnensystems gerät und dieses dann wieder verlässt.

Die Lage der Bahn im Raum wird durch die drei Winkel $i$, $\Omega$ und $\omega$ festfelegt. $i$ gibt die Neigung des Orbits gegenüber der Ekliptik an (die von der Sonne durch den Erdumlauf scheinbar am Fixsternhimmel beschriebene Bahn). $\Omega$ kennzeichnet den Winkel zwischen dem aufsteigenden Knoten (die Schnittgerade der Bahnebene mit der Ekliptik) und dem Frühlingspunkt $\gamma$, die sogenannte Knotenlänge. Schließlich beschreibt $\omega$ das Argument des Perihels P als Winkel zwischen dem aufsteigenden Knoten und der Periapsis (die Verbindungslinie zwischen Zentralmasse und dem Perihel - dem Punkt größter Annäherung an dieser).

Als letztes Bahnelement ist der Zeitpunkt $t$ zu erwähnen, zu welchem der umlaufende Körper den Perihel passiert.

Bestimmung aus kartesischen Koordinaten

Sind die Position ${\vec {r}}$ und Geschwindigkeit ${\vec {v}}$ einer umlaufenden Masse $m$ relativ zur Zentralmasse $M$ bekannt, so können daraus eindeutig die Bahnelemente abgeleitet werden. Das dazu erforderliche Vorgehen gehört zum klassischen Bestand der Himmelsmechanik und wird daher lediglich kurz skizziert. Zudem werden nur große Halbachse, Exzentrizität und Bahnneigung betrachtet, welche besonderes Augenmerk in Simulationen gefunden haben, sowie das Argument des Perihels.

Die große Halbachse folgt direkt aus den Beträgen von Orts- und Geschwindigkeitsvektor und den beiden Massen mittels der sogenannten Vis-Viva-Gleichung, welche wiederum auf der Energieerhaltung beruht:

$a={\frac {1}{{\frac {2}{r}}-{\frac {v^{2}}{G(M+m)}}}}$

Ein positiver Nenner entspricht einer Ellipse, ein negativer einer Hyperbel. Ist der Nenner gleich 0, liegt eine Parabelbahn vor.

Um mit Hilfe der großen Halbachse die Exzentrizität zu gewinnen, muss zunächst der sogenannte Bahnparameter $p$ ermittelt werden:

$p={\frac {l^{2}}{G(M+m)}}$

${\vec {l}}$ ist dabei der auf die umlaufende Masse normierte Bahndrehimpuls ${\vec {L}}$:

${\vec {l}}={\frac {\vec {L}}{m}}={\vec {r}}\times {\vec {v}}$

Je nachdem ob der Orbit eine Ellipse oder Hyperbel ist, gelten unterschiedliche Beziehungen für die Berechnung der Exzentrizität:

$e_{Ellipse}={\sqrt {1-{\frac {p}{a}}}}$

$e_{Hyperbel}={\sqrt {1+{\frac {p}{|a|}}}}$

Für die Bahnneigung wird abermals der normierte Bahndrehimpuls benötigt, und zwar die z-Komponente dessen Einheitsvektors:

$i=\arccos {\frac {l_{z}}{|{\vec {l}}|}}$

Um das Argument des Perihels zu bestimmen, wird auf eine ebenfalls normierte Darstellung des in Kapitel 4 eingeführten Runge-Lenz-Vektors ${\vec {A}}$ zurückgegriffen:

${\vec {\alpha }}={\frac {1}{m^{2}}}{\vec {A}}={\vec {v}}\times {\vec {l}}-G(M+m){\vec {e}}_{r}$

Schließlich muss man unterscheiden, ob die Bahnneigung gleich oder ungleich 0 ist. Im Allgemeinen wird die z-Komponente des Einheitsvektors von ${\vec {\alpha }}$ gebraucht:

$\omega =\arcsin {\frac {\alpha _{z}/|{\vec {\alpha }}|}{\sin {i}}}$

Für den Spezialfall einer nicht geneigten Bahn gehen die x- und y-Komponente von $\alpha$ ein:

$\omega =\arctan {\frac {\alpha _{y}}{\alpha _{x}}}$

Bahnstörungen

Grundprinzip der Berechnung

Für einen um die Sonne laufenden Planeten liegen aufgrund deren Dominanz nahezu die Bedingungen eines Zwei-Körper-Systems vor. In einem solchen sind die Bahnelemente konstant. Die von den übrigen Körpern des Sonnensystems auf einen Planeten ausgeübten Kräfte bewirken jedoch, dass diese Konstanz aufgehoben wird, die Bahnelemente im Laufe der Zeit sich allmählich ändern. Wie die nachfolgend diskutierten Beispiele zeigen, kann man diese Variationen in zwei Gruppen einteilen. Die Form der Bahn - vor allem die Exzentrizität - schwankt mehr oder weniger periodisch um einen Mittelwert, ebenso die Neigung des Orbits. Man spricht dann von (quasi)periodischen Störungen. Der Perihel entfernt sich dagegen mit jedem Umlauf immer weiter von seiner Ausgangsposition und wandert dabei im Laufe langer Zeiträume um die Sonne herum, eine oszillierende Bewegung um eine mittlere Position findet nicht statt. Eine solche Störung wird als säkular bezeichnet.

Die überwältigende Sonnengravitation erlaubt es, die Bahn eines Planeten auf Grundlage des sogenannten gestörten Zwei-Körper-Problems zu untersuchen, welches im letzten Abschnitt dieses Unterkapitels näher beleuchtet wird. Die auf einen solchen Trabanten ausgeübte Beschleunigung ${\vec {a}}$ wird als Summe zweier Komponenten dargestellt, die Sonnenbeschleunigung und eine viel schwächere Störbeschleunigung ${\vec {a}}_{S}$:

${\vec {a}}={\frac {GM}{|{\vec {r}}|^{3}}}{\vec {r}}+{\vec {a}}_{S}$

Die zeitliche Änderung sowohl von ${\vec {a}}$ als auch ${\vec {r}}$ ist durch die Bahnelemente eindeutig festgelegt. Indem man für diese wegen der Störbeschleunigung langsame Änderungen zulässt, kann man mittels verschiedener Verfahren komplizierte mathematische Formalismen entwickeln und so Gleichungen ableiten, welche die zeitliche Änderung der Bahnelemente beschreiben. Da mit geänderter Bahnform und -lage sich auch eine modifizierte Rechenvorschrift für die Störbeschleunigung ergibt, müssen auch diese Gleichungen durch Iteration gelöst werden. Nach jedem Simulationsschritt muss die Störung immer wieder aktualisiert und zur fortschreitenden Berechnung der Orbits neu vorgegeben werden.

Viele wissenschaftliche Publikationen über das Langzeitverhalten des Planetensystems stützen sich auf die Lösung des gestörten Zwei-Körper-Problems. Um die typischen Eigenschaften von Bahnstörungen zu demonstrieren, reichen die in diesem Buch vorgestellten allgemeinen Methoden wie etwa die Hermite-Polynome aber aus. Der Vorteil einer auf Störgleichungen beruhenden Lösung liegt vor allem darin, mit weit größeren zeitlichen Schrittweiten arbeiten zu können, wodurch Simulationen über sehr lange Zeiträume selbst in der Größenordnung von Milliarden Jahren möglich werden.

Für eine realistische Simulation mit einem Standardverfahren gelten zudem gewisse Einschränkungen. Um eine höchstmögliche Genauigkeit zu erreichen, verbietet es sich selbstverständlich, mehrere entfernt stehende Körper im Rahmen des Barnes-Hut-Algorithmus zu einem Objekt zusammenzufassen. Die Verwendung dynamischer Zeitschritte ist zulässig, d.h. die äußeren Planeten brauchen nicht so häufig bearbeitet zu werden wie die inneren. Man muss jedoch darauf achten, dass für das Zentralobjekt Größen wie die Beschleunigung oder der Ruck sehr gering und die daraus abgeleiteten Zeitschritte daher mit enormen Unsicherheiten behaftet sind. Um derartige Fehler zu vermeiden, wird hier empfohlen, für das Zentralobjekt die gleiche zeitliche Schrittweite anzusetzen wie für den innersten Planeten.

Alle nun dargestellten Beispiele für Bahnstörungen beruhen auf dem Hermite-Polynome-Verfahren mit hierarchischen Zeitschritten. Die Zeitskalen werden so gewählt, dass sie die Anforderung $\Delta t<0.02a/j$ erfüllen. Die relative Drift der Gesamtenergie liegt damit zwischen einem Zehn- und Hundertmilliardstel ihres Absolutbetrags pro Jahr.

Charakteristische Eigenschaften

Um die wesentlichen Eigenschaften von Bahnstörungen aufzuzeigen, wird ein stark vereinfachtes Sonnensystem vorgestellt, das zunächst nur aus Erde und Jupiter bestehen soll, wobei die Anfangsbedingungen obiger Tabelle entnommen werden. Die physikalischen Vorgänge in diesem Ensemble scheinen auf den ersten Blick sehr einfach zu sein.

Steht Jupiter der Erde am nächsten, wenn diese gerade durch den Perihel geht, so wird sie bei minimaler Entfernung zur Sonne besonders heftig vom Zentralgestirn weggezerrt. Dadurch verringert sich die Exzentrizität ihrer Bahn. Passiert die Erde zu diesem Zeitpunkt aber den Aphel, verhält es sich genau umgekehrt. Sie wird nun bei maximalem Abstand zur Sonne am stärksten von dieser fortgezogen, was eine ausgeprägtere Exzentrizität des Orbits zur Folge hat.

Befindet sich die Erde zum Zeitpunkt größter Annäherung an Jupiter besonders weit über der Ebene des Sonnensystems und jener gleichzeitig unterhalb dieser, so erfährt die Erde eine besonders starke Beschleunigung zur Ebene hin, wodurch sich die Neigung ihrer Bahn verringert. Steht in diesem Moment Jupiter aber noch höher über der Ebene als die Erde, erfolgt die Beschleunigung in umgekehrte Richtung, so dass der Erdorbit sich noch stärker neigt.

Die relativen Positionen von Erde und Jupiter zueinander ändern sich jedoch dauernd, so dass die soeben beschriebenen Szenarien einem ständigen Wechsel unterworfen sind. Hat die Erde während eines Periheldurchgangs sich Jupiter besonders genähert, so gilt dies - da Jupiter sich inzwischen auf seinem Orbt weiter bewegt hat - für die nächste Passage schon nicht mehr. Nur wenige Jahre später befindet sich die Erde während der größten Annäherung an Jupiter nahe (aber nicht exakt) des Aphels, einige Jahre danach wiederum in der Nähe (aber erneut nicht genau) des Perihels. Ebenso verhält es sich mit der Stellung von Erde und Jupiter relativ zur Ebene des Sonnensystems. Mal wird die Erde bei einer engen Begegnung mit Jupiter zu dieser hin-, mal wegbeschleunigt. Zudem ändert sich auch der minimale Abstand der beiden Planeten mit jeder Begegnung.

Die ständig wechselnden Konstellationen haben ein sehr komplexes Muster für die zeitlichen Schwankungen der Bahnelemente zur Folge, wie nachfolgend für die Exzentrizität gezeigt. Je nachdem wie Erde und Jupiter zueinander stehen, nimmt diese für einige Jahre zu und dann für einige Jahre ab. Die Änderungen von typischerweise 0.000001 bis 0.00001 pro Jahr scheinen auf den ersten Blick belanglos zu sein.

Jährliche Änderung der Exzentrizität der Erdbahn in Abhängigkeit von der Zeit für ein vereinfachtes, nur aus Erde und Jupiter bestehendes Planetensystem. Kurzfristig erkennt man einen raschen, mehr oder weniger periodischen Wechsel zwischen positiven und negativen Änderungen

Langfristig betrachtet weist die Exzentrizität jedoch ungeachtet dieses kleinskaligen Auf und Ab signifikante Trends auf. Über Jahrtausende hinweg betrachtet, lösen Epochen mit vorherrschend positiven und solche mit zumeist negativen Änderungen fortlaufend sich einander ab. Durch diese langanhaltende Dominanz jeweils eines bestimmten Trends unterliegt die Exzentrizität in der hier erörterten Simulation doch erhebliche Amplituden.

Jährliche Änderung der Exzentrizität der Erdbahn in Abhängigkeit von der Zeit für ein vereinfachtes, nur aus Erde und Jupiter bestehendes Planetensystem. Langfristig dominieren weitgehend periodisch einmal positive, einmal negative Änderungen

Exzentrizität der Erdbahn in Abhängigkeit von der Zeit für obiges Zwei-Planeten-Modell

Dieses Verhalten zeigt sich gleichermaßen für die Neigung der Erdbahn. Auch diese weist in dem vorliegenden Modell aufgrund eines wechselnden langfristigen Übergewichts jährlicher Zu- bzw. Abnahme deutliche Schwankungen auf. Anders verhält es sich hingegen mit der großen Halbachse, die im Verlauf der Simulation sich nur um wenige 0.00001 Astronomische Einheiten ändert.

Neigung der Erdbahn in Grad in Abhängigkeit von der Zeit für obiges Zwei-Planeten-Modell

Die Periheldrehung ist über lange Zeiträume betrachtet stets ein kumulativer Effekt. Das Hin- und Herpendeln des Arguments des Perihels in nachfolgendem Diagramm ist schlicht auf die zur Berechnung benutzten Arcussinus-Funktion zurückzuführen, welche nur Werte zwischen -90 und +90 Grad liefert. Hat das Argument einen Wert von +90 Grad erreicht, so läuft der Perihel anschließend keineswegs zurück, sondern behält seine Drehrichtung bei. Ebensowenig ist ein Wert von -90 Grad mit einer Umkehr der Drehrichtung verbunden. Die Zeit zwischen zwei aufeinanderfolgenden Werten von +90 bzw. -90 Grad entspricht also einem vollständigen Umlauf.

Argument des Perihels der Erdbahn in Grad in Abhängigkeit von der Zeit für obiges Zwei-Planeten-Modell

Während die Erde durch Jupiter erheblich gestört wird, bleibt jener angesichts der viel geringeren Erdmasse erwartungsgemäß fast unbehelligt. Die Exzentrizität seiner Bahn variiert nur um einige 0.0001, deren Neigung um einige 0.01 Grad. Während der Perihel der Erde während der Simulation etwa 5 komplette Umläufe um die Sonne absolviert, bewegt sich derjenige des Jupiter insgesamt nur um 3 Grad.

Mit der Hinzunahme der Venus zum Modell werden die Konstellationen der Planeten untereinander und dementsprechend die zeitlichen Schwankungen der Bahnelemente noch unübersichtlicher. Mit hoher zeitlicher Auflösung betrachtet, zeigt beispielsweise die jährliche Änderung der Exzentrizität sowohl für die Erd- als auch Venusbahn ein sehr kompliziertes Wechselspiel von positiven und negativen Werten. Dennoch schälen sich auch hier klar erkennbare Langzeitverläufe heraus, die im Vergleich zum Zwei-Planeten-Modell mit seinen sinusähnlichen Trends aber ebenfalls deutlich komplexer ausfallen.

Jährliche Änderung der Exzentrizität der Erd- und Venusbahn in Abhängigkeit von der Zeit für ein vereinfachtes, nur aus Venus, Erde und Jupiter bestehendes Planetensystem. Kurzfristig tritt jetzt ein noch komplizierterer, aber immer noch quasiperiodischen Wechsel zwischen positiven und negativen Änderungen auf

Exzentrizität der Erd- und Venusbahn in Abhängigkeit von der Zeit für obiges Drei-Planeten-Modell

Das langfristige Verhalten der Neigung von Erd- und Venusorbit in der Simulation lässt eine weitere typische Eigenschaft von Bahnstörungen erkennen. Ganz offensichtlich sind diese miteinander gekoppelt, d.h. es liegt eine Resonanz zwischen den beiden Körpern vor. Auf solche Kopplungen wird weiter unten noch genauer eingegangen.

Neigung der Erd- und Venusbahn in Grad in Abhängigkeit von der Zeit für obiges Drei-Planeten-Modell

Die große Halbachse erweist sich auch im Rahmen der Drei-Planeten-Simulation als äußerst robust, ebenso die Bahn des Jupiter.

Ergebnisse realistischer Langzeitsimulationen

Realistische Simulationen des Sonnensystems, die alle Planeten berücksichtigen und für die Erde zudem auch den Mond, zeigen angesichts der Vielzahl nun auftretender gegenseitiger Störungen noch kompliziertere Verhältnisse. Eine von Ito und Tanikawa (2002) auf Grundlage des gestörten Keplerproblems über 5 Milliarden Jahre durchgeführte Simulation lässt langfristig für die Exzentrizität nur bei Mars, Uranus und Pluto noch eine klare Periodizität erkennen. Alle übrigen Planetenbahnen weisen auch hinsichtlich der Langzeittrends eine mehr oder weniger unregelmäßig sich ändernde Exzentrizität auf. Die Schwankungen sind zudem weit stärker als in den obigen, nur aus zwei oder drei Planeten bestehenden Modellen. So finden die Autoren für die Erdbahn Variationen zwischen 0 und 0.07 auf einer Zeitskala von nur etwa 10000 Jahren. Auch die übrigen inneren Planeten sind von markanten Schwankungen der Exzentrizität betroffen. Für Mars geben Ito und Tanikawa Werte von 0 bis 0.15 an, für Merkur gar zwischen 0.01 und 0.25. Die Venus verhält sich weitgehend wie die Erde. Die Bahnen der äußeren Planeten erweisen sich als stabiler. Für Jupiter, Saturn und Uranus ermitteln die Autoren Exzentrizitäten von 0 bis 0.09. Besonders robust ist der Orbit des Neptun mit einer Exzentrizität von 0 bis 0.02. Die gegenwärtig sehr hohe Exzentrizität der Plutobahn ist kein Ausreißer. Die Langzeitsimulation liefert Resultate zwischen 0.21 und 0.27.

Wesentlich übersichtlicher verhalten sich gemäß Ito und Tanikawa langfristig die Neigungen der Planetenorbits. Hier liegt bei allen Planeten eine deutlich ausgeprägte Periodizität vor, auch wenn sich kurze und lange Zyklen überlagern. Abermals sind die Bahnen der äußeren Planeten stabiler als diejenigen der Inneren. Die ohnehin geringen Neigungen bei den vier großen Planeten von Jupiter bis Neptun variieren nur um wenige 0.1 Grad. Lediglich bei der auch heute stark geneigten Plutobahn beobachten die Autoren stärkere Schwankungen von 14 bis 17 Grad. Die Orbits der inneren Planeten unterliegen durchwegs Variationen von mehreren Grad. Am stärksten sind die Bahnen von Mars und Merkur mir Neigungen von 0 bis 7 bzw. 2 bis 10 Grad betroffen. Für Erde und Venus finden sich Bahnneigungen von 0 bis 3 Grad.

Vergleicht man die Arbeit von Ito und Tanikawa mit anderen Veröffentlichungen auf diesem Gebiet (die ebenfalls die Methode des gestörten Keplerproblems heranziehen), tun sich vor allem für den Orbit des Merkur eklatante Widersprüche auf. Eine Simulation von Boue und anderen (2012) liefert für den innersten Planeten ein geradezu dramatisches Szenario. Demgemäß soll die Exzentrizität auf einer Zeitskala von einigen Milliarden Jahren bis auf 0.8 ansteigen und so Merkur schließlich mit der Venus kollidieren können.

Zuvor hat schon Laskar (1994) auf eine mögliche Instabilität der Merkurbahn hingewiesen. Er fand bei Langzeitsimulationen nicht nur bis fast an 1 heranreichende Exzentrizitäten, sondern auch sehr starke Bahnneigungen bis über 30 Grad. Solch hohe Neigungen würden auch bei einem sehr exzentrischen Orbit eine Kollision des Merkur mit der Venus eventuell verhindern. Stattdessen könnte dem Autor zufolge die Exzentrizität möglicherweise sogar den kritischen Wert von 1 überschreiten und damit Merkur das Planetensystem gänzlich verlassen.

Für die Erdbahn ergeben sich bessere Übereinstimmungen zwischen verschiedenen Arbeiten. Laskar und andere (2011) fanden wie Ito und Tanikawa Exzentrizitäten zwischen 0 und 0.07 und ein komplexes Muster mehrerer sich überlagernder Zyklen sehr unterschiedlicher Dauer.

Schwankungen der Bahnelemente haben - zusammen mit solchen der Neigung der Rotationsachse - eine große Bedeutung für das Klima, was schon Milankowitch (1941) erkannt hat. Gegenwärtig passiert die Erde ihren Perihel im Januar, so dass auf der Nordhalbkugel die Winter aufgrund der Nähe zur Sonne etwas wärmer, die Sommer hingegen durch die größere Entfernung zur Sonne geringfügig kühler ausfallen. Kühle Sommer aber bewirken eine geringere Schnee- und Eisschmelze und begünstigen damit das Gletscherwachstum. Eine stärkere Exzentrizität der Bahn erhöht somit (bei der gegenwärtigen Jahreszeit des Periheldurchgangs) das Risiko einer Eiszeit in der nördlichen Hemisphäre. Eine mehr kreisförmige Bahn dagegen verringert die Gefahr einer Kaltzeit. Auf der Südhalbkugel verhält es sich genau umgekehrt.

Mit der gegenwärtigen Exzentrizität von etwa 0.017 beträgt im Vergleich zu einer Kreisbahn der Gewinn an Sonnenenergie während des nördlichen Winters etwa 3 Prozent, entsprechend groß ist der Verlust während des Nordsommers. Bei einer (maximalen) Exzentrizität von 0.07 aber lägen winterlicher Gewinn und sommerlicher Verlust auf der Nordhalbkugel jeweils bei circa 14 Prozent, was eine erhebliche Dämpfung der Jahreszeiten nach sich zöge. Auf der südlichen Hemisphäre hingegen würden bei unveränderter Jahreszeit der Perihelpassage die Temperaturunterschiede zwischen Sommer und Winter deutlich verstärkt.

Chaos und Resonanzen

Unvorhersehbare Bahnen

Die enorm voneinander abweichenden Aussagen über das Langzeitverhalten des Merkurs sind keineswegs auf Nachlässigkeiten der oben zitierten Autoren zurückzuführen, auch nicht auf im Detail unterschiedliche Simulationsverfahren. Laskar (1994) zeigte auf eindrucksvolle Weise, dass hier ein grundsätzliches Problem vorliegt. Wie für jede Aufgabenstellung, in deren Rahmen eine zeitliche Entwicklung prognostiziert werden soll, müssen auch für eine Mehrkörpersimulation Startwerte vorgegeben werden. Diese Anfangswerte - hier die Positionen, Geschwindigkeiten und Massen der beteiligten Körper - sind selbstverständlich nur mit endlicher Genauigkeit bekannt. Indem Laskar z.B. die Ausgangsposition der Erde nur um wenige Meter (!) verschob - weit über der Genauigkeit, mit der sich die Positionen der Planeten tatsächlich bestimmen lassen - erhielt er völlig andere Resultate hinsichtlich des Merkur. Auch die Bahnen der übrigen Planeten erwiesen sich als anfällig gegenüber selbst kleinen Modifikationen der Anfangsbedingungen, wenn auch in weit geringerem Maße. Die Bahnelemente der Erde etwa bleiben bei solchen Experimenten im Vergleich zu denjenigen des Merkur innerhalb recht enger Grenzen, der genaue zeitliche Verlauf aber ändert sich deutlich. Dieses Phänomen der äußert empfindlichen Abhängigkeit der zeitlichen Entwicklung eines dynamischen Systems von den Ausgangsbedingungen und seine daraus resultierende Unvorhersehbarkeit wird allgemein als Chaos bezeichnet und auch in vielen anderen Bereichen der Natur vorgefunden, z.B. im Klimasystem der Erde.

Das Aufschaukeln selbst kleiner Störungen zu völlig unterschiedlichen Orbits lässt sich anhand kleiner Körper wie Asteroiden ohne großen Aufwand vorführen. Aufgrund ihrer geringen Masse werden solche vor allem von den großen Planeten schon nach relativ kurzer Zeit erheblich gestört, so dass entsprechende Simulationen selbst auf einem einfachen Laptop nur wenige Minuten in Anspruch nehmen. Im folgenden Beispiel wird ein sehr einfaches System betrachtet, das nur Jupiter und einen Asteroiden mit einer Masse von 2 Milliarden Tonnen (etwa einem 1 km großen Felsbrocken entsprechend) enthält. Der Asteroid wird in zwei Simulationen zu Beginn jeweils auf die Position (4 AE / 0 / 0) gesetzt. Seine Geschwindigkeit beträgt anfangs einmal (0 / 15 km/s / 0), das andere Mal (0 / 15.001 km/s / 0), d.h. wird um lediglich 1 m/s variiert. Für Jupiter werden die Startwerte abermals obiger Tabelle entnommen.

Trotz des nur kleinen Unterschieds hinsichtlich der Anfangsbedingungen bilden sich schon nach wenigen Jahrhunderten gänzlich verschiedenen Orbits aus. Die zeitliche Entwicklung der Bahnelemente - als Beispiel soll hier die Exzentrizität genügen - unterscheidet sich für die beiden Fälle dramatisch. Zudem ist keinerlei periodisches Verhalten mehr erkennbar, die Resultate vielmehr von einem regellos erscheinenden Verlauf mit z.T. extremen Ausschlägen gekennzeichnet. Dabei nimmt die Asteroidenbahn vorübergehend die Gestalt einer Hyperbel an. Langfristig nähert sich die Exzentrizität bei beiden Simulationen dem Wert 1, d.h. es bildet sich eine extrem langgestreckte Ellipse heraus.

Exzentrizität der Asteroidenbahn in Abhängigkeit von der Zeit für ein vereinfachtes, nur aus einem Asteroiden und Jupiter bestehendes Planetensystem. Die Anfangspositionen der Körper sind für beide Simulationen gleich, die Anfangsgeschwindigkeiten des Asteroiden unterscheiden sich um 1 m/s

Angesichts des Trends zu einer immer exzentrischeren Bahn ist es nicht verwunderlich, dass der Asteroid sich z.T. weit von der Zentralmasse entfernt. In beiden Fällen treten Maximaldistanzen von über 80 Astronomischen Einheiten auf. Die trotz fast identischer Startbedingungen völlig unterschiedlichen Bewegungen fallen bei nachfolgender Darstellung der Positionen besonders prägnant ins Auge.

Vor allem der vergrößerte Ausschnitt der Orbits lässt erkennen, dass diese trotz extrem instabiler Bahnelemente eine gewisse Struktur bilden. Bestimmte Gebiete werden von dem Asteroiden offensichtlich häufig aufgesucht, andere hingegen fast völlig gemieden.

Asteroidenbahn für obiges Modell

Ausschnitt der Asteroidenbahn für obiges Modell

Ungeachtet des Simulationsbeispiels sind Asteroiden trotz ihrer geringen Masse aber nicht überall im Sonnensystem derart Spielball der Planeten, wie der Asteroidengürtel zwischen Mars und Jupiter belegt. Die Existenz relativ stabiler Refugien lässt sich anhand des hier skizzierten einfachen Modells ebenfalls gut demonstrieren. Gibt man dem Asteroiden die Startwerte (3 AE / 0 / 0) und (0 / 17 km/s / 0), so zeigen die Bahnelemente eine periodische Zeitabhängigkeit wie für das oben behandelte Zwei-Planeten-Modell. Auch gegenüber kleinen Änderungen des Anfangszustandes ist der Orbit jetzt viel robuster. Modifiziert man die Eingangsgeschwindigkeit erneut um 1 m/s, so sind für den hier überdeckten Zeitraum die Ergebnisse kaum zu unterscheiden. Allerdings bleibt die Bahn im Vergleich zu einem Planetenorbit deutlich stärker beeinträchtigt. So variiert die Exzentrizität zwischen 0.02 und 0.09, und der Perihel bewegt sich in der Simulation im Mittel in nur etwa 10000 Jahren einmal um die Sonne.

Resonanzen

Im Falle eines Asteroiden, welcher zunächst auf einer Kreisbahn mit einem Radius von 4 AE um die Sonne läuft, ist die Nähe der Jupiterbahn daran schuld, dass er sich nicht auf einem geordneten Orbit behaupten kann. Immer wieder kommt er an Jupiter bis auf nur 1 AE und schließlich sogar noch weniger heran, so dass jener eine erhebliche Störbeschleunigung ausübt und den Asteroiden in einer im Vergleich zum Alter des Sonnensystems kurzen Zeitspanne aus seiner ursprünglichen Bahn zieht.

Betrachtet man die Verteilung der Asteroiden in ihrem Gürtel gemäß Angaben des Jet Propulsion Laboratory (2007) , so erkennt man, dass erst ab Halbachsen unter 3.3 AE diese kleinen Körper signifikant gehäuft auftreten. In diesem Fall kommen sie Jupiter - sofern die Exzentrizität nahe bei 0 liegt - nur bis auf 1.7 AE nahe. Unterhalb einer Halbachse von 2.1 AE kommen erneut kaum noch Asteroiden vor, weil entsprechende Orbits der Bahn des Mars zu nahe sind.

Verteilung der Bahnhalbachsen der Asteroiden im Hauptgürtel. Die Pfeile markieren Entfernungen, in denen sich Objekte in einer Bahnresonanz mit Jupiter befänden, wobei die erste Ziffer die Zahl der Asteroidenumläufe angibt.

Innerhalb dieses Intervalls sind diese Körper aber ebenfalls nicht gleichmäßig verteilt. Neben Bereichen mit einer sehr großen Anzahl finden sich fast leere Zonen. Ein genauer Blick zeigt, dass die Umlaufszeiten in diesen Lücken immer in einem einfachen ganzzahligen Verhältnis zur Umlaufszeit des Jupiter stehen. Man spricht dann von einer störenden Resonanz.

Verhalten sich die Umlaufszeiten von Asteroid und Jupiter wie m:n, so stehen sich die beiden nach jeweils m Umläufen des kleinen Körpers bzw. n Umläufen des großen Planeten immer wieder an der gleichen Stelle auf ihren Bahnen gegenüber. Die von Jupiter ausgehenden Störungen greifen dann wiederholt den Orbit des Asteroiden stets auf die gleiche Weise an und addieren sich so im Laufe der Zeit effektiv immer mehr auf. Auch bei einem größeren Minimalabstand zu Jupiter reicht dessen Störkraft damit aus, den Asteroiden aus seiner anfänglichen Bahn zu werfen.

Je nach den Werten von m und n können auch mehrere Orte des Asteroidenorbits immer wieder der gleichen Störung durch Jupiter unterliegen. Ist mit (A / 0) und (J / 0) die größte Annäherung zwischen Asteroid und Jupiter gegeben und das Verhältnis der Umlaufszeiten 2:1, so tritt die stärkste Störung stets nach einem Umlauf des Planeten an eben diesen Positionen auf. Bei einem Verhältnis 3:1 erscheint (unter Annahme von Kreisbahnen) nach einem halben Umlauf des Jupiter an den Stellen (-A / 0) und (-J / 0) ein weiteres Abstandsminimum.

Je nach Lage dieser Minima können die dort auftretenden Störungen sich aber auch gegenseitig kompensieren und damit im Gegenteil eine besonders stabile Bahn bedingen. So findet sich bei einem Umlaufverhältnis von 7:2 mit Jupiter nicht ein Minimum, sondern stattdessen das stärkste Häufigkeitsmaximum des Asteroidengürtels überhaupt. Es liegt dann eine stabile Resonanz vor.

Bereits seit langem ist bekannt, dass ein Asteroid sich stabil im gleichen Orbit wie ein Planet aufhalten kann, wenn der Winkelabstand (die Schenkel sind jeweils durch die Verbindungslinie zur Zentralmasse definiert) der beiden untereinander genau 60 Grad beträgt. Es handelt sich hier um einen der wenigen exakt lösbaren Spezialfälle des Dreikörperproblems, die dadurch beschriebenen stabilen Positionen werden Lagrange-Punkte genannt. Für fast alle Planeten wurden an den entsprechenden Stellen mittlerweile Asteroiden gefunden, welche als Trojaner bezeichnet werden. Aufgrund seiner großen Masse besitzt Jupiter die meisten derartigen Begleiter auf seiner Bahn.

Resonanzen treten nicht nur zwischen Planeten und Asteroiden auf, sondern auch zwischen den Planeten selbst. Als Beispiel sei nochmals das aus Venus, Erde und Jupiter bestehende Modell herangezogen. Wie oben schon skizziert, verhalten sich die Neigungen der Erd- und Venusbahn in dieser Simulation exakt gegenläufig zueinander. Nimmt die Neigung des eines Orbits zu, so geht diejenige des anderen zurück und umgekehrt. Weniger klar ist das Bild für die Exzentrizitäten, doch hinterlässt auch dieses den Eindruck einer wechselseitigen Abhängigkeit. Betrachtet man die Umlaufszeiten der beiden Planeten, so scheint in der Tat eine 13:8 Resonanz zu existieren.

Auch bei andere Planeten lassen die Umlaufszeiten Resonanzen vermuten. Die Bahnperioden von Jupiter und Saturn verhalten sich nahezu wie 5:2, diejenigen von Neptun und dem Zwergplaneten Pluto fast wie 3:2.

Analyse (quasi)periodischer Bewegungen

Zerlegung in harmonische Schwingungen

Obige Beispiele lassen erkennen, dass man die zeitlichen Veränderungen der Bahnelemente als eine Überlagerung periodischer Funktionen beschreiben kann. So folgt die Exzentrizität der Erdbahn im Zwei-Planeten-Modell offensichtlich weitgehend einer Sinusschwingung, welche eine Periode von etwa 200000 Jahren und eine Amplitude von ungefähr 0.01 aufweist. Die jährliche Änderung der Exzentrizität zeigt klar auf, dass zusätzlich viel raschere Schwingungen vorhanden sind, die nur wenige Jahre dauern und deren Amplituden in der Größenordnung 0.00001 liegen.

Das Drei-Planeten-Modell lässt auch langfristig mehrere sich überlagernde Zyklen erahnen. Die Exzentrizitäten der Erd- und Venusbahn beschreiben selbst aus diesem Blickwinkel keine einfachen Sinusschwingungen, sondern weisen zusätzliche Variationen auf, die sich über etwa 10000 Jahre erstrecken und Amplituden auf einer Skala von 0.001 aufweisen.

Wie der Mathematiker Fourier zu Beginn des 19. Jh. zeigte, kann eine Funktion $f(t)$ in der Tat als eine Überlagerung harmonischer Schwingungen aufgefasst werden, wobei c eine Konstante darstellt.

$f(t)=c+\sum _{k=1}^{\infty }A_{k}\cos(k\omega t+\phi _{k})$

Jede Schwingung hat eine Amplitude $A_{k}$, eine Phase $\phi _{k}$ und eine Kreisfrequenz $k\omega$. Den ersten Summanden bezeichnet man als Grundschwingung, die übrigen als Oberschwingungen.

$f(t)$ muss nicht unbedingt wie hier für die Planetenbahnen eine Zeitfunktion sein. Wellen beispielsweise zeigen nicht nur zeitlich, sondern auch räumlich ein periodisches Verhalten, das mit Fouriers Methode ebenso untersucht werden kann.

Die Darstellbarkeit einer Funktion $f(t)$ als eine Summe harmonischer Schwingungen bedeutet nicht, dass eine solche physikalisch tatsächlich durch solche Oszillationen repräsentiert wird. Das quasiperiodische Verhalten der Bahnelemente entsteht dadurch, dass bestimmte Relativpositionen der Planeten untereinander unterschiedlich oft (fast) regelmäßig wiederkehren.

Reale harmonische Schwingungen findet man hingegen bei der soeben erwähnten Wellenausbreitung und damit z.B. bei klassischen Musikinstrumenten. Spielt man einen Grundton mit einer bestimmten Kreisfrequenz $\omega$, so erklingen in der Tat auch sogenannte Obertöne höherer Frequenz. Die Amplitudenverhältnisse all dieser Töne untereinander bestimmen die Klangfarbe des Instruments. Im Idealfall entstehen beim Musizieren nur Töne mit einem ganzzahligen Vielfachen von $\omega$. Treten auch andere Kreisfrequenzen auf, so wird der Klang oft als unharmonisch empfunden. Bei einigen Instrumenten wie z.B. dem Klavier, tragen aber gerade solche unharmonischen Obertöne wesentlich zur Klangfarbe bei.

Für die praktische Auswertung von $f(t)$ ist es einfacher, trigonometrische Funktionen ohne eine Phase zu verwenden. Mit Hilfe der Beziehung $\cos(k\omega t+\phi _{k})=\cos(k\omega t)cos(\phi _{k})-\sin(k\omega t)\sin(\phi _{k})$ gelangt man unmittelbar zu solch einer Form.

$f(t)=c+\sum _{k=1}^{\infty }(a_{k}\cos(k\omega t)+b_{k}\sin(k\omega t))$

Die Koeffizienten $a_{k}$ und $b_{k}$ lassen sich folgendermaßen in Amplituden und Phasen umrechnen:

$a_{k}=A_{k}\cos(\phi _{k})$

$b_{k}=-A_{k}\sin(\phi _{k})$

$A_{k}={\sqrt {a_{k}^{2}+b_{k}^{2}}}$

Gemäß Fourier existieren für jede Grundschwingung unendlich viele Oberschwingungen. In der Praxis jedoch ist die zu untersuchende Funktion nur mit einer endlichen zeitlichen bzw. räumlichen Auflösung bekannt, z.B. die simulierten Bahnelemente eines Planeten nur in diskreten Zeitabständen $\Delta t$. Die höchste Frequenz, die bei solch eingeschränkter Kenntnis der Zeitfunktion noch ermittelt werden kann, entspricht gemäß dem Abtasttheorem von Nyquist und Shannon dem Intervall $2\Delta t$. Befindet sich eine derart rasche Schwingung zum Zeitpunkt 0 am Minimum, so steht sie zur Zeit $\Delta t$ am Maximum. Zum Zeitpunkt $2\Delta t$ ist das nächste Minimum erreicht. Ein Zeitintervall $2\Delta t$ reicht so gerade noch aus, pro Periode 3 Messwerte zu erfassen. Versucht man die Amplituden und Phasen noch schnellerer Oszillationen zu ermitteln, liefert die Fourier-Analyse zu höheren Frequenzen hin immer mehr Artefakte, da mit nur 2 verfügbaren Messpunkten pro Schwingung deren Amplitude und Phase nicht mehr eindeutig festgelegt werden können.

Um die Koeffizienten $a_{k}$ und $b_{k}$ zu bestimmen, normiert man das Zeitintervall $T$, für welches $f(t)$ bekannt ist, auf $2\pi$. Sind für die Zeitfunktion insgesamt $N$ äquidistante Werte bekannt, welche nacheinander von 0 bis $N-1$ durchnummeriert werden, so gilt:

$a_{k}={\frac {2}{N}}\sum _{n=0}^{N-1}f_{n}\cos \left(k{\frac {2\pi n}{N}}\right)$

$b_{k}={\frac {2}{N}}\sum _{n=0}^{N-1}f_{n}\sin \left(k{\frac {2\pi n}{N}}\right)$

Mit einem auf $2\pi$ normierten Messintervall ist die Kreisfrequenz jeder Schwingung direkt durch $k$ gegeben. Macht man die Normierung wieder rückgängig, so gilt $\omega =2\pi k/T$ bzw. für die Periode der Schwingung $P=T/k$.

$k$ kann Werte von 1 bis $N/2$ annehmen. Da für jede Schwingung Summen mit $N$ Gliedern auszuwerten sind, ist auf den ersten Blick der Rechenaufwand proportional zu $N^{2}$. Spezielle Algorithmen, die bestimmte Eigenschaften der trigonometrischen Funktionen ausnutzen, gestatten jedoch eine Berechnung mit einem Aufwand proportional zu $N\log N$. Man kann dies erreichen, indem man die Summanden geschickt umsortiert.

Jeder solche ist ein Produkt aus einem Wert der Zeitfunktion und einer trigonometrischen Funktion. Alle Kombinationen aus $k$ und $n$, für welche die Division $kn/N$ den gleichen Rest liefert, weisen dabei wegen des periodischen Verhaltens der trigonometrischen Funktionen gleiche Sinus- bzw. Cosinuswerte auf. Wenn man die Summanden nach eben diesen Divisionsresten sortiert, kann man alle zu einem bestimmten Rest gehörenden $f_{n}$ addieren und muss dann nur diese Summe mit der trigonometrischen Funktion multiplizieren. Die Verringerung des Rechenaufwands geschieht also durch eine Reduzierung der Anzahl der Multiplikationen. Besonders effizient lassen sich die Summanden sortieren, wenn die Anzahl der Messwerte $N$ für die Zeitfunktion eine ganzzahlige Potenz von 2 darstellt.

Für das Umsortieren sind mehrere Verfahren bekannt, welche als schnelle Fourier-Transformationen bezeichnet werden. Solche werden von den meisten gängigen Plot-Programmen unterstützt, so dass hier nicht weiter darauf eingegangen wird.

Einschub für Fortgeschrittene: Komplexe Fourierreihe und Fourierintegral

Um eine Funktion $f(t)$ als eine Summe von Schwingungen darzustellen, verwendet man in der Praxis anstatt der beiden trigonometrischen Funktionen Simus und Cosinus oft die komplexe Exponentialfunktion:

$f(t)=\sum _{k=0}^{\infty }a_{k}e^{ik\omega t}$

Sowohl die Amplitude $A_{k}$ als auch die Phase $\phi _{k}$ einer jeden Schwingung stecken nun in einem einzigen, jedoch ebenfalls komplexen Koeffizienten $a_{k}$:

$A_{k}={\sqrt {(Re(a_{k}))^{2}+(Im(a_{k}))^{2}}}$

$Re(a_{k})=A_{k}\cos \phi _{k}$

$Im(a_{k})=A_{k}\sin \phi _{k}$

Für die Auswertung der Koeffizienten normiert man das von $N$ äquidistanten Messpunkten $f_{n}$ überdeckte Zeitintervall $T$ abermals auf $2\pi$. Anstelle des Sinus und Cosinus genügt wiederum die komplexe Exponentialfunktion:

$a_{k}={\frac {2}{N}}\sum _{n=0}^{N-1}f_{n}e^{-ik{\frac {2\pi n}{N}}}$

Die komplexe Exponentialfunktion weist die gleichen periodischen Eigenschaften auf wie Sinus und Cosinus. Dementsprechend lassen sich auch für die komplexe Darstellung Summanden, für die $kn/N$ den gleichen Rest aufweist, zusammenfassen und damit die Berechnung der $a_{k}$ erheblich optimieren.

Im Idealfall einer beliebig feinen zeitlichen (oder je nach Art des Problems auch räumlichen) Auflösung von $f(t)$ tritt an die Stelle einer Summe über harmonische Schwingungen ein Integral über ein kontinuierliches, d.h. auch nicht harmonische Anteile enthaltendes Spektrum:

$f(t)=\int _{-\infty }^{\infty }a(\omega )e^{i\omega t}d\omega$

Da die Kreisfrequenz nun beliebige Werte und nicht nur ein ganzzahliges Vielfaches einer Grundfrequenz annehmen kann, wird kein Zähler $k$ mehr benötigt. Der zu einer Frequenz jeweils gehörige Koeffizient folgt aus der Vorschrift:

$a(\omega )={\frac {1}{2\pi }}\int _{-\infty }^{\infty }f(t)e^{-i\omega t}dt$

Amplitude und Phase folgen aus den Real- und Imaginärteilen von $a(\omega )$ auf gleiche Weise wie für diskrete $a_{k}$.

Anwendung auf Simulationsbeispiele

Das Fourier-Verfahren zur Ermittlung der Grund- und Oberschwingungen einer Zeitfunktion wird nun auf den zeitlichen Verlauf der Exzentrizität angewandt. Bei einer über 500000 Jahre sich erstreckenden Simulation liefert es Amplituden für Schwingungskomponenten mit Perioden von 500000 / $k$ Jahren.

Betrachtet man die Erdbahn im Rahmen des Zwei-Planeten-Modells, so ist die Schwingung mit $k$ = 3 am stärksten. Sowohl die Amplitude von etwa 0.01 als auch die Periode von 167000 Jahren entspricht dem optischen Erscheinungsbild der Zeitfunktion. Es liegt aber keine reine Sinusschwingung vor, da vor allem die Komponenten in der Umgebung der Hauptschwingung durchaus noch vergleichbar starke Amplituden von bis zu 0.003 aufweisen. Für sämtliche Oszillationen mit einem $k$ zwischen 1 und 6 finden sich Amplituden von zumindest 0.001.

Mit zunehmendem $k$, d.h. abnehmender Periode, nimmt die Amplitude zumeist monoton ab. Im Bereich von Perioden der Größenordnung 1-10 Jahre treten jedoch markante Spitzen auf, deren Amplituden allerdings höchstens bis an eine Größenordnung von 0.00001 heranreichen. Diese Komponenten sind für die zahlreichen, jedoch geringfügigen Schwankungen der Exzentrizität verantwortlich, die sich bei einer Betrachtung mit hoher Zeitauflösung zeigen.

Schwingungskomponenten der zeitlichen Änderung der Exzentrizität der Erdbahn für obiges Zwei-Planeten-Modell

Im Drei-Planeten-Modell findet sich für die Erde die Hauptkomponente erneut bei $k$ = 3, für die Venus aber bei $k$ = 4 entsprechend einer Periode von 125000 Jahren. Die Amplitude liegt jeweils um 0.01. Vor allem bei der Venus, aber auch noch für die Erde, ist das Maximum im Vergleich zum Zwei-Planeten-Modell aber weniger scharf ausgeprägt. Für die Erde haben jetzt alle Komponenten bis $k$ = 7 eine Amplitude von zumindest 0.001, für die Venus gilt dies bis zu $k$ = 8. Bei Perioden von einigen 10000 Jahren treten für die Erde zudem noch einige signifikante Maxima mit Amplituden von einigen 0.0001 auf. All dies stimmt gut mit dem Bild eines auch auf langer Zeitskala komplizierten Verlaufs der Exzentrizität überein, der kaum noch Ähnlichkeit mit einer einfachen Sinusschwingung hat.

Wie für das Zwei-Planeten-Modell nimmt die Amplitude zu kurzen Perioden hin stark ab. Die dort sichtbaren Spitzen sind für das Drei-Planeten-Modell aber wesentlich zahlreicher, was den optischen Eindruck eines äußerst komplexen zeitlichen Verlaufs bei hoher Auflösung bestätigt.

Schwingungskomponenten der zeitlichen Änderung der Exzentrizität der Erd- und Venusbahn für obiges Drei-Planeten-Modell

Langzeitsimulationen unter Einbeziehung aller Planeten offenbaren im zeitlichen Verlauf der Bahnelemente die Existenz einer Vielzahl markanter Schwingungen mit z.T. noch viel längeren Perioden als hier im Rahmen der vereinfachten Modelle angedeutet werden kann. So geben Laskar und andere (2011) für die Exzentrizität des Erdorbits Zyklen von 95000, 124000, 405000, 688000, 1 Million und 2.4 Millionen Jahren an. All diese Komponenten weisen nur wenig unterschiedliche Amplituden im Bereich 0.001 bis 0.01 auf. Dies hat zur Folge, dass auch der langfristige Verlauf der Exzentrizität ein schwer durchschaubares Muster annimmt, wobei diese wie schon erwähnt Werte zwischen 0 und 0.07 annimmt.

Analytische Behandlung des gestörten Zwei-Körper-Problems (für Fortgeschrittene)

Eine Lösung des gestörten Zwei-Körper-Problems erfordert fundierte Kenntnisse der Differential- und Integralrechnung, so dass dieser Abschnitt an Fortgeschrittene gerichtet ist. Es existieren hierfür etliche Verfahren, welche im Rahmen dieses Buches jedoch nicht alle vorgestellt werden können. Das nun als Beispiel detailliert erörterte Vorgehen von Wisdom und Holman (1991) wurde ausgewählt, weil es dem quasiperiodischen Verhalten der Bahnelemente besonders prägnant Rechnung trägt.

Hamilton-Funktion und Hamiltonsche Gleichungen

Die Lösung vieler Probleme der Mechanik beruht auf der Hamilton-Funktion. Diese gibt die Gesamtenergie eines Systems von $N$ Massenpunkten als Funktion dessen Phasenraumes, d.h. in Abhängigkeit von den Impulskoordinaten ${\vec {p}}_{0},{\vec {p}}_{1}...{\vec {p}}_{N-1}$ und Ortskoordinaten ${\vec {q}}_{0},{\vec {q}}_{1}...{\vec {q}}_{3N-1}$ der einzelnen Mitglieder des Ensembles an.

$H=H({\vec {p}}_{0}...{\vec {p}}_{N-1},{\vec {q}}_{0}...{\vec {q}}_{N-1})$

Allgemein kann die Hamilton-Funktion auch zeitabhängig sein. Für die Planetenbewegungen sind jedoch selbst schon die von den nächsten Sternen ausgeübten Gravitationskräfte belanglos, so dass das Sonnensystem als abgeschlossen und damit die Gesamtenergie als konstant betrachtet werden darf.

Die große Bedeutung der Hamilton-Funktion liegt darin, dass sie gemäß der Hamiltonschen Gleichungen Impuls- und Ortskoordinaten miteinander verknüpft:

${\frac {\partial H}{\partial p_{i}}}={\dot {q}}_{i}$

${\frac {\partial H}{\partial q_{i}}}=-{\dot {p}}_{i}$

Die Ableitung der Gesamtenergie nach der Impulskoordinate eines Massenpunktes liefert dessen entsprechende Geschwindigkeitskoordinate, die Ableitung von $H$ nach einer Ortskoordinate die in der entsprechenden Richtung auf den Körper einwirkende Kraft (und damit Beschleunigung). Die zweite Gleichung stellt eine umgekehrte Sichtweise des Arbeitsprinzips dar. Bewegt sich ein Massenpunkt eine gewisse Strecke in einem Kraftfeld, so gewinnt (oder verliert) er potentielle Energie. Dementsprechend beschreibt die Änderung der potentiellen Energie entlang dieser Strecke die unterdessen ausgeübte Kraft.

Verwendet man gewöhnliche Orts- und Impulskoordinaten, so lautet die Hamilton-Funktion eines durch die Schwerkraft bestimmten Mehrkörpersystems (wobei $m_{i}$ die Masse des Mitglieds $i$ und $r_{ij}$ der Abstand zwischen zwei Körpern $i$ und $j$ bedeuten):

$H=\sum _{i=0}^{N-1}{\frac {p_{i}^{2}}{2m_{i}}}-\sum _{i<j}{\frac {Gm_{i}m_{j}}{r_{ij}}}$

Der erste Summand gibt die kinetische, der zweite die potentielle Energie an. Da die Bewegung jedes Planeten durch die Sonne und nicht durch die jeweils anderen Planeten dominiert wird, ist eine Darstellung der folgenden Form wünschenswert, wobei $\mu _{i}=m_{0}m_{i}/(m_{0}+m_{i})$ die reduzierte Masse des Trabanten, $m_{0}$ die Sonnenmasse sowie $p_{i}$ und $r_{i}$ Relativimpuls und -abstand zwischen Sonne und Planet angeben:

$H=\sum _{i=1}^{N-1}\left({\frac {p_{i}^{2}}{2\mu _{i}}}-{\frac {Gm_{0}m_{i}}{r_{i}}}\right)+H_{S}$

Man will also den Energiebeitrag jedes Planeten wie im ungestörten Zweikörperproblem beschreiben, wobei die von diesen untereinander definierte potentielle Energie nur eine kleine Störung $H_{S}$ der Gesamtenergie bewirkt. Unter Verwendung eines einfachen Koordinatensystems wie Schwerpunkt- oder heliozentrisches System lässt sich dies jedoch nicht bewerkstelligen.

Hamilton-Funktion mit Jacobi-Koordinaten

Die Trennung der potentiellen Energie in einen ungestörten Zwei-Körper- und einen Störanteil gelingt mittels sogenannter Jacobi-Koordinaten, welche durch folgende Iteration generiert werden. Ausgangspunkt ist die Zentralmasse $m_{0}$, welche in den Ursprung - definiert wiederum durch den Schwerpunkt des ganzen Ensembles - gesetzt wird. Die Position des 1. Planeten wird relativ zu derjenigen von $m_{0}$ angegeben, die Position des 2. Planeten relativ zum Schwerpunkt, welcher durch die Zentralmasse und den 1. Planeten gebildet wird. Allgemein wird jeder hinzutretende Planet relativ zum Schwerpunkt aller bisher eingeordneten Körper betrachtet:

$({\vec {x}}_{i})_{J}={\vec {x}}_{i}-\sum _{j=0}^{i-1}m_{j}{\vec {x}}_{j}/\sum _{j=0}^{i-1}m_{j}$

In welcher Reihenfolge diese Umrechnung erfolgt, ist an für sich egal, doch in der Praxis beginnt man zumeist mit dem innersten Planeten und geht von diesem aus sukzessive zum äußersten.

Die Geschwindigkeiten werden wie die Ortskoordinaten behandelt. Um den Impuls korrekt in das Jacobi-Koordinatensystem zu übertragen, ist aber zusätzlich eine Massentransformation erforderlich. Die Jacobi-Masse der Zentralmasse entspricht der Gesamtmasse des Systems. Für die Planeten gilt:

$(m_{i})_{J}=m_{i}\sum _{j=0}^{i-1}m_{j}/\sum _{j=0}^{i}m_{j}$

Damit ist analog zum gewöhnlichen Schwerpunktsystem $({\vec {p}}_{i})_{J}=(m_{i})_{J}({\vec {v}}_{i})_{J}$. Das Einsetzen der neuen Koordinaten in die Gesamtenergie liefert:

$H={\frac {(p_{0})_{J}^{2}}{2(m_{0})_{J}}}+\sum _{i=1}^{N-1}{\frac {(p_{i})_{J}^{2}}{2(m_{i})_{J}}}-\sum _{i=1}^{N-1}{\frac {Gm_{i}m_{0}}{r_{i0}}}-\sum _{0<i<j}{\frac {Gm_{i}m_{j}}{r_{ij}}}$

Der erste Term gibt schlicht die kräftefreie Bewegung des Schwerpunkts an und wird so in der Folge nicht weiter benötigt. Die Gesamtenergie liegt jetzt fast schon in der angestrebten Weise vor, doch sind für die beiden letzten Terme die Massen und insbesondere die Abstände zwischen diesen noch in gewöhnlichen statt Jacobi-Koordinaten gegeben. Durch Hinzufügen und gleichzeitiges Subtrahieren des Glieds $\sum _{i=1}^{N-1}Gm_{i}m_{0}/(r_{i})_{J}$ gewinnt man aber den Ausdruck:

$H=\sum _{i=1}^{N-1}\left({\frac {(p_{i})_{J}^{2}}{2(m_{i})_{J}}}-{\frac {Gm_{i}m_{0}}{(r_{i})_{J}}}\right)+\sum _{i=1}^{N-1}Gm_{i}m_{0}\left({\frac {1}{(r_{i})_{J}}}-{\frac {1}{r_{i0}}}\right)-\sum _{0<i<j}{\frac {Gm_{i}m_{j}}{r_{ij}}}$

Das Ziel der Darstellung der potentiellen Energie als Summe von ungestörter und gestörter Komponente ist damit erreicht, doch letztere erscheint noch sehr kompliziert. Wegen der Dominanz der Zentralmasse ist der Unterschied zwischen gewöhnlichen und Jacobi-Koordinaten aber gering. Dies gestattet es, den heliozentrischen Abstand $r_{i0}$ in eine Taylorreihe um den in Jacobi-Koordinaten gegebenen Abstand $(r_{i})_{J}$ zu entwickeln und diese Reihe nach dem linearen Glied abzubrechen. Unter Vernachlässigung aller höheren Terme erhält man schließlich:

$H=\sum _{i=1}^{N-1}\left({\frac {(p_{i})_{J}^{2}}{2(m_{i})_{J}}}-{\frac {Gm_{i}m_{0}}{(r_{i})_{J}}}\right)-\sum _{0<i<j}\left({\frac {Gm_{i}m_{j}}{(r_{ij})_{J}}}-{\frac {Gm_{i}m_{j}({\vec {x}}_{i})_{J}\cdot ({\vec {x}}_{j})_{J}}{(r_{i})_{J}^{3}}}\right)$

Modellierung der Störfunktion durch periodische Stöße

Der entscheidende Gedanke von Wisdom und Holman besteht nun darin, den Störanteil der Hamilton-Funktion nicht permanent, sondern periodisch wirken zu lassen. Die Autoren stellen somit die Störung der Bahn eines Planeten durch die jeweils anderen Mitglieder des Sonnensystems durch eine regelmäßige Folge kleiner Stöße da, die auf den betroffenen Körper einwirken. Mathematisch gescheht dies durch Multiplikation des Störanteils mit der periodischen Diracschen Delta-Funktion:

$H=H_{K}+2\pi \delta _{2\pi }(\Omega t)H_{S}$

$\delta _{2\pi }(\Omega t)$ ist streng genommen keine Funktion, sondern eine sogenannte Distribution, welche durch Grenzwertbildung aus einer Funktionenreihe hervorgeht:

$\delta _{2\pi }(\Omega t)=\lim _{N\to \infty }{\frac {1}{2\pi }}\sum _{n=-N}^{N}\cos(n\Omega t)$

Nachfolgendes Diagramm zeigt die periodische Diracsche Deltafunktion für den Sonderfall $\Omega =1$, angenähert durch Berücksichtigung der Summanden mit $n$ = -1..1, -3..3 und -5..5.

Die periodische Diracsche Delta-Funktion als Grenzwert der Summe $\lim _{N\to \infty }{\frac {1}{2\pi }}\sum _{n=-N}^{N}\cos(nx)$

An den Stellen 0, 2$\pi$, 4$\pi$ usw. liefert jeder Summand einen Beitrag von 1, somit divergiert dort $\delta _{2\pi }(t)$ gegen $+\infty$. An allen anderen Stellen treten mit zunehmendem $n$ immer mehr destruktive Interferenzen zwischen den Summanden auf, so dass die periodische Diracsche Delta-Funktion gegen 0 strebt. Das über eine Periode genommene Integral über $\delta _{2\pi }(t)$ bleibt gleichwohl endlich mit einem Gesamtwert von 1. Für jeden Summanden mit von 0 verschiedenem $n$ verschwindet es nämlich, für den Summanden mit $n=0$ liefert es den Wert 1. All diese Eigenschaften garantieren, dass über eine Periode gemittelt der Betrag der modifizierten Hamilton-Funktion gegenüber der originalen unverändert bleibt.

Durch die Multiplikation der Störenergie $H_{S}$ mit der periodischen Diracschen Delta-Funktion wird auch das gestörte Kepler-Problem analytisch lösbar. Der freie Parameter $\Omega$ gibt hierbei die Frequenz an, mit welcher im Modell ein Planet durch die jeweils anderen gestoßen wird, dementsprechend $1/\Omega$ die zeitliche Schrittweite, mit der ein solcher simuliert wird. Je größer die Stoßfrequenz gewählt wird, umso besser entspricht das Modell der Realität einer kontinuierlich einwirkenden Störung, umso größer ist jedoch dann auch der Rechenaufwand angesichts geringer zeitlicher Abstände zwischen aufeinanderfolgenden Stößen. Wie auch bei anderen Lösungsverfahren, dürfen für das Stoßmodell im Vergleich zu einer rein numerischen Integration die Zeitschritte viel gröber gewählt werden. Gemäß Wisdom and Holman reicht es aus, die Stoßfrequenz gleich der (ungestörten) Bahnperiode eines Planeten zu setzen. Hingegen wurden für die auf dem Hermite-Polynome-Verfahren beruhenden Beispiele in diesem Kapitel etwa 300 Zeitschritte pro Bahnperiode aufgewandt.