Quantenmechanik/ Boson-Fermion

Identische Teilchen

Ein System aus mehreren Teilchen hat quantenmechanisch nicht etwa mehrere Wellenfunktionen, sondern eine gemeinsame im Tensorprodukt der Einteilchen-Hilbert-Räume. Nur wenn für Teilchen A bei x und B bei y keine Wechselwirkung $H_{w}(x,y)$ im Hamilton-Operator auftritt, ist ein Separations-Ansatz $\psi _{A}(x)\psi _{B}(y)$ erfolgreich. Allgemein gibt es Linearkombinationen von solchen Produkten, also verschränkte Teilchen: $\Psi (x,y)=\psi _{A}(x)\psi _{B}(y)+\phi _{A}(x)\phi _{B}(y)+...$

Teilchen wie Protonen, Elektronen oder ganze Atome haben nicht das geringste Merkmal, mit dem man sie individuell verfolgen kann. Gleiche Teilchen sind ununterscheidbar. Sie haben keine Pfade wie in der klassischen Physik. Sind Teilchen A und B identisch wie die zwei Elektronen des Helium-Atoms, dann stellt sich schon ohne Wechselwirkung eine Frage.

Die Wellen $\psi (x)\phi (y)$ und $\psi (y)\phi (x)$ sind degeneriert, es könnte also sein, dass ein ganzes Kontinuum von Funktionen

\Psi (x,y)=\alpha \psi (x)\phi (y)+\beta \psi (y)\phi (x);\;|\alpha |^{2}+|\beta |^{2}=1

die selbe Physik beschreibt. Der Permutations-Operator $P:\;(P\Psi )(x,y)=\Psi (y,x)$ hätte dann keine definierten Eigenwerte?

Oder aber P hat eine einfache Darstellung, $(P\Psi )=\pm \Psi$ , das heißt die Zweiteilchen-Funktion ist entweder symmetrisch oder antisymmetrisch beim Vertausch der Koordinaten, also mit Eigenwert 1 oder -1 zu P. Auf jeden Fall hätte P einen Eigenwert auf dem Einheitskreis, so dass die observablen Betragsquadrate invariant sind beim Paarvertausch?

Die Antwort der Natur ist kategorisch und wurde von Pauli formuliert. Die Mehrteilchen-Zustände N identischer Teilchen mit halbzahligem Spin, wie der Elektronen, haben als Darstellung der Permutationsgruppe den Faktor -1 bei Paarvertauschung. Teilchen mit dieser Antisymmetrie sind Fermionen. Teilchen mit ganzzahligem Spin, die Bosonen, haben den P-Eigenwert 1, sie haben also vertauschungs-invariante N-Teilchen-Zustände.

Unmittelbare Konsequenz ist das Pauli-Verbot: keine zwei Fermionen passen auf den selben Quantenzustand: $\phi (x)\phi (y)=-\phi (x)\phi (y)=0$ . Mittelbare Konsequenz ist der Aufbau der Atome der Kernladung Z in Schalen. Die Elektronen füllen nacheinander die aufsteigende Reihe der Energie-Niveaus, die bei der Berechnung eines Elektrons im Coulomb-Potenzial herauskommen.

Die Ursache des Spin-Statistik-Prinzips ist in der relativistischen Quantenfeldtheorie zu suchen. Die identischen Teilchen sind alle die elementaren Anregungen des selben Quantenfeldes. Der Spin ergibt sich aus einer Klassifizierung der Darstellungen der Lorentz-Gruppe. Die Vielteilchen-Zustände sollen einen Hilbert-Raum aufspannen (positives Skalarprodukt) und einen stabilen Grundzustand haben (beschränktes Spektrum der Energie). Nur die experimentell untermauerte Korrespondenz von Spin und Vertauschungs-Symmetrie ist dann widerspruchsfrei möglich.

Statistik

Angenommen, es gebe $i=1...I$ Plätze oder Einzelniveaus für einzelne Teilchen. Jedem Teilchen gehöre ein Vektorraum der Dimension I mit einer Basis von Einteilchen-Zuständen $\{B(i)\}$ .

Naiv besteht ein Basissystem der Zustände von $n=0...N$ verfolgbaren Teilchen aus allen Abbildungen $\{1...N\}\to \{1...I\}$ . Es gibt $I^{N}$ Vektoren in der Basis des gewöhnlichen Tensorproduktes $T(N)$ . Sie beschreibt das System:

\{B(i_{1})\otimes B(i_{2})\otimes ...\otimes B(i_{N})\}\quad

erzeugt

T(N)

.

Die Dimension ist die Zahl der linear unabhängigen Zuständsvektoren.

Mit Fermionen sind für jedes i nur zwei Besetzungen erlaubt, 0 oder 1. Falls $I>N$ , sind die unabhängigen Zustände die Menge aller Teilmengen der Größe N, insgesamt ${\binom {I}{N}}$ Zustände. Soviel Möglichkeiten gibt es, N Plätze aus den I verfügbaren auszuwählen. Der Vektorraum des Systems wird aufgespannt von allen total antisymmetrischen Tensorprodukten, ist also ein Teilraum $A(N)\subset T(N)$ , definiert wie folgt.

