Wie ThreeBody rechnet

Zum Einstieg: die ganze Idee in drei Schritten

Wer mit Formeln wenig anfangen kann, braucht trotzdem nicht weit zu lesen: Das ganze Programm besteht aus drei Gedanken, und alles Weitere auf dieser Seite ist nur die Frage, wie man sie sauber aufschreibt.

Erstens: Zwei Massen ziehen sich an

Newtons Gravitationsgesetz von 1687. Zwei Körper mit den Massen m₁ und m₂ im Abstand r ziehen einander mit der Kraft

an. Mehr Masse, mehr Anziehung — und der Abstand geht quadratisch ein: Doppelter Abstand heißt ein Viertel der Kraft, zehnfacher Abstand ein Hundertstel. G ist eine Naturkonstante, die nur dafür sorgt, dass am Ende die richtigen Einheiten herauskommen.

Was uns eigentlich interessiert, ist aber nicht die Kraft, sondern die Beschleunigung — wie stark ein Körper seine Bewegung ändert. Die folgt aus F = m·a, und hier passiert etwas Bemerkenswertes:

Die eigene Masse kürzt sich weg. Wie stark ein Körper angezogen wird, hängt nur davon ab, wer ihn anzieht — nicht davon, wie schwer er selbst ist. Deshalb fallen Feder und Hammer auf dem Mond gleich schnell, und deshalb kann diese Simulation Staubkörnern die Masse null geben und sie trotzdem richtig fliegen lassen (Abschnitt 4).

Ein Beispiel zum Mitrechnen: Die Sonne hat 1,989·10³⁰ kg, die Erde ist 1,496·10¹¹ m entfernt. Damit fällt die Erde mit a ≈ 0,006 m/s² auf die Sonne zu — sechs Millimeter pro Sekunde schneller, jede Sekunde. Dass sie trotzdem nicht ankommt, liegt daran, dass sie sich gleichzeitig mit 30 km/s zur Seite bewegt. Genau dieses Nebeneinander von „fällt“ und „fliegt seitlich vorbei“ ist eine Umlaufbahn.

Zweitens: Kräfte addieren sich

Im Sonnensystem gibt es nicht zwei Körper, sondern zehn. Jeder zieht an jedem. Die Beschleunigung eines Körpers ist dann einfach die Summe über alle anderen — Kräfte addieren sich wie Pfeile: Länge und Richtung zählen, zwei entgegengesetzte Anziehungen heben sich teilweise auf.

Das ist der Grund, warum r in den späteren Formeln fett gedruckt ist: Es steht nicht für eine Zahl, sondern für die drei Zahlen x, y und z, also für einen Pfeil im Raum. Und es erklärt, warum dort r³ im Nenner steht, wo eben noch r² stand: Der Zähler enthält den Verbindungspfeil, und der ist selbst schon r lang. Einmal r geht also dafür drauf, aus dem Pfeil eine reine Richtung zu machen.

Genau hier endet übrigens die Schulmathematik. Für zwei Körper lässt sich die Bahn hinschreiben — das sind Keplers Ellipsen. Ab drei Körpern gibt es keine solche Formel mehr, und zwar nicht, weil sie noch niemand gefunden hätte, sondern weil es sie nicht gibt. Übrig bleibt: ausrechnen, Schritt für Schritt.

Drittens: Zweimal aufsummieren

Beschleunigung ist die Änderung der Geschwindigkeit, Geschwindigkeit die Änderung des Ortes. Wer die Beschleunigung kennt, kommt also in zwei Schritten zur Bahn — indem er sie zweimal aufsummiert. Für ein kleines Zeitstückchen Δt heißt das schlicht:

Damit ist die Simulation vollständig. Sie ist eine Schleife über diese drei Schritte, und in Pseudocode passt sie in acht Zeilen:

wiederhole sehr oft:

    für jeden Körper i:                     ← Schritt 1 und 2
        a[i] = 0
        für jeden anderen Körper j:
            a[i] += G · m[j] · (r[j] − r[i]) / |r[j] − r[i]|³

    für jeden Körper i:                     ← Schritt 3
        v[i] += a[i] · Δt
        r[i] += v[i] · Δt

Im Sonnensystem ist Δt = 3,5 Stunden. Die Erde rückt in einem solchen Schritt um rund 375 000 km weiter, etwa ein Zweitausendfünfhundertstel ihrer Bahn; für einen Umlauf sind es 2500 Schleifendurchläufe, und Ihr Browser rechnet einige hundert davon pro Sekunde.

1Das Modell

Simuliert wird das newtonsche N-Körper-Problem im Raum: N Punktmassen, die einander mit dem Gravitationsgesetz anziehen, und sonst nichts. Für jeden Körper i gilt

Das ist ein Hamiltonsches System. Mit den Impulsen p_i = m_i ṙ_i lautet die Energie

Dass sich H sauber in einen nur von den Impulsen und einen nur von den Orten abhängigen Teil zerlegen lässt, ist keine Formalie — genau darauf beruht das ganze Integrationsverfahren in Abschnitt 3.

Was das Modell nicht enthält

Beliebig viele Körper

Gerechnet wird durch Direktsummation, also jedes Paar einzeln. Der Aufwand wächst quadratisch mit der Zahl der Körper; Baum- oder Multipolverfahren, mit denen große Rechnungen Millionen von Teilchen schaffen, gibt es hier nicht. Bei höchstens 42 Körpern lohnen sie nicht.

Relativität

Gerechnet wird newtonsch. Zwei Zusätze treten hinzu und ändern daran nichts Grundsätzliches: die Rückwirkung der Gravitationsabstrahlung (Abschnitt 6), als Reibungsterm aufgesetzt und nicht hergeleitet, und ein post-newtonscher Korrekturterm erster Ordnung (Abschnitt 7), der das Perizentrum wandern lässt. Der zweite ist die Näherung für eine beherrschende Masse, nicht die vollständige Bewegungsgleichung; wo es keine solche Masse gibt, steht er deshalb nicht zur Wahl. Alles Übrige fehlt weiterhin: keine Lichtlaufzeiten, keine Lichtablenkung, keine zweite Ordnung.

