Quantenmechanik/ Spinoren

Aus Wikibooks

Spin[Bearbeiten]

Einleitung[Bearbeiten]

Die Schrödingergleichung entsteht dadurch, dass in der Energie-Impuls- Beziehung des Punktteilchens diese beiden Größen zu Ableitungs- Operatoren werden. Nämlich genau die Operatoren, die aus jeder ebenen Welle als Eigenwerte die de-Broglie-Korrespondenz hervorlocken: .

Die Quantenmechanik hat aber mehr auf Lager als diese Wellenversion der klassischen Bahngleichung. Erstens hat das Elektron ein magnetisches Moment und eine Quantisierung des Drehimpulses in zwei Stufen. Schrödingerwellen beherrschen nur ungerade Quantisierung von Bahndrehimpulsen. Zweitens taucht bei vielen gleichartigen Teilchen die Zusatzregel auf, dass die gemeinsame Wellenfunktion entweder total symmetrisch (Bosonen) oder total antisymmetrisch (Fermionen) sein muss.

Beide neue Fakten gehören zur relativistischen Physik. Aus der Diracgleichung entspringt zwanglos der Spin und das korrekte magnetische Moment des Elektrons. Die Vielteilchentheorie wurde aufgebohrt zur Quantenfeldtheorie, die die Erzeugung und Vernichtung der Teilchen sowie die Antiteilchen beschreibt. Und eine fundamentale Spin-Statistik- Regel herleitet: Teilchen mit halbzahligem Spin sind Fermionen, solche mit ganzzahligem Spin sind Bosonen.

Der Spin des Elektrons ist ein neuer Freiheitsgrad über die Schrödingergleichung hinaus. Wie gesehen, gehört ein Spektrum der Drehimpuls-Algebra mit 2 Werten zu der kleinstmöglichen Darstellung, nämlich Vektoren in mit Operatoren der Form Pauli-Matrix mal . Soweit mit Streuexperimenten belegbar, ist das Elektron punktförmig. Der Spin-Freiheitsgrad hat keine messbare räumliche Ausdehnung. Er wird dadurch modelliert, dass an jedem Punkt x die Wellenfunktion einen Wert im Darstellungsraum hat, an Stelle eines Skalarwertes. Die Werte heißen Spinoren, weil sie wie Vektoren und Tensoren bei linearen Transformationen der Ortskoordinaten ebenfalls richtig umgebogen werden müssen.

Der Raum der (Pauli-)Spinoren hat eine Basis aus 2 Zuständen und verkörpert so die kleinstmöglichen Quantensysteme, oft Qbits genannt. Alle hermiteschen Operatoren auf Spinoren sind reelle Linearkombinationen aus der Einheitsmatrix und den drei Pauli-Matrizen.

Jedem Vektor ist mit bijektiv eine spurfrei-hermitesche Matrix zugeordnet. Wenn die Norm 1 hat, ist physikalisch der Operator, der den Spin in Richtung misst.


Rechnen mit Pauli-Matrizen[Bearbeiten]

.
  • zyklisch
  • Die Pauli-Matrizen sind gleichzeitig hermitesch und unitär.
  • Sie antikommutieren miteinander und haben die Eigenwerte {1,-1}.
  • Die Matrizen sind eine Lie-Algebra-Basis
ohne den Faktor i (Mathe-Konvention): .
  • Die haben Quadrate -1, sind antiunitär, antikommutieren.
  • Zusammen mit sind sie eine Basis der Quaternionen-Algebra.

Ist eine Matrix M diagonalisierbar mit einer Ähnlichkeitstransformation, lässt sich so manche Funktion f(x) mit reellem bzw. komplexem Definitions- und Wertebereich zu einer matrixwertigen Matrixfunktion f(M) erweitern. Man wende sie einfach in der Diagonalform auf alle Eigenwerte an. Zum Beispiel Exponential- Funktion und Quadratwurzel haben dann gewohnte Eigenschaften.

Beispiel, Exponential der Matrizen mit als Reihenentwicklung berechnet. Es funktioniert wie mit :

Die Familien sind unitäre Ein-Parameter-Gruppen.

Magnetisches Feld[Bearbeiten]

Spinor- und Vektor-Wellenfunktionen[Bearbeiten]

Die Schrödingerwelle ist ein komplexer Skalar. Ein Objekt kann zusätzlich Eigenschaften haben, die mit räumlicher Orientierung zusammenhängen, doch näherungsweise keine Ausdehnung besitzen. Dann zieht man eine Wellenfunktion heran, deren Werte Spinoren, Vektoren oder Tensoren sind. Die Operatoren / Matrizen auf diesen Wertebereichen bescheren dem Objekt neue Quantenzahlen und neue Terme im Hamilton-Operator, womit sie an Nachbarteilchen und an externe Kraftfelder ankoppeln.

Angenommen, dass das System unter den Drehungen des Koordinatensystems invariant ist. Der Hamilton-Operator, ein Skalar also, sollte nach der Drehung gleich aussehen, wenn diese Symmetrie gilt. Ein Raumpunkt x wird zu verändert. Weil die nicht-skalaren Wellenfunktionen orientierte Werte haben, müssen diese auch verformt werden, also wird gemäß umgeschrieben. ist eine Transformation, im linearen Fall eine Matrix, die auf dem Wertebereich der Wellen operiert.

Die Paarung von mit soll Gruppeneigenschaften beibehalten. Zu , zwei Transformationen hintereinander, soll gehören.

Gibt es eine Abbildung , so ist sie eine Darstellung der Drehgruppe, ein Gruppenhomomorphismus. Es kann auch umgekehrt ein Homomorphismus bestehen; die Drehgruppe ist eine Darstellung einer vielleicht größeren Gruppe, die auf dem Wellenwert-Raum operiert.

Genau dies tritt auf bei Spinoren und allen Darstellungen mit halbzahligem Spin der Drehimpuls-Algebra. Die Gruppe SU(2) auf dem Spinor-Raum, alle Speziellen (Determinante=1) Unitären komplexen 2-mal-2-Matrizen, überlagert die Drehgruppe SO(3) zwei mal. Nach einer Raumdrehung hat ein Spinor eine halbe Runde hinter sich, ist aber proportional zum Startwert mit Vorzeichen -1. Das ist neutral für die Quantenmechanik, denn diese akzeptiert Strahldarstellungen auf dem Hilbertraum für ihre Symmetrien. Warum? Wenn ein Zustandsvektor einen konstanten Faktor vom Betrag 1 bekommt, bleiben alle physikalisch messbaren Dinge gleich. Gegeben eine Observable, also ein hermitescher Operator A, ein Eigenwert a und sein Eigenvektor . Dann geht die Wellenfunktion nur als Prognose-Wahrscheinlichkeit ein, in der Form , als Betragsquadrat.


Motivation des Spinor-Raums[Bearbeiten]

In einem Gedankenexperiment schicken wir Silberatome durch drei Stern-Gerlach-Maschinen, die am magnetischen Moment in den drei Richtungen x,y,z angreifen. Irgendwie biegen wir den Strahl so, dass er senkrecht genug zu den Magnetfeldern ist. Solche Felder haben einen longitudinalen Gradienten und können damit Spin-Eigenzustände räumlich trennen. In allen Reihenfolgen xy, yz, xz ... teilt der erste Stern-Gerlach den Strahl in zwei Teile, der zweite jeden davon weiter in zwei gleiche Teile. Nur wenn zwei Maschinen direkt hintereinander die gleiche Richtung haben, teilt die zweite nichts weiter auf.

Gesucht wird eine minimale Darstellung von Operatoren , die jeweils das Eigenwert-Paar {+1,-1} haben und die hintereinander geschaltet wie die Stern-Gerlach-Kombinationen funktionieren, also die modellhaft Projektions-Messungen auf Zustandsvektoren ausführen.

Der Zustandsraum hat mindestens die Dimension 2. wird zuerst gewählt, als die

Matrix ,
mit Eigenvektoren u:=(1,0), v:=(0,1).

soll u oder v in zwei gleiche Teile zerlegen,