Sei $P_{N}$ die Gruppe aller N! Permutationen von $\{1...N\}$ . Die Gruppe hat zwei elementare Homomorphismen zu den ganzen Zahlen,

den trivialen, ${\text{one}}(p)=1\quad \forall p\in P_{N}$ und
das Permutations-Vorzeichen $\operatorname {sgn}(p)=-1$ genau dann wenn

p aus einer ungeraden Zahl von Nachbar-Vertauschungen gebaut ist, 1 sonst.

Alle Paarvertauschungen, nicht nur die Transpositionen von Nachbarn, haben negatives $sgn()$ . Eine beliebige Permutation kann in ein Produkt aus Zyklen zerlegt werden. Ein Zyklus der Größe $G$ hat das Vorzeichen $-(-1)^{G}$ .

Der Antisymmetrie-Operator A auf $T(N)$ ist die lineare Fortsetzung von:

A(B(i_{1})\otimes B(i_{2})\otimes ...\otimes B(i_{N}))=\sum _{p\in P_{N}}{\frac {\operatorname {sgn}(p)}{N!}}\;[B(i_{p1})\otimes B(i_{p2})\otimes ...\otimes B(i_{pN})]

Der Vektorraum des Fermionsystems ist die Bildmenge $A(N):=A(T(N))$ . A ist eine Projektion, $A^{2}=A$ . Als Projektion hat A die Tendenz, Vektoren zu "verkürzen" und erhält daher nicht die Normen.

In der Ortsdarstellung geht die Vielteilchen-Welle aus Produkten hervor. Man wende den Operator A an auf $W(x_{1},...,x_{N})=B_{1}(x_{1})B_{2}(x_{2})...B_{N}(x_{N})$ .

(AW)(x_{1},...,x_{N})=\sum _{p\in P_{N}}\operatorname {sgn}(p)/N!\;[B_{1}(x_{p1})B_{2}(x_{p2})...B_{N}(x_{pN})]

Eine für orthogonale Produkte norm-erhaltende Variante hat folgende Form:

{\bar {A}}\,W={\sqrt {N!}}(AW)(x_{1},...,x_{N})={\frac {1}{\sqrt {N!}}}\;\det {\begin{vmatrix}B_{1}(x_{1})&\ldots &B_{N}(x_{1})\\\vdots &\ldots &\vdots \\B_{1}(x_{N})&\ldots &B_{N}(x_{N})\end{vmatrix}}

.

Diese Schreibweise einer N-Fermion-Wellenfunktion ist bekannt under dem Namen Slater-Determinante. Alle beteiligten Einteilchen-Wellen müssen linear unabhängig sein, damit sie nicht verschwindet. Die Slater-Determinanten eines Einteilchen-Orthonormalsystems bilden ein Vielteilchen-Orthonormalsystem.

Die Argumente $x_{n}$ können allgemeiner Tupel sein aus Ortskoordinaten, Spin-Indizes und anderen Quantenzahlen. Bei N Elektronen können eventuell nur die Spin-Faktoren der Wellenfunktion antisymmetrisch sein; dann werden die Orts-Faktoren symmetrisch für die geforderte totale Antisymmetrie. Alle möglichen Kombinationen der Symmetrien mit Tupeln können mit der Darstellungstheorie der Permutationsgruppen (Young-Tableaus etc) klassifiziert werden.

Bei Bosonen darf man N "Ziehungen mit Zurücklegen" aus den I Einzelniveaus machen, was der Zahl von Möglichkeiten ${\binom {N+I-1}{N}}$ entspricht.

Die immer wieder erquickliche Umformulierung beweist diese Kombinatorik:

Die N Kugeln sollen auf I Fächer $i=1...I$ verteilt werden. Die Fächer kriegen (I-1) Trennwände, so dass jede Verteilung genau einer Folge aus (N+I-1) Dingen entspricht, entweder Kugel oder Wand. Gesucht ist die Zahl aller Möglichkeiten, die (I-1) Wände zu setzen. Das ergibt ${\binom {N+I-1}{I-1}}={\binom {N+I-1}{N}}$ .

Der Boson-Vektorraum ist der Unterraum der total symmetrischen Tensorprodukte in $S(N)\subset T(N)$ , er ist die Bildmenge des Symmetrierung-Operators S. Operator $S:T(N)\to S(N)$ ist definiert wie folgt (linear zu erweitern):

S(B(i_{1})\otimes B(i_{2})\otimes ...\otimes B(i_{N}))=\sum _{p\in P_{N}}{\frac {{\text{one}}(p)}{N!}}\;[B(i_{p1})\otimes B(i_{p2})\otimes ...\otimes B(i_{pN})]

Auch der Operator S ist eine Projektion. Eine norm-erhaltende Version für Produkte orthogonaler Funktionen wäre: ${\bar {S}}=({\sqrt {N!}}\,S)$ .