Ausdehnung

Alle Körper sind Punktmassen. Ihr Radius spielt nur an zwei Stellen eine Rolle: er entscheidet, wann zwei Körper verschmelzen, und wann ein Testteilchen verschluckt wird. Gezeiten, Abplattung, Rotation: nichts davon.

Alles Nichtgravitative

Kein Gasdruck, keine Reibung, kein Strahlungsdruck, keine Magnetfelder. Die „Akkretionsscheibe“ der Schwarzen Löcher ist deshalb kein hydrodynamisches Modell, sondern ein Ring aus Probekörpern.

2Einheiten, und warum sie sich lohnen

Jedes Szenario rechnet in seinem eigenen Satz aus Länge L, Masse M und Zeit T, nicht in SI. Umgerechnet werden dabei nicht die Zahlen, sondern die Konstanten:

In AE, Sonnenmassen und Jahren wird daraus G = 4π² ≈ 39,478 — das dritte Keplersche Gesetz in Reinform. Der Prüflauf rechnet genau das nach.

Der Grund ist die Fließkommaarithmetik. Ein double hat 53 Bit Mantisse, also eine relative Auflösung von rund 10⁻¹⁶; der absolute Rundungsfehler wächst mit der Größe der Zahl. In SI wäre eine Bahnhalbachse 10¹¹, ihre dritte Potenz 10³³, das Produkt Gm 10²⁰ — mit Absolutfehlern in derselben Größenordnung mal 10⁻¹⁶. Bleiben alle Zahlen dagegen in der Nähe von 1, bleiben es auch die Fehler. Die Idee stammt aus der Java-Fassung von 2016 und ist der billigste Genauigkeitsgewinn im ganzen Programm.

3Zeitintegration

3.1 Warum symplektisch

Ein N-Körper-System ist ein Hamiltonsches System, und sein Fluss hat Eigenschaften, die ein beliebiges Integrationsverfahren zerstört: Er erhält das Phasenraumvolumen, er ist zeitumkehrbar, und er erhält die Energie. Ein Runge-Kutta-Verfahren vierter Ordnung ist pro Schritt genauer als das hier benutzte, verliert aber genau diese Struktur — sein Energiefehler wächst säkular, also stetig in eine Richtung. Nach einer Million Schritten ist der Planet dann nicht mehr auf seiner Bahn, sondern in einer Spirale.

Symplektische Verfahren machen den umgekehrten Handel. Sie lösen die gestellte Aufgabe nicht exakt, aber sie lösen exakt eine benachbarte: Zu einem symplektischen Verfahren der Ordnung p existiert ein modifiziertes Hamiltonsches System

dessen Bahnen die berechneten Punkte über sehr lange Zeiten treffen. Die wahre Energie H weicht davon um 𝒪(h^p) ab — aber sie schwankt um diesen Betrag, statt wegzulaufen. Genau das sieht man im Ausleser „Energiefehler“: eine Zahl, die auf und ab geht und nicht wächst. Das ist der einzige Grund, warum eine Simulation im Browser über Millionen von Schritten überhaupt etwas Sinnvolles zeigen kann.

3.2 Velocity-Verlet

Der Grundbaustein ist der Verlet-Schritt in seiner Kick-Drift-Kick-Form. Weil H = T(p) + U(r) separabel ist, sind beide Teilflüsse einzeln exakt lösbar: unter T allein fliegen alle Körper geradeaus („Drift“), unter U allein ändern sich nur die Geschwindigkeiten („Kick“). Ihre symmetrische Verkettung ist

Das ist zweiter Ordnung, zeitumkehrbar und kostet eine Kraftauswertung je Schritt — die Beschleunigung am Ende ist die vom Anfang des nächsten. Entscheidend für die Erhaltungssätze: jede Paarkraft wirkt im selben Kick auf beide Partner, mit umgekehrtem Vorzeichen. Der Gesamtimpuls bleibt dadurch bis auf Rundung exakt erhalten, unabhängig von der Schrittweite.

3.3 Yoshida, vierte Ordnung

Aus drei Verlet-Schritten mit passend gewählten, teils negativen Zeitanteilen entsteht ein Verfahren vierter Ordnung (Yoshida 1990). Mit Φ^V_h für den Verlet-Schritt (5) ist

Zahlenwerte: w₁ ≈ 1,351207, w₀ ≈ −1,702414, und 2w₁ + w₀ = 1, wie es sein muss. Der mittlere Schritt läuft also in der Zeit rückwärts — die Fehler dritter Ordnung der beiden äußeren Schritte heben sich damit gegen ihn weg. Der Preis sind drei Kraftauswertungen statt einer; er lohnt sich, weil die Schrittweite dafür deutlich größer sein darf. Beide Verfahren stehen in der Oberfläche zur Wahl, und der Unterschied im Energiefehler ist leicht zu sehen.

3.4 Wie die Schrittweite zustande kommt

Pro Bild soll die Simulation immer gleich weit kommen, nämlich h₀·n₀·Tempo. Wie viele Teilschritte dafür gerechnet werden, hängt davon ab, wie schnell das Gerät ist:

λ ∈ [0,15, 1] ist die Drosselung: Braucht ein Bild länger als 11 ms, sinkt sie, unter 6,5 ms steigt sie wieder.

Die Simulation wird also nicht langsamer, wenn der Rechner nicht mitkommt, sondern ungenauer — sichtbar am Energiefehler und an der Anzeige „Schritte“ neben der Bildrate. Das ist eine bewusste Entscheidung: eine ruckelnde Bahn sieht kaputt aus, eine leicht ungenaue nicht.

In fünf Szenarien kommt eine Kopplung an die Physik dazu. Für jedes Paar lässt sich eine dynamische Zeitskala angeben, die bei nahen Begegnungen klein wird:

η = 0,002 beim Pythagoräischen Problem, 0,01 bei den drei Sonnen, 0,012 bei Sgr A*, 0,015 im Sternhaufen, 0,02 im ausführlichen Sonnensystem.

Im Sonnensystem ist (8) nicht wegen einer Begegnung da, sondern wegen der Monde. Das Minimum bilden dort nie zwei Planeten, sondern dauerhaft ein Paar: Io und Jupiter, 422.000 km auseinander, 1,77 Tage für einen Umlauf. Daraus folgt eine dynamische Zeit von 6,8 Stunden und mit η = 0,02 eine Schrittweite von 8,1 Minuten — für die ganze Rechnung, auch für Neptun, der damit gut zehn Millionen Schritte je Umlauf bekommt statt der nötigen paar hundert. Ein Verfahren mit eigener Schrittweite je Körper würde das auflösen; hier ist es der Preis dafür, dass alle Körper in einem Takt laufen. Deshalb hat die einfache Fassung des Sonnensystems keine Monde, und deshalb fehlen auch der ausführlichen die schnellsten — Phobos umkreist Mars in acht Stunden und würde die Schrittweite noch einmal fünfteln.

Ohne (8) läuft Burraus Aufgabe schlicht falsch aus: Bei zu grober Schrittweite bleiben alle drei Körper beieinander, statt dass die Masse 3 hinausgeworfen wird. Die Aufgabe war deshalb bis in die 1960er Jahre ein Prüfstein für Integratoren.

4Masselose Testteilchen

Asteroidengürtel, Akkretionsscheibe und Staubring bestehen aus Teilchen mit Masse null. Sie spüren die Schwerkraft, üben aber selbst keine aus:

Die Kopplung geht nur in eine Richtung. Deshalb wächst der Aufwand für M Teilchen nur wie N·M statt wie (N+M)², und deshalb stört der Asteroidengürtel die Planeten nicht. Die Kirkwood-Lücken bei 2,50 und 2,83 AE entstehen hier folglich nicht von selbst — in der einfachen Fassung des Sonnensystems sind sie beim Ziehen der Startwerte ausgespart, echt ist nur ihre Lage. In der ausführlichen Fassung sind die Teilchen keine Ziehung, sondern 2848 wirkliche Asteroiden mit ihren gemessenen Bahnen; die Lücken stehen dort von vornherein in den Daten. Herausgerechnet werden sie auch dann nicht: Ein masseloses Teilchen spürt zwar die Resonanz mit Jupiter, aber die Millionen Jahre, die das Ausräumen gebraucht hat, laufen hier nicht ab.

Der Aufwand liegt trotzdem fast vollständig beim Staub. Im Sonnensystem stehen 45 Körperpaaren 14 000 Teilchen-Körper-Paare gegenüber, also 99,7 % der Rechenarbeit. Wer diese Simulation schneller machen will, muss genau hier ansetzen.

4.1 Ein eigener Schritt für den Staub

Bis vor Kurzem liefen die Teilchen im Takt der Körper mit: Sie wurden in jedem der drei Verlet-Teilschritte von (6) mitgezogen, bekamen also drei Kraftauswertungen je vollem Schritt. Beides ist unnötig.

Unnötig ist erstens die vierte Ordnung. Kein Ergebnis hängt von den Teilchen ab; sie sind Anschauung, nicht Messgröße. Ein gewöhnlicher Verlet-Schritt genügt. Unnötig ist zweitens die kleine Schrittweite: Die Körper zwingen sich ihre Schrittweite gegenseitig auf (das schnellste Paar bestimmt sie für alle), der Staub steht außerhalb dieser Kopplung. Im Sonnensystem ist h₀ = 4·10⁻⁴ Jahre genau die Zahl, die Merkur verlangt — der Asteroidengürtel bei 2,7 AE verträgt das Siebzehnfache.

Der Staub läuft deshalb mit der eigenen Schrittweite H = k·h, wobei k ganzzahlig ist. Damit fallen Anfang und Ende jedes Staubschritts immer auf eine Position der Körper, und es muss nichts interpoliert werden. Gerechnet wird Verlet, aber vollständig am Ende des Intervalls:

Dass die ersten beiden Zeilen erst am Ende des Intervalls ausgewertet werden, ändert am Ergebnis nichts — in ihnen stehen nur gespeicherte Größen und die verstrichene Zeit. Der Gewinn ist praktischer Natur: H muss nicht vorhergesagt werden. Ändert die Bildschleife ihre Schrittweite mittendrin, etwa weil die Drosselung λ anspringt, wird H einfach als tatsächlich verstrichene Zeit abgelesen. Es gibt keinen Fall, in dem Körper- und Staubuhr auseinanderlaufen.

4.2 Wie groß darf k sein?

Maßgeblich ist dieselbe dynamische Zeitskala wie in (8), diesmal über alle Paare aus Teilchen und Körper:

η = 1/100 bedeutet rund 2π/η ≈ 630 Schritte je Umlauf am innersten Rand der Verteilung — für masselose Deko-Teilchen reichlich. Das Minimum wird ohne zusätzliche Kosten gebildet: Es fällt in derselben Schleife an, die ohnehin jedes Paar durchläuft, und kostet dort drei Multiplikationen.

Die beiden Schranken sind wichtiger als die Formel dazwischen. Nach unten: Gerät ein Teilchen einem Körper nahe, sinkt k von selbst, notfalls bis auf 1 — die Ablenkung wird dann so fein gerechnet wie die Körper selbst. Nach oben ist k ≤ n, die Zahl der Teilschritte je Bild: Der Staub soll nie länger als ein Bild lang auf der Stelle stehen, sonst ruckelt er sichtbar, statt bloß gröber gerechnet zu werden.

4.3 Was es bringt

Verglichen werden 300 Bilder eines Szenarios bei Tempo 1× — vorher gegen nachher, gleiche Startwerte, gleiche simulierte Zeit. Die Zahl der Kraftauswertungen hängt nur vom Verfahren ab, nicht vom Gerät:

Das Lagrange-Dreieck zeigt den ungünstigen Fall: Sein Staubring liegt eng um drei Sterne, (11) drückt k dort immer wieder auf 1 herunter. Selbst dann bleiben 5,7× übrig — der Faktor 3 aus dem Verzicht auf die Yoshida-Verkettung trägt auch bei k = 1.

Bezahlt wird mit Genauigkeit, die niemand braucht. Im ruhenden Zentralfeld, wo die Lösung bekannt ist, hält ein Teilchen auf einer Kreisbahn seinen Radius über 14 Jahre:

4·10⁻⁵ relativ sind bei einer Bahn von 3 AE und einem Bildausschnitt von 26 AE Breite etwa ein Hundertstel Pixel. Die vorherigen 10⁻¹² waren Genauigkeit, die niemand je zu sehen bekam.

4.4 Was dabei liegen bleibt

• Ein Stride für alle. (11) nimmt das Minimum über sämtliche Teilchen. Ein einziges, das einem Körper nahe kommt, zwingt die übrigen 1999 auf denselben feinen Schritt. Der nächste Ausbauschritt wären Blockstufen — Teilchen in Klassen h, 2h, 4h, … einteilen, wie es Aarseth-artige Sternhaufencodes für ihre Körper tun.
• Das Einfangen wird seltener geprüft. Ob ein Teilchen in einem Körper verschwindet, wird jetzt einmal je Staubschritt geprüft, nicht mehr je Körperschritt. Ein Teilchen kann einen sehr kleinen Körper dadurch theoretisch überspringen. Praktisch verhindert (11) das: Wo ein Körper nahe ist, ist k ohnehin klein.
• Der Staub hinkt bis zu einem Bild hinterher. Zwischen zwei Staubschritten stehen die Teilchen still und springen dann. Bei k ≤ n ist der Sprung genau die Strecke, die sie in einem Bild ohnehin zurücklegen — sichtbar ist davon nichts.

5Verschmelzungen

Berühren sich zwei Körper — d ≤ R₁ + R₂ — werden sie zu einem. Der Stoß ist vollkommen unelastisch, Masse und Impuls bleiben erhalten:

Der neue Radius hängt davon ab, was verschmilzt. Bei Schwarzen Löchern ist der Schwarzschildradius r_s = 2GM/c² proportional zur Masse, die Radien addieren sich also. Alles andere behält sein Volumen:

Beim Stoß geht die Relativbewegung verloren, die Energie springt. Über einen solchen Sprung lässt sich keine Fehlerbilanz fortführen, deshalb beginnt sie danach neu: Der Sollwert wird auf die Energie unmittelbar nach der Verschmelzung gesetzt, und der Ausleser sagt ab dann „Energiefehler seit der letzten Verschmelzung“. Der verlorene Betrag wird als Ereignis mitgeschrieben und nicht stillschweigend unterschlagen.

6Gravitationswellen

Ein enges Doppelsystem verliert Energie an Gravitationswellen. Die Quadrupolformel für eine Kreisbahn (Peters & Mathews 1963) gibt die abgestrahlte Leistung

Eingebaut wird sie nicht als post-newtonscher Term, sondern als Reibungskraft entlang der Relativgeschwindigkeit v = v₂ − v₁ des schwersten Paares:

Die zweite Gleichung ist der Punkt: Die Konstruktion entzieht dem System exakt die Leistung (14), und weil die beiden Kräfte entgegengesetzt gleich sind, bleibt der Gesamtimpuls erhalten. Das genügt, um die richtige Einspiralzeit zu bekommen. Für eine Kreisbahn gibt Peters (1964)

Für 36 + 29 M☉ bei 2000 km Anfangsabstand sind das 4,771 s. Die Simulation verschmilzt nach 4,774 s — 0,06 % daneben, und das ist der einzige Zweck der Übung.

7Relativistische Korrektur

Der Schalter „Relativistische Korrektur“ schlägt auf die newtonsche Beschleunigung (1) einen Zusatzterm erster post-newtonscher Ordnung auf — 1PN, das erste Glied einer Entwicklung nach (v/c)². Er ist der Grund, aus dem Merkurs Perihel wandert und die Bahn von S2 um das Schwarze Loch im Galaktischen Zentrum keine geschlossene Ellipse ist. Für einen Körper am Ort r mit der Geschwindigkeit v, beides relativ zur schwersten Masse M des Systems, lautet er

r = ‖r‖, v = ‖v‖. Für c → ∞ verschwindet der Zusatz, und es bleibt genau (1) stehen.

Das ist die Schwarzschild-Näherung: Jeder Körper außer dem schwersten läuft als Testteilchen im Feld dieser einen Masse. Die vollständigen post-newtonschen Bewegungsgleichungen — nach Einstein, Infeld und Hoffmann — sind das ausdrücklich nicht. Sie wären eine Dreifachsumme über alle Körper und ließen sich gegen keinen bekannten Wert mehr sauber nachrechnen; (17) lässt sich prüfen, und deshalb steht (17) hier.

Aus derselben Näherung folgt etwas, das wie ein Versehen aussieht und keines ist: Die Zentralmasse bekommt keine Gegenkraft. Ein Testteilchen übt keine aus. Der Gesamtimpuls, den die Bauart des Verlet-Kicks (5) sonst bis auf Rundung exakt erhält, ist unter dem Zusatzterm damit nicht mehr streng erhalten. Eine Gegenkraft nachzutragen würde daran nichts verbessern: Wo sie ins Gewicht fiele, wäre bereits (17) die falsche Formel.

7.1 Was dabei herauskommt

Der Term liefert die Drehung des Perizentrums — des Punktes, an dem der Körper der Zentralmasse am nächsten kommt — genau in der Größe, die die Allgemeine Relativitätstheorie vorhersagt:

je Umlauf, mit der großen Halbachse a und der Exzentrizität e. Im Nenner steht der Halbparameter a(1 − e²); eng und exzentrisch heißt also viel Drehung.

Merkur bekommt daraus 42,98″ je Jahrhundert; integriert man Sonne und Merkur allein über hundert Jahre, kommen 42,92″ heraus. Gemessen sind es 43″. Bei S2 um Sgr A* liefert die Integration 11,999′ je Umlauf — die GRAVITY Collaboration hat den Vorschub an dieser Bahn 2020 nachgewiesen und dabei 1,10 ± 0,19 mal den vorhergesagten Wert von 12,1′ gefunden.

In genau diesen beiden Szenarien — Sonnensystem und Sgr A* — ist der Term deshalb von Anfang an eingeschaltet. Beide geben gemessene Bahnen wieder, und beide sind die Aufstellung, für die (17) gemacht ist: eine schwere Masse, um die herum alles andere leicht ist. Ohne den Term liefe die Rechnung um einen Betrag daneben, der seit 1859 bekannt ist. Abschalten lässt er sich, und die Wahl bleibt stehen.

Anderswo steht er dagegen gar nicht zur Wahl. „Beherrschend“ ist eine Masse, wenn die schwerste neben ihr höchstens ein Hundertstel von ihr wiegt: Genau dieses Verhältnis ist der Fehler, den (18) mit M statt M + m macht, und ein Hundertstel hält ihn unter einem Prozent. Sonne gegen Jupiter sind 1:1047, Sgr A* gegen die masselos mitlaufenden S-Sterne unendlich. Die beiden Schwarzen Löcher mit 36 gegen 29 Sonnenmassen sind es nicht, die drei Sonnen und der Sternhaufen erst recht nicht — und ein im Editor ins Sonnensystem geworfenes Loch von zehn Sonnenmassen nimmt der Sonne die Vorherrschaft ebenfalls. Dort gäbe (17) nur noch die Größenordnung der Korrektur her; eine Zahl, die wie Physik aussieht und keine ist, ist schlechter als gar keine, und deshalb verschwindet mit dem Bezugskörper auch das Kästchen.

In der Anzeige steht dafür die Zeile „Perizentrum“. Sie rechnet diese Zahl, sie misst sie nicht, und das ist die ehrlichere Angabe: Merkurs Perihel dreht sich um 574″ je Jahrhundert, und nur 43″ davon stammen aus (17). Am Bild abgelesen wäre die Drehung deshalb überwiegend die newtonsche Störung der übrigen Planeten — und man müsste für sie ein Jahrhundert simulieren. (18) ist dagegen genau die Aussage des eingeschalteten Terms, ausgewertet an den Bahnelementen des Augenblicks. Angezeigt wird der Körper, bei dem sie am größten ausfällt, dazu die seit dem Start aufgelaufene Summe. Ungebundene Bahnen bleiben dabei außen vor — S111 im Galaktischen Zentrum etwa: Ein Perizentrum, das nicht wiederkommt, dreht sich nicht „je Umlauf“.

7.2 Warum kein neuer Integrator nötig war

(17) hängt von v ab, nicht nur von r. Damit ist der Schluss-Kick des Verlet-Schritts (5) implizit — die Geschwindigkeit, die herauskommen soll, steht schon auf der rechten Seite:

a_N ist die newtonsche Beschleunigung aus (1), a_1PN der Zusatz aus (17).

Gelöst wird das als Fixpunkt: rechte Seite auswerten, das Ergebnis wieder einsetzen, fertig. Die Iteration konvergiert mit demselben Faktor, um den der Zusatz kleiner ist als das newtonsche Feld, also mit (v/c)² — 10⁻⁸ im Sonnensystem, 10⁻⁴ bei S2. Zwei Durchläufe liegen damit überall unter der Rundung des Zahlenformats; ein Abbruchkriterium wäre teurer als der Durchlauf, den es einspart.

Der eigentliche Kunstgriff ist, was dabei nicht mititeriert. a_N ist die Paarschleife aus (1) — der teure Teil, der mit dem Quadrat der Körperzahl wächst —, und sie hängt allein von den Orten ab. Nach dem Drift steht sie fest. Sie wird einmal ausgewertet und gemerkt; die Schleife darüber fasst nur den Zusatzterm an, und der ist eine Rechnung je Körper gegen die Zentralmasse, wächst also bloß linear. Die Zahl der teuren Kraftauswertungen je Schritt bleibt damit genau dieselbe wie ohne den Schalter — die Iteration kostet fast nichts. Yoshida (6) erbt die ganze Konstruktion umsonst, weil ein Yoshida-Schritt aus drei Verlet-Schritten besteht.

Runge-Kutta vierter Ordnung wäre der naheliegende Ausweg — ein explizites Verfahren, dem geschwindigkeitsabhängige Kräfte nichts ausmachen — und wurde verworfen. Es ist nicht symplektisch: Die Energie driftet säkular, statt um einen festen Wert zu schwanken (3.1). Das würde ausgerechnet die Kennzahl ruinieren, an der diese Seite ihre Qualität misst, und zwar in jedem Szenario — um eines Terms willen, der in den meisten davon in der achten Stelle steht.

Umsonst ist die Iteration trotzdem nicht zu haben. Ein Kick, dessen Beschleunigung von der Geschwindigkeit abhängt, kommt nicht mehr aus einem Hamiltonian, und der Schritt ist damit nicht mehr streng symplektisch — dieselbe Einschränkung, die schon für den Reibungsterm (15) gilt. Im ausführlichen Sonnensystem kostet das knapp eine Größenordnung Energiefehler: 6,3·10⁻¹⁴ mit ausgeschaltetem Schalter, 6,6·10⁻¹³ mit eingeschaltetem. Das bleibt gut neun Zehnerpotenzen unter dem Wert, ab dem sich die Anzeige färbt.

7.3 Die Energie

