Home| Inhalt| Theorie| Experimente| Versionen| PDF
| de

2.8.9 Numerische Berechnung der Feldstärke beliebig bewegter Punktladungen

2.8.9.1 Herleitung des Algorithmus

In den vorangegangenen Abschnitten wurde der Quantinodruck für einige Spezialfälle analytisch untersucht und berechnet. Für allgemeine Situationen, wie z.B. bei den besonders interessanten Mehrteilchenwechselwirkungen, ist eine symbolische Lösung des Integrals (2.8.3.15) in aller Regel unmöglich. Aus diesem Grund bedarf es einer Methodik, um den Quantinodruck, welcher die retardierende elektromagnetische Feldstärke einer beliebig bewegten Punktladung vollständig beschreibt, numerisch auszuwerten. Die Chance, dass es gelingt ein brauchbares und effizientes Verfahren zu finden, ist nicht schlecht. Immerhin ist der Quantinodruck durch Gleichung (2.8.3.15) bereits vollständig gegeben und muss nicht erst aus einem Differentialgleichungssystem abgleitet werden.

Zur elektrischen Feldstärke gelangt man mit Hilfe der Gleichung (2.8.5.3):
$$\vec{F} = \frac{\sigma_e\,q_d\,q_s}{e^2}\,\vec{\mathcal{P}}_e(\vec{r},\vec{v}).$$ (2.8.9.1.1)
Hierbei ist $\vec{F}$ die Kraft, die eine Quellladung $q_s$ auf eine Zielladung $q_d$ ausübt. Da die elektrische Feldstärke $\vec{E}$ durch
$$\vec{E} := \frac{\vec{F}}{q_d}$$ (2.8.9.1.2)
definiert ist, folgt
$$\vec{E} = \frac{\sigma_e\,q_s}{e^2}\,\vec{\mathcal{P}}_e(\vec{r},\vec{v}).$$ (2.8.9.1.3)
Mit Hilfe der Gleichung (2.8.5.4) lässt sich die unbekannte Konstante $\sigma_e$ ersetzen. Man gelangt zu
$$\vec{E} = \frac{q_s\,\mu_0}{a_c\,m_{pho}}\,\vec{\mathcal{P}}_e(\vec{r},\vec{v}).$$ (2.8.9.1.4)
Ein Vergleich mit Gleichung (2.8.3.15) zeigt, dass die vom Zahlenwert her unbekannten Konstanten $a_c$ und $m_{pho}$ nicht benötigt werden, da sie sich herauskürzen.

Da sich Quantinodruck und elektrische Feldstärke nur durch einen konstanten Faktor unterscheiden, ist es letztlich ohne Bedeutung, welche der beiden Größen berechnet wird. Einen Algorithmus zur numerischen Berechnung findet man, indem man die Ausdrücke (2.8.3.15), (2.8.3.11) und (2.8.3.12) einige Momente betrachtet. Man erkennt, dass der auf ein Objekt mit der Bahnkurve $\vec{r}_d(t)$ wirkende Quantinodruck nur von dessen Ort zum Zeitpunkt $t$, sowie von dessen Geschwindigkeit $\dot{\vec{r}}_d(t)$ abhängt. Mit anderen Worten: Die Historie des Empfängers ist irrelevant. Aus diesem Grund kann man $\vec{r}_d(t)$ einfach mit $\vec{r}$ und $\dot{\vec{r}}_d(t)$ mit $\vec{v}$ abkürzen. Für die aussendende Ladung mit der Bahnkurve $\vec{r}_s(\tau)$ gilt jedoch das Gegenteil. Hier ist die Historie von entscheidender Bedeutung und man benötigt alle Orte $\vec{r}_s(\tau)$, sowie Geschwindigkeiten $\dot{\vec{r}}_s(\tau)$ für alle $\tau$ vom "Anbeginn der Zeit" bis zur Gegenwart $\tau = t$.