Eigenvektoren:
Matrix: .

Nun braucht nochmal das gleiche Verhalten bezüglich sowohl als auch . Muss dafür die Dimension der Vektoren höher gehen? Nein, aber die Größe des Zahlenkörpers! Mit komplexen Zahlen gelingt es.

Eigenvektoren:
Matrix: .

Daher ist die kleinstmögliche Darstellung ein komplexer und die observablen Operatoren sind proportional zu Pauli-Matrizen. Die komplexen Zahlen erweisen sich praktisch als eine Notwendigkeit zur Beschreibung der Quantenwelt.

Mit der Einheitsmatrix I sind die Pauli-Matrizen die Standardbasis für den Raum aller hermiteschen Matrizen, also speziell aller quantenmechanischen Operatoren im Zwei-Niveau-System.

Die Pauli-Matrizen sind hermitesch, spurfrei und befolgen die Regeln:

Die spurfreien hermiteschen Matrizen sind ein dreidimensionaler reeller Vektorraum

Die allgemeinen hermiteschen Matrizen sind ein Raum der Dimension 4:

Ist ein Einheitsvektor

dann hat seine hermitesche Matrix die Eigenwerte {+1,-1}.

Die Exponentialfunktion dieser speziellen Matrizen folgt aus der Reihe, in der alle geraden Potenzen wegen die Einheitsmatrix ergeben.

Die Matrizen U haben die Determinante 1, denn mit  :

Außerdem sind sie unitär, denn mit allen erwähnten Rechenregeln

Die Menge der Matrizen U ist die Lie-Gruppe SU(2) der unitären Matrizen mit Determinante 1. Die Menge der Matrizen

ist die zugehörige Lie-Algebra su(2).

Adjungierte Darstellung von SU(2)[Bearbeiten]

Bei unitären Ähnlichkeitstransformationen bleiben die Eigenschaften 'hermitesch' und 'spurfrei' erhalten. Mit der Parametrierung als Dreiervektoren gilt daher

Die Gruppe SU(2) der unitären 2x2-Matrizen mit Determinante 1 ist eine der einfachen Lie-Gruppen, ist kein Produkt von nichttrivialen kontinuierlichen Untergruppen. Elemente ihrer Lie-Algebra, also die infinitesimalen Transformationen, die Ableitungen am Einselement I von Ein-Parameter-Untergruppen, sind genau die spurfreien anti-hermiteschen Matrizen iG, U(t)=exp(iGt). Die Physik hat den immerwährenden Reflex, einen Faktor i abzuspalten, so dass nur hermitesche Matrizen G vorkommen.

Matrizen U ∈ SU(2) induzieren eine Abbildung , wie folgt auszurechnen.

Nach einigen Manipulationen bleibt eine Kombination der Pauli-Matrizen

Mit den unverwüstlichen Rechenregeln:

Eine dreidimensionale Drehung um die Achse, die durch den Einheitsvektor definiert wird, mit dem Winkel φ, wirft einen Vektor auf eine Linearkombination aus und sieht so aus:

Ist auf der posititven z-Achse und auf der positiven x-Achse, dann dreht ein positiver Winkel gegen den Uhrzeigersinn zur positiven y-Achse hin.

Damit definiert das berechnete B eine Drehung in drei Dimensionen.

Die Abbildung D, mit Parameter φ und Einheitsvektor a:

ist ein Gruppen-Homomorphismus von SU(2) auf SO(3). Er ist zwei-zu-eins, denn um π verschobene Winkel φ ergeben dieselbe Drehung in Raum. SU(2) wird als die Überlagerungsgruppe von SO(3) bezeichnet. U(a,φ) ∈SU(2) transformiert Spinor-Zustandsvektoren. Seine zugeordnete adungierte Darstellung R(a,-2φ) transformiert 3D-Vektor-Operatoren (=hermitesche 2x2-Matrizen). Eine physikalische Vektor-Observable muss zwei volle Drehungen machen, wenn ein Spinor nur eine Runde macht. Das ist kein Drama, denn in beobachtbaren Größen kommen die Spinoren stets paarweise vor, genauer gesagt hermitesch-sesquilinear in Matrixelementen von Operatoren.

Die halben 'Drehungen' in SU(2), abgebildet auf eine ganze Drehung gleich Identität im Ortsraum, sind die Matrizen (-I) — ein einfacher Vorzeichenwechsel. Die zwei isomorphen Gruppen sind SO(3) und die Quotientengruppe SU(2)/{±I}. Die recht triviale Untergruppe {±I} ist das Zentrum von SU(2), die Menge der Elemente, die mit allen anderen der Gruppe vertauschen.

Einige Abhandlungen definieren sowas wie mehrwertige Funktionen, die wahrscheinlich Unbehagen in der Welt der strengen Mathematik auslösen. Hier wäre eine solche die Relation, mit Einheitsvektoren a und Winkeln β:

Ein anderer Begriff läuft auch um, die Strahldarstellung. Das ist eine Zuordnung von einer Symmetriegruppe zu Paaren aus Elementen einer Matrixgruppe und Phasenfaktoren. Derart trägt man dem Umstand Rechnung, dass die Quantenzustände Strahlen von Vektoren sind, also eindimensionale Teilräume eines Hilbertraums.

Jedenfalls soll im Rest des Kapitels mit D meistens die folgende Abbildung bezeichnet werden:

Diagonalisierung hermitescher spurfreier Matrizen[Bearbeiten]

Eine Matrix

kann mit einer unitären Ähnlichkeitstransformation auf die diagonale Matrix abgebildet werden und hat daher die zwei Eigenwerte ±1. Es genügt, eine Rotation R(a,-2φ) zu finden, die den Vektor auf den Einheitsvektor der z-Achse wirft. Die Achse der Drehung steht senkrecht auf der Ebene, in der sich der Start- und Zielvektor befinden. Man kann direkt die Eigenvektoren von ausrechnen:

In Polarkoordinaten

Eigenwerte

Mit lösen wir die reellen Gleichungen

(u,w) sind typisch (cos,sin) des halben Winkels, denn mit Abkürzungen:

Zur Verschönerung wird noch mit der halben Phase von -φ gedreht und wir erhalten ein Paar von normierten Eigenspinoren  :

Aus den komplex konjugierten Zeilen dieser Eigenvektoren ergibt sich eine unitäre Matrix D, die als Ähnlichkeitstransformation verwandelt.

Folgende Projektions-Matrizen im Spinorraum haben als Bildmengen die zwei Eigenspinor-Strahlen der Matrix M:

Bewegungsgleichung der Systeme von Dimension 2[Bearbeiten]

Die Spinoren sind ein Prototyp für alle Quantensysteme in zwei Dimensionen, also praktisch alle Systeme mit bloß zwei Energie-Niveaus. Dazu gehören auch die Qubits, Bauelemente der Quantencomputer. Die Schrödinger-Gleichung i(d/dt)ψ=Hψ kann allgemein für alle hermiteschen Matrizen gelöst werden. Sowohl für konstantes H wie für periodische Störungen. Es ist immer möglich, das Niveau des Potenzials im Hintergrund so zu verschieben, dass H spurfrei wird. Anders gesagt, der Hamilton-Operator H ist eine reelle Linearkombination der Pauli-Matrizen.

Der Referenz-Hamiltonoperator koppelt den Spin des Elektrons an ein Magnetfeld und hat die Gestalt

Dieses H fällt nicht vom Himmel, sondern purzelt als Näherung aus der relativistischen Dirac-Gleichung heraus. Obwohl die spinorwertigen Wellenfunktionen, die jene Gleichung lösen, streng genommen nicht als eine Verallgemeinerung der Schrödingerschen Zustände durchgehen. Sieht man sich an, wie ein Vektorpotenzial in die Schrödinger-Gleichung einkoppelt, dann erscheint als Multiplikator des Magnetfelds insgesamt die Linearkombination aus dem Operator des Bahndrehimpulses und dem Spin-Operator. Mit dem gyromagnetischen Faktor g=2, wie Dirac herausfand. Die "wahre" Theorie der Elektronen und Photonen, die Quantenelektrodynamik, schaffte es später, die kleine Abweichung (g-2) mit hoher Genauigkeit zu berechnen.