Wer einen Term zur Beschleunigung addiert, muss auch die Erhaltungsgröße nachziehen — sonst zeigte der Energiefehler (Abschnitt 8) einen systematischen Versatz an, obwohl nichts kaputt ist. Je Körper und Einheitsmasse ist die erhaltene Größe

wieder mit r und v relativ zur Zentralmasse. Die ersten beiden Summanden sind die newtonsche Energie, die in (2) ohnehin steht; die eckige Klammer ist der Zusatz.

Dass er zählt, sieht man am besten daran, was ohne ihn dastünde: im Sonnensystem 3,6·10⁻¹⁰ statt 1,3·10⁻¹¹ (Fassung „simpel“) und 2,0·10⁻¹⁰ statt 6,6·10⁻¹³ („detailliert“), also das 28- beziehungsweise 302-fache. Nicht die Kleinheit des Terms hält die Zahl klein, sondern dass er in der Bilanz steht. Wer den Schalter im Lauf umlegt, sieht das Gütemaß deshalb auch nicht springen: Der Sollwert zieht um genau den Betrag mit, um den sich der Energieausdruck ändert. Der bis dahin aufgelaufene Verfahrensfehler bleibt stehen — verschwiege man ihn, wäre die Anzeige wertlos.

(20) gilt für ein Testteilchen, und das ist keine Nachlässigkeit, sondern dieselbe Näherung wie in (17). Die allgemeine 1PN-Energie zweier Körper hat einen vierten Summanden ½(GM/r)ṙ², davor aber das Massenverhältnis ν = m₁m₂/(m₁ + m₂)² — und das ist null, solange der Begleiter als Testteilchen läuft. Der Summand fällt hier also weg.

Nachprüfen lässt sich das nicht an einem einzelnen Zahlenwert, sondern nur an der Ordnung des Restes. Mit (20) schrumpft die Schwankung der Gesamtenergie auf einer Bahn mit e = 0,4 wie 1/c⁴ und auf der Kreisbahn wie 1/c⁶ — genau um die Ordnung, die als nächste fehlt (2PN). Nimmt man den vierten Summanden dazu, bleibt stattdessen ein Rest von 1/c²-Ordnung stehen, so groß wie der Korrekturterm selbst. Der Unterschied ist ein Faktor 16 gegen 4 bei jeder Verdopplung von c, und er ist der Unterschied zwischen einer Näherung und einem Rechenfehler.

7.4 Wo der Schalter steht

Überall dort, wo das Szenario mit einer endlichen Lichtgeschwindigkeit rechnet: im Sonnensystem in beiden Fassungen, bei Sgr A* und den S-Sternen, bei den Schwarzen Löchern und auf der leeren Bühne in ihren beiden maßstäblichen Fassungen. Bei „Drei Sonnen“ und im Sternhaufen gibt es ihn nicht — dort ist keine Masse dominierend, und (17) wäre eine Zahl ohne Aussage. Die dimensionslosen Modellfälle (Figur-Acht, pythagoräisches Problem, Lagrange-Dreieck) haben ohnehin kein c; für sie ist die Frage gar nicht gestellt. Der Staub bekommt den Term überall dort mit, wo die Körper ihn bekommen — ein Staubkorn ist der Fall, für den die Näherung eigentlich gemacht ist —, allerdings ohne die Iteration aus (19): Er läuft mit seiner eigenen, gröberen Schrittweite (Abschnitt 4), und nichts hängt von ihm ab.

8Was der Energiefehler misst

Der Ausleser „Energiefehler“ ist das Gütemaß der laufenden Rechnung. Er vergleicht die aktuelle Gesamtenergie (2) nicht mit ihrem Startwert, sondern mit dem Wert, den sie haben sollte:

Abgestrahlte Energie wird also herausgerechnet, und nach einer Verschmelzung wird die Bilanz neu aufgesetzt. Übrig bleibt der reine Verfahrensfehler. Größenordnungen im laufenden Betrieb: 10⁻¹⁴ beim Sonnensystem und beim Lagrange-Dreieck, 10⁻⁶ beim Einspiral der Schwarzen Löcher — dort ist die Reibung (15) selbst die Fehlerquelle. Ab 10⁻³ färbt sich die Anzeige.

Drei Dinge stehen ausdrücklich nicht in dieser Bilanz: die masselosen Testteilchen (sie tragen keine Energie), der Impuls (er ist durch die Bauart von (5) ohnehin erhalten) und der Drehimpuls L = Σ m_i (x_iv_y,i − y_iv_x,i), der zwar berechnet werden kann, aber nirgends angezeigt wird.

9Startwerte

Die gemessenen Szenarien bekommen ihre Startwerte als Bahnelemente und müssen daraus Orte und Geschwindigkeiten machen. Mit dem Bahnparameter p = a(1 − e²) und der wahren Anomalie ν ist im Perizentrumssystem

Anschließend wird um die Länge des Perizentrums ϖ gedreht. Wo ein Körper auf seiner Bahn steht, folgt aus der Zeit: Die mittlere Anomalie M = n(t − t_peri) wächst gleichförmig, die wahre Anomalie nicht. Dazwischen steht die Keplergleichung, die nur numerisch aufzulösen ist — hier mit dem Newton-Verfahren, Abbruch bei 10⁻¹⁴:

Der ungebundene Fall ist nicht der Vollständigkeit halber da: S111 im Galaktischen Zentrum hat nach Gillessen et al. die Exzentrizität 1,09 und verlässt das Zentrum tatsächlich. Bei sehr langgestreckten Bahnen — die S-Sterne gehen bis e = 0,98 — startet die Iteration bei π statt bei M, sonst springt sie weg.

Die Lage der Bahn im Raum

Nach (22) liegt die Bahn noch in ihrer eigenen Ebene. Wie diese Ebene im Raum steht, sagen drei weitere Winkel, und sie werden als drei Drehungen angebracht:

Ω Knotenlänge, i Neigung, ω Argument des Perizentrums.