Dass der Quantinodruck einer Ladung von der gesamten Historie der Ladung abhängt, ist von weitreichender Bedeutung. Letztlich entstehen dadurch faszinierende physikalische Effekte, wie z.B. die Fähigkeit eines Dipols mit seiner eigenen Welle zu interferieren. Sollte der Leser dabei an die Quantenmechanik denken, so liegt er hierbei vollkommen richtig. Auch das zweite Newtonsche Gesetz $Kraft = Masse \cdot Beschleunigung$ hat letztlich seine Ursache in dieser Abhängigkeit der elektromagnetischen Kraft von der Historie der aussendenden Ladung.

Die Abhängigkeit des Quantinodrucks von der Vergangenheit der Quellladung ist aber auch eine praktische Schwierigkeit, da man auf diese Weise gezwungen wird, unendlich viele Anfangsbedingungen zu berücksichtigen. Glücklicherweise wird der Einfluss der Historie umso kleiner, je weiter diese Vergangenheit zurück liegt. Dies wird zum einen durch den Term $t-\tau$ unter dem Bruchstrich im Integral der Quantinodichte (2.8.3.10) bewirkt. Zum anderen tritt auch im Parameter der Funktion $\Gamma$ ein $t-\tau$ unterhalb des Bruchstriches auf. Und da die Funktionswerte der Funktion $\Gamma(w)$ im Wesentlichen linear zur Geschwindigkeit $w$ sind, bewirkt dies für größer werdendes $t-\tau$ eine kleiner werdende Gewichtung durch $\Gamma$. In der Praxis reicht es daher aus, eine gewisse Historie zu berücksichtigen und alles zuvor Gewesene zu vernachlässigen, indem man das Integral nach unten begrenzt.

Berechnet man die Bahnkurve eines Objektes unter Einfluss einer Kraft numerisch (z.B. mit dem Runge-Kutta-Verfahren), so erhält man eine Folge von Orten
$$\vec{r}_i = \{\vec{r}_0, \vec{r}_1, \vec{r}_2, ...,\vec{r}_n\}$$ (2.8.9.1.5)
und Geschwindigkeiten
$$\vec{v}_i = \{\vec{v}_0, \vec{v}_1, \vec{v}_2, ...,\vec{v}_n\}$$ (2.8.9.1.6)
für die Zeitpunkte
$$\tau_i = \{\tau_0, \tau_1, \tau_2, ...,\tau_n\}.$$ (2.8.9.1.7)
Eine naheliegende Grundidee für eine numerische Lösung besteht darin, die Orte $\vec{r}_s(\tau)$ und Geschwindigkeiten $\dot{\vec{r}}_s(\tau)$ der aussendenden Ladung stückweise linear zu interpolieren und als Stützstellen der Interpolation die Werte (2.8.9.1.5), (2.8.9.1.6) und (2.8.9.1.7) zu verwenden. Mit diesem Ansatz folgen die Näherungen
$$\vec{r}_s(\tau) \approx \sum\limits_{j=0}^{n-1} \intfunc_{\tau_{j}}^{\tau_{j+1}}(\tau)\left(\vec{A}_j + \tau\,\vec{B}_j\right)$$ (2.8.9.1.8)
und
$$\dot{\vec{r}}_s(\tau) \approx \sum\limits_{j=0}^{n-1} \intfunc_{\tau_{j}}^{\tau_{j+1}}(\tau)\left(\vec{C}_j + \tau\,\vec{D}_j\right),$$ (2.8.9.1.9)
wobei zur Verkürzung der Schreibweise die Konstanten
$$\vec{A}_j := \frac{\vec{r}_{j}\,\tau_{j+1}-\vec{r}_{j+1}\,\tau_j}{\tau_{j+1}-\tau_j},\quad \vec{B}_j := \frac{\vec{r}_{j+1}-\vec{r}_{j}}{\tau_{j+1}-\tau_j}$$ (2.8.9.1.10)
und
$$\vec{C}_j := \frac{\vec{v}_{j}\,\tau_{j+1}-\vec{v}_{j+1}\,\tau_j}{\tau_{j+1}-\tau_j},\quad \vec{D}_j := \frac{\vec{v}_{j+1}-\vec{v}_{j}}{\tau_{j+1}-\tau_j}$$ (2.8.9.1.11)
eingeführt wurden. Die Interpolationen (2.8.9.1.8) und (2.8.9.1.9) lassen sich nun in die Definitionen (2.8.3.11) und (2.8.3.12) einsetzen. Man erhält die Näherungen
$$\vec{w}(\tau) \approx \sum\limits_{j=0}^{n-1} \intfunc_{\tau_{j}}^{\tau_{j+1}}(\tau)\,\vec{w}_j(\tau)$$ (2.8.9.1.12)
mit
$$\vec{w}_j(\tau) := \frac{\vec{r} - \vec{A}_j -\tau\,\vec{B}_j}{\tau_n-\tau} - \vec{C}_j - \tau\,\vec{D}_j$$ (2.8.9.1.13)
und
$$\vec{u}(\tau) \approx \sum\limits_{j=0}^{n-1} \intfunc_{\tau_{j}}^{\tau_{j+1}}(\tau)\,\vec{u}_j(\tau)$$ (2.8.9.1.14)
mit
$$\vec{u}_j(\tau) := \frac{\vec{r} - \vec{A}_j -\tau\,\vec{B}_j}{\tau_n-\tau} - \vec{v}.$$ (2.8.9.1.15)
Setzt man dieses in die Gleichung (2.8.3.15) ein, so zerfällt das Integral in eine Summe von Teilintegralen, die sich jeweils über einen einzelnen Zeitschritt erstrecken. Man erhält
$$\vec{\mathcal{P}}_e = \frac{a_c\,m_{pho}}{8\,\pi} \,\sum\limits_{i=0}^{n-1}\,\int\limits_{\tau_i}^{\tau_{i+1}}\,\vec{K}_i(\tau)\,\d{\tau}$$ (2.8.9.1.16)
mit
$$\vec{K}_i(\tau) := \frac{u_i(\tau)^2\,\vec{w}_i(\tau)\,\sgn\left(\vec{w}_i(\tau)\,\vec{u}_i(\tau)\right)\,\Gamma\left(w_i(\tau)\right) \,\intfunc_{0}^{c}\left(u_i(\tau)\right)}{(\tau_n-\tau)^3\,w_i(\tau)^3}.$$ (2.8.9.1.17)