Begibt man sich ins Heisenberg-Bild, dann ist der Vektor-Operator des Spins S(t) eine Funktion der Zeit und wird dargestellt nicht durch die konstanten Pauli-Matrizen, sondern durch ein Triplett von zeitveränderlichen Matrizen σ(t), die aber zu jedem Zeitpunkt die Vertauschungsrelationen erfüllen:

Die Heisenberg-Gleichung gibt an, wie dieser Spin sich entwickelt:

Als Dreier-Vektor-Gleichung geschrieben:

In der klassischen Mechanik kommt diese Gleichung bei Kreiseln vor, für gewöhnliche Vektoren. Ihre Lösung ist eine Präzession um die Richtung des Feldes

In Worten: die Komponente von σ parallel zur Feldrichtung bleibt konstant, während die zwei Komponenten senkrecht dazu eine Kreisbewegung mit der Winkelgeschwindigkeit ω vollführen. σ ist ein Vektor-Operator, allgemeiner als bloß ein Vektor. Die Adjektive parallel, senkrecht beziehen sich hier auf die Vektoren, die ein beliebig ausgewähltes Matrixelement des Operators abwirft: ⟨ψ|σ(t)|φ⟩.

Das vektorielle Produkt lässt sich als der 'infinitesimale' Erzeuger von Drehungen auffassen, denn mit der oben eingeführten Notation

Den Erwartungswert des Spin-Operators auf einem normierten, im Heisenberg-Bild konstanten Zustands-Spinor ψ bezeichnen wir als die Polarisation, p(t)=⟨ψ|σ(t)|ψ⟩.

Die Polarisation ist ein reellwertiger zeitabhängiger Vektor und er bewegt sich mit Präzession:

Es gilt |p(t)|=1, wie gleich gezeigt wird.

Was bedeutet es, wenn von einem Vektor-Operator die Rede ist? Wenn die Spinoren mit einer speziellen unitären Transformation U bedacht werden, dann sollen sich die Erwartungswerte der Komponenten eines Vektor-Operators mit einer zugehörigen Drehmatrix transformieren.

Der Operator der Polarisation ist so ein Vektor-Operator. Mit der besprochenen Beziehung der SU(2) zur Drehgruppe sieht es so aus:

Nach obiger Konvention ist hier . Setzen wir für b die drei Einheitsvektoren der Achsen ein.

Wenn bedacht wird, dass der Homomorphismus D mit dem anderen Vorzeichen definiert war (D:) und dass bei Dreh-Matrizen Transponierte gleich Inverse sind, folgt genau, dass die Polarisation die Anforderung an einen Vektor-Operator erfüllt.

Wir unterscheiden gern zwei Typen von Transformationen.

• Ein passiver Koordinatenwechsel im Spinor-Raum, bei dem alle Operatoren und Spinoren gleichzeitig so verbogen werden, dass die Erwartungswerte gleich bleiben:

• Eine Drehbewegung im Ortsraum, die von vektorwertigen Observablen und ihren Erwartungswerten mitgemacht wird. Um das zu erreichen, erhält der Spinor-Raum eine unitäre Strahl-Darstellung der Drehgruppe, oder besser umgekehrt der Ortsraum hat eine orthogonale Darstellung der Gruppe SU(2):

Passiver Art ist der Wechsel zwischen Heisenberg- und Schrödinger-Bild. Erwartungswerte haben die gleiche Zeitentwicklung im Heisenberg-Bild wie im Schrödinger-Bild, zu dem wir jetzt zurück kehren. Daher

Die Polarisation ist ein Einheitsvektor, wenn ψ einen normierten Spinor und σ das Triplett der Pauli-Matrizen bezeichnet.

Beweis. Eine hermitesche 2x2-Matrix hat die Spurformel

impliziert also die Zerlegung der Zustandsmatrix

Solche Zustands-Operatoren ganz allgemein sind idempotent, denn

Hier als Matrizengleichung

Aus dem Koeffizientenvergleich folgt (p·p)=1.

Finden wir stationäre Lösungen von .

Die Eigenwerte sind α=±ω. Hat der Vektor ω die Polarkoordinaten (r,θ,φ), dann sind zwei Eigenvektoren von , genau wie oben zur Diagonalisierung ausgerechnet, die folgenden:

Ist die ausgezeichnete Richtung die z-Achse, φ=0, θ=0, dann fallen natürlich die Eigenspinoren von ab:

Haben wir zum Zeitpunkt t=0 einen Eigenzustand zum Spin-Operator in der Ausrichtung (θ,φ) bezüglich der z-Achse, dann:

Die Zustandsmatrix hat die Form

mit dem Einheitsvektor der Polarisation

Selbstverständlich kommt im Schrödinger-Bild dieselbe Bewegung der Polarisation heraus wie im Heisenberg-Bild, die Präzession um die Feldrichtung nämlich.

Die Zustandmatrix ist ein Spezialfall der allgemeinen Konstruktion zur Beschreibung inkohärent gemischter, also nicht "phasenreiner" Quantensysteme, die under dem Namen Zustandsoperator läuft. Das Tupel oder die Folge ist eine Wahrscheinlichkeitsverteilung. Mitunter wird ρ in alten Büchern als Statistischer Operator und in noch älteren als Dichte-Matrix oder Dichte-Operator gehandelt.

to-do: Thema bekommt eine Unterseite, im Zusammenhang mit der Dekohärenz.

Die Polarisation, Erwartungswert des Operators σ, charakterisiert eindeutig einen normierten Spinor-Zustand ψ. Denn die Menge der Polarisationen ist eine zweidimensionale Mannigfaltigkeit, die Fläche der Einheitskugel. Sie wird als Bloch-Kugel bezeichnet. Die normierten -Spinoren haben drei freie Parameter. Aber da ihre Phase in Erwartungswerten herausfällt, sind die Zustände eine zweidimensionale Menge von Äquivalenzklassen. Eine Zustandsmatrix Z definiert den Spinzustand genauso vollständig und hat eine umkehrbare Zuordnung zum Polarisationsvektor,

Nur das zweidimensionale Quantensystem hat die Besonderheit, dass man die Dynamik der Zustände anschaulich durch die Bewegung eines Punktes auf der Kugel darstellen kann.

Spin im periodisch zeitveränderlichen Feld[Bearbeiten]

Ein Modell, das gut durchgerechnet werden kann, überlagert ein konstantes Magnetfeld in z-Richtung mit einem zweiten Feld, dessen Richtung mit konstanter Kreisfrequenz eine Präzession um die z-Achse vollführt.

Die Frequenz ω repräsentiert den hauptsächlichen Energie-Abstand von zwei Niveaus und ist relativ groß im Vergleich zu α welches eine kleine 'Störung' einbringt. Allerdings kommt in die Nähe von ω die Frequenz γ womit die Störung-sprich-Anregung am System rüttelt. Werden beide Hochfrequenzen gleich, kann sich eine Resonanz abspielen.

Die Tricks bestehen darin, zu einer Gleichung von zeitveränderlichen Spinoren zu wechseln, für die das statische wegfällt, und dann noch in ein rotierendes xyz-Koordinatensystem zu gehen (Spinoren und Vektoren drehen verzahnt). Am Ende wird das Problem so auf den vorher behandelten Fall des konstanten Feldes abgebildet.

Zuerst die Umrechnung auf eine Art von Wechselwirkungs-Bild. Transformiert wird mit dem unitären . Die Operatoren U(t) und der statische vertauschen.

Die Komponenten des transformierten Operators H' können mit

und sonstigen Regeln der Pauli-Matrizen eingedampft werden.

Die Feldrichtung von H' rotiert auf demselben Kegel wie der Operator der Anregung, aber H' dreht sich mit δ viel langsamer.