Ein rückläufiger Umlauf braucht dabei keine Sonderbehandlung: Er folgt aus einer Neigung über 90° von selbst. Halley steht mit seinen 162° also tatsächlich auf dem Kopf, zieht seine Bahn gegen die Planeten und taucht dabei fast zehn AE unter deren Ebene. Die 40 S-Sterne stehen zufällig im Raum, 22 von ihnen rückläufig — was in einer ebenen Rechnung ein Durcheinander gewesen wäre, ist im Raum genau das, was gemessen wurde.

Nachgeprüft wird das auf einem anderen Weg als dem, auf dem es entstanden ist: Aus dem fertigen Zustandsvektor werden über den Drehimpuls h = r × v, den Knotenvektor und den Exzentrizitätsvektor alle sechs Elemente zurückgerechnet. Sie stimmen auf 2·10⁻¹⁴.

Wenn die Quellen zu verschiedenen Zeiten gelten

Die knapp dreitausend Asteroiden der ausführlichen Fassung des Sonnensystems kommen aus einer anderen Quelle als die Planeten, und deren Elemente gelten zu einer anderen Epoche — gut fünf Monate später. Ihre Orte dürfen deshalb nicht einfach übernommen werden. Zusammengeführt wird über die mittlere Anomalie, die als einzige der sechs Größen gleichförmig wächst:

Das ist eine rein keplersche Fortschreibung: Sie tut so, als zöge in diesen fünf Monaten nur die Sonne. Für den Anfangszustand eines Gürtelasteroiden reicht das — danach rechnet die Simulation ohnehin mit allen Kräften, und ein einzelner Asteroid ist als masseloses Testteilchen kein Messobjekt, sondern einer von dreitausend Punkten, die zusammen die Form des Gürtels ergeben. Der Unterschied zwischen den Zeitskalen UTC und TDB, gut siebzig Sekunden, bleibt unberücksichtigt; er verschöbe die mittlere Anomalie um zwei Zehntausendstel Grad.

Die Monde bringen eine zweite Art von Zusammenführung mit: Ihre Elemente beziehen sich nicht auf die Sonne, sondern auf den Mittelpunkt ihres Planeten. In (22) gehört dann die Summe aus Planeten- und Mondmasse in den Zähler, und der fertige Zustandsvektor wird auf den des Planeten addiert. Dass dabei wirklich der Planet steht und nicht die Sonne, wird gegen die von der Quelle mitgelieferte Umlaufzeit geprüft — mit der Sonnenmasse käme für den Erdmond ein Jahr statt eines Monats heraus.

Eine Umrechnung noch, weil sie hübsch ist: Gillessen et al. geben die großen Halbachsen in Bogensekunden am Himmel an. In der Entfernung d (in Parsec) entspricht eine Bogensekunde gerade d Astronomischen Einheiten — bei 8,32 kpc also 8320 AE.

Selbst gesetzte Startwerte

Startwerte müssen nicht aus einer Quelle kommen: Der Editor erlaubt es, Körper einzusetzen, zu entfernen und zu verändern — ein Schwarzes Loch quer durch das Sonnensystem etwa. Für das Verfahren ändert sich dadurch nichts, für seine Voraussetzungen aber zweierlei.

Die Schrittweite ist nicht mehr die gewählte. Jede Aufstellung dieser Seite hat eine Schrittweite, die zu ihr passt (Abschnitt 3.4). Ein eingeworfener Körper kann einem anderen beliebig nahe kommen, und bei einem Vorbeiflug auf ein Hundertstel AE hilft es nicht, die Schrittweite um den Faktor vier zu verkleinern: Die dynamische Zeit der Begegnung liegt dann Größenordnungen unter dem Schritt, und der Fehler ist nicht mehr eine Frage der Ordnung des Verfahrens, sondern der Auflösung. Eine bearbeitete Aufstellung koppelt ihre Schrittweite deshalb immer an die kürzeste dynamische Zeit, auch wo die unbearbeitete das nicht tut. Beim genannten Vorbeiflug macht das den Unterschied zwischen einem Energiefehler von 2·10⁻³ und einem von 5·10⁻⁸.

Der Schwerpunkt ruht nicht mehr. Alle Szenarien werden nach dem Aufbau so verschoben, dass der Schwerpunkt im Ursprung liegt und der Gesamtimpuls null ist — sonst wandert das ganze System aus dem Bild. Eine bearbeitete Aufstellung wird nicht zentriert: Sonst verschöbe jede Eingabe alle anderen Körper, und der eingetippte Ort wäre nicht der, der danach im Feld steht. Der Schwerpunkt trägt dann im Allgemeinen einen Impuls und wandert; weil die Kamera ihm folgt, fällt es nicht auf. Die Erhaltungsgrößen gelten unverändert, sie haben nur nicht mehr die bequemen Werte null.

Am Energiefehler in der Anzeige ändert das nichts: Er misst weiter, was Abschnitt 8 beschreibt. Er ist gerade bei selbst gebauten Aufstellungen die wichtigste Zahl auf dem Bildschirm — wer zwei Körper aus der Ruhe übereinander setzt, sieht dort sofort, dass das Ergebnis nichts bedeutet.

10Nachgerechnet

Jede Aussage dieser Seite lässt sich gegen etwas Bekanntes prüfen — gegen eine analytische Formel, ein Ergebnis aus der Literatur oder einen Messwert. Genau das passiert bei jeder Änderung, und was die Prüfung nicht besteht, wird nicht veröffentlicht.

Die Umlaufzeiten der S-Sterne, die sich aus Halbachse und Zentralmasse ergeben, liegen ein bis zwei Prozent neben den in derselben Arbeit angegebenen. Das ist kein Rechenfehler: Die Größen stammen aus einer gemeinsamen Anpassung mit zusammenhängenden Unsicherheiten, und die Entfernung allein ist nur auf zwei Prozent genau. Wo ein Stern steht, folgt deshalb aus seiner gemessenen Umlaufzeit, nicht aus der Halbachse.