Die Lösung der Integrale des Typs
$$\int\limits_{\tau_i}^{\tau_{i+1}}\,\vec{K}_i(\tau)\,\d{\tau},$$ (2.8.9.1.18)
kann mit üblichen numerischen Integrationsverfahren (z.B. Simpsonregel) erfolgen. Dabei ist jedoch darauf zu achten, das der Integrand $\vec{K}_i(\tau)$ Unstetigkeiten und Singularitäten enthalten kann, wobei ersteres zum Einen daher rührt, dass der Term $\intfunc_{0}^{c}\left(u_i(\tau)\right)$ verschwindet, wenn der Wert von $u_i(\tau)$ die Lichtgeschwindigkeit $c$ überschreitet. Zum anderen kann auch die Vorzeichenfunktion $\sgn$ zu Sprüngen im Funktionswert von $\vec{K}_i(\tau)$ führen. Singularitäten können hingegen auftreten, wenn der Term $w_i(\tau)$ verschwindet und die Funktion $\Gamma$ einen linearen Term enthält. Einen konstanten Term enthält sie grundsätzlich nicht, wie bereits gezeigt wurde. Weiterhin würde eine Singularität auftreten, wenn der Term $\tau_n-\tau$ Null wird. Dieses entspricht aber der Frage nach dem Wert des Quantinodrucks am Ort der Ladung selbst. Für diesen Ort ist die Gesetzmäßigkeit (2.8.3.15) aber sowieso nicht gültig, da die Modellannahme einer kontinuierlichen Quantinodichte in unmittelbarer Nähe der emittierenden Ladung nicht erfüllt ist und hier der Effekt der Quantisierung dominiert.

Singularitäten sind bei der Berechnung des Integrals also eher die Ausnahme. Unstetigkeiten treten jedoch regelmäßig auf und müssen bei der Berechnung systematisch berücksichtigt werden, da numerische Quadratur üblicherweise einen stetigen Integranden voraussetzt. Unstetigkeiten sind jedoch kein Problem, wenn bekannt ist, wo sie sich befinden, da es in diesem Fall möglich ist, die Integrationsintervalle so zu wählen, dass der Integrand keine Unstetigkeiten enthält.