Die Transformation induziert natürlich eine Transformation aller Erwartungswerte der Vektor-Operatoren, besonders der Polarisation. Sie machen eine Drehung mit konstanter Geschwindigkeit ω um die z-Achse; das folgt aus dem Homomorphismus D von SU(2) auf SO(3) .

Eine weitere Transformation kann nun einen Operator H" herbei schaffen, der nicht mehr explizit von der Zeit abhängt. Man nehme im Spinor-Raum die SU(2)-Transformation U(t), die als Drehung im Ortsraum dafür sorgt, dass

Alle physikalischen Vektor-Operatoren machen nach dieser Transformation eine Drehung D(U(t)) um die z-Achse mit der Kreisfrequenz δ. Wird also ein vektorieller Erwartungswert wie die Polarisation mit H" ausgerechnet, bekommt man die gewünschte Zeitentwicklung unter dem vorherigen Operator H', indem man die Rotation D(U(t)) nachschaltet.

Für den Operator H" ist bekannt, wie die Bewegung eines Vektors der Polarisation aussieht. Er kreist mit der Frequenz α um den Richtungsvektor des Feldes, das hier in der xz-Ebene liegt und im häufigsten Spezialfall cos(β)=0 auf der x-Achse. Dieser Richtungsvektor macht aber seinerseits eine Präzession um die z-Achse, weil die eben eingeführte Operator-Transformation auf allen Dreiervektoren eine Rotation um die z-Achse mit der Kreisfrequenz δ induziert.

Kombiniert mit der Auswirkung der ersten Transformation im Raum der Vektor-Operatoren, rotiert der Axialvektor um die z-Achse mit der Hochfrequenz ω+δ=γ und die Polarisation um diesen herum mit der Niederfrequenz α. Die hohe Frequenz kommt also herein, wenn man vom Wechselwirkungs-Bild zum Schrödinger-Bild zurück geht. Sie hat keine messbaren Effekte und der Phasenfaktor, den sie produziert, wird weiter unten im Grafik-Skript neutralisiert.

Ergebnis. Ein reiner Spinor-Zustand entspricht eindeutig einem Polarisationsvektor auf der Einheitskugel. Jede Polarisation des Spin-(1/2)-Systems im konstanten plus drehenden Magnetfeld vollführt eine doppelte Präzession. Die Richtung des konstanten Feldes wird als die z-Achse gewählt; die zwei Energie-Eigenwerte des Spins haben den Abstand ℏω. Ein Hilfsvektor, sozusagen eine sekundäre Achse, kreist mit einer Frequenz δ um die z-Achse auf einem Kegel mit dem Winkel β. Seine Kreisfrequenz ist δ=(γ—ω), die Differenz von Anregungs-Frequenz und ungestörter Eigenfrequenz. Ein Erwartungswert der Polarisation rotiert um die kreisende Achse herum mit der vergleichsweise kleinen Frequenz α, die der Energie zum Absolutbetrag des Störfeldes entspricht. Wird δ auf Null getrimmt (Resonanz) und steht das Drehfeld senkrecht zum statischen Feld, dann kann eine Polarisation in z-Richtung vollständig umgeklappt werden. Die Geschwindigkeit des Prozesses wird durch die Feldstärke ~α kontrolliert und dieser Vorgang kann eventuell als niederfrequente Schwingung gemessen werden.

Bewegung der Polarisationsvektoren[Bearbeiten]

Statt der Bewegungsgleichung eines Spinors kann gleichwertig die seiner Zustandsmatrix und/oder seines Polarisationsvektors zum Einsatz kommen.

Die Matrizen Z,H sind hermitesch, deshalb gibt es drei Gleichungen.

Umrechnung zwischen Polarisation und Zustandsmatrix:

Daraus folgt die Zeitentwicklung der Polarisation.

Die Hamilton-Matrix definiert einen Vektor, der zeitabhängig sein darf:

Die Matrix der Präzession mit einem Vektor a der Achse, Betrag = Winkelgeschwindigkeit, sei definiert als die gleichförmige Rotation

Die Präzession löst die Differenzialgleichung für kreisende Vektoren b(t):

Kurz, ist die Matrix der Abbildung .

Außerdem ist die Einheitsmatrix.

Die Bewegung einer Polarisation b(t) im rotierenden Feld sieht nun keine konstante Achse, sondern eine, die ihrerseits rotiert:

Auszunutzen: Das Kreuzprodukt ist kovariant unter allen Rotationen, :

Die Vektorfunktion

hat eine Zeitableitung laut Produktregel

Es folgt mit dem vorhin angeschriebenen zweiten Term

Der transformierte Vector d hat also eine konstante Achse und daher eine bekannte Lösung seiner Bewegungsgleichung:

Aus der Definition von d folgt die gesuchte Trajektorie:

Interpretation. Die Vektorrechnung läuft auf dasselbe Verhalten hinaus wie die vorher diskutierte Spinor-Rechnung. Das anregende Magnetfeld rotiert um eine Achse c und es stellt für die Polarisation eines Spins selbst eine Drehachse a dar. Beide Axial-Vektoren sind etwa gleich lang nahe der Resonanz, also die Vektordifferenz (a(0)-c) ist relativ klein. Die Polarisation bewegt sich mit einer doppelten Präzession. (Die Differenz entspricht dem obigen Wechselwirkungs-Teil des Hamilton-Operators, der statische Teil in Richtung der Achse c wurde abgespalten.) Zuerst wird 'langsam' um die kurze Differenz-Achse gedreht, das Resultat wird dann synchron zum Drehfeld um die Hauptachse c rotiert. Die Komponente der Polarisation in Richtung der Hauptachse macht die schnelle Drehung nicht mit, sondern oszilliert mit der langsamen Frequenz. Die Effekte sind manchmal als Rabi-Oszillationen zu beobachten, am langsamsten und mit größter Amplitude auf der Resonanz.

Zeitentwicklung der Spin-Zustände, direkt berechnet[Bearbeiten]

Variablentransformation auf konstante Koeffizienten. Ansatz

Nach dem Einsetzen verschwindet mit a=iγ/2 die explizite Zeit.

Gesucht wird oszillierendes Zeitverhalten mit Faktor exp(izt/2). Die Matrix der Form hat Eigenwerte z gemäß

0=(x-z)(-x-z)-y²=z²—(x²+y²). Das bedeutet hier:
bestimmt das Verhältnis

Damit ergibt sich eine Basis von zwei nicht-normierten Lösungen der homogenen Differenzialgleichung, in den ursprünglichen Variablen:

In den allgemeinen Lösungen, Linearkombinationen dieses Paares, erscheinen als Zeitfaktoren die Frequenzen (Ω±γ)/2. Kleine Frequenzen bis Null gibt es für gewisse Kombinationen von α,βγ.

Stichwort Resonanz, Rabi-Oszillationen. Das Thema kommt weiter unten zur Sprache, wo Spin-Systeme gedämpfte Schwingungen ausführen.

Warum kommt in allen lehrsamen Büchern und Texten — auch diesem hier — immerzu das rotierende Feld dran und nicht ein einfacheres, das in einer Dimension schwingt? Weil das rotierende Feld eine leichter berechenbare Approximation des oszillierenden ist. Das soll nun erläutert werden.

In den Übungen der vorigen Abschnitte war H von der Form

Zum Vergleich soll der Koppelterm mal nur reell sein:

Ein Standard-Ansatz eliminiert die Diagonale einer Hamilton-Matrix, indem man zu einem "Wechselwirkungs-Bild" übergeht.

In den neuen Variablen ist die Koppelmatrix zwar nulldiagonal, aber oszillierend zeitabhängig. Kein Problem scheinbar für unser Beispiel, der Phasenfaktor addiert sich nur zu dem schon vorgegebenen.

Bei der kreisenden Anregung sieht es so aus:

Die explizite Zeitabhängigkeit steckt in der Differenzfrequenz (α—γ). Normalerweise ist die Physik der Sache interessant im Bereich der Resonanz: die Frequenz der Anregung liegt nahe an der Frequenz, die dem Energiesprung zwischen den Eigenzuständen des Diagonalteils enspricht. Die Differenz ist wesentlich kleiner als die Resonanzfrequenz.