Der Grundzustand mit vielen Bosonen ist der, wo alle auf dem selben niedrigsten Energie-Niveau hocken: Bose-Einstein-Kondensation. Er hat mehr Gewicht bei insgesamt ${\binom {N+I-1}{N}}$ im Vergleich zu den $I^{N}$ klassisch markierbaren möglichen Zustands-Dimensionen.

Die Multi-Fermion- und Multi-Boson-Zustände liegen in Eigenräumen zum Eigenwert 1 der Operatoren A bzw. S. Der Hamilton-Operator und alle Operatoren von Observablen sollten diesen Zustandsraum nicht verlassen können und sollten mit allen Permutations-Operatoren von Teilchenpaaren (Transpositionen) kommutieren. Es ist auch nicht möglich, dass ein System aus gleichen Teilchen sein Verhalten dynamisch ändert, von Symmetrie zu Antisymmetrie oder zurück.

(In der spekulativen Physik gibt es Super-Symmetrien mit Operatoren, die Bosonen und Fermionen ineinander umwandeln! Bisher wurde nichts dergleichen in der Natur gefunden.)

Auf der direkten Summe aller Boson-Vektorräume zu allen Teilchenzahlen

\Sigma _{S}=\mathbb {C} \oplus S(1)\oplus S(2)\oplus S(3)\oplus \cdots

gibt es ein Produkt. Es wird auf Monomen von Produktzuständen von N bzw. M Teilchen als $U_{N}\bullet V_{M}=S_{N+M}(U_{N}\otimes V_{M})$ definiert, bilinear fortzusetzen auf alle Polynome. Das Produkt ist assoziativ und kommutativ; der Boson-Raum ist damit eine Kommutative Algebra. Die erste Dimension von komplexen Zahlen wird hinzugenommen und ist physikalisch der Zustand des Vakuums, das null Teilchen enthält. Damit hat die Algebra ein Einselement.

Auf der direkten Summe der Fermion-Vektorräume zu allen Teilchenzahlen

\Sigma _{A}=\mathbb {C} \oplus A(1)\oplus A(2)\oplus A(3)\oplus \cdots

ist das entsprechende Produkt auf Monomen von N bzw. M Teilchen als $U_{N}\wedge V_{M}=A_{N+M}(U_{N}\otimes V_{M})$ definiert, bilinear fortzusetzen auf alle Polynome. Das Produkt ist assoziativ und alternierend. Der Fermion-Raum ist damit eine Grassmann-Algebra. Die Vertauschung funktioniert mit der Vorzeichenregel:

U_{N}\wedge V_{M}=(-1)^{NM}V_{M}\wedge U_{N}

.

Zwei Monome mit ungeraden Einteilchen-Faktoren antikommutieren, alle anderen Paare kommutieren.

Die beiden Algebren haben einen Grad bei Monomen, die Zahl der elementaren Faktoren (wie bei gewöhnlichen Polynomen der Homogenitäts-Grad). Sie werden graduierte Algebren genannt. Der Grad ist additiv bei Multiplikation von homogenen Elementen. Viele Operatoren auf den Algebren variieren den Grad um 0,+1,-1 Einheit. Manche agieren auf Produkten mit einer Leibniz-Regel, ähnlich zur Ableitung vom Produkt zweier Funktionen; sie heißen dann Derivationen bzw. Anti-Derivationen.

Exponential-Verteilung

In den einfachsten klassischen Modellen der statistischen Mechanik sind die Energie-Zustände der Teilchen exponentiell verteilt. Die relative Häufigkeit der Teilchen mit Energie $E_{i}\;{\text{ist}}\;W(i)\propto \exp(-\beta E_{i})$ .

Quantenmechanisch gibt es weit weniger 'Freiheitsgrade', wie schon aus der Abzählung der Basen hervorging. Pauli verbietet den Fermionen die Doppelbelegung eines Niveaus, und auch bei Bosonen zählen gleiche Partikelbelegungen der Niveaus nur als ein einziger Zustand.

Ein Multi-Boson-Zustand ist eine Folge ${\bar {n}}=\{n(i)\}$ ganzer Zahlen; seine Häufigkeit ist ${\bar {W}}({\bar {n}})\propto \exp(-\beta \sum n(i)E_{i})$ . Interessant ist die relative Häufigkeit von Niveau $i,W(i)$ , wofür der Erwartungswert von $n(i)$ aus der Verteilung ${\bar {W}}({\bar {n}})$ zu berechnen ist.

W(i)\propto (\sum _{\bar {n}}n(i){\bar {W}}({\bar {n}}))/(\sum _{\bar {n}}{\bar {W}}({\bar {n}}))

Faktoren $exp(-\beta \sum _{j;j\neq i}n(j)E(j))$ können aus allen ${\bar {W}}$ gestrichen werden. Es bleibt mit $\alpha ^{n(i)}:=\exp(-\beta n(i)E_{i})$ :

W(i)=(\sum _{n}n\alpha ^{n})/(\sum _{n}\alpha ^{n})

.