Der erste Problemterm im Integrand $\vec{K}_i(\tau)$ ist $\intfunc_{0}^{c}\left(u_i(\tau)\right)$. Dieser ist nur für $u_i(\tau) = \Vert\vec{u}_i(\tau)\Vert = c$ unstetig. Wegen Gleichung (2.8.9.1.14) folgt
$$\left(\frac{\vec{r} - \vec{A}_i -\tau\,\vec{B}_i}{\tau_n-\tau} - \vec{v}\right)\cdot\left(\frac{\vec{r} - \vec{A}_i -\tau\,\vec{B}_i}{\tau_n-\tau} - \vec{v}\right) = c^2.$$ (2.8.9.1.19)
Nach etwas Umformen erhält man die quadratische Gleichung
$$\left\Vert\vec{E}_i + \tau\,\vec{F}_i\right\Vert^2 = c^2 (\tau_n-\tau)^2$$ (2.8.9.1.20)
mit
$$\vec{E}_i := \vec{r} - \vec{A}_i - \vec{v}\,\tau_n\quad\mathrm{und}\quad\vec{F}_i := \vec{v} -\vec{B}_i.$$ (2.8.9.1.21)
Die beiden Lösungen der Gleichung (2.8.9.1.20) lauten
$$\begin{split}\tau_{1/2} = & \frac{1}{c^2 - F_i^2} \left(c^2\,\tau_n + \vec{E}_i\,\vec{F}_i \pm \right. \\ & \left. \sqrt{\left(c^2\,\tau_n + \vec{E}_i\,\vec{F}_i \right)^2 + \left(c^2 - F_i^2\right)\,\left(E_i^2 - c^2\,\tau_n^2\right)} \right),\end{split}$$ (2.8.9.1.22)
wobei eine Lösung nur dann von Interesse ist, wenn diese in der Vergangenheit liegt. Dieses ist aber nur für $\tau_2$ gegeben, wie man sich überlegen kann. Aus diesem Grund ist es sehr leicht möglich, vor jeder numerischen Integration in der Summe (2.8.9.1.16) zu überprüfen, ob $u_{i}(\tau_i)$ oder $u_{i}(\tau_{i+1})$ größer ist als $c$. Sind beide kleiner, so kann die Integration über das komplette Intervall erfolgen. Ist einer der beiden größer $c$, so ist das Intervall auf den Bereich zu begrenzen, in dem $w_{i}(\tau) \leq c$ gilt. Sind sogar beide größer, so kann die Integration komplett entfallen, da die Intervallfunktion im gesamten Integrationsgebiet Null liefert.

Die zweite Quelle für Unstetigkeit ist der Term $\sgn\left(\vec{w}_i(\tau)\,\vec{u}_i(\tau)\right)$. Das bedeutet, dass für jedes Integrationsintervall zu prüfen ist, ob die Gleichung
$$\left(\frac{\vec{r} - \vec{A}_i -\tau\,\vec{B}_i}{\tau_n-\tau} - \vec{C}_i - \tau\,\vec{D}_i\right)\cdot\left(\frac{\vec{r} - \vec{A}_i -\tau\,\vec{B}_i}{\tau_n-\tau} - \vec{v}\right) = 0$$ (2.8.9.1.23)
eine oder mehrere Lösungen enthält, die innerhalb der Integrationsgrenzen liegen. Durch Multiplikation beider Seiten mit $(\tau_n-\tau)^2$ und etwas Umsortieren folgt die kubische Gleichung
$$\left(\vec{G}_i + \tau\,\vec{H}_i + \tau^2\,\vec{D}_i\right)\cdot\left(\vec{E}_i + \tau\,\vec{F}_i\right) = 0$$ (2.8.9.1.24)
mit
$$\vec{G}_i := \vec{r} - \vec{A}_i - \vec{C}_i\,\tau_n\quad\mathrm{und}\quad\vec{H}_i := \vec{C}_i - \vec{B}_i-\vec{D}_i\,\tau_n.$$ (2.8.9.1.25)
Die drei Nullstellen ließen sich beispielsweise mit den Cardanischen Formeln berechnen. Sollten tatsächlich Lösungen existieren, die innerhalb der Integrationsgrenzen liegen, so könnte man das Integral in entsprechende Teilintegrale zerlegen. Da die Berechnung der Nullstellen aufwendig ist, bietet es sich an, den Integrationsfehler zu kontrollieren, indem man jedes Integral einmal für das komplette Intervall und einmal als Summe zweier gleich großer Teilintervalle berechnet. Beides sollte, sofern die Funktion im Intervall hinreichend gutartig ist, also keine Unstetigkeiten enthält, zum beinahe gleichen Ergebnis führen. Falls nicht, wird die Integration weiter rekursiv in Teilintervalle zerlegt, bis der Fehler hinreichend klein ist. Diese Vorgehensweise ist auch in Hinblick auf den generellen numerischen Fehler von Vorteil, also auch dann, wenn keine Unstetigkeiten vorliegen.

