Quantenmechanik/ Naeherungen

Zurück zu Drehimpuls-Kopplung | Hoch zum Quantenmechanik Anker | Vor zu Numerik-Beispiel

Näherungsverfahren

Meist ist die Schrödinger-Gleichung nicht in geschlossener Form lösbar, wenn Energien von gebundenen Zuständen bei Molekülen usw. gesucht sind. Dann kommen Approximationsverfahren zum Einsatz, die mehr oder weniger mit der Brechstange sich von bekannten Lösungen mit einfachen Hamilton-Operatoren vorarbeiten zu den ungefähren Lösungen mit komplizierteren Operatoren. Kurze Einführung zur Störungsrechnung und Variationsrechnung. (Nur viel weitergehende Abhandlungen könnten sich mit dem Namen „Störungstheorie“ schmücken.)

Algorithmus Störungsrechnung

Angefangen sei mit einem Hamilton-Operator $H_{0}$ und dessem orthonormalen vollständigen System von Eigenfunktionen $F_{0i};\;H_{0}F_{i}=E_{0i}F_{0i}$. Dann wird mit dem Faktor t eine Störung dazu geschaltet $H(t)=H_{0}+tH_{1}$. Polynome in t dienen hier zum Sortieren nach der Ordnung. Die Eigenvektoren und Eigenwerte variieren bis zur zweiten Ordnung.

$(H_{0}+tH_{1})(F_{0i}+tF_{1i}+t^{2}F_{2i})=(E_{0i}+tE_{1i}+t^{2}E_{2i})(F_{0i}+tF_{1i}+t^{2}F_{2i})$

Die Terme von Ordnung t liefern: $H_{1}F_{0i}+H_{0}F_{1i}=E_{1i}F_{0i}+E_{0i}F_{1i}$.

Das Skalarprodukt mit $\langle F_{0i}|$ macht ein $E_{0i}$ mit dem hermiteschen $H_{0}$. Es bleibt übrig $\langle F_{0i}|H_{1}F_{0i}\rangle =E_{1i}$ weil die $F_{0i}$ normiert sind. Ergebnis: In erster Ordnung ist der korrigierte Energie-Eigenwert gleich dem alten plus dem Erwartungswert des Stör-Operators $H_{1}$ mit dem alten Eigenvektor.

Das Skalarprodukt mit einem anderen ungestörten $\langle F_{0j}|$ bringt:

$\langle F_{0j}|H_{1}F_{0i}\rangle =(E_{0i}-E_{0j})\langle F_{0j}|F_{1i}\rangle ;\quad i\neq j$

Bei Entartung $E_{0i}=E_{0j}$ ist das grundfalsch! Denn die Matrixelemente von $H_{1}$ verschwinden nicht. Um doch die Konsistenz der Näherung zu sichern, macht man erst die hermitesche Matrix $\langle F_{0j}|H_{1}F_{0i}\rangle$ aus all den $F_{0}$, die den entarteten Teilraum zum selben Energie-Eigenwert $E_{0}$ aufspannen. Man diagonalisiert diese Matrix, was eine neue orthonormale Teilbasis mit $\langle G_{0j}|H_{1}G_{0i}\rangle =\delta _{ij}h_{i}$ hervorbringt. Mit den $G_{0i}$ wird jede gefährliche Gleichung weggeschaltet und man rechnet mit diesen weiter.

Ohne Entartung gilt: $<F_{0j},F_{1i}>=\langle F_{0j}|H_{1}F_{1i}\rangle /(E_{i}-E_{j})$. Nun soll die Bedingung ausgewertet werden, dass $F_{0i}+tF_{1i}$ normiert sei:

$1=\langle F_{0i}+tF_{1i}|F_{0i}+tF_{1i}\rangle =1+t\langle F_{0i}|F_{1i}\rangle +t\langle F_{1i}|F_{0i}\rangle +O(t^{2})$

woraus folgt: $0=Re\langle F_{0i}|F_{1i}\rangle$. Jetzt nutzt man aus, dass die komplexe Phase von $F_{0i}$ beliebig wählbar ist, ohne die Orthonormalität und die erste Energiekorrektur zu verändern. Die Phase sei so gedreht, dass auch $Im\langle F_{0i}|F_{1i}\rangle =0$ gilt. Mit $\langle F_{0i}|F_{1i}\rangle =0$ ist dann erreicht, dass solche Paare orthogonal sind. Mit dem Orthonormalsystem zu $H_{0}$ folgt dann die Näherungsformel für $F_{1i}$:

$F_{1i}=\sum _{j}\langle F_{0j}|F_{1i}\rangle F_{0j}=\sum _{j,j\neq i}F_{0j}<F_{0j}|H_{1}F_{0i}>/(E_{i}-E_{j})$

Dies ist die Entwicklung der Vektor-Störung in der ungestörten Basis.

Die Korrekturen der Ordnung $t^{2}$, genannt $F_{2i}$ und $E_{2i}$, sind nun fällig.

$H_{0}F_{2i}+H_{1}F_{1i}=E_{0i}F_{2i}+E_{1i}F_{1i}+E_{2i}F_{0i}$

Hierin nach erster Ordnung:

$E_{1i}=\langle F_{0i}|H_{1}F_{0i}\rangle ;\;F_{1i}=\sum _{j\neq i}F_{0}j\langle F_{0j}|H_{1}F_{0i}\rangle /(E_{i}-E_{j})$

Das Skalarprodukt der zweiten Ordnung mit $\langle F_{0i}|$ macht $\langle F_{0i}|H_{0}F_{2i}\rangle$, was wegen Hermitizität mit $E_{0i}\langle F_{0i}|F_{2i}\rangle$ herausfällt. Daher folgt:

$\langle F_{0i}|H_{1}F_{1i}\rangle =E_{2i}+E_{1i}\langle F_{0i}|F_{1i}\rangle$

Der letzte Term wurde oben mittels Phasendrehung zu Null wegdiskutiert. Resultat nach Einsetzen der Entwicklung von $F_{1i}$:

$E_{2i}=\sum _{j\neq i}|\langle F_{0j}|H_{1}F_{0i}\rangle |^{2}/(E_{i}-E_{j})$

Während die erste Energie-Korrektur nur den Eigenzustand selbst berührt, bezieht die zweite bevorzugt alle Nachbarzustände mit ein. Wenn $F_{0i}$ der Grundzustand war, $E_{j}>E_{i}\;\forall j$, ist diese zweite Korrektur negativ.

In einem Zwei-Niveau-System sind die zweiten Korrekturen gegenläufig:

$\delta E_{1}=\langle F_{1}H_{1}F_{1}\rangle +|\langle F_{2}H_{1}F_{1}\rangle |^{2}/(E_{1}-E_{2})$

$\delta E_{2}=\langle F_{2}H_{1}F_{2}\rangle -|\langle F_{2}H_{1}F_{1}\rangle |^{2}/(E_{1}-E_{2})$

Die unendlichen Summen der zweiten Korrektur geben nicht selten Konvergenz-Probleme auf, so etwas wie Infrarot- oder Ultraviolett-Katastrophen.

Algorithmus Variationsrechnung

Diese Berechnung ist vor allem gebrächlich für Grundzustände von Molekülen. Die kleinste Energie ist das absolute Minimum aller Erwartungswerte des Hamiltonoperators. Beweis. Sei ein beliebiger Zustand im orthonormalen Eigensystem des Operators H entwickelt,

$F=\sum _{n}F_{n}\langle F_{n}|F\rangle ;\;\langle F_{n}|F_{m}\rangle =\delta _{nm};\;HF_{n}=E_{n}F_{n}$

dann ist sein Erwartungswert

$\langle FHF\rangle /\langle F|F\rangle =[\sum _{n}E_{n}|\langle F_{n}|F\rangle |^{2}]/[\sum _{n}|\langle F_{n}|F\rangle |^{2}]$

größer oder gleich $E_{0}$ wo $E_{n}\geq E_{0}\;\forall n$.

Angemerkt sei noch, dass die Eigenwerte von H Extrema oder stationäre Punkte von Erwartungswerten $E=\langle FHF\rangle /\langle F|F\rangle$ sind, wenn F beliebige Kurven F(z) beschreibt. Man setze dF/dz = f(z). Wenn H F= E F (H hermitesch), verschwindet dieses f:

$E(z)={\frac {\langle F(z)|H|F(z)\rangle }{\langle F(z)|F(z)\rangle }}$

${\frac {d}{dz}}E(z)={\frac {1}{\langle F|F\rangle ^{2}}}(\langle fHF\rangle +\langle FHf\rangle )\langle F|F\rangle -\langle FHF\rangle (\langle f|F\rangle +\langle F|f\rangle )$.

Die Strategie des Variations-Algorithmus ist, eine Familie von Test-Wellenfunktionen F(p) zu erzeugen mit n freien Parametern. Es wird das Minimum der Erwartungswerte von $H,\langle F(p)HF(p)\rangle /\langle F(p)|F(p)\rangle$ im Parameter-Raum gesucht. Damit wird bei richtiger Wahl der Grundzustand gut angenähert.

In der einfachsten Variante sind die Testfunktionen auf einer Hyperebene; sie hängen linear von n Koordinaten ab. Die Berechnung des Minimums führt (mit einem Lagrange-Faktor der Normierung?) auf Eigenwert-Gleichungen. Man sucht die kürzste Achse eines Ellipsoids in n Dimensionen. Es geht auch ohne orthonormierte Vektoren, wie ein Beispiel unten zeigt.

In anderen Varianten wird sowohl ein hermitesches H wie auch Vektor F längs einer Kurve beispielsweise so variiert $H=H(z),\;F=F(z),$ dass F immer Eigenvektor bleibt: $H(z)F(z)=E(z)F(z)$.

Wenn auch noch die Normierungs-Bedingung $\langle F(z)|F(z)\rangle =1$ greift, dann folgt die Hellmann-Feynman-Gleichung

${\frac {d}{dz}}E(z)=\langle {\frac {d}{dz}}F(z)H(z)F(z)\rangle +\langle F(z){\frac {d}{dz}}H(z)F(z)\rangle +\langle F(z)H(z){\frac {d}{dz}}F(z)\rangle$

$=\langle F(z){\frac {d}{dz}}H(z)F(z)\rangle +E(z)[\langle {\frac {d}{dz}}F(z)F(z)\rangle +\langle F(z){\frac {d}{dz}}F(z)\rangle ]$

$=\langle F(z){\frac {d}{dz}}H(z)F(z)\rangle$

denn die Klammer = $[{\frac {d}{dz}}\langle F(z)|F(z)\rangle ]=0$ nach Voraussetzung.

Es können auch andere als Grundzustände berechnet werden, wenn man Quantenzahlen direkt einbaut. Um zum Beispiel Zustände im Zentralpotenzial mit Drehimpuls $\ell$ zu erzwingen, bekommt der Satz Testfunktionen die Winkelanteile $Y_{\ell m}(\theta ,\phi )$ verpasst.

Anharmonischer Oszillator

Der Potenzial-Topf $V(x)\propto x^{2}$ eines Oszillators sei verformt mit einem unsymmetrischen Anteil in der dritten Potenz der Ortskordinate.

$H_{0}={\frac {p^{2}}{2m}}+{\frac {1}{2}}m\omega ^{2}x^{2}$

$H_{1}=s\;\hbar \omega \left({\frac {x}{x_{0}}}\right)^{3}$

Die Störung hat drei Zutaten: den dimensionslosen Parameter s, die Referenz-Energie $\hbar \omega$ und Referenz-Länge $x_{0}$.

Eine klassisches Pendel mit der Dynamik $H_{0}+H_{1}$ beschreibt eine verformte Sinus-Schwingung, hat also Oberwellen im Frequenz-Spektrum. Auch die Grundfrequenz fällt als Funktion der Amplitude ab. Quantenmechanisch speisen wir den $x^{3}$ Anteil in die Störungsrechnung ein und schätzen die Korrekturen an der Folge der Niveaus, die für den ungestörten harmonischen Oszillator den gleichen Abstand haben.

Formeln des Harmonischen Oszillators mit der Basis |n> und Leiter-Operatoren:

Normiert: $\langle n|m\rangle =\delta _{nm};\;a^{\dagger }|n\rangle ={\sqrt {n+1}}|n+1\rangle \;a|n\rangle ={\sqrt {n}}|n-1\rangle$

$N|n\rangle =a^{\dagger }a|n\rangle =n|n\rangle ;\;[a,a^{\dagger }]=1;\;H_{0}=\hbar \omega (N+{\tfrac {1}{2}})$

$a={\sqrt {\frac {m\omega }{2\hbar }}}(x+{\frac {i}{m\omega }});\quad a^{\dagger }={\sqrt {\frac {m\omega }{2\hbar }}}(x-{\frac {i}{m\omega }})$