11Grenzen

• Ein flacher Bildschirm. Gerechnet wird im Raum (außer in den drei ebenen Modellfällen), gezeigt wird eine Parallelprojektion, die sich drehen lässt. Anders als bei einer perspektivischen Darstellung bleibt dabei der Maßstab über das ganze Bild derselbe, was den Maßstabsbalken erst gültig macht; dafür fehlt der Größeneindruck, den ein Fluchtpunkt gibt. Tiefe steckt nur in der Punktgröße des Staubs und darin, was vor was gezeichnet wird.
• Relativität nur in erster Ordnung, und nur um eine Masse. Mehr als der Reibungsterm (15) und der 1PN-Zusatz (17) ist nicht zu haben — letzterer läuft im Sonnensystem und bei Sgr A* von Anfang an mit. Merkurs Perihel wandert damit richtig, und die Bahnen der S-Sterne drehen sich in der gemessenen Größe. Alles Weitere fehlt: keine Lichtlaufzeiten — die Anziehung wirkt augenblicklich —, keine Lichtablenkung, keine zweite post-newtonsche Ordnung, keine EIH-Gleichungen, keine Gegenkraft auf die Zentralmasse. Wo es keine beherrschende Masse gibt — bei „Drei Sonnen“, im Sternhaufen und auch bei den beiden Schwarzen Löchern —, gibt es den Schalter gar nicht erst. Und der schwarze Punkt in der Mitte bleibt ein newtonscher Massenpunkt mit einem Radius, der zufällig 2GM/c² heißt.
• Direktsummation. Der Aufwand wächst wie N²; Baum- oder Multipolverfahren gibt es nicht. Bei höchstens 42 Körpern lohnt sich das nicht — der Aufwand steckt im Staub, und der ist linear. Der Übergang von zwei auf drei Dimensionen hat ihn übrigens nur um ein Viertel erhöht und nicht um die Hälfte: Wurzel und Division fallen je Paar einmal an, gleichgültig wie viele Komponenten der Abstandsvektor hat.
• Ein Stride für allen Staub, siehe 4.4. Blockstufen wären der nächste Schritt.
• Gemeinsame Schrittweite für alle Körper. Merkur bestimmt sie im Sonnensystem für Neptun mit, der dadurch mit rund 400 000 Schritten je Umlauf gerechnet wird statt mit den nötigen 600. Das kostet hier fast nichts (45 von 14 045 Paaren), wäre in einem großen Haufen aber der entscheidende Punkt. Der elegante Ausweg für ein System mit dominanter Zentralmasse wäre ein gemischtes Verfahren nach Wisdom und Holman: die Keplerbahn exakt lösen und nur die Störungen numerisch integrieren. Für die chaotischen Szenarien ohne Zentralkörper hilft das nicht.
• Rundungsfehler. Sie sind der Grund, warum die Figur-Acht nach vielen Umläufen ihre Form verliert und das Lagrange-Dreieck von selbst auseinanderbricht. Bei der Acht ist das ein Schönheitsfehler, beim Dreieck genau der Punkt: Die Lösung ist exakt, aber instabil.

Literatur

1. Yoshida, H. 1990, Construction of higher order symplectic integrators, Phys. Lett. A 150, 262.
2. Hairer, E., Lubich, C., Wanner, G. 2006, Geometric Numerical Integration, 2. Aufl., Springer — Rückwärtsfehleranalyse und modifizierter Hamiltonian, Gl. (4).
3. Aarseth, S. 2003, Gravitational N-Body Simulations, Cambridge University Press — individuelle Schrittweiten und Blockstufen.
4. Wisdom, J., Holman, M. 1991, Symplectic maps for the n-body problem, AJ 102, 1528.
5. Peters, P. C., Mathews, J. 1963, Phys. Rev. 131, 435 — Gl. (14).
6. Peters, P. C. 1964, Phys. Rev. 136, B1224 — Gl. (16).
7. Abbott, B. P. et al. 2016, Phys. Rev. Lett. 116, 061102 — GW150914, Massen des Doppelsystems.
8. Einstein, A., Infeld, L., Hoffmann, B. 1938, The Gravitational Equations and the Problem of Motion, Ann. Math. 39, 65 — die vollständigen post-newtonschen Bewegungsgleichungen, an deren Stelle hier (17) steht.
9. Weinberg, S. 1972, Gravitation and Cosmology, Wiley — Bahn eines Testteilchens im Schwarzschildfeld in erster post-newtonscher Ordnung, Gl. (17) und (18).
10. Blanchet, L. 2014, Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries, Living Rev. Relativity 17, 2, doi:10.12942/lrr-2014-2 — die allgemeine 1PN-Energie zweier Körper, aus der (20) im Grenzfall ν → 0 hervorgeht.
11. Will, C. M. 2014, The Confrontation between General Relativity and Experiment, Living Rev. Relativity 17, 4, doi:10.12942/lrr-2014-4 — die gemessene Periheldrehung Merkurs und ihr Anteil an den beobachteten 574″ je Jahrhundert.
12. GRAVITY Collaboration 2020, Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole, A&A 636, L5, doi:10.1051/0004-6361/202037813 — 12,1′ je Umlauf, gemessen als das 1,10 ± 0,19-fache der Vorhersage.
13. Gillessen, S. et al. 2017, ApJ 837, 30, doi:10.3847/1538-4357/aa5c41 — Bahnen der S-Sterne, Masse und Entfernung von Sgr A*.
14. Chenciner, A., Montgomery, R. 2000, A remarkable periodic solution of the three-body problem in the case of equal masses, Ann. Math. 152, 881.
15. Burrau, C. 1913, Astron. Nachr. 195, 113; Szebehely, V., Peters, C. F. 1967, AJ 72, 876 — das pythagoräische Problem und seine Lösung.
16. JPL Horizons, NASA/JPL — Bahnelemente, GM und Radien des Sonnensystems zur Epoche 2026-01-01.