Bei der Anregung mit linearer Schwingung kommt so eine Gleichung an:

Entsprechend für die andere Spinor-Komponente und auch etwa für den Polarisations-Vektor, welcher Realteil und Imaginärteil des Koppelterms in Anspruch nimmt. Die Frequenz α—γ ist klein in praktischen Fällen, die Frequenz α+γ dagegen extrem groß. Man geht soweit, die Hochfrequenz-Terme zu vernachlässigen. Messbare Polarisation ist träge genug und der Mittelwert der Terme über viele Perioden dieser doppelten Resonanzfrequenz geht gegen Null. Das Problem des oszillierenden Feldes wird also reduziert auf das des Drehfeldes. Die Näherung mit der rotierenden Anregung wird also gern benutzt und führt schnell zu Lösungen, deren Kurvendiskussion weniger Mühe macht. Die Bloch-Gleichungen sind die Version des Modells für die Variablen der Zustandsmatrix. Es gibt Rabi-Oszillationen zu entdecken.

Die Approximation des rotierenden Feldes benutzt man bei der Diskussion von Absorption und Emission von Licht durch Systeme mit zwei Niveaus. Das klassisch behandelte Strahlungsfeld regt das Atom periodisch an. Die Bloch-Gleichungen bekommen Relaxations-Terme, welche die Spektrallinien verbreitern. Durch spontane Emission, durch Kollisionen, den Doppler-Effekt, Sättigungs-Effekt usw. Das halbklassische Modell funktioniert vielfach recht gut, obwohl es die theoretisch korrekte Quantenelektrodynamik absolut ignoriert.

Relaxation von Spin-Systemen[Bearbeiten]

In praktischen Experimenten hat man es zu tun mit einer großen Gesamtheit von Spin-Freiheitsgraden, häufig einer Menge von Protonen. Sie sollen sich so wenig gegenseitig beeinflussen, dass die inkohärente Mischung von Quantensystemen ein ausreichendes Modell ist. Mit einer geeigneten resonanten Anregung können (fast) alle Spins im konstanten Magnetfeld auf den Zustand der höheren Energie gebracht werden. Jedes Proton hat eine relativ ungeordnete Umgebung von Atomen, die mit temperatur-bewegt zeitlich veränderlichen Feldern einwirken. Je nach chemischer Umgebung wird die mittlere Polarisation mehr oder weniger schnell vom Anfangswert auf der Einheitskugel auf Werte im Inneren der Kugel absinken. Verschiedene Zeitkonstanten dieser Relaxation können gemessen werden und liefern Informationen über die Natur des Stoffes, der die Protonen enthält. Die Technik der magnetischen Kernspin-Resonanz vermag solche Dynamik räumlich aufgelöst zu sondieren und liefert dreidimensionale Bilder, die besonders der medizinischen Diagnose einen wichtigen Schub gaben.

Die Zustandmatrix des Gemisches ist eine gewichtete Summe von reinen Zustandsmatrizen. Für ein vollkommen klassisches Gemenge ohne Quanten-Interferenzen ist sie diagonal, die nicht-diagonalen komplexen Informationen über Phasen verschwinden. Alle Zustandsmatrizen, ob rein oder gemittelt, haben die Spur 1. Ihre diagonalen Einträge sind als dimensionslose 'Wahrscheinlichkeiten' zu interpretieren.

Die Relaxation wird durch zwei Zeitkonstanten beschrieben. Die erste oder 'longitudinale' treibt die Wahrscheinlichkeitsverteilung in die Richtung einer statischen Gleichgewichtsverteilung {g,(1-g)}. Die zweite oder 'transversale' treibt die nichtdiagonalen Einträge gegen Null. Longitudinal bezieht sich auf die Richtung des statischen magnetischen Feldes im Hintergrund, in dem die Indizes 1,2 zu den zwei Eigenzuständen der Energie gehören. Das variable Feld soll senkrecht zum statischen stehen und kann im Experiment ein- und ausgeschaltet werden. Auch das statische Feld kann bei Gelegenheit umgepolt werden.

Die Differenzialgleichung der gemischten Z-Matrix besteht aus den Termen der Quantenmechanik plus den Termen der Relaxation, die eine anwachsende Durchmischung oder Unordnung beschreiben sollen. Wenn der Punkt die totale Zeitableitung andeutet, "∂" den quantendynamischen unitären Teil, und g den Bruchteil der 1-Zustände im (thermischen) Gleichgewicht, dann hat das Modell die Gleichung:

Mit den Gleichungen der Polarisation interpretiert, fallen die Komponenten exponentiell ab mit der transversalen Zeitkonstanten und dem widerfährt Gleiches mit der longitudinalen Zeit . Wegen reicht es aus, das Paar 11,12 von Elementen der Z-Matrix zu diskutieren. Eine reelle und eine komplexe Zahl. Sie enthalten die drei reellen Variablen der Polarisation.

Der Endzustand dieses dissipativen (irreversiblen) Prozesses ist . Das Modell enthält zwei einfache klassische Markov-Prozesse, die über den nicht-dissipativen Hamilton-Operator gekoppelt werden. Das Gleichgewicht wird mit gedämpften Oszillationen erreicht. Im Gegensatz dazu schwingt jeder reine Quantenzustand im statischen wie im periodischen Feld unendlich lange. Die Interpretation hat zwei Aspekte: Ensemble-Mittelung über Spins, die alle zufällig leicht verschiedene Hamilton-Operatoren sehen. Oder Zeit-Mittelung für einen Spinzustand, dessen Hamilton-Operator immer wieder stochastisch verformt wird, etwa weil im Wärmebad andere Moleküle ihn anstoßen.

Eine Modellrechnung soll die Differenzialgleichung des Zwei-Niveau- Systems für alle Anfangsbedingungen lösen. Die Rechenschritte werden (hoffentlich) genau (genug) erklärt; am Ende kommt statt einer seitenlangen Formel ein Skript, das die Lösungen zeichnet. Auch seitenlang, aber vielleicht praktischer. Daneben wird der Spezialfall gleicher Zeitkonstanten diskutiert, der einfachere Ausdrücke hergibt.

Interessant ist etwa die Anfangsbedingung . Die ganze Population der Spins befindet sich im tieferen Energieniveau 2 des diagonalen Anteils von H, also ohne Anregung bei β=0. Der zeitvariable Term von H sei vergleichsweise klein und gefragt wird nach der Dynamik und dem stationären Endpunkt der relativen Häufigkeit vom angeregten Zustand 1, wenn zur Zeit Null der periodische Koppelterm eingeschaltet wird. Das Modell benutzt die Approximation des rotierenden Feldes und die Relaxation soll ohne Anregung das ganze Spin-Ensemble auf eine Gleichgewichtsverteilung (g,1-g) ziehen.

Die Gleichungen für Z sind inhomogen linear. Eine spezielle Lösung, möglichst stationär mit oszillierenden Werten, wird zuerst gesucht.

Beide Gleichungen werden mit a=0, γ+b=0 konsistent zeitkonstant.

Nach Real- und Imaginärteil aufgelöst, B =: x+iy:

Man kann eventuell den Nullpunkt der Zeit so setzen, dass reell wird, denn die einzige Zeitabhängigkeit der Lösung ist die Phasendrehung im Gleichtakt mit der Anregung .

Der Wert A ist also die Gleichgewichts-Konzentration des angeregten Niveaus, offenbar dimensionslos, gleich der statischen Gleichgewichts-Konzentration g wenn das oszillierende Feld β verschwindet. Wird die Anregung β extrem groß, konvergiert der Faktor y gegen (g-1/2)·(η/β), also A gegen (1/2). Unabhängig vom statischen Gleichgewicht enden die zwei Referenz-Niveaus mit gleicher Konzentration.

Als Funktion der Frequenz γ der Anregung kommt eine Resonanz auf am Punkt γ=α. Verschiedene Amplituden der stationären Lösung haben im Nenner eine typisch resonante Formel der Art (Δ²+(γ—α)²).