Operatoren $a,\;a^{\dagger }$ sind dimensionslos. x hat die Referenzlänge $x_{0}$:

$x={\sqrt {\frac {\hbar }{2m\omega }}}(a^{\dagger }+a)=x_{0}\cdot (a^{\dagger }+a)$

Mit $q=s\;\hbar \omega$ und $x/x_{0}=a^{\dagger }+a$:

$H_{1}=q[(a^{\dagger })^{3}+a^{3}+3Na^{\dagger }+3(N+1)a];\;N:=a^{\dagger }a$

$H_{1}$ ändert die Quantenzahl n um {-3,-1,1,3}. Er hat folgende Matrixelemente, die von den 4 Termen in ihrer Reihenfolge abstammen, sowie Energiedifferenzen:

$\langle n+3|H_{1}|n\rangle =q\cdot {\sqrt {(n+1)(n+2)(n+3)}};\quad E_{n}-E_{n+3}=-3\hbar \omega$

$\langle n-3|H_{1}|n\rangle =q\cdot {\sqrt {(n(n-1)(n-2)}};\quad E_{n}-E_{n-3}=3\hbar \omega$

$\langle n+1|H_{1}|n\rangle =q\cdot 3(n+1){\sqrt {n+1}};\quad E_{n}-E_{n+1}=-\hbar \omega$

$\langle n-1|H_{1}|n\rangle =q\cdot 3n{\sqrt {n}};\quad E_{n}-E_{n-1}=\hbar \omega$

Es kommen nur Energiekorrekturen von zweiter Ordnung $s^{2}$ vor. Es folgt nach fleißigem Ausrechnen mit der Energie-Formel der Störungsrechnung:

$E_{n}=\hbar \omega [(n+{\frac {1}{2}})-s^{2}(30n^{2}+30n+11)]$

$E_{n}-E_{n-1}=\hbar \omega [1-s^{2}60n]$

Die Niveaus haben einen ungleichen Abstand, der enger wird mit steigendem n. In der Tat zeigen die Vibrations-Spektren von zweiatomigen polaren Molekülen Absorptions-Linien, die solchen Reihen von Energieniveaus entsprechen. Es gibt dafür Modelle des Potenzials vom Typ Lennard-Jones: $V(r)/V_{0}=(r/r_{0})^{2n}-2(r/r_{0})^{n}$.

Helium-Atom Grundzustand

In der vereinfachten Schrödinger-Gleichung gibt es zwei Elektronen an Orten ${\vec {r}}_{a},{\vec {r}}_{b}$, ein ruhendes Zentrum mit anziehendem Potenzial $-Ze^{2}/r$ und das Abstoßungs-Potenzial der zwei Elektronen $U({\vec {r}}_{a},{\vec {r}}_{b})=+e^{2}/|{\vec {r}}_{a}-{\vec {r}}_{b}|$.

Mit $H({\vec {\nabla }},{\vec {r}})=-{\frac {\hbar ^{2}}{2m}}{\vec {\nabla }}^{2}-2{\frac {e^{2}}{r}}$:

$[H({\vec {\nabla }}_{a},{\vec {r}}_{a})+H({\vec {\nabla }}_{b},{\vec {r}}_{b})+U({\vec {r}}_{a},{\vec {r}}_{b})]\Psi ({\vec {r}}_{a},{\vec {r}}_{b})=E\;\Psi ({\vec {r}}_{a},{\vec {r}}_{b})$

Der gesamte Hamilton-Operator ist invariant unter der Vertauschung von ${\vec {r}}_{a},{\vec {r}}_{b}$. Also will man Lösungen, die auch Eigenfunktionen des hermiteschen Vertausch-Operators sind. Konkret, symmetrische und antisymmetrische Funktionen,

$\Psi _{s}({\vec {r}}_{a},{\vec {r}}_{b})=\Psi _{s}({\vec {r}}_{b},{\vec {r}}_{a});\quad \Psi _{a}({\vec {r}}_{a},{\vec {r}}_{b})=-\Psi _{a}({\vec {r}}_{b},{\vec {r}}_{a})$.

Wie anderswo im Buch steht, benötigen Elektronenpaare antisymmetrische Wellen, aber auch Spin-Freiheitsgrade. Die Forderung ist erfüllt, wenn bei symmetrischem Orts-Anteil die Spin-Faktoren antisymmetrisch sind und umgekehrt. Im Grundzustand kombinieren sich die Spins zum antisymmetrischen Singulett-Zustand mit Spin Null. Daher nimmt dieses Beispiel ein symmetrisches Paar von Ortsfunktionen an.

In nullter Ordnung setzt man U zu Null und nimmt simple Produkte von bekannten Wasserstoff-Atom-Lösungen:

$\Psi ({\vec {r}}_{a},{\vec {r}}_{b})=g({\vec {r}}_{a})\cdot h({\vec {r}}_{b});\;H({\vec {\nabla }},{\vec {r}})h=E_{a}h;\;H({\vec {\nabla }},{\vec {r}})g=E_{b}g;\;E=E_{a}+E_{b}$.

Das niedrigste Niveau 1s,1s hat mit $g({\vec {r}})=h({\vec {r}})={\sqrt {\frac {Z^{3}}{\pi a^{3}}}}\exp(-{\frac {Zr}{a}})$ eine symmetrische Funktion. $E_{a}=E_{b}=-Z^{2}E_{0}=-4E_{0};\;E_{0}=$ 13,6 eV; a ist der Bohr-Radius.

Die erste Korrektur zur Energie ist dann nach dem Störungs-Algorithmus:

$\delta E=\int d{\vec {r}}_{a}d{\vec {r}}_{b}U({\vec {r}}_{a},{\vec {r}}_{b})|g({\vec {r}}_{a})g({\vec {r}}_{b})|^{2}={\frac {5}{4}}{\frac {e^{2}}{a}}$ = 34 eV.

Das ergibt -75 eV für die totale Ionisierung des Heliums, verglichen mit experimentell -79 eV.

Danach wird eine Variation errechnet, die eine bessere Näherung von -77,5 eV ergibt. Die Funktion g erhält einen freien Parameter s, der die Abschirmung der Zentralladung durch das Partner-Elektron simuliert. $g({\vec {r}},x):=g({\vec {r}})$ mit Z ersetzt durch (Z-x). Der Parameter x wird eingestellt auf das Minimum des Erwartungswerts E(x) des Gesamt-H-Operators.

Ergebnis: $E(x={\tfrac {5}{16}})=-(Z-{\tfrac {5}{16}})^{2}E_{0}=$ -77,5 eV

Einzelheiten der Rechenschritte:

Funktionen $f(t)=t^{n}e^{-t}$ haben das Integral

$\int _{0}^{\infty }t^{n}e^{-t}dt=n!\quad (=\Gamma (n+1)$ wenn n keine ganze Zahl ist).

$\Gamma (1)=\Gamma (2)=1;\;\Gamma ({\tfrac {1}{2}})={\sqrt {\pi }}$.

Es folgt $\int _{0}^{\infty }t^{n}e^{-t/b}dt=b^{n+1}n!$

Die Funktion $f({\vec {x}})=e^{-r/b};\;r=|{\vec {x}}|$ hat das Volumenintegral $\int f({\vec {x}})\;d{\vec {x}}=4\pi \int _{0}^{\infty }r^{2}e^{-r/b}\;dr=8\pi b^{3}$.

Die Funktion $g=e^{-r/a}$ hat die Quadratintegral-Norm $\int g^{2}({\vec {x}})\;d{\vec {x}}=\pi a^{3}$ (setze b=a/2 in voriger Formel).

Diese Funktion g ist der Wasserstoff-Grundzustand mit a= Bohr-Radius, normiert also: $\phi ({\vec {x}})=(\pi a^{3})^{-1/2}e^{-r/a}$.

Es gibt eine Entwickung von $v(x,y)=1/|{\vec {x}}-{\vec {y}}|$ in Kugelfunktionen:

Sei $r_{+}=\max(|{\vec {x}}|,|{\vec {y}}|);\;r_{-}=\min(|{\vec {x}}|,|{\vec {y}}|);\;\theta ({\hat {x}},{\hat {y}})$ der Winkel zwischen beiden Richtungen.

${\frac {1}{|{\vec {x}}-{\vec {y}}|}}={\frac {1}{r_{+}}}\left[\sum _{\ell =0}^{\infty }({\frac {r_{-}}{r_{+}}})^{\ell }P_{\ell }(\cos \theta ({\hat {x}},{\hat {y}}))\right]$

Das Legendre-Polynom $P_{\ell }$ kann weiter in $Y_{\ell m}$ mit absoluten Winkeln zerlegt werden, wenn es sein muss:

$P_{\ell }(\cos \theta )={\frac {4\pi }{2\ell +1}}\sum _{m=-\ell }^{\ell }Y_{\ell m}^{*}({\hat {x}})Y_{\ell m}({\hat {y}})$

Die Abstoßung zweier 1s Elektronenwolken ist also ein Matrixelement mit dem Integral vom Typ:

$I=\int d{\vec {x}}\;d{\vec {y}}\;v({\vec {x}},{\vec {y}})e^{-r/b}e^{-s/b}$

wo $r=|{\vec {x}}|,s=|{\vec {y}}|$. Bei der Entwicklung des Coulomb-Terms in $P_{\ell }$ bleibt nur $\ell =0$ übrig, denn die Wellen-Integranden hängen nur von Radien ab, $P_{0}$ ist konstant 1 und alle übrigen Integrale verschwinden, da sie orthogonal zu diesem $P_{0}$ sind. Die Winkelintegration über zwei Kugeln bringt den Faktor $16\pi ^{2}$ mit. Die Integration über die zwei Radien sieht daher so aus:

$I(b)=16\pi ^{2}\int _{0}^{\infty }r^{2}e^{-r/b}[\int _{0}^{r}(s^{2}/r)e^{-s/b}\;ds+\int _{r}^{\infty }se^{-s/b}\;ds]\;dr$

Alle Integrale sind einfach aber zeitraubend für Langsamrechner, Ergebnis: $I(b)=16\pi ^{2}{\frac {5}{4}}b^{5}$.

Im Folgenden werden die Energien bezogen auf den Wasserstoff-Grundzustand $(-E_{0})$, sowie Längen auf den Bohr-Radius $a_{0}$. Mit einer Kernladung (Ze) hat der Ein-Elektron-Zustand die Energie $K(Z)=-Z^{2}E_{0}$ und Radius $a(Z)=a_{0}/Z$. Der Faktor $e^{2}$ in Coulomb-Termen kann als $(2a_{0}E_{0})$ kodiert werden.

Der Helium-Grundzustand hat Energie 2K(Z) von Kernbindung der zwei Elektronen, plus der Elektron-Elektron-Abstoßung A(Z), auzuwerten als das Matrixelement $\int dxdy\langle \phi (x)\phi (y)|V(x,y)|\phi (x)\phi (y)\rangle$. Das führt auf obiges Integral mal $e^{2}$, geteilt durch die Normierung:

$A(Z)=e^{2}I(a(Z)/2){\frac {e^{2}}{(\pi a(Z)^{3})^{2}}}=e^{2}{\frac {16\pi ^{2}}{\pi ^{2}a(Z)^{6}}}{\frac {5}{4}}{\frac {a(Z)^{5}}{32}}={\frac {e^{2}}{a(Z)}}{\frac {5}{8}}=E_{0}{\frac {a_{0}}{a(Z)}}{\frac {5}{4}}={\frac {5}{4}}E_{0}Z$.

Damit das geschätzte Energieniveau: $E(Z)/E_{0}=-2Z^{2}+{\tfrac {5}{4}}Z$.

Die zwei Elektronen schirmen eines fürs andere die Kernladung teilweise ab. Daher ist zu erwarten, dass eine Variation der 1s-Wellen zu einer kleineren Ladung Y < Z, das heißt zu größerem Radius a(Y), einen besseren, niedrigeren Erwartungswert für den Grundzustand aufweist. Die Welle zu Y ist aber keine Eigenfunktion mehr zur wahren Kernladung Z, daher muss der Erwartungswert von $[({\vec {p}})^{2}/(2m)-Ze^{2}/r]$ für $\exp(-r/a(Y))$ neu berechnet werden. Energien werden wieder geteilt durch die Referenz-Energie E_0.

Es gilt: $\int d{\vec {x}}\;(e^{-r/a})^{2}/r=(1/a)\int d{\vec {x}}\;(e^{-r/a})^{2}$. Es folgt der Coulombterm: $\langle -Ze^{2}/r\rangle /E_{0}=-2Za_{0}/a$.

In den kinetischen Term geht ein die Mittelung über den Laplace-Operator

$\Delta (e^{-r/a})=[f''(r)+2f'(r)/r]=[{\frac {1}{a^{2}}}-{\frac {2}{ar}}]e^{-r/a}$

Dies macht mit dem Integral von vorhin = $-1/a^{2}$. Geteilt durch $E_{0}$ kommt dimensionslos heraus: $\langle {\vec {p}}^{2}/2m\rangle /E_{0}=ma_{0}^{2}/a^{2}$.

Ersetzt man noch a=a(Y) durch das dimensionslose $Y,\;Y\cdot a(Y)=a_{0}$, ergibt sich das einfache Modell für die Bindungsenergie zum Kern, pro Elektron: $E(Y)/E_{0}=Y^{2}-2YZ$. Andererseits gilt von vorhin der Kalkül der Abstoßung: $A(Y)/E_{0}={\tfrac {5}{4}}Y$.

Die Variation besteht nun darin, Y = Z-x mit variablem x bis zum Minimum der Gesamtenergie herunterzufahren.

$f(x)=(2E(Y)+A(Y))/E_{0}=2(Z-x)^{2}-4(Z-x)Z+{\tfrac {5}{4}}(Z-x)$

Ergebnis fürs Mimimum: $f'(x)=0:\;x=5/16;\;f(x)=-2(Z-x)^{2}$.

Die Variationsrechnung verbessert also die Energie der ersten Störungs-Näherung. Diese bleibt aber immer noch über der experimentellen, ganz wie es theoretisch erwartet wird.

Wasserstoff-Bindung

Das einfachste System, das eine Vorstufe für ein Molekül darstellt, ist das Ion aus zwei Protonen und einem Elektron, das $H_{2}^{+}$ - Ion. Mit der Variationsmethode wird versucht zu zeigen, dass es einen optimalen Abstand R zwischen den Protonen gibt, wo sie vom gemeinsamen Elektron gebunden werden. Das heißt, verglichen mit dem Grundzustand von Wasserstoff-Atom plus weit entferntem Proton, ist da der Energie-Erwartungswert niedriger, wenn die Wellenfunktion geeignet überlagert wird aus denen, die beim einen oder anderen Kern konzentriert sind. Schon die einfache beste Linearkombination aus den zwei 1s-Zuständen mit Bohr-Radius führt zu einem Energie-Minimum, wenn R variiert wird. Eine verbesserte Rechnung mit zwei variablen Parametern, Abstand R und Radius der Elektronenwolken, ergibt schon den experimentell richtigen Wert des Protonen-Abstands.

Das angehängte Python-Skript macht eine rudimentäre Grafik-Datei, wo das Energie-Minimum (mit der besseren Variation) klar beim Proton-Abstand von zwei Bohr-Radien liegt. Die Superposition aus zwei Wellen produziert zwei Eigenzustände, den bindenden und den anti-bindenden, der auch gezeigt ist.

In dieser Näherung sind die Protonen unbeweglich. Nur das eine Elektron kommt in die Wellenfunktion. Man denke sich einen trägen Faktor für die Protonen, deren Freiheitsgrade abgespalten wurden von der theoretisch korrekten Dreiteilchen-Wellenfunktion.

Die Versuchsfunktion ist die Superposition von zwei Wasserstoff-1s-Wellen $\propto e^{-r/a}$, je an einem der zwei Kerne. Im ersten Anlauf wird nur der Abstand der Protonen variiert und der niedrigste Energie-Erwartungswert gesucht. Im verbesserten Schema wird noch der Bohr-Radius der 1s-Wellen geschrumpft, bei jedem Abstand optimal. Denn wenn die Protonen sich nahe kommen, sollte man vom Elektron einen ähnlichen Zustand wie an einem Helium-Kern erwarten, also mit kleinerem Radius.

Durchgeführt wird die Variationsrechnung mit nicht-orthonormierter Basis. Gegeben, linear unabhängige $|f_{i}\rangle$ und ihre Skalarpordukt-Matrix $S_{ik}=\langle fi|fk\rangle$.

Gesucht sind Linearkombinationen $\sum _{i}c_{i}|f_{i}\rangle$, die Eigenvektoren des symmetrischen (oder hermiteschen) H-operators sein sollen. Der kleinste der Eigenwerte wird der gesuchte Minimal-Erwartungswert E sein.

$H(\sum c_{k}f_{k})=E(\sum c_{k}f_{k})$. Skalarprodukte mit den f_i ergeben:

$\sum c_{k}\langle fi|H|f_{k}\rangle =E\sum _{i}S_{ik}\;c_{k}$.

Mit Matrix $H_{ik}:=\langle f_{i}|H|f_{k}\rangle$ ist also eine verallgemeinerte Eigenwertgleichung zu lösen: $\sum _{k}(H_{ik}-ES_{ik})c_{k}=0$.

Die Determinante der Matrix ([H]-E[S]) ist ein Polynom in E, dessen Nullstellen die gesuchten Energiewerte sind. Die aufwändige Etappe ist die Integralrechnung, um die Matrizen [S] und [H] herzustellen.

Im vorliegenden Fall sind H-Operator und Testfunktionen reell, die Matrizen symmetrisch, sogar $H_{11}=H_{22},S_{11}=S_{22}$.

In Dimension 2 wird die Eigenwertgleichung einfach:

${\text{Det}}{\begin{pmatrix}(H_{11}-ES_{11})&(H_{12}-ES_{12})\\(H_{12}-ES_{12})&(H_{11}-ES_{11}\end{pmatrix}}=0$

$(H_{11}-ES_{11})^{2}=(H_{12}-ES_{12})^{2}$

$E_{a}=(H_{11}-H_{12})/(S_{11}-S_{12});\quad E_{b}=(H_{11}+H_{12})/(S_{11}+S_{12})$

Etwas Integralrechnung

Von den kartesischen Koordinaten der Welle $\Psi (x,y,z)$ wird übergegangen zu elliptischen Koordinaten $(u,v,w)$ mit folgender Definition. Die Protonen sitzen auf der z-Achse bei Punkten $A=(0,0,R/2)$ und $B=(0,0,-R/2)$.

Nach folgendem Schema bekommt der Punkt $X=(x,y,z)$ neue Koordinaten $(u,v,w)$:

• $w=\arctan(y/x)$ sei der Winkel nach Projektion auf die xy-Ebene
• $r_{a}=|X-A|;\;r_{b}=|X-B|$ die Abstände zu den zwei Zentren A,B
• $u=(r_{a}+r_{b})/R;\;v=(r_{a}-r_{b})/R$ ; wo R der Abstand der zwei Zentren ist.

Die Erfindung von $u,v,w$ hat Vorteile bei der Integralrechnung, wenn Faktoren im Integranden um Punkt A zentralsymmetrisch sind und andere um Punkt B. Was sind die Wertebereiche, wie sieht das Integrationsmaß aus?

• $u:1\leq u<\infty$ denn $(r_{a}+r_{b})\geq R$
• $v:-1\leq v\leq 1$ denn $|r_{a}-r_{b}|\leq R$
• $w:0\leq w\leq 2\pi$
• $dx\;dy\;dz=(R/2)^{3}(u^{2}-v^{2})\;du\;dv\;dw$

(Ausrechnung noch im Mathe-Kapitel nachzuliefern).

Integranden vom Typ $f=(r_{a})^{n}\exp(-r_{a}/q_{a})(r_{b})^{m}\exp(-r_{b}/q_{b})$ werden nach Einsetzen $r_{a}=R/2(u+v);\;r_{b}=R/2(u-v)$ zu Summen von Termen des gleichen Typs in u und v. Jeder Term kann entkoppelt über u und v integriert werden.

$I=\int d^{3}x(r_{a})^{n}\exp(-r_{a}/q_{a})(r_{b})^{m}\exp(-r_{b}/q_{b})=\int _{1}^{\infty }du\int _{-1}^{1}dv\int _{0}^{2\pi }dw\;(R/8)^{3}(u^{2}-v^{2})f(x)=$

$2\pi (R/2)^{n+m+3}\int _{1}^{\infty }du\int _{-1}^{1}dv(u+v)^{n+1}(u-v)^{m+1}\exp(-[(u+v)R/(2q_{a})+(u-v)R/(2q_{b})])$

Die drei Integranden des Beispiels sind von der angesprochenen Sorte:

$S(R,q):=\int dxdydz\;\exp(-r_{a}/q)\;\exp(-r_{b}/q)$

$T(R,q):=\int dxdydz\;\exp(-r_{a}/q)\;(1/r_{b})$

$U(R,q):=\int dxdydz\;\exp(-r_{a}/q)/r_{a}\;\exp(-r_{b}/q)$.

Nach leidiger, langer Rechnerei voller einfacher partieller Integrationen:

$S(R,q)=2\pi (R/2)^{3}\int _{1}^{\infty }du\int _{-1}^{1}dv\;(u^{2}-v^{2})\exp(-Ru/q)$

$S(R,q)=(\pi R^{3})[1/3(q/R)+(q/R)^{2}+(q/R)^{3}]\exp(-R/q)$

$T(R,q)=2\pi (R/2)^{2}\int _{1}^{\infty }du\int _{-1}^{1}dv\;(u+v)\exp(-R(u+v)/(2q))$

$T(R,q)=\pi R^{2}[z^{3}-(z^{2}+z^{3})exp(-R/q)]\quad$ mit $z:=2q/R$

$U(R,q)=2\pi (R/2)^{2}\int _{1}^{\infty }du\int _{-1}^{1}dv\;(u+v)\exp(-Ru/q)$

$U(R,q)=\pi R^{2}((q/R)+(q/R)^{2})\exp(-R/q)$

Fortsetzung ($H_{2}^{+}$)-Ion

Der Hamilton-Operator enthält das Potenzial von zwei Kernen, die das Elektron anziehen, sowie die Abstoßung der beiden Protonen.

$H=-\hbar ^{2}\Delta /(2m)+e^{2}[-1/r_{1}-1/r_{2}+1/R]$

Versuchsfunktion: $f({\vec {x}})=f(r)=\exp(-r/a);\;\int d^{3}{\vec {x}}\;f^{2}(r)=\int d^{3}{\vec {x}}\exp(-2r/a)=\pi a^{3}$.

$\Delta (e^{-r/a})=[f''(r)+2f'(r)/r]=({\frac {1}{a^{2}}}-{\frac {2}{ar}})e^{-r/a}$

Sei ein Proton am Nullpunkt. Hat die Testfunktion eine stationäre Lösung?

$[{\frac {\hbar ^{2}}{mar}}-{\frac {\hbar ^{2}}{2ma^{2}}}-{\frac {e^{2}}{r}}]f=-E\;f$

Natürlich ist f eine Eigenfunktion, genau wenn a der Bohr-Radius ist; $e^{2}={\frac {\hbar ^{2}}{ma}};E={\frac {\hbar ^{2}}{2ma^{2}}}$.

Ab hier wird dimensionslos, relativ zur Wasserstoff-Energie E gerechnet. Mit $e^{2}/E=2a$ und $(\hbar ^{2}/(2m))/E=a^{2}$ ist alles im Bohr-Radius verpackt:

$H^{'}=(H/E)=-a^{2}\Delta +2a[-1/r_{1}-1/r_{2}+1/R]$

Die Testfunktionen $e^{-r/q}$ werden mit allgemeinem Radius q um a herum variiert.

Man nehme $\{f_{1}=e^{-r_{1}/q};f_{2}=e^{-r_{2}/q}\}$ als Basis der Testfunktionen. Dann hat $\Delta$ die Form wie oben, wenn a durch q ersetzt wird. Also wird

$H^{'}f_{1}=[{\frac {2a^{2}}{qr_{1}}}-{\frac {a^{2}}{q^{2}}}+2a(-1/r_{1}-1/r_{2}+1/R)]f_{1}$

$H^{'}f_{2}=[{\frac {2a^{2}}{qr_{2}}}-{\frac {a^{2}}{q^{2}}}+2a(-1/r_{1}-1/r_{2}+1/R)]f_{2}$

Es ist $N(q)=\pi q^{3}$ das Normquadrat und $S(R,q)$ das Skalarprodukt von $f_{1},f_{2}$. Damit ist die Matrix [S] der Variation fertig.

$S_{11}=S_{22}=N(q);\;S_{12}=S_{21}=S(R,q)$

Die Matrixelemente des Operators H' benutzen noch den Coulomb-Erwartungswert

$C(q)=\int d^{3}x(1/r)e^{-2r/q}=\pi q^{2}$.

Mit all den definierten Ausdrücken $N(q),C(q),S(R,q),T(R,q),U(R,q)$ wird als Funktion von 3 Parametern $(R,a,q)$:

$H_{11}^{'}={\frac {2a^{2}}{q}}C(q)-{\frac {a^{2}}{q^{2}}}N(q)-2aC(q)-2aT(q/2)+2(a/R)N(q)$

$H_{21}^{'}={\frac {2a^{2}}{q}}U(R,q)-{\frac {a^{2}}{q^{2}}}S(R,q)-2aU(R,q)-2aU(R,q)+2(a/R)S(R,q)$

Das Problem ist dimensionsbefreit: werden R,a,q mit demselben Faktor z skaliert, bleibt die Eigenwertgleichung gleich ($z^{3}$ in allen Matrixelementen). Man darf den Bohr-Radius a=1 setzen und bekommt die Energie als Funktion des Proton-Abstands in Einheiten von a. Das Skript setzt a=1, variiert erst R von 1 bis 4 mit q=a; die Protonen durchfahren 1 bis 4 Bohr-Radien Abstand. Dann noch einmal R von 1 bis 4, wo bei jedem Schritt der Wellenradius q probeweise von a auf 0,25a herunterfährt. Die Werte $E_{a}$ und $E_{b}$ nach obigen Formeln werden gespeichert, eine schnelle unordentliche Grafik entsteht.

Beispiel Stark-Effekt am Wasserstoffatom

Es soll berechnet werden, wie sich die Energie-Niveaus des H-Atoms ändern, wenn ein äußeres homogenes, konstantes elektrisches Feld einwirkt. Das Stör-Potenzial ist $H_{1}=qFz$, F bezeichnet die elektrische Feldstaärke. Mit dem Bohr-Radius $a:H_{1}=qFa\;(z/a)$. Gebraucht wird also die Matrix eines dimensionslosen Operators $(z/a)$ zwischen allen Eigenzusänden des Atoms. Diese haben eine Basis aus den folgenden Wellenfunktionen:

$|n,l,m\rangle =\exp(im\phi )P_{lm}(\theta )r^{l}\exp(-r/na)R_{nl}(r/a);\quad P_{l,-m}=P_{lm}$

$n=1,2,3,\dots ;\;l=0,1,\dots ,n-1;\;m=-l,-l+1,\dots ,l$

$R_{nl}(z):=L_{n-l-1}^{2l+1}(2z/n)$ in Zugeordneten Laguerre-Polynomen.

Die entarteten Energie-Eigenwerte bilden die Balmer-Folge $E_{n}=E_{1}/n^{2}$.

Das Integral in Kugelkoordinaten: $\int d{\vec {x}}=\int _{0}^{2\pi }d\phi \int _{0}^{\pi }\sin \theta \;d\theta \int _{0}^{\infty }r^{2}\;dr$

Der Stör-Operator in Kugelkoordinaten: $(z/a)=(r/a)\cos \theta$.

Operator $H_{1}$ vertauscht mit Drehimpuls-Komponente $L_{3}$, deren Eigenwerte die Quantenzahlen m der Tripel (nlm) der Wasserstoff-Niveaus sind.

Diagonal-Elemente von (z/a) verschwinden, da der Integrand über den Winkel $\theta$ dann $(\sin \theta \,\cos(\theta ))$ mal ein Polynomquadrat $Q^{2}(\cos \theta ,\sin \theta )$ ist, also antisymmetrisch um die Mitte $\theta =\pi /2$ herum. Allgemeiner hat das Produkt zweier $P_{lm}(s,c)$ bei gleichem m nur Terme von gerader Potenz, wenn die Differenz der Quantenzahl l gerade ist. Auch dann verschwinden die Matrixelemente von (z/a).

Die erste Ordung ist daher nicht-diagonal und besteht nur aus den Eigenwerten der Stör-Matrix innerhalb der entarteten Mengen {(n,m)=const, l variabel}. Ist ein Unterraum wie (n,m)= (1,0) oder (2,1) oder (3,-2) nur eindimensional, entfällt der Stark-Effekt in erster Ordnung.

Für m=0, n=3 gibt es ohne $H_{1}$ drei entartete Zustände $\{l=0...n-1;\;|m|\leq l\}$. Jeder Unterraum m=const ist invariant unter $H_{1}$ und verschiedene m haben Matrixelemente Null mit $H_{1}$. Wenn $H_{1}$ die Entartung aufhebt, dann errechnet sich die erste Ordnung der Aufspaltung durch Diagonalisieren der Matrix, zum Beispiel der 3x3-Matrix von $H_{1}$ auf dem Unterraum $\{(n,m)=(3,0);\;l=0,1,2\}$.

Das lästige Ausrechnen der Matrixelemente und Normen wurde automatisiert. Ein Python-Skript baut erst die unnormierten Radialfunktionen $R_{nl}(r)$ und Winkel-Polynome $P_{lm}(s,c)=P_{lm}(\sin \theta ,\cos \theta )$ aus deren Rekursions- Rezepten. Es beherrscht die Arithmetik mit rationalen Zahlen (Paare ganzer Zahlen) und deren Quadratwurzeln (Quadrupel), etwas Manipulation von entsprechenden Polynomen, sowie die nötige Integralrechnung in den Koordinaten Radius und Winkel \theta. Die Normquadrate und Normen der Basisfunktionen werden ausgerechnet. Die Teil-Matrizen des Stör-Operators (z/a) für relevante Paare (n,m), mit normierten Wellen, sind dann das Endergebnis des Skripts. Hier einige:

• Polynome der Radialfunktionen $R_{nl}(r/a)$, Radial-Normen für a=1.

${\begin{array}{ll}R_{10}=1&{\sqrt {\frac {1}{2}}}\\R_{20}=1-{\frac {1}{2}}\,r^{1}&{\sqrt {2}}\\R_{21}=1&2{\sqrt {6}}\\R_{30}=1-{\frac {2}{3}}\,r^{1}+{\frac {2}{27}}\,r^{2}&3{\sqrt {\frac {1}{2}}}\\R_{31}=1-{\frac {1}{6}}\,r^{1}&{\frac {27}{4}}\\R_{32}=1&{\frac {81}{2}}{\sqrt {5}}\\R_{40}=1-{\frac {3}{4}}\,r^{1}+{\frac {2}{16}}\,r^{2}-{\frac {1}{192}}\,r^{3}&2{\sqrt {2}}\\R_{41}=1-{\frac {1}{4}}\,r^{1}+{\frac {1}{80}}\,r^{2}&8{\sqrt {\frac {6}{5}}}\\R_{42}=1-{\frac {1}{12}}\,r^{1}&32{\sqrt {10}}\\R_{43}=1&384{\sqrt {70}}\\\end{array}}$

• $P_{lm}(\theta )$ als Polynome von (s,c) mit $s=\sin \theta ,\;c=\cos \theta$. Normen.

${\begin{array}{ll}P_{00}=1\,&{\sqrt {2}}\\P_{10}=1\,c&{\sqrt {\frac {2}{3}}}\\P_{11}=-1\,s&2{\sqrt {\frac {1}{3}}}\\P_{20}=-{\frac {1}{2}}\,+{\frac {3}{2}}\,c^{2}&{\sqrt {\frac {2}{5}}}\\P_{21}=-3\,sc&2{\sqrt {\frac {3}{5}}}\\P_{22}=3\,s^{2}&4{\sqrt {\frac {3}{5}}}\\P_{30}=-{\frac {3}{2}}\,c+{\frac {5}{2}}\,c^{3}&{\sqrt {\frac {2}{7}}}\\P_{31}={\frac {3}{2}}\,s-{\frac {15}{2}}\,sc^{2}&2{\sqrt {\frac {6}{7}}}\\P_{32}=15\,s^{2}c&4{\sqrt {\frac {15}{7}}}\\P_{33}=-15\,s^{3}&12{\sqrt {\frac {10}{7}}}\\P_{40}={\frac {3}{8}}\,-{\frac {15}{4}}\,c^{2}+{\frac {35}{8}}\,c^{4}&{\frac {1}{3}}{\sqrt {2}}\\P_{41}={\frac {15}{2}}\,sc-{\frac {35}{2}}\,sc^{3}&{\frac {2}{3}}{\sqrt {10}}\\P_{42}=-{\frac {15}{2}}\,s^{2}+{\frac {105}{2}}\,s^{2}c^{2}&4{\sqrt {5}}\\P_{43}=-105\,s^{3}c&4{\sqrt {70}}\\P_{44}=105\,s^{4}&16{\sqrt {35}}\\\end{array}}$

Die Matrizen M des (z/a)-Operators, $\langle n,m,l|(z/a)|n,m,l'\rangle$, enthalten die gesuchten Energie- Verschiebungen x als Eigenwerte, also als Nullstellen der Polynome $det(M-x\mathbf {1} )=0$, die so genannte Säkular-Gleichung.

Matrix n=2 m=0 Liste l= 0 1

${\begin{pmatrix}0&-3\\-3&0\end{pmatrix}}$

Matrix n=3 m=0 Liste l= 0 1 2

${\begin{pmatrix}0&-3{\sqrt {6}}&0\\-3{\sqrt {6}}&0&-3{\sqrt {3}}\\0&-3{\sqrt {3}}&0\end{pmatrix}}$

Matrix n=3 m=1 Liste l= 1 2

${\begin{pmatrix}0&{\frac {-9}{2}}\\{\frac {-9}{2}}&0\end{pmatrix}}$

Matrix n=4 m=0 Liste l= 0 1 2 3

${\begin{pmatrix}0&-6{\sqrt {5}}&0&0\\-6{\sqrt {5}}&0&{\frac {-24}{5}}{\sqrt {5}}&0\\0&{\frac {-24}{5}}{\sqrt {5}}&0&{\frac {-18}{5}}{\sqrt {5}}\\0&0&{\frac {-18}{5}}{\sqrt {5}}&0\end{pmatrix}}$

Matrix n=4 m=1 Liste l= 1 2 3

${\begin{pmatrix}0&{\frac {-36}{5}}{\sqrt {\frac {5}{3}}}&0\\{\frac {-36}{5}}{\sqrt {\frac {5}{3}}}&0&{\frac {-24}{5}}{\sqrt {\frac {5}{2}}}\\0&{\frac {-24}{5}}{\sqrt {\frac {5}{2}}}&0\end{pmatrix}}$

Matrix n=4 m=2 Liste l= 2 3

${\begin{pmatrix}0&-6\\-6&0\end{pmatrix}}$

Für die Niveaus n=2,m=0 sind Eigenwerte $\pm 3$ der Energie-Verschiebung möglich, gerechnet in der Einheit (qFa). Da $n=2,m=\pm 1$ unverändert bleiben, spaltet sich n=2 insgesamt in drei Niveaus mit Abstand 3 Einheiten.

Für $n=3,m=\pm 1$ gibt es Verschiebungen von $\pm 9/2$ Einheiten, keine für $n=3,m=\pm 2$. Und für n=3,m=0 erscheinen als Eigenwerte einer 3x3-Matrix, $\pm 9$ Einheiten. Insgesamt kann n=3 fünf Niveaus liefern, im Abstand von 4,5 Einheiten.

Für n=4 sieht es nach 7 Niveaus aus; die Eigenwerte der drei berechneten Matrizen liegen bei $\pm 6,\;\pm 12,\;\pm 18$ Einheiten.

Für die Korrektur zweiter Ordnung etwa an (n,m)=(0,0) bestimmt man die Matrixelemente von $H_{1}$ zwischen (0,0,0) und allen anderen (n,l,0) und summiert deren Betragsquadrate, geteilt durch die Energiedifferenz nach der Balmer-Serie; ganz nach Vorschrift. Das numerische Resultat ist sehr klein, während die Aufspaltung erster Ordnung, also zum Beispiel für n=3 in fünf mögliche Niveaus bei hohen elektrischen Feldern durchaus im Spektrum zu sehen ist.

Es gibt eine Rechenmethode in Parabelkoordinaten, die das Wasserstoff-Spektrum statt in den Quantenzahlen (n,l,m) in alternativen Tripeln $(k_{1},k_{2},m)$ ausdrückt. Man fängt ganz von vorne an und erhält eine andere Eigenfunktions-Basis. In einer solchen Matrixdarstellung ist $H_{1}$ praktisch diagonal, seine Erwartungswerte sind proportional zu $(k_{1}-k_{2})$, wobei $n=k_{1}+k_{2}+m+1$.

Python-Skript

# hydrogen atom eigenfunctions and stark effect perturbation

def hasfact(x,z) : return (x>0) and ((x % z)==0)
def div(p,q) : return int(p/q)

def division(p,q) : # to print integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return s*p,q

''' sphericals
- integ(0,pi) dt sin(t) Pl'm'(cos t) Plm(cos t) = 
-  2/(2n+1) (l+m)!/(l-m)! delta(l,l') delta(m,m')

z= cos(t), x=sin(t). shorten z=c, x=s. 
l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
Pl=Pl0. recurs Pl(z):= (1/2^l l!) d^l (zz-1)^l
Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1
'''

def getprimes(n) : # want n prime numbers
  p=[2]; i=1; plast=2
  while i<=n :
    q=plast+1; found=False
    while not found:
      j=0
      while ((q%p[j])!=0) and ((p[j]*p[j])<=q) :  j +=1
      if (q%p[j] != 0) : p += [q]; plast=q; found=True; i+=1
      else : q+=1
  return p

def greatestdiv(x) : # greatest divisor, skip 0s in list
  n=len(x); y=[0]*n; w=[0]*n; h=0; ymin=-1; gcd=1
  pn=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]; ln= len(pn) # small primes
  for i in range(n) :
    z=abs(x[i])
    if z>0 : ymin= z if(ymin<0) else min(ymin,z); w[h]=z; h+=1
  w=w[:h]; ip=0
  while (ip<ln)and(pn[ip]<=ymin) :
    z=pn[ip]; ok=True
    while ok :
      for i in range(h) : ok= ok and ((w[i]%z)==0)
      if ok : 
        for i in range(h) : w[i]= div(w[i],z)
        ymin= div(ymin,z); gcd *= z
    ip +=1 
  return gcd

def smallestmult(x) : # smallest multiple of nonzero numbers
  n=len(x); y=abs(x[0])
  for i in range(1,n) : 
    z=abs(x[i]); y = y if (z==0) else div(y*z,greatestdiv([y,z]))
  return y

def divide(p,q) : # improve integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  if p==0 : q=1
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return [s*p,q]

def polynscale(u,p,q) : # scale factor on polynomial
  n=len(u); z=[[]]*n
  for i in range(n) : z[i]= divide(p*u[i][0],q*u[i][1])
  return z

def mult(a,b) : # multiply rationals or integers, result rational
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return [pa*pb,qa*qb]

def add(a,b) : # add rational or int
  ta=str(type(a)); inta= (ta.find('int')>0)
  tb=str(type(b)); intb= (tb.find('int')>0)
  pa,qa= (a,1) if inta else (a[0],a[1])
  pb,qb= (b,1) if intb else (b[0],b[1])
  return divide(pa*qb+pb*qa,qa*qb)

def polynadd(u,v,a,b) : # rational-coeff polynomial added, with int? factors
  n=len(u); m=len(v); k=max(n,m); w=[[]]*k
  for i in range(k): 
    uu= mult(a,u[i]) if (i<n) else [0,1]
    vv= mult(b,v[i]) if (i<m) else [0,1]
    q=uu[1]*vv[1]; p=uu[0]*vv[1]+vv[0]*uu[1]; w[i]= divide(p,q)
  return w

def polynmul(u,v) : # multiply polynomials
  n=len(u); m=len(v); w=[[]]*(n+m-1)
  for i in range(n) :
    for k in range(m) :
      a=u[i][0];b=u[i][1]; c=v[k][0];d=v[k][1]; z= divide(a*c,b*d); j=i+k
      if w[j]==[] : w[j]=z # first use of j
      else : p=z[0]*w[j][1]+z[1]*w[j][0]; q=z[1]*w[j][1]; w[j]=divide(p,q)
  return w

def polyndiff(u) : # derivative of polynomial
  n=len(u); w=[[]]*(n-1)
  for i in range(1,n) : w[i-1]= [i*u[i][0],u[i][1]]
  return w

def polynpurge(u) : # omit leading zero terms and scale to integers
  # return new polyn and scalefactor applied
  n=len(u); j=n-1; k=0 
  while (u[j][0]==0)and (j>=0) : j -=1
  while (u[k][0]==0)and (k<n) : k+=1
  j +=1; v= u[:j]; p=[0]*j; q=[0]*j; sign= (-1) if (u[k][0]<0) else 1
  for i in range(j) : p[i]= v[i][0]; q[i]=v[i][1]
  z=smallestmult(q) 
  for i in range(j) : r= divide(p[i]*z,q[i]); p[i]=r[0]
  y=greatestdiv(p); 
  for i in range(j) : r=divide(p[i],y); v[i]=[sign*r[0],1]
  return v, divide(z,y)

def present(s,q,mode,noisy) :
  # mode 0 no x, mode 1 complementary x, mode -n const x^n, mode>0 l=mode 
  n=len(q); dmax=1; j= (n-1) if (mode==1) else 0; s+='('
  if mode>0 : j=mode # new
  v= '' if (mode >=0) else ('x'+str(-mode))
  for i in range(n) : dmax= q[i][1] if (q[i][1]>dmax) else dmax
  for i in range(n) :
    a=q[i][0]; b= q[i][1]; bad= (dmax % b) != 0; sig=''
    if bad : print('present bug exit'); exit()
    t=('x'+str(j)) if (j>0) else ''; u=('z'+str(i)) if (i>0) else ''
    if a!=0 : cf= div(dmax,b)*a; sig='+' if (cf>0) else ''
    if a!=0 : s+= ' '+sig+str(cf)+v+t+u
    j -=1
  s+=')/'+str(dmax)
  if noisy : print(s)
  return (s+'\n') 

def plmore(pl,l,noisy) :
  #Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1
  # option: homogenize with products (xx+zz) when deficient
  ql=pl;pwr=1;sign=-1; pi=''; qlist=[[l,0,[]+pl]]
  for i in range(1,l+1) :
    ql= polyndiff(ql); pwr += 1; ql=polynscale(ql,sign,1); qlist+=[[l,i,ql]]
    pi += present('Plm['+str(l)+str(i)+']=', ql, -i, noisy)
    # power i of x included, z^k must get (l-(k+i))/2 factors (1+zz). 
  return pi,qlist
 
def plzero(n,noisy) : # rational polys of z
  #  l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
  #  l Pl(z,x)= (2l-1) z P(l-1)(z,x) - (l-1)(zz+[1~=xx])P(l-2) 
  pl=[[]]*20; pl[0]= [[1,1]]; z= [[0,1],[1,1]]; un=[[1,1],[0,1],[1,1]]
  pli='Plm[0,0]=(1)\n'; 
  pllist=[[0,0,[[1,1]] ]] # l=0 m=0 constant
  for i in range(1,n) :
    y= polynmul(pl[i-1],z); y= polynscale(y, 2*i-1,1)
    if (i-2)>= 0 : w= polynscale(pl[i-2],i-1,1); y=polynadd(y,w,1,-1)
    pl[i]= polynscale(y,1,i); pli += present('Plm['+str(i)+'0]=',pl[i],0,noisy)
    pi,qlist = plmore(pl[i],i,noisy); pllist += qlist 
    pli+=pi
  return pl,pli,pllist

'''
  I(m,n)= integ 0 to pi of sin^m cos^n [ n even ]
  - m += 2 implies  I(m+2,n) = (m+1)/(m+n+2) *I(m,n)
  - n += 2 implies  I(m,n+2) = (n+1)/(m+n+2) *I(m,n)
  - start functions : 1, sin, cos, sincos: 0,0 1,0 0,1, 1,1 
  0,0: pi. 1,0: 2 0,1 : 0  1,1 : 0
'''
def isico0pi(m,n) : # integral 0 to pi of sin^m(x) cos^n(x)
  if (n%2 != 0) : return [0,1] # ugly trick: negative val encode pi factor
  nn=0; mm= 0 if (m%2 == 0) else 1; p=-1 if (mm==0) else 2; q=1 
  while nn<n : p *= (nn+1); q *= (mm+nn+2); nn += 2
  p,q= tuple(divide(p,q))
  while mm<m : p *= (mm+1); q *= (mm+nn+2); mm += 2
  return divide(p,q)

def csproduct(n,a,b,c, m,d,e,f) : # product of 2-var polynomials p(x,z)
  # length n and m. a,d= lists of coeffs, b,e powers of x, c,f powers of z.
  # returns list of coeffs and list of power pairs they target
  # isico0pi uses these power pairs, results to be summed up then
  # variant: a,d are pairs=rationals.
  def ratpr(a,b) : return [a[0]*b[0],a[1]*b[1]]
  def ratsum(a,b) :
    if a==[] : a=[0,1]
    p=a[0]*b[1]+a[1]*b[0];q=a[1]*b[1]; return divide(p,q)
  def getindex(ix,p,q) : # where to store x^p y^q
    i=0; n=len(ix); found=False
    while (not found) and (i<n) and (ix[i]!=[]) :
      found= (ix[i][0]==p)and(ix[i][1]==q); i= (i+1) if (not found) else i
    if not found : ix[i]=[p,q] # new exponent combination
    return i 
  t=str(type(a[0])); ratnum= (t.find('list')>0)
  x= ([[]]*(n+m)) if ratnum else ([0]*(n+m)) 
  ix=[[]]*(n+m); jmax=0
  for i in range(n) :
    for k in range(m) :
      yy=b[i]+e[k]; zz=c[i]+f[k]; j=getindex(ix,yy,zz) # side effect on ix
      if ratnum : x[j]= ratsum(x[j],ratpr(a[i],d[k]))
      else : x[j] += a[i]*d[k]
      if j>jmax : jmax=j
  return x[:(jmax+1)],ix[:(jmax+1)]

def intprod(pa,ia,la, pb,ib,lb) : # integral of product of 2 plm(x,z) polyns
  # pa(x,z)= sum(k,pa[k]*x^ia*z^k, pb same. nonhomog, const x power=i+j
  # status ok for plzero
  na=len(pa); cfa=[0]*na; pxa=[0]*na; pza=[0]*na
  for i in range(na) : cfa[i]=pa[i]; pxa[i]=ia; pza[i]=i
  nb=len(pb); cfb=[0]*nb; pxb=[0]*nb; pzb=[0]*nb
  for i in range(nb) : cfb[i]=pb[i]; pxb[i]=ib; pzb[i]=i
  cf,pw= csproduct(na,cfa,pxa,pza, nb,cfb,pxb,pzb); nt=len(cf); sum=[]
  # print('intprod cf,pw='+str([cf,pw]))
  for j in range(nt) :
    pw[j][0] +=1 # sin factor integral measure
    z= isico0pi(pw[j][0],pw[j][1]); a=z[0]; b=z[1]
    z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
    sum= polynadd(sum,z, 1, cf[j])
  # print('pa,sum='+str([pa,sum]))
  return sum

'''
hydrogen eigenfuncs
 |n,l,m> = \exp(im\phi)P_{lm}(\theta) r^l \exp(-r/na) Rad{nl}(r/na)
 Rad{nl}(z) := L_{n-l-1}^{2l+1}(2z) associated Laguerre polynomial
  int_0^\infty\;dr\; r^m exp(-r/na)) = (na)^(m+1) m!
 unnormalized (nlm) is exp(im\phi)Plm(s,c) L()(2s) exp(-s)*r^l; s=r/na
   here L(2s) = polyn(r/n), stored as Rnl(r). Bohrradius a=1.
 radial funcs for different n same l must be orthogonal, checked.
 L(q-p,p)= p-th deriv of L(q) poly of deg q.

perturbation z: inspect n=1,2,3,4 m=0..n-1 submatrices.
  matrix elements with z=r \cos\theta
'''

def sqroot(pq) : # root-of-rationals are size-4 int lists
  def hasfact(u,v) : return (u>=v) and ((u % v)==0)
  pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
  p,q= tuple(pq) ; rp=1; rq=1
  if (p<0)or(q<=0) : print('sqroot bad arg. exit'); exit()
  if p==0 : q=1; ln=0
  for i in range(ln) :
    z=pn[i]; y=z*z
    while hasfact(p,y) : rp *= z; p=div(p,y) #;print('rp,i,p,y'+str([rp,i,p,y]))
    while hasfact(q,y) : rq *= z; q=div(q,y) #;print('rq,i,q,y'+str([rq,i,q,y]))
  return [p,q,rp,rq] # meaning: sqrt(p/q)*(rp/rq)

def normalize(x,na,nb) : # x is rational, na nb are square-roots, 4-component
  if len(x)==1 : x=x[0] # some x are inside a one-element list
  if len(x) != 2 : print('wrong type of argument.exit.'); exit()
  xrp,xrq=tuple(x); sign,xrp= (1,xrp) if (xrp>0) else (-1,-xrp)
  ap,aq,arp,arq=tuple(na); bp,bq,brp,brq=tuple(nb)
  xp=aq*bq; xq=ap*bp; xrp *= (arq*brq); xrq *= (arp*brp)
  yp,yq,yrp,yrq= tuple(sqroot([xp,xq]))
  return divide(yp,yq) + divide(sign*xrp*yrp,xrq*yrq)

def mulroot(x,y) : # two square-root rationals
  xp,xq,xrp,xrq= tuple(x); yp,yq,yrp,yrq= tuple(y)
  zp,zq,zrp,zrq= tuple(sqroot([xp*yp,xq*yq])); zrp*= xrp*yrp; zrq *= xrq*yrq
  return divide(zp,zq) + divide(zrp,zrq)

def stroot(x) : # root-number to string
  p,q,rp,rq= tuple(x); s='0'
  if (p!=0)and(rp!=0) :
    t= str(rp) if(rq==1)   else (str(rp)+'/'+str(rq))
    u= str(p)  if (q==1)   else (str(p)+'/'+str(q))
    s= t       if (u=='1') else (t+'*sqrt('+u+')')
    if (t=='1')and(u != '1') : s= 'sqrt('+u+')'
  return s

def lxroot(x) : # root-number to latex-string
  p,q,rp,rq= tuple(x); s='0'
  if (p!=0)and(rp!=0) :
    t= str(rp) if(rq==1)   else ('\\frac{'+str(rp)+'}{'+str(rq)+'}')
    u= str(p)  if (q==1)   else ('\\frac{'+str(p)+'}{'+str(q)+'}')
    s= t       if (u=='1') else (t+'\\sqrt{'+u+'}')
    if (t=='1')and(u != '1') : s= '\\sqrt{'+u+'}'
  return s

def lxmatrix(m) : # m is 2-dim array of strings
  nrow= len(m); ncol=len(m[0]) # print('nrow,ncol='+str([nrow,ncol])) 
  s=':<math>\\begin{pmatrix}\n'
  for i in range(nrow) :
    for k in range(ncol) : 
      t= m[i][k]; u= (' & '+t) if (k>0) else t; s +=u
    s += ('\\\\\n') if (i<(nrow-1)) else '\n'
  s +='\\end{pmatrix}</math>\n'; return s

def lxradials(rnl) : # list with normsquare int r^(2l+2) R^2(r) exp(-2r/n)
  n=len(rnl); s='<math>\\begin{array}{ll}\n'
  for i in range(n) :
    nn,nl,p,nm,nmr = tuple(rnl[i]); t='R_{'+str(nn)+str(nl)+'}='; d=len(p)
    for k in range(d) :
      x,y= tuple(p[k]); x,sx= (x,'+') if (x>=0) else (-x,'-')
      u= str(x) if(y==1)  else ('\\frac{'+str(x)+'}{'+str(y)+'}')
      if (k==0)and(sx=='+') : sx=''
      t += sx+u; t+= ('\\,r^{'+str(k)+'}') if (k>0) else ''
    s += t+' & '+lxroot(nmr)+'\\\\\n'
  s +='\\end{array}</math>\n'; return s

def lxsphericals(plm) : # format Plm(s,c)= s^|m| c^k, s,c= sin,cos (theta)
  n=len(plm); s='<math>\\begin{array}{ll}\n'
  for i in range(n) :
    nl,nm,p,nr,nrr = tuple(plm[i]); t='P_{'+str(nl)+str(nm)+'}=';d=len(p);pfx=0
    for k in range(d) :
      x,y= tuple(p[k]); x,sx= (x,'+') if (x>=0) else (-x,'-')
      u= str(x) if(y==1)  else ('\\frac{'+str(x)+'}{'+str(y)+'}')
      if (pfx==0)and(sx=='+') : sx=''
      sc= ('s^'+str(nm)) if(nm>1) else ('s' if (nm==1) else '') 
      sc+=('c^'+str(k)) if (k>1) else ('c' if(k==1) else '')
      if x != 0 : t += sx+u +'\\,'+sc; pfx=1 # pfx=1 show prefix even if +
    s += t+' & '+lxroot(nrr)+'\\\\\n'
  s +='\\end{array}</math>\n'; return s

def radialpolys(lmax=3,mode=0) : # coeffs of radial hydrogen atom polyns
  # related to associated Laguerre L(n-l-1;2l+1). L(1-1,1) is poly of deg zero
  pols=[]  # n=1,2,3,4, l=0..n-1, polynomials in (r/a). 
  for n in range(1,lmax+1) :
    for l in range(n) :
      j=2*l+2; k=0; p=1; q=1; t=''; pol=[]
      while p != 0 :
        p,q= division(p,q); pol += [[p,q]]
        t += ' '+ (str(p) if(q==1) else (str(p)+'/'+str(q)))
        q *= (k+1)*(k+j); p *= (-2*n+2*k+j); k+=1 
      #print('n='+str(n)+' l='+str(l)+' '+t)
      j=len(pol); q=1 # pol has args ~ r/(an).  multiply divisors by n^h
      for h in range(j) : pol[h][1] *= q; q *= n
      pols += [[n,l,pol,[0,1],[]]]; # print(str([n,l,pol]))
  return pols

def radialprd(pols,h,ia,ib,mode) :
  na,nb,la,lb= pols[ia][0],pols[ib][0], pols[ia][1],pols[ib][1] 
  a=1; b= [na*nb, na+nb] # exp [-r/(na a) - r/(nb a) ]
  pa=pols[ia][2]; pb=pols[ib][2]; q=polynmul(pa,pb); sum=[0,1]
  #print('pa='+str(pa)+' pb='+str(pb)) 
  #print('b='+str(b)+' q='+str(q))
  for i in range(len(q)) :
    p=la+lb+i+h+2 # f=[b[0],b[1]] # p true power of r, add q[i]*b^(p+1) p!?
    f=[1,1] #  optional const prefactor b ?
    for k in range(1,p+1) : f= mult(f,b); f[0] *= k # = b^p p!
    z=mult(q[i],f); sum=add(sum,z) 
    #print('i,p,k='+str([i,p,k])+' q[i]='+str(q[i])+' f='+str(f)+' z='+str(z))
  if (h==0)and(na==nb)and(la==lb) :
    pols[ia][3]= sum; pols[ia][4]= sqroot(sum) # append normsquare
  if mode==0 : print('radial h nl nl sum',str([h,na,la,nb,lb,sum]))
  return sum

def radialprod(pols, h, na,la, nb,lb,mode) : # radial matrix of operator r^h
  # polyn for na has coeffs for r'= r/(na), already divided
  # radial integral adds factor r^(2+la+lb)  mix should be (r/na)^la?
  ia=0; ib=0; m=len(pols); # print('h,na,la,nb,lb='+str([h,na,la,nb,lb]))
  while (ia<m) and ((pols[ia][0] != na)or(pols[ia][1] != la)) : ia +=1
  while (ib<m) and ((pols[ib][0] != nb)or(pols[ib][1] != lb)) : ib +=1
  #print('ia,ib='+str(ia)+','+str(ib))
  return radialprd(pols,h, ia,ib,mode)

def normsquared(pols,nmax,mode) :
  for n in range(nmax+1) :
    for l in range(n) : sum= radialprod(pols,0, n,l, n,l,mode)

def crossradial(pols,na,nb,mode) :
  bugs=0 # different n, same l must be orthogonal
  for la in range(na) :
    for lb in range(nb) :
      sum= radialprod(pols,0, na,la, nb,lb,mode); z=sum[0]
      if (na!=nb)and(la==lb)and(z!=0): print('BUG'); bugs +=1
  print('BUGS for ('+str(na)+','+str(nb)+') ='+str(bugs))

def getradials(nmax=4,mode=0) : # radial polyns
  pols=radialpolys(nmax,mode)
  normsquared(pols,nmax,mode)
  np=len(pols); pmx=[[]]*np # symm perturbation matrix
  for i in range(np) : pmx[i] =[[]]*np #; print(str(pols[i]))
  for i in range(np) : # later: divide each by sqrt(normsquares)
    for k in range(i+1) : 
      pmx[i][k]= radialprd(pols, 1,i,k,mode); pmx[k][i]= pmx[i][k]
  #for i in range(np) : print(str(pols[i][4]))
  return pols, pmx

def integplms(pa,ia, pb,ib, fc) :
  # integral of product of 2 plm(x,z) polyn, times sin * (cos power fc)
  # pa(x,z)= sum(k,pa[k]*x^ia*z^k, pb same. x=sin, z=cos. const x-power
  na=len(pa); cfa=[0]*na; pxa=[0]*na; pza=[0]*na
  for i in range(na) : cfa[i]=pa[i]; pxa[i]=ia; pza[i]=i
  nb=len(pb); cfb=[0]*nb; pxb=[0]*nb; pzb=[0]*nb
  for i in range(nb) : cfb[i]=pb[i]; pxb[i]=ib; pzb[i]=i
  cf,pw= csproduct(na,cfa,pxa,pza, nb,cfb,pxb,pzb); nt=len(cf); sum=[]
  # print('intprod cf,pw='+str([cf,pw]))
  for j in range(nt) :
    pw[j][0] +=1 # sin factor = integral measure
    pw[j][1] += fc # cos factor for z operator
    z= isico0pi(pw[j][0],pw[j][1]); a=z[0]; b=z[1]
    z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
    sum= polynadd(sum,z, 1, cf[j])
  # print('pa,sum='+str([pa,sum]))
  return sum

def getsphericals(lmax=4,mode=0) :
  # return list of Plm and matrix for z operator. mode==0 debug output
  plz,plri,pllist= plzero(lmax+1,mode==0); ipl=len(pllist); mpl=[[]]*ipl
  # plz contains the Legendre polyn and associates, alt. pllist is used
  for i in range(ipl) : # recursion, nonnormalized Plm, and their integrals
    l=pllist[i][0]; m=pllist[i][1]; plm= pllist[i][2]
    y= integplms(plm,m, plm,m, 0); ry= sqroot(y[0]) # explicit integral
    if mode==0 : print(str(pllist[i])+' y='+str(y)+' ry='+str(ry))
    pllist[i].append(y); pllist[i].append(ry) # store with pllist
    mpl[i]= [[]]*ipl # matrix element mpl[i][k] of z=r cos(theta) operator 
  for i in range(ipl) :
    la=pllist[i][0]; ma=pllist[i][1]; pa= pllist[i][2]
    for k in range(i+1) :
      lb=pllist[k][0]; mb=pllist[k][1]; pb= pllist[k][2]
      y= integplms(pa,ma, pb,mb, 1); mpl[i][k]=y; mpl[k][i]=y
  return pllist, mpl

def matdisplay(m) : # m is 2-dim array of strings
  nrow= len(m); ncol=len(m[0]); wcol=[0]*ncol; sp=' '*80; s=''
  for i in range(nrow) :
    for k in range(ncol) : w= len(m[i][k]); wcol[k]= max(w,wcol[k])
  for i in range(nrow) :
    for k in range(ncol) : 
      t=m[i][k]; w=len(t); t= sp[:(wcol[k]+2-w)] +t; s +=t
    print(s); s=''

def perturbation(mode=0) : # stark effect on hydrogen degenerate subspaces
  # mode=1 latex-type  output, mode=0 debug screen output
  # for all pairs (n,m) need the matrix of all l with Rnl*Plm elements
  # rnl entries have entries [n,l, polynomial, sqnorm,norm]; ditto plm
  rnl,rmat= getradials(4,mode) # radial func Rnl,norms,their r-matrix elements
  plm,smat= getsphericals(4,mode) # Plm func,norms,cos(theta) matrix elements
  nr=len(rnl); ns=len(plm); rad=[0]*32;lrad=[0]*32;sph=[0]*32;lsph=[0]*32
  pair=[[]]*32; noisy=True; noisy=False 
  ##           perturbation matrix degenerate level n,m
  def setupmatrix(n,m,pair,noisy) : # pair[i]= [l-value, irad,iang]
    p=len(pair); sl=''; mat=[[]]*p; lxm=[[]]*p # terminal and latex matrix
    for i in range(p) : sl +=' '+str(pair[i][0])
    print('Matrix  n='+str(n)+'  m='+str(m)+' Liste l= '+sl)
    for k in range(p) : mat[k]= ['']*p; lxm[k]= ['']*p
    for i in range(p) :
      irad=pair[i][1]; iang=pair[i][2]; irnl=rnl[irad]; iplm=plm[iang]
      #print(' i='+str(i)+ str(irnl)+ str(iplm)) 
      irnorm=irnl[4]; ipnorm=iplm[4]
      for k in range(i+1) :
        krad=pair[k][1]; kang=pair[k][2]; krnl=rnl[krad]; kplm=plm[kang]
        krnorm=krnl[4]; kpnorm=kplm[4]; r=rmat[irad][krad];a=smat[iang][kang]
        if noisy and (a[0][0] != 0) and (mode==0) :
          print('  ik,irnorm,ipnorm,krnorm,kpnorm '+str(i)+str(k)+' '
            +str(irnorm)+str(ipnorm)+str(krnorm)+str(kpnorm))
          print(' mat-el r,a='+str(r)+str(a))
        # matrixelement is: (r/irnorm/krnorm)*(a/ipnorm/kpnorm)
        r= normalize(r,irnorm,krnorm); a=normalize(a,ipnorm,kpnorm)
        ra=mulroot(r,a); z= stroot(ra); r=stroot(r); a=stroot(a) 
        mat[i][k]=z; mat[k][i]=z
        z=lxroot(ra); lxm[i][k]=z; lxm[k][i]=z
        # print('  mat,rad,ang:  '+z+'  , '+r+','+a) # angular=0 on diagonals
    return mat,lxm
  ##
  for n in range(1,5) : # find all l that match radial(n,l) and angular(l,m) 
    h=0
    for i in range(nr) : 
      if rnl[i][0]==n : rad[h]= i;lrad[h]=rnl[i][1]; h+=1
    # print('n='+str(n)+' lrad='+str(lrad[:h]))
    for m in range(5) : # find all l from Plm and Rnl 
      g=0; ix=0 # pair[ix]=[lvalue, index to rnl, index to plm]
      for i in range(ns) : 
        if plm[i][1]==m : sph[g]=i; lsph[g]=plm[i][0]; g+=1
      # print('   m='+str(m)+' lsph='+str(lsph[:g])) 
      for j in range(h) :
        ll=lrad[j]
        for k in range(g) : 
          if (ll==lsph[k]) : pair[ix]=[ll, rad[j],sph[k]]; ix+=1
      if ix>1 : # matrix, not a single nondegen level
        mat,lxm= setupmatrix(n,m,pair[:ix],noisy); noisy=False
        if mode==0 : matdisplay(mat)
        if mode==1 : lm= lxmatrix(lxm); print(lm)
  if mode==1 :
    print(lxradials(rnl)) 
    print(lxsphericals(plm)) 
  
#perturbation(0) # debug
#print('*** Wikibook-Latex Output ***')
perturbation(1)

Die WKB-Methode

Die Wentzel-Kramers-Brillouin-Methode versucht, näherungsweise die stationäre Schrödingergleichung in einer Dimension zu lösen. Die Energie sollte genügend weit vom Grundzustand entfernt sein. In jeder Umgebung gibt es eine Basis von zwei Wellen, da die Differenzialgleichung zweiter Ordnung und reell ist. Es werden diese Fälle um einen Punkt x herum behandelt:

• E weit über dem Potenzial V(x), ergibt ein Paar oszillierender Funktionen
• E weit unter V(x), eine exponentiell fallende und eine steigende Funktion
• E=V(x), das heißt, x ist ein Umkehrpunkt für klassische Bahnen.

Unter der Annahme, dass V(x) genügend weit linear um x herum ist, leisten am Umkehrpunkt dann Paare aus Besselfunktionen die nötige Approximation. Danach werden Linearkombinationen aller Segmente so bestimmt, dass sie stetig-differenzierbar aufeinander treffen und global als Lösung durchgehen.

Mit $k(x)={\sqrt {{\frac {2m}{\hbar ^{2}}}|(E-V(x))|}}$ wird die Schrödinger-Gleichung

• $u''(x)+k^{2}(x)u(x)=0\quad$ in einem Intervall wo $E>V(x)$
• $u''(x)-k^{2}(x)u(x)=0\quad$ in einem Intervall wo $E<V(x)$

Die Näherung postuliert nun, dass die Amplitude von u wesentlich langsamer variiert als die Phase. Im Fall E>V um Punkt $x_{0}$ wird angesetzt:

$u(x)=a(x)\exp(\pm i\int _{x_{0}}^{x}dx'\;k(x'))$.

Einsetzen in die Gleichung liefert dies für a(x) :

$a''\pm 2ik(x)a'\pm ik'(x)a=0$

Wird hier $a''$ vernachlässigt, folgt $(a'/a)=-{\tfrac {1}{2}}(k'/k)$, was eine Lösung $a(x)=1/{\sqrt {k(x)}}$ hat. Damit wird die WKB-Näherung für E>V:

$u(x)\propto {\frac {1}{\sqrt {k(x)}}}\exp[\pm i\int k(x')dx']$

Mit $a=1/{\sqrt {k}}$ wird $a''={\frac {-k''}{2k^{3/2}}}+{\frac {3k'^{2}}{4k^{5/2}}}$. Die Bedingung $|a''|\ll |k'a|$ kann dann übersetzt werden in folgende Langsamkeits-Kriterien für k(x): $|k''/k'|\ll k$ und $|k'/k|\ll k$.

Im Gebiet E<V, worin die Funktion u per Tunneleffekt exponentiell eindringt, ergibt sich nach dem gleichen Muster die Approximation:

$u(x)\propto {\frac {1}{\sqrt {k(x)}}}\exp[\pm \int k(x')dx']$

Wie geht man vor, um die Umkehrpunkte V(w)=E zu überbrücken?

An einem Punkt w mit E=V(w) sei $(V(x)-E)\propto (x-w)$ mit V>E für x>w. Also $k(x)=k_{0}{\sqrt {x-w}}$.

Mit der Variablen $y=\int _{w}^{x}k(x')dx'={\frac {2k_{0}}{3}}(x-w)^{3/2}$ hat u(y) die Schrödingergleichung

${\frac {d^{2}u}{dy^{2}}}+{\frac {1}{3y}}{\frac {du}{dy}}-u=0$

Auf der anderen Seite x<w, E<V, ist $k(x)=k_{0}{\sqrt {w-x}}$.

Dort setze $y=\int _{x}^{w}k(x')dx'={\frac {2k_{0}}{3}}(w-x)^{3/2}$. Die Schrödingergleichung wird zu

${\frac {d^{2}u}{dy^{2}}}+{\frac {1}{3y}}{\frac {du}{dy}}+u=0$

Die Differenzialgleichungen für u(y) am Umkehrpunkt haben die Form

$u''+{\frac {a}{y}}u'+bu=0\quad$ wo $\quad a=1/3;\quad b=\pm 1$

Der Trick ist, die Gleichung für Besselfunktionen z(t) da herauszuschälen:

${\ddot {z}}+{\frac {z}{t}}+(1-{\frac {s^{2}}{t^{2}}})z=0;\quad z(t)=J_{\pm s}(t)$

Dazu hier die funktionierende Variablen-Substitution mit Parametern p,q:

$t=qy;\quad u(y)=({\frac {y}{q}})^{p}z(qx);\quad y={\frac {t}{q}};\quad z'=q{\dot {z}};z''=q^{2}{\ddot {z}}$

Einsetzen in die Gleichung für u und Auflösen nach z und t führt zu:

${\ddot {z}}+{\frac {2p+a}{t}}{\dot {z}}+[{\frac {ap+p(p-1)}{t^{2}}}+{\frac {b}{q^{2}}}]z=0$

Mit der geforderten Bessel-Gleichung ergeben sich q,p und s aus a und b:

$q^{2}=b;\quad 2p+a=1;\quad -s^{2}=p(p-1)+ap$

Die zwei Wurzeln für s ergeben zwei linear unabhängige Besselfunktionen. q kann imaginär sein, s ebenfalls? Die Familie von Besselfunktionen $J_{s}(t)$ ist analytisch in beiden Argumenten s und t, t kann imaginäre Werte haben.

$a={\frac {1}{3}};\;b=1\;\Rightarrow \;q=1;\;p={\frac {1}{3}}\;\;s=\pm {\frac {1}{3}}$

$a={\frac {1}{3}};\;b=-1\;\Rightarrow \;q=i;\;p={\frac {1}{3}}\;\;s=\pm {\frac {1}{3}}$

Auf der Tunnel-Seite E<V des Umkehrpunkts gibt es das Paar von Lösungen:

$u(y)=y^{1/3}I_{\pm 1/3}(y)\quad$ mit $I_{s}(t):=(i)^{-s}J_{s}(it);\;i^{-s}=\exp(-is\pi /2)$.

Es wird angenommen, dass u(y) bis in einen Bereich $y\gg 1$ hinein gut ist, so dass dort die Exponential-Lösung angeflanscht wird. Es gibt dort folgende Asymptotik:

$y^{1/3}I_{1/3}(y)\sim [\exp(y)(1+O(1/y))+\exp(-y-5i\pi /6)(1+O(1/y))]/{\sqrt {2\pi y^{1/3}}}$

$y^{1/3}I_{-1/3}(y)\sim [\exp(y)(1+O(1/y))+\exp(-y-i\pi /6)(1+O(1/y))]/{\sqrt {2\pi y^{1/3}}}$

Nun war $y\propto (x-w)^{3/2}$, daher $1/{\sqrt {y^{1/3}}}\propto 1/{\sqrt {k(x)}}$.

Die Asymptote bekommt glatten Anschluss an die WKB-Lösung, wenn sie exponentiell abfällt, daher: $u\propto y^{1/3}(I_{1/3}(y)-I_{-1/3}(y))$.

Auf der oszillierenden Seite E>V: $u(y)=y^{1/3}J_{\pm 1/3}(y)$. Das Grenzverhalten für kleine x auf beiden Seiten des Punktes w wird gebraucht, um die stetig differenzierbare Fortsetzung zu erreichen:

$y^{1/3}I_{1/3}(y)\to {\frac {y^{2/3}}{2^{1/3}\Gamma (4/3)}}={\frac {k_{0}^{2/3}}{2^{1/3}\Gamma (4/3)}}(x-w)$

$y^{1/3}I_{-1/3}(y)\to {\frac {2^{1/3}}{\Gamma (2/3)}}$

$y^{1/3}J_{1/3}(y)\to {\frac {y^{2/3}}{2^{1/3}\Gamma (4/3)}}={\frac {k_{0}^{2/3}}{2^{1/3}\Gamma (4/3)}}(w-x)$

$y^{1/3}J_{-1/3}(y)\to {\frac {2^{1/3}}{\Gamma (2/3)}}$

Der $I_{+}$-Teil hat glatten Anschluss an $-J_{+}$ und der $I_{-}$-Teil an $J_{-}$. Also ergibt sich als Fortsetzung von der Tunnelseite in (E>V)-Seite:

$u(y)\propto y^{1/3}(J_{1/3}(y)+J_{-1/3}(y))$

Auch hier hoffentlich bis in einen Bereich $y\gg 1$ mit der Asymptotik:

$y^{1/3}[J_{1/3}(y)+J_{-1/3}(y)]\to {\sqrt {\frac {2}{\pi y^{1/3}}}}[\cos(y-(\pi /4)-(\pi /6))+\cos(y-(\pi /4)+(\pi /6))]$

$\propto y^{-1/6}\cos(y-(\pi /4))\propto k(x)^{-1/2}\;\cos(\int _{x}^{w}k(x')dx'-\pi /4)$

Andererseits gilt auf offener Strecke (E>V) zwischen Umkehrpunkten nach WKB eine Linearkombination aus $u\propto k^{-1/2}\exp(\pm i\int k(x')dx')$. Diese kann nun fixiert werden bezüglich des nächsten Umkehrpunktes w>x:

$u\propto k(x)^{-1/2}\;\cos(\int _{x}^{w}k(x')dx'-\pi /4)$.

In einem Potentialtopf zwischen zwei Umkehrpunkten a<x<b muss symmetrisch die andere Randbedingung verwertet werden:

$u\propto k(x)^{-1/2}\;\cos(\int _{a}^{x}k(x')dx'-\pi /4)$, also muss gelten

$\cos(\int _{x}^{b}k(x')dx'-\pi /4)\propto \cos(\int _{a}^{x}k(x')dx'-\pi /4)$

Mit positivem Faktor ergäbe sich daraus

$\int _{x}^{b}k(x')dx'-\pi /4=\int _{a}^{x}k(x')dx'-\pi /4+n\pi$.

Aber als Funktion von x variieren hier beide Seiten im Gegensinn, was diese Wahl ausschließt. Vielmehr stimmt der Faktor -1 und man kann die Integrale sofort vereinen zu einer Forderung an die Funktion k(x):

$\int _{a}^{b}k(x')dx'=(n+1/2)\pi ,\quad$ mit nicht-negativem ganzen n.

Funktion k(x) ist die Rate, mit der die Phase der Welle sich ändert. Die summierte Phase muss, von Randeffekten abgesehen, ein Vielfaches von $\pi$ sein, damit u(x) als stehende Welle in eine Potenzialmulde passt.

Die WKB-Methode findet also quantisierte Niveaus. Sie produziert sogar für den harmonischen Oszillator exakt die richtigen Energiewerte, während für die l=0 Wellen im Coulomb-Potenzial eine verschobene Energie-Folge herauskommt, nämlich anstatt $1/n^{2}$ etwa $1/(n+3/4)^{2}$. Dies unterstreicht nur, dass die Methode nicht zu nahe am Grundzustand des Systems anzuwenden ist.

Tunneleffekt mit WKB

Zur Behandlung eines Potenzial-Hügels statt eines Tals wird noch die Anschluss- Bedingung an Umkehrpunkten gebraucht, die sich exponentiell ansteigend in den Tunnelbereich fortsetzt. Dann wird folgende Aufgabe angegangen: Eine freie Welle kommt von links, dringt in die Potenzial-Mauer V(x) ein und wird zum Teil rechts zu einer durchgehenden Welle, zum anderen Teil reflektiert. Die Linearkombinationen aus je zwei WKB-Wellen in den drei Bereichen sind mit den Übergangs-Formeln stetig zu verbinden. Gesucht ist die Amplitude für die Transmission, also auch die Tunnel-Wahrscheinlichkeit als Quadrat davon.

Mit Hilfe von asymptotischen Formeln für Besselfunktionen, analog zu oben, begründet man folgenden Regelsatz für die Verbindung von WKB-Wellen:

• V(x)<E für x<a, Umkehrpunkt a, V(x)>E für x>a :

${\frac {2}{\sqrt {k}}}\cos(\int _{x}^{a}k\,dx'-\pi /4)\longleftrightarrow {\frac {1}{\sqrt {k}}}\exp(-\int _{a}^{x}k\,dx')$

${\frac {1}{\sqrt {k}}}\sin(\int _{x}^{a}k\,dx'-\pi /4)\longleftrightarrow -{\frac {1}{\sqrt {k}}}\exp(\int _{a}^{x}k\,dx')$

• V(x)>E für x<b, Umkehrpunkt b, V(x)<E für x>b :

${\frac {1}{\sqrt {k}}}\exp(-\int _{x}^{b}k\,dx')\longleftrightarrow {\frac {2}{\sqrt {k}}}\cos(\int _{b}^{x}k\,dx'-\pi /4)$

$-{\frac {1}{\sqrt {k}}}\exp(\int _{x}^{b}k\,dx')\longleftrightarrow {\frac {1}{\sqrt {k}}}\sin(\int _{b}^{x}k\,dx'-\pi /4)$

Ansatz der Wellenfunktion für einen Potential-Berg V(x)>E im Intervall a<x<b, setze Variable x immer als Integral-Obergrenze:

$(x<a):\quad u_{1}={\frac {A}{\sqrt {k(x)}}}\exp(i\int _{a}^{x}k\,dx')+{\frac {B}{\sqrt {k(x)}}}\exp(-i\int _{a}^{x}k\,dx')$

$(a<x<b):\quad u_{2}={\frac {C}{\sqrt {k(x)}}}\exp(-\int _{a}^{x}k\,dx')+{\frac {D}{\sqrt {k(x)}}}\exp(\int _{a}^{x}k\,dx')$

$(b<x):\quad u_{3}={\frac {F}{\sqrt {k(x)}}}\exp(i\int _{b}^{x}k\,dx')+{\frac {G}{\sqrt {k(x)}}}\exp(-i\int _{b}^{x}k\,dx')$

Um die Anpassungs-Regeln einzusetzen, werden erst $\exp(\pm iz)$ in $\cos(z-\pi /4),\sin(z-\pi /4)$ umcodiert und bei Bedarf Integral-Grenzen umgeklappt:

$u_{1}=[A(1-i)+B(1+i)]{\frac {\cos(\int _{x}^{a}k\,dx')-\pi /4}{\sqrt {2k}}}+[A(-1-i)+B(-1+i)]{\frac {\sin(\int _{x}^{a}k\,dx')-\pi /4}{{\sqrt {(}}2k)}}$

$u_{3}=[F(-1+i)+G(-1-i)]{\frac {\cos(\int _{b}^{x}k\,dx')-\pi /4}{\sqrt {2k}}}+[F(1+i)+G(1-i)]{\frac {\sin(\int _{b}^{x}k\,dx')-\pi /4}{{\sqrt {(}}2k)}}$

$u_{2}=(C/I)\exp(\int _{x}^{b}k\,dx()+DI\exp(-\int _{x}^{b}k\,dx')$

mit $I=\int _{a}^{b}k\,dx'$. I misst die Opazität und ist typisch eine große Zahl bei dicker,hoher Barriere (V(x)-E). Die Verbindung geschieht dann nach dem Dreisatz: Soll $2\cos(z_{1}-\pi /4)$ sich mit $\exp(-z_{2})$ treffen, dann muss der erste Faktor in $u_{1}$ gleich 2C sein. Alles zusammen ergibt 4 lineare Gleichungen mit A,B,C,D,E,F. Aufgelöst nach A,B ergibt sich mit $u:=(2I+(2I)^{-1});\;v:=(2I-(2I)^{-1})$ :

$A={\tfrac {1}{2}}(uF+ivG)$

$B={\tfrac {1}{2}}((-i)vF+uG)$

Das erste Matrixelement (u/2) ergibt die Durchlass-Wahrscheinlichkeit der Welle, oder auch den Durchlass-Koeffizienten, wenn G zu Null gesetzt wird (keine einlaufende Welle von rechts).

$T=(|F|/|A|)^{2}=4/u^{2}\approx 1/I^{2}$ weil typischerweise $I\gg 1$.

Alternative Behandlung von Umkehrpunkten Kurz erwähnt sei ein Verformungs-Verfahren für Potenziale mit Umkehrpunkten. Sein Anwendungsbereich ist 'quasiklassisch' wie WKB. Das Prinzip: Sowohl der Definitions- wie der den Wertebereich der zeitunabhängigen Schrödinger-Gleichung sollen nichtlinear so transformiert werden, dass in den neuen Koordinaten eine äquivalent Schrödinger-Gleichung mit einem leichter lösbaren Potenzial-Term auftaucht, $\psi (x)\mapsto \phi (s)$.

Aus $\psi ''(x)+v^{2}(x)\psi (x)=0$ (v reell oder imaginär in Intervallen) wird $\phi ''+p^{2}(s)\phi (s)=0$. Das gewählte $p^{2}(s)$ soll dieselbe Folge von Umkehrpunkten und postiven/negativen Intervallen haben wie $v^{2}(s)$. Die Transformation habe die Form $s=s(x)$; es sei $T(x)\phi (s(x))=\psi (x)$ mit einer reellen Amplitudenmodulation T(x).

$\psi ''=T''\psi +(2T's'+Ts'')(d\psi /ds)+Ts'^{2}(d^{2}\psi /ds^{2})$

Der mittlere Term wird wegtransformiert, wenn T(x) als Potenz von $s'(x)$ gewählt wird: $T(x)=s'(x)^{-1/2}$. Es bleibt diese Gleichung mit $p^{2}(s)$:

$\psi ''=[T''-Ts'^{2}p^{2}(s(x))]\phi (s(x))$

Andereseits folgt aus der Original-Schrödinger-Gleichung für das Paar {T,s}:

$T''(x)+[v^{2}(x)-s'^{2}(x)p^{2}(s(x))]T(x)=0$.

Jetzt wird die quasiklassische Näherung gemacht, $T''$ zu vernachlässigen. Damit wird $s'(x)p(s(x))=v(x)$, woraus sich die Funktion s(x) als Obergrenze des folgenden Integrals indirekt ergibt. $(x_{0},s_{0})$ ist dabei ein Umkehrpunkt:

$\int _{s_{0}}^{s(x)}p(t)dt=\int _{x_{0}}^{x}v(z)dz$.

Die Approximations-Lösung der Schrödinger-Gleichung folgt dem Rezept:

• wähle ein gutmütiges Potenzial $p^{2}(s)$ kompatibel mit $v^{2}(x)$
• berechne aus p() und v() die Funktionen $s(x),s'(x)$, eventuell numerisch
• prüfe mit $T=1/{\sqrt {s'}}$, ob $T''(x)$ klein genug ist
• finde Lösungen der Schrödinger-Gleichung, $\phi (s)$ zum Potenzial $p^{2}(s)$
• $\psi (x)=\phi (s(x))/{\sqrt {s'(x)}}$

Zeitabhängige Störungsrechnung

Das Ziel einer zeitabhängigen Rechnung ist es, mehr herauszufinden als die stationären Zustände. Nämlich die Dynamik der Prozesse. Geschwindigkeiten der Übergänge zwischen Zuständen, Regeln für erlaubte und verbotene Prozesse, Verteilungen von Häufigkeiten, Energien und Winkeln von Spaltprodukten bei Explosionen und Kollisionen, usw.

Anhang Python-Skript

#h2plus.py  basicmodel of hydrogen ion bond, variation method

from math import exp,sqrt
def div(a,b) : return int(a/b)

# h2plus energy level approx
def erglvl2(r,q,a) : # r=proton distance, q=wave radius, a=bohr radius
  # all factors pi are suppresed, nq= S11=S22, s=S12=S21
  # to return energies with respect to hydrogen, factor 2a
  nq= q*q*q; cq= q*q; x=q/r; z=2*q/r; y=exp(-r/q)
  s= r*r*r*y*( x/3.0 + x*x + x*x*x)
  u= r*r * y * (x+x*x)
  q2=q/2.0; z=2.0*q2/r; y2=exp(-r/q2); t=r*r*((z*z+z*z*z)*(-y2) + z*z*z)
  f=(1.0/r-a/(2.0*q*q))
  h11= f*nq + (a/q-1)*cq - t
  h12= f*s + (a/q-2.0)*u
  # solve (h11-nq u)^2=(h12-s u)^2: 
  uplus= 2*a*(h11-h12)/(nq-s); uminu= 2*a*(h11+h12)/(nq+s)
  return uplus,uminu, [nq,s,t,u, cq]

def findmini(x,y) : # parabola interpolation y=a+bx+cx^2
  x1,x2,x3=tuple(x); y1,y2,y3=tuple(y) # find y=a+bx+cxx
  z1= (y2-y1)*(x3*x3-x2*x2); z2= (y3-y2)*(x2*x2-x1*x1)
  w1= (x2-x1)*(x3*x3-x2*x2); w2= (x3-x2)*(x2*x2-x1*x1) # b(w1-w2)=z1-z2
  b=(z1-z2)/(w1-w2); c= ((y2-y1)-b*(x2-x1))/(x2*x2-x1*x1)
  a= y1 - b*x1 - c*x1*x1; xmin=-b/(2*c); ymin=a-b*b/(4.0*c)
  return xmin,ymin

def scanerglevel(r,a) :
  noisy=False; fmax=1.2; df=0.04; refine=True; tbl=[]
  for k in range(101) : # r from 1 to 4 step 0.03
    emin=1e9; imin=-1; r= 1.0+ 0.03*k; y=0
    u0,v0,bug= erglvl2(r,a,a)
    if noisy : print('erglevel scan r='+str(r)); 
    for i in range(20) :
      f=fmax- df*i; q=a*f 
      u,v,bug= erglvl2(r,q,a); w=min(u,v)
      if noisy : print('i='+str(i)+' e='+str(w))
      emin,imin= (w,i) if (w<emin) else (emin,imin)
    if refine :
      x=fmax-df*imin; q1=a*x; q0=a*(x+df); q2=a*(x-df)
      u,v,bug= erglvl2(r,q0,a); y0=min(u,v)
      u,v,bug= erglvl2(r,q1,a); y1=min(u,v)
      u,v,bug= erglvl2(r,q2,a); y2=min(u,v)
      x,y= findmini([q0,q1,q2],[y0,y1,y2]) 
      u1,v1,bug= erglvl2(r,x,a); tbl +=[[r,u0,v0,u1,v1]]
    print('erg r,i,emin ='+str([r,imin,emin,y])) 
  return tbl

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 myplot(tbl) :
  n=len(tbl); xmin=1e9;xmax=-xmin; ymin=xmin;ymax=-ymin; k=len(tbl[0])
  print('n='+str(n)+' k='+str(k))
  xys=[[]]*(k-1)
  for j in range(k-1) : xys[j]= [0.0]*(2*n)
  for i in range(n) :
    xx= tbl[i][0]; xmin=min(xmin,xx); xmax=max(xmax,xx)
    for j in range(1,k) :
      z=tbl[i][j]; ymin=min(ymin,z);ymax=max(ymax,z)
      xys[j-1][2*i]=xx; xys[j-1][2*i+1]= z
  print('[xmin,xmax,ymin,ymax]='+str([xmin,xmax,ymin,ymax]))
  m=len(xys); print('len xys='+str(m))
  for i in range(m) : print(' i='+str(i)+' len='+str(len(xys[i])))
  xmin=1;xmax=4;ymin=-1.5;ymax=0;yref=-1
  xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin,
        xmin,yref, xmax,yref, xmax,ymin,3,ymin,3,ymax,2,ymax,2,ymin] # frame
  xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale start user units, min in pixels
  ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#000','#000','#00a','#077','#0a0','#770','#a00']
  fn='h2plus';  xys= [xyf] + xys
  text='H2+ Ion Energy levels. Base line = Hydrogen atom ground state.';
  sd=svgdump(); buf=sd.plotbegin(800,600,1)
  buf += sd.labl(text,50,500); i=0
  text='Verticals: 1,2,3,4 Bohr radius';
  buf += sd.labl(text,50,530); i=0
  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()

tbl= scanerglevel(2,1); myplot(tbl)
print('File h2plus.svg written.')

Hoch zum Seitenanfang