f(\alpha )=\sum _{n}\alpha ^{n}=1/(1-\alpha );\quad \alpha f'(\alpha )=\sum _{n}n\alpha ^{n}=\alpha /(1-\alpha )^{2}\quad \Rightarrow

W(i)=\alpha /(1-\alpha )=1/[\exp(\beta E_{i})-1]

Ein Multi-Fermion-Zustand ist eine Folge $\{n(i)\}$ die nur aus den Werten 0 und 1 besteht. Die Wahrscheinlichkeit für ein Teilchen auf dem Niveau i wird dann ganz analog mit Summen über {0,1} statt $\{0...\infty \}$ :

W(i)=\alpha /(1+\alpha )=1/[\exp(\beta E_{i})+1]

??? Normierung? Diese Verteilung wird nahezu Null für E_i weit oberhalb der thermischen Energie k_BT=1/\beta, nahezu konstant für Energien weit darunter? negative?

Zusammenfassung

Klassische Verteilung $\quad W(i)\propto \exp(-\beta E_{i})$ .
Fermi-Dirac-Verteilung $\quad W(i)\propto 1/[\exp(\beta E_{i})+1]$
Bose-Einstein-Verteilung $\quad W(i)\propto 1/[\exp(\beta E_{i})-1]$

Periodensystem

Radiale Dichte der Elektronenladung von Edelgas-Atomen. Zeigt den Aufbau des Elektronensystems in Schalen

Nach dem Pauli-Prinzip bauen sich die Atome um einen Kern der Ladung Z auf. In erster Näherung belegt das erste Elektron das niedrigste Energieniveau von den Lösungen im Coulomb-Potenzial. Die weiteren Elektronen belegen höhere Plätze. Sie sehen ein verformtes Potenzial wegen der Abschirmung des Kerns durch alle anderen Elektronenwolken.

Die einfachste Näherung wäre, für alle Z Elektronen das selbe mittlere radiale Potenzial $V(r)$ einzusetzen. Der vereinfachte Hamilton-Operator wäre

H=\sum _{i=1}^{Z}\left[({\vec {p}}_{i}^{2}/2m)+V({\vec {x}}_{i})\right]

.

Die Wellenfunktionen kämen aus einer Basis der Slater-Determinanten von wasserstoff-ähnlichen Lösungen.

Mit besseren Methoden und Computer-Hilfe konnte die Wellenmechanik nach Schrödinger und Pauli die Grundzustände aller Atome des Periodensystems naturgetreu berechnen. Die Entartung aller (n,l) zu gleichem n wird durch das kollektive Potenzial aufgehoben. Die Reihenfolge, in der die Paare (n,l) das höchste belegte Niveau stellen, konnte genau bestätigt werden.

Die folgende Tabelle zeigt die Ordnung der Niveaus mit Hauptquantenzahl 1-6 und Drehimpuls 0-3 (spdf), wie sie mit steigender Energie gefüllt werden. Bei gegebener Hauptquantenzahl haben die Drehimpulse s,p,d,f, zusammen mit den Paaren vom Elektronenspin, je 2,6,10,14 Zustände mit fast gleicher Energie. Jede Reihe beginnt mit einem Alkali-Atom und endet mit einem Edelgas. Vor dem Edelgas steht jeweils ein Halogen. Das Edelgas mit seiner stabilen Elektronenschale ist chemisch inert. Das Alkalimetall gibt reaktionsfreudig sein vereinzeltes Elektron ab. Genauso begierig nimmt das Halogen ein Elektron auf, um die Edelgaskonfiguration zu erreichen. Die theoretische und numerische Chemie ist nichts anderes als angewandte Quantenmechanik.

Niveaus	Anzahl	Elemente
1s	2	H He
2s 2p	2+6	Li ... F Ne
3s 3p	2+6	Na ... Cl Ar
4s 3d 4p	2+10+6	K ... Br Kr
5s 4d 5p	2+10+6	Rb ... I Xe
6s 4f 5d 6p	2+14+10+6	Cs ... At Rn

Die Konfigurationen der Elektronen der Alkalimetalle Natrium und Kalium und des Halogen-Atoms Chlor sehen schematisch aus wie folgt.

Thomas-Fermi-Modell

Dies ist ein vereinfachtes Modell für das kollektive Zentralpotenzial V(r), das von den Elektronen um einen Kern der Ladung Z gefühlt wird. Gefordert:

Asymptotisch für $r\to 0:rV(r)\to Ze$
Für $r>0:\Delta V=(1/r)({\frac {d^{2}}{dr^{2}}}r)V(r)=4\pi e\rho (r)$
$\rho (r)=\sum _{i=1}^{Z}|\psi _{i}|^{2}$ (Elektron-Dichte im Grundzustand)
$Z=4\pi \int _{0}^{\infty }\rho (r)r^{2}\;dr$

V ist am Zentrum von der Kernladung bestimmt und anderswo von der Ladungsdichte der gesamten Z Elektron-Wellen-Faktoren. Diese stecken in einer Slater-Determinante, die den Grundzustand annähert.