Ausgedrückt als Polarisations-Vektor, bewegen sich seine Komponenten in der xy-Ebene auf einem Kreis synchron mit dem Feld der Anregung. Das Verhältnis vom Betrag der Polarisation zur Feldstärke β ist die Suszeptibilität, die hier eine nicht-lineare Form annimmt. Das System geht in Sättigung.

Im Wechselwirkungs-Bild definiert man die Zustandmatrix als

,

relevant hier also

Mit δ:=γ—α gibt es dann einen Parameter weniger; es kommt zu folgendem Paar von Bloch'schen Gleichungen.

Es müssen noch alle Lösungen der homogenen Bloch-Gleichungen her. Dazu bemühen wir die Version in den ρ-Variablen und dürfen nach dem Vergessen der konstanten Terme mit g und 1 noch eine andere Phasendrehung ausüben:

In den Variablen u,w sind nur noch konstante Koeffizienten da.

Gesucht wird eine Basis von Eigenlösungen:

Das charakteristische Polynom sollte eine reelle Wurzel haben, sowie ein komplex-konjugiertes Paar k=(a±ib). Aus dessen Exponential müssen danach reelle Lösungen kombiniert werden.

Die seltene Gelegenheit, eine kubische Gleichung lösen zu müssen. Hier der halbe Algorithmus, aus Formelsammlungen entliehen:

k³+bk²+ck+d=0, k=y-b/3 ⇒ y³+3py+2q=0.
p=(-b²/3+c)/3, q=(2b²/27-bc/3+d)/2.

Das Vorzeichen von z:=q²+p³ bestimmt, wie es weiter geht. Ist es positiv, gibt es eine reelle und zwei komplex konjugierte Lösungen. Das wird hier angenommen; negative Realteile sind erwartet. (Die Prozedur bei drei reellen Nullstellen sei hier verschwiegen.)

Zu jedem der drei findet man einen Eigenvektor

Es kommen zwei verschiedene negative Realteile vor in den drei k-Werten, . Sowie ein imaginärer Wert . sind konjugierte Paare. Die homogenen Übergangslösungen ('Transienten') fallen also exponentiell mit der Zeit ab.

Allerdings: Weil reell sein muss, ist eine physikalische Basis homogener Lösungen der Satz der drei folgenden Tripel:

(Eigentlich wären sie schöner geschrieben als Spaltenvektoren?) Sie haben ein Zeitverhalten mit den drei Faktoren

Die inhomogene Lösung hat die Anfangswerte (a,b,c) mit

Die allgemeine Lösung ist eine Überlagerung der speziellen, zum Feld synchronen, inhomogenen Lösung und der gedämpften Oszillationen, also der drei Transienten, welche die homogenen Gleichungen auswerfen. Die Anfangsbedingungen für zur Zeit t=0 ergeben die drei reellen Koeffizienten dafür.

Die reelle 11-Komponente von ist zeitkonstant. Die von hat exponentiellen Abfall, die zwei anderen fallen exponentiell mit einer gedämpften Schwingung. Zusätzlich ist die komplexe 12-Komponente hochfrequent moduliert, bei allen vier y-Termen mit exp(-iγt).

Nehmen wir die folgende Anfangsbedingungen an,

Die Startwerte (a,b+if) der synchronen Lösung wurden schon ausgerechnet. Gleichungen für (c,d,e) zur Zeit t=0:

In der ersten Gleichung steht (e mal Null) weil .

Die Werte A,B,C, a,b,f sind bekannt, das Tripel c,d,e wird aus den drei linearen Gleichungen ermittelt. Damit ist das Anfangswertproblem vollständig gelöst.

(Leider verwirrend an dieser Stelle, die Rechnerei. Die transienten Lösungen kommen aus einer Komplexifizierung, die nicht dieselbe ist wie die des Koppelterms der Zustandsmatrix. Die (u,v) werden hier aus zwei Realteilen oder zwei Imaginärteilen gebaut.)

Die gedämpften Schwingungen, mit denen die Zustandsmatrix oder äquivalent der Polarisations-Vektor nach dem Umschalten einer Anregung reagiert, sind die beobachtbaren Rabi-Oszillationen. Sie entstehen, wenn an einer großen Anzahl gleichartiger atomarer Systeme die mittlere Polarisation auf die eine oder andere Art gemessen werden kann, wobei die anregende Frequenz in der Nähe der Resonanz liegt, also nahe bei der Übergangsfrequenz zwischen zwei Niveaus. Die absolute Stärke der Anregung definiert quantenmechanisch die Frequenz — wesentlich kleiner als die Resonanzfrequenz — mit der die Polarisation des Spin-Gemenges sich einschwingt auf einen letztlich synchronisierten Zustand zum antreibenden Feld.

Spezialfall mit einfacheren Formeln[Bearbeiten]

In einer abgespeckten Version kommt die Physik der periodischen Anregung deutlicher zum Vorschein. Keine Kubischen Gleichungen, wenn beide Zeitkonstanten der Relaxation gleich sind.

g=1; η=ν; δ:=γ—α.

Zuerst die eingeschwungene synchrone Lösung:

Wird die Stärke β der periodischen Anregung der dominierende Parameter, dann wandert die stationäre Konzentration des Zustandes 1 vom (thermischen?) Ruhepunkt g auf den Grenzwert (1/2) zu. Als Funktion der Frequenz γ geht der Effekt durch ein Maximum bei der Resonanz 0=δ=(γ-α) und beschreibt die übliche Lorentzkurve. Deren Verbreiterung wird von der Relaxationszeit und von der Amplitude β der Anregung bestimmt. Es gibt einen Effekt der "Sättigung".

Das charakteristische Polynom für die homogene Gleichung zerfällt.

Es hat die reelle Lösung k=—ν und mit (k+ν)²=—β²—δ² das konjugierte Paar

Die Zeitkonstante des exponentiellen Abfalls von Transienten ist also in allen Fällen .

Die Frequenz der Rabi-Schwingungen hat ihr Minimum, gleich β, auf dem Punkt der Resonanz δ=0.

Basis der Transienten: Zu jedem der drei gab es allgemein einen Eigenvektor mit

Hier:
Für Index 2 also:

Aber die Eigenfunktion zum reellen bekommt nun ein etwas anderes Tripel verpasst; ohne durch Null zu teilen sieht sie so aus:

Die relevanten Teile der Funktion zu ergeben die zwei anderen speziellen Lösungen:

Die Tripel hier sind jeweils Komponenten für . Die Linearkombination mit folgenden Koeffizienten (c,d,e) ist der transiente Teil der Lösung

Zusammengefasst. Die Zustandsmatrix einer Spin-Gesamtheit entspricht eindeutig dem Polarisationsvektor gemäß

Die Dynamik in einem oszillierenden Kraftfeld wird approximativ mit einem Drehfeld beschrieben — extrem hochfrequente Terme fallen dabei unter den Tisch.

Die Bloch-Gleichungen des Zustands im rotierenden Feld (Hochfrequenz γ)
mit fünf Parametern
und den Anfangsbedingungen
haben diese Lösung mit der 'kleinen' Frequenz κ und der Dämpfung ν:

Skript: gedämpfte Rabi-Oszillationen[Bearbeiten]

Das Skript in Python hat etwa 300 Zeilen, hängt von keiner anderen Code-Bibliothek ab und malt als Datei 'example1.svg' eine Lösung der Bloch-Gleichungen. Die Prozedur folgt genau dem, was in diesem Kapitel allgemein mit zwei Zeitkonstanten dargestellt wurde. Wie schon angesagt, für die grafische Darstellung wird ein hochfrequenter Phasenfaktor im nichtdiagonalen Element der Zustandsmatrix heraus dividiert.

# example plot for Bloch equation

from math import sqrt,pi,sin,cos,exp,log, atan2
from random import random as rnd # rnd() is in interval [0,1)

def cc(z) : return z.real - 1j*z.imag # conjugate
def expc(z) : # complex exp
  u=exp(z.real); a=z.imag; c=cos(a); s=sin(a); i=0+1j
  return u*(c+i*s)