Eine Beispielimplementierung des Algorithmus in C++ findet sich hier. Im nachfolgenden Abschnitten wird die Verwendung anhand einiger Beispiele demonstriert und diskutiert.

2.8.9.2 Verifizierung anhand einer ruhenden Ladung

Der im vorangegangenen Abschnitt hergeleitete Algorithmus und dessen Implementierung wird zunächst anhand einfacher Spezialfälle mit bekannten analytischen Lösungen verifiziert. Der einfachste Fall ist das elektrische Feld einer ruhenden Einheitsladung. Abbildung 2.8.9.2.1 zeigt das numerisch berechnete Feld einer solchen Ladung in einem Gebiet von $2 \times 2\,\mu m^2$. Daneben ist in Abbildung 2.8.9.2.2 die hundertfach verstärkte Differenz zum Feld dargestellt, welches mit dem Coulombgesetz ermittelt wurde. Der numerische Fehler kann durch eine andere Wahl der Puffer-Größe (Parameter in Zeile 51) weiter verringert werden. Da dadurch allerdings der Rechen- und Speicheraufwand ansteigt, ist dieses nur in begrenztem Umfang möglich.

Abbildung 2.8.9.2.1: Numerisch berechnetes Feld einer ruhenden positiven Einheitsladung.
Abbildung 2.8.9.2.2: Das hundertfach verstärkte Residuum.

Der Glue-Code der Simulation (also ohne diesen Code) kann hier heruntergeladen werden. Die Simulation erfolgte unter Linux, wobei zur Darstellung der Feldlinien die Bibliothek Cairo verwendet wurde. Es muss im Übrigen erwähnt werden, dass die numerische Berechnung der Felder nur für einen hinreichend gefüllten Puffer gute Ergebnisse liefert. Das heißt, es müssen erst einige Zeitschritte vergehen, bevor das Integral der Formel (2.8.3.15) genügend weit in die Vergangenheit reicht.

2.8.9.3 Das Feld einer schnell bewegten Ladung

Abbildung 2.8.9.3.1: Numerisch berechnetes Feld einer positiven Einheitsladung aus Sicht eines sich mit der Geschwindigkeit $0.99\,c$ in x-Richtung bewegenden Empfängers.
Eine wesentlich interessantere Fragestellung als die eben untersuchte ist, welche Form das Feld einer Einheitsladung aus Sicht eines Empfängers hat, der sich beinahe mit Lichtgeschwindigkeit bewegt. In der maxwellschen Elektrodynamik bleibt diese Frage offen. Im Rahmen der Quantinotheorie kann sie allerdings sogar analytisch beantwortet werden, da mit Formel (2.8.4.13) für den unbeschleunigten Fall ein Weg zur Berechnung des Feldes für beliebige Relativgeschwindigkeiten kleiner $c$ gegeben ist. Abbildung 2.8.9.3.1 zeigt das Feld einer ruhenden Einheitsladung aus Sicht eines sich mit der Geschwindigkeit $0.99\,c$ nach rechts bewegenden Empfängers. Das gleiche Feld erhält man, wenn man eine analytische Berechnung durchführt.