Nun kommt eine halb-klassische Idee hinzu: die Dichte von gebundenen Quantenzuständen im Phasenraum. In der WKB-Näherung sind die Energie-Niveaus des eindimensionalen Modells gequantelt nach der Regel:

\int k(x)dx=(N+{\tfrac {1}{2}})\pi ;\quad (\hbar k(x))^{2}/(2m)=E-V(x)

.

Die Ebene der Punkte $(q,p)=(x,\hbar k)$ heißt der Phasenraum, seine allgemeine Version für N Freiheitsgrade hat die Dimension 2N.

Mit $p=\hbar k;2\pi \hbar =h$ kann die Quantelung so beschrieben werden: Das klassische System im Potenzialtopf oszilliert zwischen zwei Umkehrpunkten $x=a,\;x=b$ und beschreibt geschlossene Kurven im Phasenraum $(p,q)=(\hbar k,x)$ . Zuerst mit $+k(x)$ für x von a bis b, dann zurück mit $-k(x)$ für x von b nach a. Das Integral über eine Runde ist die umrundete Fläche in der $(p,q)$ -Ebene.

I=\oint p(q)dq=2\hbar \int k\;dx=h(N+{\tfrac {1}{2}})

.

Einer Stufe der Energie E entspricht ein Zuwachs von "Volumen im Phasenraum" $\oint pdq$ vom Wert h. In drei Dimensionen wächst der Phasenraum um den Wert $h^{3}$ zwischen Energieniveaus. Die Quantenmechanik schreibt also den gebundenen Fermionen wegen des Pauli-Verbots vor, mit einer konstanten Packungsdichte $(2/h^{3})$ den Phasenraum aufzufüllen, angefangen mit dem niedrigsten Energie-Niveau. Der Faktor 2 zählt die 2 Dimensionen des Spins.

Nach dem Thomas-Fermi-Modell sind die Energiezustände im 6-dimensionalen Phasenraum $({\vec {x}},{\vec {p}})$ lokalisiert mit der Dichte $D({\vec {x}},{\vec {p}})=(2/h^{3})$ falls die Energie $E({\vec {x}},{\vec {p}})=p^{2}/(2m)-eV({\vec {x}})$ kleiner ist als das maximal besetzte Niveau $E_{0}$ . Man darf $E_{0}=0$ festlegen, denn eine Verschiebung ändert nicht das Verhalten von $(rV(r))$ am Nullpunkt. Der Faktor 2 berücksichtigt die zwei Spin-Richtungen der Elektronen. Die Elektronendichte im Raum:

\rho ({\vec {x}})=\int _{\{E<0\}}d{\vec {p}}\;D({\vec {x}},{\vec {p}})=

(2/h^{3})\;(4\pi )\int _{\{p^{2}<(2meV)\}}p^{2}\;dp={\frac {8\pi }{3h^{3}}}\,(2meV)^{3/2}

.

Die Dichte $\rho$ hängt also von $V(r)$ ab und ist positiv wenn $V>0$ , Null sonst. Es gilt also diese Differenzialgleichung für $V(r)$ :

(1/r)(d^{2}/dr^{2})(rV(r))=4\pi e\rho (V(r))

.

Es werden $(V,r)$ durch dimensionslose Variablen $(u,x)$ ersetzt:

r=(b/{\sqrt {Z}})x\;{\text{mit}}\;b={\frac {1}{2}}\left({\frac {3\pi }{4}}\right)^{2/3}{\frac {\hbar ^{2}}{me^{2}}}

V=(Ze/r)u\;

mit dem Grenzwert

u(x=0)=0

wegen

rV(r)\to 0\;{\text{für}}\;r\to 0

.

Damit: $\rho ={\frac {Z^{2}}{4\pi h^{3}}}(u/x)^{3/2}$ im Gebiet $\{r:\;u(r)>0\}$ .
Die Differenzialgleichung für u lautet in der Kugel wo $\rho$ positiv ist:

u''=(d^{2}u/dx^{2})={\sqrt {u^{3}/x}};\;u(0)=1;\ u\geq 0\;{\text{für}}\;x\leq x_{0}

.

Wann kommt die erste Nullstelle $x_{0}$ von u(x)? Ausgerechnet wird dazu eine Normierung von u, die daraus folgt dass $\rho$ sich zu Z summiert:

1=\int _{0}^{x_{0}}{\sqrt {xu^{3}}}\;dx=\int _{0}^{x_{0}}xu''\;dx=(xu'-u)|_{0}^{x_{0}}=x_{0}u'(x_{0})+1

.

Daraus folgt $u'(x_{0})=0$ an der Nullstelle. Das weist darauf hin, dass $x_{0}$ bei Unendlich liegt. $u''={\sqrt {u^{3}/x}}$ ist also zu lösen mit den zwei Randbedingungen $u(0)=1;\;u(\infty )=0$ . Ein Skript im Anhang integriert die Gleichung numerisch und macht den Graphen als SVG-Datei.

Aus der Funktion u() folgen die radialen Potenziale $V(r)$ und Ladungsdichten $\rho (r)$ , die unter der Annahme $Z\gg 1$ großer Atome angenähert beschreiben, wie dicht die Elektronen sich drängeln dürfen. Man definiert einen Atomradius $R(\alpha )$ als die Kugel, in der der Bruchteil $(1-\alpha )$ der Elektronenladung ist. Typischerweise $\alpha =(1/Z)$ , also die Kugel für alle bis auf ein Elektron. Gleichung: $(1-\alpha )Z=Z-1=4\pi \int _{0}^{R}\rho r^{2}\;dr$ . Mit dimensionslosem $X(\alpha );\;R(\alpha )=Z^{-1/2}b\,X(\alpha )$ , ergibt sich X als numerische Lösung von $u(X)-Xu'(X)=\alpha$ . Ergebnis: R variiert recht langsam mit steigendem Z. Der Phasenraum wird schneller in der Impuls- als in der Orts-Dimension aufgeblasen.

Das Skript integriert numerisch die Gleichung $u''(x)={\sqrt {u^{3}/x}}$ mit der Probeschuss-Technik wie folgt, ohne Anspruch auf hohe Präzision:

Intervall $x\in [0,1]$ wird in 100 gleiche Schritte d=0.01 geteilt.

Der typische Rechenschritt bestimmt $u[i+1]$ aus $(u[i+1]+u[i-1]-2u[i])/d^{2}={\sqrt {u[i]^{3}/x[i]}}$

Intervall $x\in [1,200]$ wird als geometrische Folge in 1% große Schritte geteilt.

Hier ist u zur Funktion $v(y);\;y=\ln(x)$ umgeschrieben und die Rechenschritte sind die Differenz-Version von $(v''-v')=(v\exp(y))^{3/2}$ .
Es wird $u[0]=1$ festgesetzt und $u[1]\approx 0.98$ mit Dichotomie variiert, bis am Ende $u(x=200)\sim 0$ herauskommt, von Rundungsfehlern abgesehen. Die SVG-Grafik zeigt die Funktion mit logarithmischer x-Achse.

Anhang: Skript

Skript fürs Thomas-Fermi-Modell. Schreibt eine Datei thomferm.svg.

'''
Thomas-Fermi potential model.
solve x(f")^2 = f^3 on interval [0,inf), f(0)=1, f(1)=0
implemented direct shoot algorithm with log section:
f(x)=F(y); y=log(x); x=exp(y); f'=F'/x; f"= F"/x^2-F'/x^2=(F"-F')/(xx)
F"-F'= xx sqrt(F^3/x)= (F*x)^(3/2) = (F*exp(y))^(3/2)= z
(F(2)+F(0)-2F(1))/(dd)-(F(2)-F(0)/(2d) = z
(1-d/2)F(2)/dd = z + (2F(1)-F(0))/dd - F(0)/(2d) 
(1-d/2)F(2)= ddz+ 2F(1)- (1+d/2) F(0)
'''

from math import sqrt,exp,log
def div(p,q) : return int(p/q)

###### plot

class svgdump :
  def plotbegin(sd, w,h,f) : # args: width height zoomfactor
    W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
    s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
     '<!DOCTYPE svg>\n'+
     '<svg xmlns="http://www.w3.org/2000/svg"\n'+
     '  xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
     '  width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
    return s

  def plotend(sd) : return '</svg>\n'
  def move(sd,x,y) : return 'M'+str(x)+','+str(y)
  def line(sd,x,y) : return 'L'+str(x)+','+str(y)
      
  def labl(sd,s,x,y, size=14,fill='#000') :
    f=' font-family="Arial"'
    t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
    return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')

  def path(sd,w,s,lc='#000') : # path w=with lc=linecolor
    t='<path fill="none" stroke="'+lc+'" stroke-width="'+str(w)+'" '
    t+= ' stroke-linecap="round" ' 
    return t + 'd="\n'+s+'" />\n'

# end class svgdump

def aprx(x) : # 3-digit approximation
  y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
  return z/1000.0

def curve(sd,scale,xy, fill='#000', width=1) :
  (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2) 
  for i in range(n) : 
    x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
    s+= sd.move(x,y) if(i==0) else sd.line(x,y)
  return sd.path(width,s,lc=fill)

def label(sd,scale,txt,tx,ty) : # tx,ty plot coordinates not pixels
  (xs,fx,xmin, ys,fy,ymin) = scale 
  x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
  return sd.labl(txt,x,y)

###########

def shootlog(y,f,i,k) : # x=exp(y) version
  d=y[i]-y[i-1]; z=f[i]; j=i; p=1.0-0.5*d; q=1.0+0.5*d
  while (j<k) and (z>=0) and (z<10) : # stop if out of bounds
    t=exp(y[j])*z;  z=sqrt(t*t*t)
    u= d*d* z + 2*f[j]-q*f[j-1]
    f[j+1]= u/p; y[j+1]=y[j]+d; j+=1; z=f[j]
  return j,z  # print(str(f[j+1]))

def shootlin(x,f,i,k) : # i=1, f[0] and f[1] known, x[0]=0, x[1]=d.
  d=x[i]-x[i-1]; j=i; z=f[j]
  while (j<k) and (z>=0) and (z<10) :
    u= d*d* sqrt(z*z*z/x[j])
    f[j+1]= u+2*z-f[j-1]; x[j+1]=x[j]+d; j+=1; z=f[j]
  return j,z    # print(str(f[j+1]))

def logplot(fn,x,f,y,g) : # log x axis, y axis 0 to 1.
  txt=['Function y\"= sqrt(y^3/x), y(0)=1.',
    'xscale log 0.01...100','yscale 0...1']
  n=len(x); m=len(y); ok=(g[1]==f[n-1]); k=0; xy=[0.0]*(2*n+2*m)
  if not ok : print('logplot no join'); exit()
  for i in range(1,n) : xy[k]=log(x[i]);xy[k+1]=f[i]; k+=2
  for i in range(2,m) : xy[k]=y[i]; xy[k+1]=g[i]; k+=2
  xmin=log(0.01); xmax=log(100); ymin=0; ymax=1.0; xy=xy[:k]
  xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin ] # frame
  xg=[1,2,4,6,8]; fg=[0.01,0.1,1,10,100]; z=0; 
  lbl=['0.01','0.1','1','10','100']
  for i in range(len(fg)) :
    for k in range(len(xg)) :
      f=fg[i]; g=f*xg[k]; h=log(g)
      if (i==0)and(k==0) : xp=xmin; z=0
      elif (h<=xmax)and(z==0) : xyf += [xp,ymin,h,ymin,h,ymax];xp=h; z=1
      elif (h<=xmax)and(z==1) : xyf += [xp,ymax,h,ymax,h,ymin];xp=h; z=0
  xys=[xyf,xy]
  xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale user units, min in pixels
  ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#777','#000','#00a','#077','#0a0','#770','#a00']
  sd=svgdump(); buf=sd.plotbegin(800,600,1); i=0; dy=(ymax-ymin)/20.0
  for j in range(len(txt)) :
    buf += label(sd,scale, txt[j], (xmin+xmax)*0.5,ymax-dy-j*dy)
  for j in range(len(lbl)) :
    buf += label(sd,scale, lbl[j], log(fg[j]),ymax)
  for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

def shooter() : # dichotomic shots to solve equation
  # adjust f[1] to hit target
  # logscale y,g(y) after initial x,f(x) linear section 
  n=50; n1=n+1; x=[0]*n1; f=[0]*n1; x[0]=0;x[1]=0.01; f[0]=1; fhi=1; flo=0.5
  m=300; m1=m+1; y=[0]*m1; g=[0]*m1; logsection=True
  targ=1e-8; ff=0.5*(fhi+flo); f[1]=ff; zhi=10; zlo=-1; 
  i=0; busy=True;zbefore=1
  while busy and (i<100) :
    j,z= shootlin(x,f,1,n); print(str([f[1],j,z]))
    if logsection and (z>0)and(z<1) :
      g[0]=f[j-1]; g[1]=f[j]; y[0]=log(x[j-1]); y[1]=log(x[j])
      k,z= shootlog(y,g,1,m); xmax=exp(y[k]); print('log: '+str([xmax,k,z]))
    if z>targ : fhi=ff;zhi=z 
    else : flo=ff;zlo=z
    ff=0.5*(fhi+flo);
    if (zlo>0)and(zlo<zhi)and(zhi<1) : ff= flo+(targ-zlo)*(fhi-flo)/(zhi-zlo)
    f[1]=ff; i+=1; busy= (abs(z-zbefore)>1e-8); zbefore=z
  print('Done at i='+str(i))
  logplot('thomferm', x,f,y,g)
  atomradius('bla', x,f,y,g)

def atomradius(fn,x,f,y,g) : # solve x(a): u(x)-xu'(x)=a where a=0.01..0.1
  # u-xu'= u[i]-xx[i]*(u[i+1]-u[i-1])/(xx[i+1]-xx[i-1])
  n=len(x); m=len(y); k=0; xx=[0.0]*(n+m);u=[0.0]*(n+m); v=[0.0]*(n+m)
  for i in range(1,n) : xx[k]=x[i];u[k+1]=f[i]; k+=1
  for i in range(2,m) : xx[k]=exp(y[i]); u[k+1]=g[i]; k+=1
  vmax,jmax,vmin,jmin= -1e6,0,1e6,0
  for j in range(1,k-1) :
    v[j]= u[j]-xx[j]*(u[j+1]-u[j-1])/(xx[j+1]-xx[j-1])
    if v[j]>vmax: vmax,jmax= v[j],j
    if v[j]<vmin: vmin,jmin= v[j],j
  if False : #debug
    print('v='+str([v[1],v[50],v[100],v[k-2]]))
    print('vmax,xmax='+str([vmax,xx[jmax]]))
    print('vmin,xmin='+str([vmin,xx[jmin]]))
  print('Percent electron charge; Radius of ball')
  for p in range(1,21) : # p=percentage
    a=0.01*p ; j=jmax
    while v[j]>a : j+=1
    xa=xx[j-1];xb=xx[j]; va=v[j-1];vb=v[j]
    xp=xa+(a-va)*(xb-xa)/(vb-va); print('p,xp='+str([100-p,xp]))  

### following stuff not working

def geterror(f,y,err, n) :
  d=1.0/n
  for i in range(1,n) :
    x=y[i]; z=1.0-x; a=f[i]; f3=a*a*a
    f1=(f[i+1]-f[i-1])/(2.0*d); f2=(f[i+1]+f[i-1]-2*f[i])/(d*d)
    b=f2-2*z*f1; c=x*z*z*b*b; err[i]= f3-c

def shootstep(y,f,i,n) : # get f[i+1] for zero error.
  d=1.0/n; x=y[i]; z=1.0-x; a=f[i]; f3=a*a*a;
  # init bracket problem minimum find
  def err(i,ff) :
    f1=(ff-f[i-1])/(2.0*d); f2=(ff+f[i-1]-2*f[i])/(d*d)
    b=f2-2*z*f1; c=x*z*z*b*b; return f3-c
  # shoot is initialized with f[0] and f[1], the first slope
  def bugsweep(i,fh,eh) : # nearly degen pair of zeroes!
    df=1e-2; fmin=fh;emin=1e9
    for k in range(25) :
      ff=fh+(k-10)*df; ef=err(i,ff);
      fmin,emin = (ff,ef) if (abs(ef)<emin) else (fmin,emin)
      print('bug k,ff,ef='+str([k-10,ff,ef]))
    df=1e-3; fs=fmin; ok=False; fl=fmin;el=emin
    for k in range(25) :
      ff=fs+(k-10)*df; ef=err(i,ff)
      print('bug2 k,ff,ef='+str([k-10,ff,ef]))
      if (not ok) and ((ef*eh)<0) : ok=True; fl=ff;el=ef
    return fl,el
  #
  def descent(i,fh) :
    # search for bracket for f[i+1], fh<0 search a point >0
    eps=0.001; df= max(abs(eps*fh),eps); fl=fh-df; j=0; nf=0
    eh= err(i,fh); el=err(i,fl); found=(el*eh)<0; q=el/eh
    if (not found) and (q>1) : df=-df; fl=fh-df; el=err(i,fl)
    #print('i,eh,el,fh,fl='+str([i, eh,el, fh,fl]))
    while (j<100)and(not found) :
      q=el/eh
      #if q>1 : df=-0.5*df; fl -= df; nf+=1; print('flip j,eh,el='+str([j,eh,el]))
      if q>0.9 : df= 4*df; fl -=df
      elif q>0.5 : df= 2*df; fl -=df
      else : fl -= df
      el=err(i,fl); found=(el*eh)<0; j= 100 if (nf>4) else (j+1)
    if not found : fx,ex= bugsweep(i,fh,eh); 
    # print('i,j,eh,el,fh,fl='+str([i,j, eh,el, fh,fl]))
    return fl,el,j
  #
  fh=a; eh= err(i,fh); ok=True; k=0
  fl,el,j =descent(i,fh)
  while (k<20)and(ok) : # this part working, slow convergence
    if (eh*el)>=0 : ok=False;print('shoot exit i,k='+str(i)+','+str(k))
    ff= (fh*el-fl*eh) / (el-eh); ef=err(i,ff)
    if ef*eh>0 : fh=ff; eh=ef
    else : fl=ff; el=ef
    k +=1
  if ok : f[i+1]=ff; print('i,j,k, ff,ef='+str([i,j,k,ff,ef]))
  return ok

def smooth(x) :
  n=len(x); n1=n-1; y=[0]*n
  for i in range(1,n1) : y[i]=(2.0*x[i]+x[i-1]+x[i+1])/4.0
  return y

# def downerr(y,f,g,y,ef,eg, n) :
# tweak f[i],g[i] so that err at [i-1] doesnt worsen 

def explore(n) :
  d=1.0/n;n1=n+1; f=[0]*(n+1); g=[0]*n1; y=[0]*n1; f[0]=1; f[n]=0; 
  for i in range(n+1) : f[i]= 1.0-d*i; y[i]=d*i; g[i]=f[i] - 0.7*y[i]*f[i]
  print('start end', str([f[0],g[0],f[n],g[n]]))
  errf=[0]*n1; errg=[0]*n1; h=[0]*n1; errh=[0]*n1; h[0]=1
  geterror(f,y,errf, n); ef=smooth(errf); ef=errf
  geterror(g,y,errg, n); eg=smooth(errg); eg=errg
  for i in range(1,n) : h[i]=(f[i]*eg[i]-g[i]*ef[i]) / (eg[i]-ef[i])
  for i in range(1,n) : h[i]=(f[i]+g[i])/2.0
  geterror(h,y,errh,n) # secant fails, with eg,ef worse than average
  #for i in range(1,n) :
  #  print(str([i,errf[i],errg[i],errh[i]]))
  for i in range(1,n) : ok= shootstep(y,g,i,n)

# explore(20)
shooter()

Hoch zum Seitenanfang