def abs2(z) : a,b= z.real,z.imag; return(a*a+b*b)

class relaxation : # damped state-matrix oscillation. symbol self==rn 

  def cubiceqn(rn, b,c,d) : # solve xxx+bxx+cx+d=0, real coeffs b,c,d
    def root3(x) :
      (z,s) = (x,1) if x>0 else (-x,-1); return s*exp(log(z)/3.0)
    # 
    p= (-b*b/3.0+c) /3.0; q=(2*b*b*b/27.0-b*c/3.0+d) /2.0; z=q*q+p*p*p
    if z>0 : # yyy+3py+2q=0
      rz=sqrt(z); u=root3(rz-q); v=root3(-rz-q); w=abs(u-v)*sqrt(3)/2.0
      y1=u+v; y2=-y1/2.0 + 1j*w; y3=-y1/2.0 - 1j*w
      f=b/3.0; x1=y1-f; x2=y2-f; x3=y3-f; x=[x1,x2,x3]; err=0
      for i in range(3) : t= d+x[i]*(c+x[i]*(b+x[i])); err+=abs2(t)
      if err>1e-12 : print('cubic x1,x2,err='+str([x1,x2,err])); exit()
    else : print('cubiceqn: No complex root.'); exit()
    return x # x[0] real, x[1],x[2] conjugate

  def linear3(rn, a,b) : # solve ax=b
    x=[0,0,0]; y=[0,0,0]; err=[0,0,0]; errs=0
    det=( a[0][0]*(a[1][1]*a[2][2]-a[2][1]*a[1][2])
       -a[0][1]*(a[1][0]*a[2][2]-a[2][0]*a[1][2])
       +a[0][2]*(a[1][0]*a[2][1]-a[2][0]*a[1][1]) +0.0)
    x[0]=( b[0]*(a[1][1]*a[2][2]-a[2][1]*a[1][2])
       -a[0][1]*(b[1]*a[2][2]-b[2]*a[1][2])
       +a[0][2]*(b[1]*a[2][1]-b[2]*a[1][1])) /det
    x[1]=( a[0][0]*(b[1]*a[2][2]-b[2]*a[1][2])
      -b[0]*(a[1][0]*a[2][2]-a[2][0]*a[1][2])
      +a[0][2]*(a[1][0]*b[2]-a[2][0]*b[1])) /det
    x[2]=( a[0][0]*(a[1][1]*b[2]-a[2][1]*b[1])
      -a[0][1]*(a[1][0]*b[2]-a[2][0]*b[1])
      +b[0]*(a[1][0]*a[2][1]-a[2][0]*a[1][1])) /det
    for i in range(3) :
      for k in range(3) : y[i]+= a[i][k]*x[k]
      err[i]= y[i]-b[i]; errs+= abs(err[i])
    if errs>1e-12 : print('linear3 err='+str(err)); exit()
    return x

  def initpars(rn, pars) :
    if pars==[] : # testing, random setup
      g=rnd(); alpha=1.0; x=rnd(); y=rnd()
      t1=10*(1+rnd()); t2=15*(1+rnd())
      gamma= alpha*(1+(x-0.5)*0.1); beta=alpha*y*0.1
    else : g,t1,t2, alpha,beta,gamma, k =tuple(pars)
    eta=1.0/t1; nu=1.0/t2; delta= gamma-alpha
    b=eta+2*nu; nd=nu*nu+delta*delta; c=nd+2*eta*nu+beta*beta
    d= eta*nd +nu*beta*beta 
    z= rn.cubiceqn(b,c,d) # 'frequencies' of damped transients 
    return [g,t1,t2, alpha,beta,gamma, z]

  def initgeneral(rn, aa,bb,cc,pars=[]) : # return coeffs (c,d,e)
    pars= rn.initpars(pars) # random parameter setup
    g,t1,t2, alpha,beta,gamma, k =tuple(pars); delta=gamma-alpha
    i=1j; s=[0]*3; r=[0]*3; eta=1.0/t1; nu=1.0/t2
    f=beta*(g-0.5)/(beta*beta/eta+nu+delta*delta/nu)
    a= g - beta*f/eta # synchr solution; z11(0)=a, z12(0)=b+if
    b= -f*delta/nu
    if t1==t2 : # simpler
      k2= beta*beta+delta*delta; kappa=sqrt(k2)
      e= -(beta/kappa)*(cc-f)
      d= beta*((aa-a)*beta+(bb-b)*delta)/k2
      c= (aa-a)-d; x=[c,d,e]
    else :
      s1=-(eta+k[0])/beta;  r1= -delta* s1/(nu+k[0])
      s2=-(eta+k[1])/beta;  r2= -delta* s2/(nu+k[1])
      mat=[[1,1,0], [r1,r2.real,r2.imag], [s1,s2.real,s2.imag]]
      y=[aa-a, bb-b, cc-f]; x= rn.linear3(mat,y) # transient coeffs
    return pars, x 
 
  def transient(rn,t, pars, c,d,e, ang=0) : # transient superpositions
    g,t1,t2, alpha,beta,gamma, k =tuple(pars); delta=gamma-alpha
    i=1j; s=[0]*3; r=[0]*3; eta=1.0/t1; nu=1.0/t2; z=[0,0,0]; jmin=0 
    if t1==t2 : r[0]=-beta/delta; s[0]=0; jmin=1 # special, eta==mu
    for j in range(jmin,2) :
      s[j]=-(eta+k[j])/beta;  r[j]=-(delta* s[j])/(nu+k[j])
    p1=exp(k[0]*t);  q1=[p1,r[0]*p1,s[0]*p1]
    p2=expc(k[1]*t); q2=[p2,r[1]*p2,s[1]*p2]
    for h in range(3) : z[h]= c*q1[h] + d*q2[h].real + e*q2[h].imag 
    u=z[0]; w= (z[1]+i*z[2])*expc(-i*(ang+gamma*t)) 
    return [u,1-u,w]

  def stationary(rn,t, pars, ang=0) : # stationary solution @ time t
    g,t1,t2, alpha,beta,gamma, aux =tuple(pars); i=1j 
    eta=1.0/t1; nu=1.0/t2; p=gamma-alpha; q= beta*beta/eta+nu+p*p/nu
    y= beta*(g-0.5)/q;  x= -p*y/nu;  u=g-beta*y/eta; v=1-u
    w= (x+i*y)*expc(-i*(ang+gamma*t))
    return [u,v,w]

  def general(rn,t, pars, c,d,e, ang=0) :
    s=rn.stationary(t, pars, ang)
    r=rn.transient(t, pars, c,d,e, ang)
    return [s[0]+r[0], 1-s[0]-r[0], s[2]+r[2]]

  def blocheqn(rn,z,t, pars, ang=0) : #  dz(t)/dt of general z=[u,v,w]
    u,v,w=tuple(z) # Z11, Z22, Z12
    g,t1,t2, alpha,beta,gamma, aux =tuple(pars); i=1j 
    xt=expc(i*(ang+gamma*t)); yt=cc(xt)
    qu= (beta/2.0)*(yt*cc(w)-xt*w)  # oscillating quantum parts
    qw= alpha*w + (beta/2.0)*yt*(1-2*u)
    du = -(u-g)/t1 -i*qu # relaxation parts added
    dv = -du; dw = -w/t2 -i*qw
    return [du,dv,dw]

  def testgeneral(rn, aa,bb,cc) : # check the Bloch equn
    pars,z = rn.initgeneral(aa,bb,cc); c,d,e= tuple(z)
    y=rn.general(0,pars,c,d,e) # want y[0]=aa; y[2]= bb+i*cc
    t=rnd()*10.0; z0=rn.general(t,pars, c,d,e); dt=1e-6
    dz=rn.blocheqn(z0,t,pars); du,dv,dw=tuple(dz); u0,v0,w0=tuple(z0)
    # compare derivative and difference approximations
    d1=dt; z1=rn.general(t+d1, pars, c,d,e); u1,v1,w1=tuple(z1)
    d2=0.1*dt; z2=rn.general(t+d2,pars, c,d,e); u2,v2,w2=tuple(z2)
    cu=(u1-u0)/d1 - du; cw=(w1-w0)/d1 - dw; err1=abs2(cu)+abs2(cw)
    cu=(u2-u0)/d2 - du; cw=(w2-w0)/d2 - dw; err2=abs2(cu)+abs2(cw)
    print('Gen cde='+str([c,d,e])+' errs='+str([err1,err2]))
    err=abs2(y[0]-aa)+abs2(y[2]-(bb+1j*cc)) # init condition ok?
    print('Init abc='+str([aa,bb,cc])+' err='+str(err)) 