Wie zu sehen ist, verliert das Feld durch die hohe Geschwindigkeit des Empfängers seine kreisrunde Form und wird stattdessen linsenförmig. Dieser Effekt verstärkt sich, je weiter die Geschwindigkeit zunimmt. Im Grenzfall, also wenn der Empfänger die Lichtgeschwindigkeit erreicht, verschwindet das Feld völlig. Dieses hat entscheidende Auswirkungen, denn es zeigt, dass es grundsätzlich unmöglich ist, eine elektrische Ladung mittels einer ruhenden Anordnung, wie zum Beispiel einem Teilchenbeschleuniger, auf Überlichtgeschwindigkeit zu bringen. Der Effekt wird im Abschnitt 2.7.7 anschaulich und ausführlich beschrieben.

Der Sourcecode der Simulation kann von hier herunter geladen werden. Bei der Berechnung ist zu beachten, dass die Puffergröße sehr groß gewählt werden muss. Der Grund besteht bildlich gesprochen darin, dass der Raum zunächst mit langsamen Quantinos gefüllt werden muss.

2.8.9.4 Das Feld eines strahlenden Dipols

Abbildung 2.8.9.4.1: Der Dipol beginnt zum Zeitpunkt 0 mit einer Frequenz von einem Peta-Hertz (UV-Licht) zu schwingen.
In den vorangegangenen beiden Beispielen wurden statische Felder berechnet. Der Algorithmus ist jedoch nicht auf solche begrenzt, sondern eignet sich ganz im Gegenteil vor allem für dynamische Analysen. Abbildung 2.8.9.4.1 zeigt die Ausbreitung des elektrischen Feldes eines aus zwei Einheitsladungen bestehenden Dipols, der bis zum Zeitpunkt $t=0$ ruht. Da die Einheitsladungen entgegengesetzte Vorzeichen haben, neutralisieren sie sich bis zum Start der Berechnung perfekt. Von der Wirkung her kann der umgebende Raum daher bis zum Zeitpunkt $t=0$ als quantinofrei angesehen werden. Durch die Schwingung kommt es für $t>0$ zu Dichteschwankungen, die sich dann in Form eines sich mit Lichtgeschwindigkeit ausbreitenden elektrischen Feldes bemerkbar machen. In der Animation gut zu sehen ist, wie sich das Feld von Wellenfront zu Wellenfront umpolt. Interessant ist ebenfalls, dass sich in der Nähe des Dipols das vom statischen Dipol her bekannte typische Dipolfeld zeigt. Dies führt zu Longitudinalwellen (x-Achse), als auch zu Transversalwellen (y-Achse). Da keine weiteren Dipole vorhanden sind, breiten sich beide Wellentypen gleichmäßig in alle Richtungen aus. Die Amplitude der Schwingung verringert sich aus diesem Grund proportional zum Quadrat des Abstandes.

Auch zu diesem Feld existiert eine analytische Näherungslösung. Ein Vergleich mit Abbildung 2.8.7.1 zeigt, dass das Ergebnis der numerischen Berechnung gut mit dem theoretisch ermittelten Verhalten übereinstimmt. Der Source-Code der Berechnung ist hier erhältlich.

2.8.9.5 Das Feld beim Bohrschen Atommodell

Abbildung 2.8.9.5.1: Das Feld beim Bohrschen Atommodell. Die Zeichenfläche hat eine Größe von 200 mal 200 nm.
Das Feld des Hertzschen Dipols ist bereits gut aus der klassischen Elektrodynamik bekannt. Weniger bekannt dürfte das Feld sein, welches man im Bohrschen Atommodell erhalten würde, wenn dieses denn korrekt wäre. Da die Anordnung als starrer, rotierender Dipol angesehen werden kann, die Ausbreitung der Kraft aber durch die Lichtgeschwindigkeit begrenzt wird, sollte das Feld diesmal spiralförmig erscheinen. Abbildung 2.8.9.5.1 zeigt das Ergebnis der numerischen Berechnung. Das Elektron umkreist das Proton im Abstand von ca. $5.3\,10^{-11}\,m$. Wie erwartet, ist das resultierende Feld spiralförmig. Die Wellenlänge der abgegebenen Strahlung entspricht der Umlaufzeit des Elektrons (etwa $0.15\,fs$) mal der Lichtgeschwindigkeit, also ca. $45\,nm$, was weicher Röntgenstrahlung entspricht.

Der Source-Code der Berechnung findet sich hier.