### graph output, svg plot
def div(p,q) : return int(p/q)

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 top125(y) :
  z=1e-9; y= max(abs(y),z)
  while z<y : z=10*z
  if (z/5.0)>y : z=z/5.0
  elif (z/2.0)>y : z=z/2.0
  return z

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 ball(sd, x,y,r,col) :
    return ('<circle fill="'+col+'" cx="'+str(x)+'" cy="'+str(y)+
      '" r="'+str(r)+'" />\n')
      
  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',fc='none') :
    # path w=with lc=linecolor fc=fillcol
    t='<path fill="'+fc+'" stroke="'+lc+'" stroke-width="'+str(w)+'" '
    t+= ' stroke-linecap="round" ' 
    return t + 'd="\n'+s+'" />\n'

  def curve(sd,scale,xy, fill='#000', width=1, bkgr='none') :
    (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,fc=bkgr)

  def grids(sd,scale, gx, gy, fill='#000', width=1) :
    (xs,fx,xmin, ys,fy,ymin) = scale; u=''
    (x0,dx,x1, ya,yb) = gx; (y0,dy,y1, xa,xb)= gy # in user units
    ya= aprx(ymin+fy*(ya-ys)); yb= aprx(ymin+fy*(yb-ys))
    xa= aprx(xmin+fx*(xa-xs)); xb= aprx(xmin+fx*(xb-xs))
    while x0<(x1-0.1*dx) : 
      x= aprx(xmin+fx*(x0-xs))
      s= sd.move(x,ya) + sd.line(x,yb); x0+=dx;
      u+= sd.path(width,s,lc=fill,fc='none')
    while y0<(y1-0.1*dy) : 
      y= aprx(ymin+fy*(y0-ys))
      s= sd.move(xa,y) + sd.line(xb,y); y0+=dy
      u+= sd.path(width,s,lc=fill,fc='none')
    return u

  def frame(sd,scale, xmin,ymin, xmax,ymax, dx,dy) :
    xyf=[xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
    bf= sd.curve(scale, xyf, fill='#777',bkgr='#ffd',width=3) 
    if (dx>0)and(dy>0)and(xmax>xmin)and(ymax>ymin) :
      gx=(xmin+dx,dx,xmax, ymin,ymax); gy=(ymin+dy,dy,ymax, xmin,xmax)
      bf+= sd.grids(scale, gx, gy, fill='#ccc', width=2) 
    return bf

  def curveset(sd,scale,xyz, fill='#000', width=1, bkgr='none') :
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(xyz) 
    #print('curveset n='+str(n))
    for i in range(n) : 
      px,py,pz= tuple(xyz[i])
      x= aprx(xmin+fx*(px-xs)); y= aprx(ymin+fy*(py-ys))
      s+= sd.move(x,y) if(pz==0) else sd.line(x,y)
    return sd.path(width,s,lc=fill,fc=bkgr)

  def ballset(sd,scale,plt, fill='#000', r=5) :
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(plt)
    for i in range(n) : 
      x= aprx(xmin+fx*(plt[i][0]-xs)); y= aprx(ymin+fy*(plt[i][1]-ys))
      s+= sd.ball(x,y,r,fill)
    return s

  def ballsrf(sd,scale,plt) : # vary radius and fillcolor
    (xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(plt)
    for i in range(n) : 
      x= aprx(xmin+fx*(plt[i][0]-xs)); y= aprx(ymin+fy*(plt[i][1]-ys))
      s+= sd.ball(x,y,plt[i][2],plt[i][3])
    return s

  def label(sd,scale,txt,tx,ty) : # tx,ty user 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 infobox(sd,s,x,y) :
    t=s.split('\n'); k=len(t); u=''
    for i in range(k) : u+= sd.labl(t[i],x,y+i*26,size=24)
    return u

# end class svgdump

def myplot(fn, xa,ya,xb,yb,dx,dy, pball, pcurve, s, style=0) :
  # dots(x,y,r,color) in pball array. window size 3:2 
  winx=1000;winy=667; dwx=20;dwy=20 # winx=1000;winy=1000;dw=25
  xmin=xa;ymin=ya;xmax=xb;ymax=yb 
  xyf=[xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
  xs,fx,xpix= xmin, (winx-2*dwx)/(xmax-xmin+0.0),dwx
  # scale user units, xpix,ypix is pixel-pos at user-min.
  ys,fy,ypix= ymin,-(winy-2*dwy)/(ymax-ymin+0.0),winy-dwy
  scale=(xs,fx,xpix, ys,fy,ypix)
  colo=['#777','#000','#00a','#0a0','#077','#770','#a00']
  wid= [3]*len(colo); ch='fedccba98765'
  sd=svgdump(); buf=sd.plotbegin(winx,winy,1)
  buf+= sd.frame(scale, xa,ya, xb,yb, dx,dy)
  buf+= sd.ballsrf(scale,pball)
  for k in range(len(pcurve)) : 
    j=(k+1) % (len(colo))
    buf+= sd.curveset(scale,pcurve[k], fill=colo[j], width=wid[j])
  if len(s)>0 : buf+= sd.infobox(s,3*winx/5,2*winy/3)
  g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()

### main

def example1(fn, mode=0) :  # plot a solution to Bloch equn.
  ### set 6 model parameters alpha,beta,gamma,t1,t2,g
  f=0.95; alpha=1; gamma= f*alpha; beta=0.2*alpha; g=0.5
  eta=gamma/100.0; tau=2*eta; t1=1.0/eta; t2=1.0/tau;  k=[0,0,0]
  if mode==1 : t2=t1; tau=eta # force t1=t2, some formulae change
  pars=[ g,t1,t2, alpha,beta,gamma, k]; rn=relaxation()
  ### set initial conitions of state matrix Z11=A, Z12=B+iC 
  z11=0.9; z12=0.25
  ### get linear coeffs c,d,e for 3 transient solutions
  pars, cde = rn.initgeneral(z11,z12,0, pars); c,d,e=tuple(cde)
  ### compute and plot the array of function values 
  np=250+1; tmax=50.0/beta; # scan some periods of LF-oscillation
  dt=tmax/(np-1.0); u=[0]*np;v=[0]*np;w=[0]*np
  pen=0; pu=[[]]*np; pv=[[]]*np; pw=[[]]*np
  for n in range(np) : # remove HF-oscillations from Z12
    t=n*dt; z= rn.general(t,pars, c,d,e); x=z[2]*expc(1j*gamma*t) 
    u[n],v[n],w[n]= z[0], x.real+0.5, x.imag+0.5 
    pu[n]=[t,u[n],pen]; pv[n]=[t,v[n],pen]; pw[n]=[t,w[n],pen]; pen=1
  umin=min(u);umax=max(u); #print('u min,max '+str([umin,umax]))
  vmin=min(v);vmax=max(v); #print('v min,max '+str([vmin,vmax]))
  wmin=min(w);wmax=max(w); #print('w min,max '+str([wmin,wmax]))
  s=('black: Z11\nblue,green: Re,Im(Z12)\nTmax='+str(tmax)+
    '\nAlpha,Beta,Gamma='+str([alpha,beta,gamma])+
    '\nG,Eta,Tau='+str([g,eta,tau])) ; 
  myplot(fn, 0,0, tmax,1, 0.1*tmax, 0.1, [], [pu,pv,pw], s)

if True :
  #rn=relaxation(); rn.testgeneral(0.5,pi,0)
  example1('examp1')
  #example1('examp2',1)