Komplexe Dynamische Systeme

Barbara Drossel, Technische Universitat Darmstadt

Wintersemester 2009/10

Inhaltsverzeichnis

1 Einfuhrung 1

2 Stochastische Prozesse 2

2.1 Markov-Ketten . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.2 Random Walk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.1 Einfuhrung . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.2 Kontinuumsbeschreibung . . . . . . . . . . . . . . . . . . 72.2.3 Die Diffusionsgleichung . . . . . . . . . . . . . . . . . . . 82.2.4 Varianten des Random Walk, die dieselbe Kontinuumsbe-

schreibung haben . . . . . . . . . . . . . . . . . . . . . . . 82.3 Master-Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3.1 Raten und Wahrscheinlichkeiten . . . . . . . . . . . . . . 102.3.2 Die Mastergleichung . . . . . . . . . . . . . . . . . . . . . 112.3.3 Exkurs: Monte-Carlo-Simulationen zur Bestimmung sta-

tistischer Eigenschaften . . . . . . . . . . . . . . . . . . . 122.4 Kramers-Moyal-Entwicklung und Fokker-Planck-Gleichung . . . . 14

2.4.1 Absorbierende Randbedingungen . . . . . . . . . . . . . . 162.5 Langevin-Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . 19

2.5.1 Brownsche Bewegung . . . . . . . . . . . . . . . . . . . . 212.5.2 Brownsches Teilchen im Kraftfeld . . . . . . . . . . . . . . 222.5.3 Brownsche Bewegung im Doppelmuldenpotenzial: Kramers

Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3 Nichtlineare Dynamik 25

3.1 Nichtlineare Dynamik in einer Dimension . . . . . . . . . . . . . 293.1.1 Phasenportraits . . . . . . . . . . . . . . . . . . . . . . . . 293.1.2 Stabilitat von Fixpunkten . . . . . . . . . . . . . . . . . . 293.1.3 Sattelknotenbifurkationen . . . . . . . . . . . . . . . . . . 313.1.4 Die transkritische Bifurkation . . . . . . . . . . . . . . . . 323.1.5 Die Heugabelbifurkation . . . . . . . . . . . . . . . . . . . 333.1.6 Bifurkationen mit mehr als einem Parameter . . . . . . . 36

3.2 Nichtlineare Dynamik in zwei Dimensionen . . . . . . . . . . . . 403.2.1 Fixpunkte . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.2.2 Bifurkationen in zwei Dimensionen . . . . . . . . . . . . . 423.2.3 Phasenportraits und topologische Uberlegungen . . . . . . 46

3.3 Seltsame Attraktoren . . . . . . . . . . . . . . . . . . . . . . . . . 503.4 Die fraktale Dimension . . . . . . . . . . . . . . . . . . . . . . . . 51

3.4.1 Kastchenzahlmethode . . . . . . . . . . . . . . . . . . . . 52

2

3.4.2 Das Spektrum der fraktalen Dimensionen . . . . . . . . . 52

4 Phasenubergange und kritische Phanomene 55

4.1 Einfuhrung in Phasenubergange . . . . . . . . . . . . . . . . . . . 554.2 Mean-Field-Theorie . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.2.1 Landau-Theorie fur Phasenubergange . . . . . . . . . . . 594.2.2 Dynamik . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.2.3 Kritische Exponenten . . . . . . . . . . . . . . . . . . . . 61

4.3 Ginzburg-Landau-Funktional und Fluktuationen in gaußscher Nahe-rung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.4 Ortsraumrenormierung und Skalengesetze . . . . . . . . . . . . . 714.4.1 Ortsraumrenormierung . . . . . . . . . . . . . . . . . . . . 714.4.2 Skalengesetze . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.5 Perkolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.5.1 Definition und Beispiele . . . . . . . . . . . . . . . . . . . 794.5.2 Kritische Exponenten und Skalenverhalten . . . . . . . . . 804.5.3 Ortsraumrenormierung . . . . . . . . . . . . . . . . . . . . 82

5 Musterbildung in raumlich ausgedehnten Systemen 86

5.1 Wellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 865.1.1 Fisher-Kolmogoroff-Gleichung . . . . . . . . . . . . . . . . 875.1.2 Storungstheorie fur die Form der Front . . . . . . . . . . . 895.1.3 Stabilitat der Front . . . . . . . . . . . . . . . . . . . . . . 90

5.2 Turing-Muster . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915.2.1 Allgemeiner Formalismus . . . . . . . . . . . . . . . . . . 915.2.2 Zwei Klassen von Systemen . . . . . . . . . . . . . . . . . 945.2.3 Das Schnakenberg-Modell . . . . . . . . . . . . . . . . . . 955.2.4 Das Gierer-Meinhardt-Modell . . . . . . . . . . . . . . . . 975.2.5 Tierfellmuster . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.3 Die komplexe Ginzburg-Landau-Gleichung . . . . . . . . . . . . . 995.3.1 Ordnungsparameterdynamik in der Nahe eines kritischen

Punktes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.3.2 Anwachsen einer diffusionsgetriebenen Instabilitat . . . . 1005.3.3 Die komplexe Ginzburg-Landau-Gleichung . . . . . . . . . 101

3

Kapitel 1

Einfuhrung

Diese Vorlesung befasst sich mit fortgeschrittenen Themen der statistischenPhysik und nichtlinearen Dynamik, die fur die aktuelle Forschung von Relevanzsind. Wahrend sich die Kursvorlesung uber statistische Physik uberwiegend aufdie Gleichgewichtsphysik konzentriert, sind die meisten realen Systeme offen,d.h. sie konnen nicht von ihrer Umgebung abgekoppelt werden und konnenkomplexe Dynamik und Strukturbildung zeigen. Heutzutage werden Methodender statistischen Physik auf so verschiedene Themen wie die Ausbildung vonVerkehrsstaus, die Fluktuationen der Borsenkurse, die Populationsdynamik inOkosystemen, die Evolution der DNA und das Entstehen von Sanddunen ange-wandt. Diese Vorlesung soll in die verschiedenen Methoden und Konzepte aufdiesem Gebiet einfuhren.

Im zweiten Kapitel werden zunachst stochastische Prozesse behandelt. Im-mer dann, wenn fur die Anderungen des Wertes einer (oder mehrerer) Varia-blen nur Wahrscheinlichkeiten angegeben werden konnen, liegt ein stochasti-scher Prozess vor. Eines der alteren und beruhmten Beispiele ist die BrownscheBewegung. Allerdings geht die Theorie stochastischer Prozesse auf das Glucks-spiel zuruck. Es werden Markov-Ketten, Master-Gleichungen, Fokker-Planck-Gleichungen und Langevin-Gleichungen vorgestellt und ein paar Beispiele vor-gerechnet.

Im dritten Kapitel behandeln wir nichtlineare Dynamik. Die mathemati-schen Gleichungen hierzu sind deterministisch. Es muss also vorausgesetzt wer-den, dass stochastische Effekte entweder unwichtig sind oder sich herausmit-teln, weil die betrachteten Variablen makroskopische Großen eines Systems ausvielen mikroskopischen Freiheitsgraden sind. Hier interessieren wir uns fur dieAbhangigkeit des dynamischen Verhaltens von den Parametern des Systems undbehandeln Fixpunkte, Grenzzyklen und ihre Bifurkationen. Die Ubungsbeispielesind meist aus der Okologie gewahlt.

Ab dem vierten Kapitel wechseln wir dann zu raumlich ausgedehnten Syste-men. Wir behandeln zunachst Systeme, die sich in der Zeit einfach verhalten(namlich die im Gleichgewicht sind), und befassen uns mit Phasenubergangenund kritischen Phanomenen. Wir betrachten die beiden klassischen Beispiele,das Ising-Modell und Perkolation. Im funften Kapitel schließlich wird das Aus-bilden von Mustern in Nichtgleichgewichtssystemen behandelt. Ein wichtigerTeil des Kapitels befasst sich mit Wellen, ein anderer mit Turing-Instabilitaten,die z.B. Tierfellmuster erzeugen.

1

Kapitel 2

Stochastische Prozesse

Die Theorie stochastischer Prozesse begann im 17. Jahrhundert in Verbindungmit dem Glucksspiel. Als Schlusselereignis gilt dabei ein Briefwechsel zwischenBlaise Pascal und Pierre de Fermat im Jahr 1654, in dem sie sich mit Wahr-scheinlichkeiten beim Wurfeln und der fairen Aufteilung eines Spielgewinns be-fassten. Vom Konzept her nicht allzuweit vom Glucksspiel entfernt ist die Fi-nanzmathematik, die im Jahr 1900 mit der Doktorarbeit “Theorie de la specula-tion” von Louis Bachelier begann. Heute sind die Anwendungen stochastischerProzesse sehr vielfaltig: uberall dort, wo viele unvorhersagbare oder nicht imeinzelnen messbare Ereignisse beteiligt sind, hat man es mit stochastischen Pro-zessen zu tun. Hierzu gehoren thermisches Rauschen, Prozesse mit Ubergangs-raten (Mutationen in biologischer Evolution, Geburtenraten und Sterberaten,chemische Reaktionsraten), biologische Motoren, Verkehrsphysik, Finanzphysik,Versicherungsphysik, Krankheitsausbreitung. Die Große, die dem stochastischenProzess unterworfen ist, nennt man eine Zufallsvariable. Beispiele fur Zufalls-variablen sind gewurfelte Zahlen, Borsenkurse, die Große einer Population, dieZahl der in den Vorlesungsstunden anwesenden Studenten, etc.

Wir behandeln im Folgenden nur Markov-Prozesse. Dies sind Prozesse, beidenen die Wahrscheinlichkeiten fur das nachste Ereignis nicht von der Vorge-schichte, sondern nur vom momentanen Zustand des Systems abhangen. Wirbeschranken uns weiterhin auf stationare Markov-Prozesse, d.h. Prozesse, beidenen die Wahrscheinlichkeiten nicht explizit von der Zeit abhangen.

Wir beginnen mit Systemen, die nur eine begrenzte Zahl von moglichenZustanden haben, zwischen denen Ubergangswahrscheinlichkeiten vorgegebensind. Die Zustandsabfolge in derartigen Systemen ist eine Markov-Kette. Alsnachstes behandeln wir das Paradebeispiel eindimensionaler stochastischer Pro-zesse, den “Random Walk”, anhand dessen wir die diskrete und kontinuierlicheFormulierung stochastischer Prozesse in Raum und Zeit einuben konnen. In denfolgenden Abschnitten werden wir uns dann mit Master-Gleichungen, Fokker-Planck-Gleichungen und Langevin-Gleichungen allgemeiner befassen.

2.1 Markov-Ketten

Markov-Ketten sind Markov-Prozesse, die diskrete Werte der Zufallsvariablenhaben und diskret in der Zeit sind. Wenn die Zufallsvariable nur endlich viele

2

Werte annehmen kann, kann man die Ubergangswahrscheinlichkeiten zwischenihnen durch eine Matrix Q darstellen, wobei das Element Qji die Ubergangs-wahrscheinlichkeit pro Zeitschritt von Zustand i nach Zustand j angibt. Alsogilt Qji ∈ [0, 1] fur alle Elemente von Q. Da die Wahrscheinlichkeit normiertsein muss, ist die Summe der Elemente einer Spalte 1, d.h.

∑

i Qij = 1. Alseinfaches Beispiel betrachten wir zwei Gefaße. Im linken Gefaß befinden sichanfangs 2 weiße Kugeln, im rechten Gefaß drei rote Kugeln. Wir greifen nunmit geschlossenen Augen in beide Gefaße, nehmen aus jedem eine Kugel herausund vertauschen die beiden Kugeln. Wenn dies beliebig oft gemacht wird, be-kommen wir einen stochastischen Prozess. Als Variable konnen wir die Zahl yder weißen Kugeln im linken Gefaß nehmen. Sie hat drei mogliche Werte, 0, 1,2. Die Matrix Q ist also eine 3× 3 Matrix. Durch Abzahlen aller Moglichkeitenkonnen wir die Elemente dieser Matrix ermitteln und finden

Q =

1/3 1/3 02/3 1/2 10 1/6 0

Der oben beschriebene Anfangszustand ist der Spaltenvektor ~r0 = (0, 0, 1)t.Wenn wir an ihn von links die Matrix Q multiplizieren, erhalten wir die Wahr-scheinlichkeiten dafur, dass nach einem Schritt im linken Gefaß 0,1 oder 2 Ku-geln sind, ~r1 = Q~r0. Nach n Schritten erhalten wir ~rn = Qn~r0. In den Ubungenwerden wir zeigen, dass sich im Grenzfall n → ∞ die stationare Verteilung~rn → r∞ = (0.3, 0.6, 0.1)t ergibt, und zwar unabhangig vom Anfangszustand.

Im Folgenden leiten wir wichtige Eigenschaften der Matrix Q und ihrer Ei-genwerte und Eigenvektoren her, die fur alle Markov-Ketten gelten. Sei ~v einEigenvektor und λ der zugehorige Eigenwert. Fur jedes naturliche n gelten fol-gende Beziehungen:

∑

i

|∑

j

(Qn)ijvj | =∑

i

|λnvi| = |λ|n∑

i

|vi| (2.1)

∑

i

|∑

j

(Qn)ijvj | ≤∑

i

(∑

j

(Qn)ij |vj |) =∑

i,j

(Qn)ij |vj | =∑

i

|vi| (2.2)

Daraus erhalten wir die folgenden Aussagen:

1. Der großte Eigenwert von Q ist 1. Dass 1 ein Eigenwert ist, folgt unmit-telbar aus der Tatsache, dass der Vektor (1, 1, . . . 1) ein linker Eigenvektorvon Q ist, weil die Summe der Elemente jeder Spalte von Q 1 ist. Ausdem Vergleich von (2.1) und (2.2) folgt sofort, dass es keinen Eigenwertmit einem Betrag großer als 1 geben kann.

2. Die Summe∑

i vi der Elemente eines Eigenvektors muss verschwinden,wenn der zugehorige Eigenwert nicht 1 ist. Dies sieht man, wenn man dieBeziehungen (2.1) und (2.2) ohne Betragsstriche schreibt. Dann wird das≤ zum Gleichheitszeichen, und die rechten Seiten der beiden Gleichungenmussen gleich sein. Dies geht nur, wenn λ = 1 ist oder

∑

i vi = 0.

3. Wenn es ein n gibt, so dass Qn lauter positive Eintrage hat, dann hatQ nur einen Eigenvektor zum Eigenwert 1, und dieser hat lauter positiveEintrage. Alle anderen Eigenwerte haben einen Betrag kleiner als 1. Umdies zu zeigen, stellen wir zuerst fest, dass das Gleichheitszeichen im ersten

3

Schritt von (2.2) fur ein System mit lauter positiven Qnij nur dann gilt,

wenn alle Komponenten von ~v positiv sind. Gabe es einen Eigenwert 6= 1mit dem Betrag 1, waren die Komponenten des zugehorigen Eigenvektorswegen dem vorigen Punkt nicht positiv, und man hatte einen Widerspruchzwischen der Ungleichung (2.2) und der Gleichung (2.1). Gabe es zweiverschiedene Eigenvektoren ~u und ~w zum Eigenwert 1, konnte man eineLinearkombination ~v = α~u + β ~w bilden, deren Elemente nicht alle positivsind. ~v ware aber ebenfalls ein Eigenvektor zum Eigenwert 1. Damit wurdeman wieder einen Widerspruch zwischen (2.1) und (2.2) erhalten.

Anschaulich bedeutet ein Qn mit lauter positiven Eintragen, dass dasSystem ergodisch ist: Wenn alle Elemente von Qn positiv sind, kann man inn Schritten von jedem Zustand zu jedem anderen kommen. Wie oft welcherZustand in einer unendlich langen Zeitreihe besucht wird, kann daher nichtvom Anfangszustand abhangen. Der Eigenvektor zum Eigenwert 1 gibtden Anteil der Zeit an, die das System in einer unendlich langen Zeitreihein jedem Zustand verbringt. Gleichzeitig gibt er fur ein Ensemble vonSystemen an, welcher Anteil des Ensembles in jedem der Zustande ist. Dasobige Beispiel mit den zwei weißen und drei roten Kugeln ist ergodisch,da Q2 nur positive Eintrage hat, wie man leicht nachprufen kann.

4. Ein ergodisches System geht fur jede Anfangsverteilung ~r0 mit der Zeitin die stationare Verteilung ~v1 uber. ~v1 ist der einzige Eigenvektor zumEigenwert 1, alle seine Elemente sind positiv. Um dies zu zeigen, schreibenwir die Anfangsverteilung als Linearkombination der stationaren Vertei-lung und der anderen Eigenvektoren von Q,

~r0 =∑

µ

cµ~vµ .

Nach n Schritten wird daraus

~rn =∑

µ

cµλnµ~vµ .

Fur n → ∞ bleibt nur der Eigenvektor zum Eigenwert 1 ubrig. Der zweit-großte Eigenwert bestimmt, wie schnell sich das System der stationarenVerteilung annahert, da die Beitrage der anderen Eigenvektoren schnellergegen Null gehen. (Wenn es keine Basis aus Eigenvektoren gibt, kann manmit Hilfe von (2.2) dennoch zeigen, dass mit jeder Iteration der Abstand∑

i |(~rn)i − (~v1)i| zur stationaren Verteilung kleiner wird.)

5. In einem ergodischen System gilt, dass die Matrix Qn fur n → ∞ in jederSpalte ~v1 stehen hat. Denn es ist

Qn = (Qn−1)Q .

Eine Spalte der Matrix auf der linken Seite berechnet sich also als Produktvon Qn−1 mit der entsprechenden Spalte von Q. Im Limes n → ∞ istaber das Produkt von Qn−1 mit jedem beliebigen Vektor, der im Prinzipeine Anfangsverteilung darstellen kann, der Eigenvektor ~v1 (siehe vorigerPunkt).

4

Wir haben uns in den letzten drei Punkten auf ergodische Systeme speziali-siert, da das anfangs eingefuhrte Beispiel ergodisch ist. Zum Schluss betrachtenwir Beispiele fur Systeme, die nicht ergodisch sind. Die erste Beispiel ist eineMatrix Q, die Blockdiagonalform hat. Dann gibt es offensichtlich mehrere Ei-genvektoren zum Eigenwert 1. Das System zerfallt in disjunkte Teilsysteme, indenen eine Zeitreihe jeweils gefangen bleibt.

Auch Systeme mit nur einem Eigenvektor zum Eigenwert 1 konnen nicht-ergodisch sein. Solche Systeme mussen transiente Zustande haben. Dies sindZustande, die mit einer nichtverschwindenden Wahrscheinlichkeit nie wiederauftreten, wenn man von ihnen startet. Sie konnen daher in der stationaren Ver-teilung nicht auftreten, und alle Elemente von ~v1, die zu transienten Zustandengehoren, sind folglich Null. Ein rekurrenter Zustand dagegen kommt mit Si-cherheit wieder vor und tritt daher auch in der stationaren Verteilung auf. Einergodisches System hat nur rekurrente Zustande. Einfache Beispiele fur Syste-me mit transienten Zustanden sind in der folgenden Abbildung gezeigt. In denKreisen stehen die verschiedenen Zustande, die einfach von 1 bis 3 durchnum-meriert sind, und an den Pfeilen stehen die Ubergangswahrscheinlichkeiten. Imersten Beispiel sind die Zustande 1 und 2 transient. Der Zustand 1 im zweitenund dritten Beispiel ist ein besonderer transienter Zustand, weil man in ihnnicht zuruckkehren kann, wenn man ihn einmal verlassen hat. Man nennt ihnGarten-Eden-Zustand. Zustand 3 im ersten Beispiel und Zustand 2 und 3 imdritten Beispiel sind Zustande, aus denen man nicht nehr herauskommt, wennman einmal drin ist. Man nennt solche Zustande absorbierende Zustande.

1

2

3

1 32

1

2 3

0.5

0.5

0.20.8

1

0.50.5

0.5

0.5

1

1 1

0.2 0.2

0.6

Aufgaben

1. Gib fur die dargestellten einfachen Beispiele jeweils die Matrix Q und dieMatrix limn→∞ Qn an.

2. Berechne fur das Beispiel mit den beiden Gefaßen die Eigenwerte undEigenvektoren von Q.

3. Die Zufallsvariable y sei die beim Wurfeln gewurfelte Zahl. Schreibe dieMatrix Q fur diesen Fall auf. Was sind die Eigenwerte von Q, und wassind die Potenzen Qn ? Interpretiere die Ergebnisse anschaulich.

4. Betrachte einen geschlossenen Ring aus N Platzen. Auf einem dieser Platzesteht ein Spielstein. Um zu entscheiden, in welcher Richtung sich der Spiel-stein bewegt, wird eine Munze geworfen. Bei “Kopf” bewegt sich der Spiel-

5

stein um einen Platz im Uhrzeigersinn, bei “Zahl” um einen Platz gegenden Uhrzeigersinn. Was ist die Matrix Q, was ist der Eigenvektor zumEigenwert 1?

5. Betrachte eine Kette aus N Platzen. Auf einem dieser Platze steht einSpielstein. Um zu entscheiden, in welcher Richtung sich der Spielsteinbewegt, wird eine Munze geworfen. Bei “Kopf” bewegt sich der Spielsteinum einen Platz nach rechts, bei “Zahl” um einen Platz nach links. Wennder Spielstein sich nicht in der gewunschten Richtung bewegen kann, weiler an einem der Randplatze sitzt, bewegt er sich in der anderen Richtung.Was ist die Matrix Q, was ist der Eigenvektor zum Eigenwert 1?

2.2 Random Walk

Von jetzt an heben wir die Beschrankung auf nur endlich viele mogliche Werteder Zufallsvariablen auf. Wir betrachten aber nur solche Zufallsvariablen, dieganzzahlige oder reelle Werte annehmen konnen. Die entsprechenden stochasti-schen Prozesse spielen sich also in einer Dimension ab. Die Verallgemeinerungauf mehr Dimensionen ist nicht schwer, aber wird in dieser Vorlesung kaumerwahnt.

Wir betrachten nur solche stochastischen Prozesse, bei denen die Ubergangeauf eine lokale Umgebung begrenzt sind, d.h. bei denen die Werte der Zufallsva-riable nach einem Schritt oder einer Zeiteinheit in der Nahe des Ausgangswertesbleiben. Das Paradebeispiel fur derartige Prozesse ist der Random Walk, denwir in diesem Teilkapitel behandeln. Die allgemeineren eindimensionalen Pro-zesse, die wir danach betrachten werden, konnen mit Hilfe des Random Walksleichter verstanden werden.

2.2.1 Einfuhrung

Man stelle sich einen Besoffenen vor, der versucht, nach Hause zu gehen. Da eretwas durcheinander ist, geht nicht jeder seiner Schritte in die richtige Richtung.Wir gehen davon aus, dass jeder seiner Schritte die Große ∆x hat und mitWahrscheinlichkeit p nach rechts geht und mit Wahrscheinlichkeit q = 1 − pnach links. Wenn p > q ist, geht er haufiger nach rechts als nach links undkommt seinem Zuhause im Laufe der Zeit naher.

Der Weg des Besoffenen ist ein Beispiel fur einen Zufallspfad (“randomwalk”) auf einem eindimensionalen Gitter mit der Gitterkonstanten ∆x. DerPfad beginnt am Platz 0 und geht in jedem Schritt mit Wahrscheinlichkeit pnach rechts und mit Wahrscheinlichkeit q nach links.

Wir fragen als erstes, wie groß die Wahrscheinlichkeit p(m, N) ist, dass derPfad nach N Schritten an der Position m ist. Um nach N Schritten bei m zusein, mussen (N+m)/2 von den N Schritten nach rechts und (N−m)/2 Schrittenach links gehen. (m + N muss eine gerade Zahl sein.) Also ist

p(m, N) =

(

NN+m

2

)

p(N+m)/2q(N−m)/2 . (2.3)

Wenn N groß ist, erhalten wir unter Verwendung des zentralen Grenzwertsatzeseine Gaußverteilung. Der Mittelwert der Gaußverteilung ist N mal die mittlere

6

Positionsanderung eines Schrittes (in Einheiten von ∆x), 〈m〉 = p− q = 2p− 1,und die Varianz der Gaußverteilung ist N mal die Varianz eines Schrittes σ2 =4pq. Also haben wir nach N Schritten

p(m, N) =2√

8πNpqe−

(m−N(p−q))2

8Npq . (2.4)

Der Vorfaktor ist so gewahlt, dass die Summe uber alle m-Werte normiert ist.Hierbei muss man beachten, dass fur (un)gerade N nur (un)gerade Werte vonm auftauchen, so dass der Abstand zweier benachbarter Werte von m die Große2 hat. Diese Schrittgroße steht im Zahler des Vorfaktors.

Aufgaben

1. Peter und Paul treffen sich jeden Abend zum Kartenspielen. Derjenige, deram Ende des Abends gewonnen hat, bekommt vom anderen 10 Euro. Beidespielen gleich gut, so dass jeder die gleiche Gewinnchance hat. Inzwischensind 10 Jahre vergangen. Was ist die Wahrscheinlichkeit, dass einer derbeiden insgesamt mehr als 5000 Euro verloren hat? Berechne die Antwortunter der Annahme, dass der zentrale Grenzwertsatz verwendet werdendarf. Warum ist die Naherung, die man hierbei macht, eine gute Naherung?

2. Zur Zeit kommen in Deutschland bei 48.5 Prozent aller Geburten Madchenzur Welt und bei 51.5 Prozent aller Geburten Jungen. Was ist die Wahr-scheinlichkeit dafur, dass bei 10000 Geburten mehr als 5500 Jungen gebo-ren werden? Berechne die Antwort unter der Annahme, dass der zentraleGrenzwertsatz verwendet werden darf. Warum ist die Naherung, die manhierbei macht, eine gute Naherung?

3. Dennis und Lotta sollen beide die Wahrscheinlichkeit berechnen, dass derBesoffene nach 1000 Schritten der Große 1m zuhause angekommen ist,wenn sein Haus 1 km vom Lokal entfernt ist. Der Besoffene ist so besoffen,dass sie p = q = 1/2 wahlen durfen. Lotta verwendet Formel (2.3) underhalt das Ergebnis (1/2)1000 ≃ 9 · 10−302. Dennis verwendet wegen dergroßen Schrittzahl Formel (2.4) und erhalt das Ergebnis ≃ 1, 8 · 10−219,also eine extrem viel hohere Wahrscheinlichkeit. Wessen Ergebnis stimmt,und warum ist das andere Ergebnis falsch? Warum ist der Unterschied derbeiden Ergebnisse in der Praxis eigentlich nicht wichtig?

2.2.2 Kontinuumsbeschreibung

Wir wollen nun zu einer Kontinuumsbeschreibung ubergehen und die Schritt-lange ∆x sowie den zeitlichen Abstand ∆t immer kleiner machen. Wir definieren

x = m∆x und t = N∆t (2.5)

und fuhren die Parameter

D = 2pq(∆x)2

∆tund v = (p − q)

∆x

∆t(2.6)

ein. Dann konnen wir obige Gaußverteilung auf die folgende Weise schreiben:

p(m∆x, N∆t) =2∆x√4πDt

e−(x−vt)2

4Dt .

7

Wir machen den Grenzubergang ∆x → 0, ∆t → 0 bei konstantem D und kon-stantem v. Dies bedeutet, dass wir die beiden “makroskopischen” Charakteri-stika des Zufallspfades, die mittlere Driftgeschwindigkeit und die Diffusionskon-stante vorgeben, aber Raum und Zeit kontinuierlich machen. Bei diesem Prozessschrumpft ∆x mit der Wurzel von ∆t, und die Parameter p und q andern sichso, dass p − q ∝ ∆x ist. Beide Parameter nahern sich also dem Wert 1/2 an.Wir erhalten dann eine Wahrscheinlichkeitsdichte

P (x, t) =1√

4πDte−

(x−vt)2

4Dt . (2.7)

2.2.3 Die Diffusionsgleichung

Wir kennen Gaußverteilungen, deren Breite wie√

t anwachst, als Losung der Dif-fusionsgleichung. Wir werden spater im Zusammenhang mit der Fokker-Planck-Gleichung zeigen, dass sich der Random Walk in der Kontinuumsbeschreibungals Diffusionsprozess formulieren lasst.

Hier geben wir die Diffusionsgleichung, deren Losung (2.7) ist, nur an:

∂P (x, t)

∂t= −v

∂P (x, t)

∂x+ D

∂2P (x, t)

∂x2. (2.8)

Gleichung (2.8) wird von der Verteilung (2.7) gelost, wenn als Anfangsverteilungein Deltapeak am Ursprung gegeben ist. Die beiden Terme auf der rechten Seitevon Gleichung (2.8) haben eine einfache anschauliche Bedeutung: Der erste Term(mit v) verschiebt die Wahrscheinlichkeitsdichte P (x, t) mit Geschwindigkeit vnach rechts (wenn v positiv ist; sonst wird die Verteilung nach links geschoben).Man nennt diesen Term den Driftterm. Der zweite Term (mit D) macht dieVerteilung breiter. Er ist der Diffusionsterm.

2.2.4 Varianten des Random Walk, die dieselbe Kontinu-

umsbeschreibung haben

Die Kontinuumsbeschreibung ist eine gute Naherung, wenn die Positionsande-rungen und Zeitintervalle, die man betrachtet, sich aus sehr vielen einzelnenSchritten zusammensetzen. Dann kann man den zentralen Grenzwertsatz ver-wenden und die Gauß-verteilte Wahrscheinlichkeitsdichte (2.7) ableiten. DasErgebnis (2.7) enthalt keine Information mehr daruber, dass die Schritte desRandom Walk alle dieselbe Große haben und in einem festen Zeittakt statt-finden. In das Ergebnis geht nur die Driftgeschwindigkeit v ein, die bestimmt,mit welcher Geschwindigkeit sich die Gaußverteilung verschiebt, und die Dif-fusionskonstante D, die bestimmt, wie schnell die Breite der Gaußverteilunganwachst.

Varianten des Random Walk, in denen verschiedene Schrittgroßen moglichsind oder in denen Schritte nicht in einem festen Takt auftreten, sollten nach demZentralen Grenzwertsatz auch auf eine Gaußverteilung der Form (2.7) fuhren.Wir haben also die folgenden vier Varianten:

1. Der ursprungliche Random Walk, wie wir ihn bisher behandelt haben.

2. Variante, die kontinuierlich in der Zeit und diskret im Ort ist: Statt dassder Besoffene in einem festen Takt seine Schritte macht, macht er nun mit

8

einer Rate µr einen Schritt nach rechts und mit einer Rate µl einen Schrittnach links. Im infinitesimalen Zeitintervall dt macht der Besoffene also mitWahrscheinlichkeit µrdt einen Schritt der Große ∆x nach rechts und mitWahrscheinlichkeit µldt einen Schritt nach links. Also betragt die in derZeit dt zuruckgelegte mittlere Entfernung vdt = (µr − µl)∆xdt, und wirerhalten v = (µr −µl)∆x. Die Varianz der Position des Besoffenen nimmtin der Zeit dt um (µr + µl)(∆x)2dt zu (nachrechnen!). Also haben wir2D = (µr + µl)(∆x)2. Wenn viele Schritte zuruckgelegt wurden, gilt alsowieder die Gauß-Verteilung (2.7). Auch hier kann man den Grenzubergang∆x → 0 bei festem v und D durchfuhren. In diesem Fall mussen µr undµl sich andern gemaß den Gleichungen

µr =2D + v∆x

2(∆x)2

und

µl =2D − v∆x

2(∆x)2.

Beide mussen also ∝ (∆x)−2 sein, und ihre Differenz muss mit 1/∆xanwachsen. Das Verhaltnis der beiden Raten nahert sich also immer mehrdem Wert 1 an.

3. Variante, die kontinuierlich im Ort und diskret in der Zeit ist: Der Besoffe-ne macht seine Schritte in einem festen Takt, hat aber eine Schrittgroßen-verteilung g(∆x). Das Vorzeichen von ∆x gibt jetzt an, ob der Schritt nachlinks oder rechts geht. Die Wahrscheinlichkeit, dass der nachste Schritt ei-ne Große im Intervall [∆x, ∆x + dx] hat, betragt also g(∆x)dx. die Funk-tion g(∆x) muss auf 1 normiert sein, also

∫

g(∆x)d(∆x) = 1. Damit derzentrale Grenzwertsatz angewendet werden darf, ist es wichtig, dass dieSchrittgroßenverteilung eine endliche Varianz hat, d.h. dass das zweiteMoment

∫

g(∆x)(∆x)2d(∆x) existiert, also das Integral nicht divergiert.Dann erhalten wir

v∆t =

∫

g(∆x)(∆x)d(∆x) = 〈∆x〉

und2D∆t = 〈∆x2〉 − 〈∆x〉2.

Einen Kontinuumslimes ∆t → 0 durchzufuhren, geht hier nicht eindeutig,außer wenn die sich die Verteilung g(∆x) durch zwei Parameter charakte-risieren lasst und man vorgibt, dass die Verteilung beim Verkleinern von∆t ihre funktionale Form beibehalt, also dass sich nur die Werte der Pa-rameter andern. Besonders wichtig ist der Fall, dass g(∆x) eine Gaußver-teilung ist. Denn dann kann man auf der einen Seite den eben erwahntenUbergang zur Kontinuumsbeschreibung durchfuhren (siehe Ubungsaufga-be). Zum anderen ergibt sich fur g(∆x) immer dann eine Gauß-Verteilung,wenn sich das Zeitintervall ∆t in Wirklichkeit aus vielen kleinen Teilin-tervallen zusammensetzt, in denen jeweils Schritte getan werden. Dieseeinzelnen Schritte werden aber nicht aufgelost oder nicht registriert, daman immer nur in Zeitabstanden ∆t nachsieht, wo der Besoffene sich be-findet.

9

4. Variante, die kontinuierliche Werte fur die Schrittgroßen und fur die Zeitab-stande zwischen den Schritten hat: Diese Variante wird durch die Raten-dichte µ(∆x) spezifiziert. µ(∆x)d(∆x)dt ist die Wahrscheinlichkeit, dassder Besoffene in der infinitesimalen Zeit dt einen Schritt im Großeninter-vall [∆x, ∆x + d(∆x)] macht. Die Wahrscheinlichkeit, dass er in der Zeitdt uberhaupt einen Schritt macht, ist also dt

∫

µ(∆x)d(∆x). Im Gegensatzzu g(∆x) ist µ(∆x) nicht normiert. Wir erhalten in diesem Fall

v =

∫

µ(∆x)(∆x)d(∆x) (2.9)

und

2Ddt = dt

∫

(∆x)2µ(∆x)d(∆x) −(

dt

∫

∆xµ(∆x)d(∆x)

)2

= dt

∫

(∆x)2µ(∆x)d(∆x) ,

weil dt infinitesimal klein ist. Also bleibt

2D =

∫

(∆x)2µ(∆x)d(∆x) . (2.10)

Aufgaben

1. Zeige explizit, dass Gleichung (2.8) von der Verteilung (2.7) gelost wird.Hast Du eine Idee, wie man die Losung von Gleichung (2.8) bekommenkann, wenn die Anfangsverteilung keine δ-Funktion ist?

2. Fuhre fur die zweite Variante des Random Walk fur den Fall, dass g(∆x)eine Gaußverteilung ist, den Kontinuumsubergang explizit durch.

3. Random Walk in drei Dimensionen: Eine Anwendung hierfur ist die Wahr-scheinlichkeitsverteilung des End-zu-End-Abstands eines Polymers. Wirstellen ein Polymer idealisiert dar als Kette von Monomeren, die alle die-selbe Lange a haben, aber in beliebige Richtung orientiert sind. Wenn dasPolymer N Monomere hat, was ist der mittlere quadratische Abstand zwi-schen dem Anfangs- und dem Endpunkt? Und was ist die Wahrscheinlich-keitsverteilung des Abstands? Was findet Ihr an dem Modell unrealistisch?Welches verbesserte Modell wurdet Ihr vorschlagen? Erwartet Ihr fur dasverbesserte Modell ein ahnliches Ergebnis?

2.3 Master-Gleichungen

2.3.1 Raten und Wahrscheinlichkeiten

Master-Gleichungen beschreiben Systeme mit kontinuierlicher Zeit, aber diskre-ten Zustanden. Statt der Ubergangswahrscheinlichkeiten (pro Zeitschritt), dieden Elementen Qji der Matrix Q entsprechen, oder die fur den Random Walkdie Werte p und q haben, gibt es nun Ubergangsraten wji von Zustand i nachZustand j. Raten sind Wahrscheinlichkeiten pro Zeiteinheit. Durch Multiplika-tion einer Rate wji mit einem infinitesimal kleinen Zeitintervall erhalt man die

10

(infinitesimal kleine) Wahrscheinlichkeit dafur, dass in diesem Zeitintervall einUbergang vom Zustand i in den Zustand j passiert.

Als erstes Beispiel betrachten wir den radioaktiven Zerfall eines Isotops mitder Zerfallskonstanten λ. Jedes Atom zerfallt mit der Rate λ. Im Zeitintervalldt zerfallt es mit der Wahrscheinlichkeit λdt. Die Wahrscheinlichkeit, dass dasAtom wahrend der Zeit t = Ndt nicht zerfallt, betragt folglich (1−λdt)N = e−λt.Betrachtet man makroskopisch viele Atome, so folgt daraus, dass die Zahl dernicht zerfallenen Atome wie n(t) = n(0)e−λt abklingt. Dies ist aquivalent zurDifferenzialgleichung n = −λn. Wenn die Zahl der Atome nicht makroskopischgroß ist, konnen wir keine deterministische Gleichung fur die Zahl der Atomeals Funktion der Zeit schreiben, sondern wir haben es mit einem stochastischenProzess mit der Variablen n zu tun. Es gibt nur Ubergange von Werten n zumnachstkleineren Wert n− 1. Die Ubergangsrate hierfur ist λn. Wenn man einensolchen Prozess am Computer simulieren will, muss man darauf achten, dassman das Zeitintervall ∆t (das jetzt naturlich nicht infinitesimal klein sein kann)so klein wahlt, dass nicht nur λ∆t klein ist, sondern auch λn∆t. Dann kannman die Wahrscheinlichkeit dafur, dass in demselben Zeitintervall zwei odermehr Atome zerfallen, vernachlassigen gegenuber der Wahrscheinlichkeit, dassein Atom zerfallt. Die Wahrscheinlichkeit dafur, dass zwei Atome im Zeitintervall∆t zerfallen, betragt 0.5n(n − 1)λ2(∆t)2.

Als zweites Beispiel betrachten wir wieder den Random Walk. Statt dassder Besoffene in einem festen Takt seine Schritte macht, macht er nun mit einerRate µr einen Schritt nach rechts und mit einer Rate µl einen Schritt nach links,wie in der zweiten Variante. Die mittlere Zeit zwischen zwei Schritten nachrechts betragt 1/µr, und die mittlere Zeit zwischen zwei Schritten nach linksbetragt 1/µl. Die mittlere Zeit zwischen zwei aufeinanderfolgenden Schritten inbeliebiger Richtung betragt 1/(µr + µl). Die Wahrscheinlichkeit dafur, dass erin einem Zeitintervall t keinen Schritt macht, betragt e−(µr+µl)t.

2.3.2 Die Mastergleichung

Die zeitliche Entwicklung der Wahrscheinlichkeit, in einem bestimmten Zustandzu sein, ist durch die Master-Gleichung gegeben:

∂p(i, t)

∂t=

∑

j 6=i

[p(j, t)wij − p(i, t)wji] (2.11)

Die stationare Verteilung andert sich nicht mit der Zeit, und sie erfullt daherdie Gleichung

0 =∑

j 6=i

[peq(j)wij − peq(i)wji] . (2.12)

Wenn nicht nur die Summe, sondern jeder einzelne Summand verschwindet, alsowenn

0 = peq(j)wij − peq(i)wji , (2.13)

dann liegt detailliertes Gleichgewicht vor. Anschaulich bedeutet dies, dass Uber-gange von i nach j genauso haufig sind wie Ubergange von j nach i. Wennman also eine Filmaufnahme eines Systems im detaillierten Gleichgewicht machtund den Film dann ruckwarts laufen lasst, sieht er genauso realistisch aus wie

11

vorwarts. Stationare Systeme, die kein detailliertes Gleichgewicht haben, sehenunter Umkehr der Zeitrichtung anders aus.

Detailliertes Gleichgewicht liegt z.B. im mikrokanonischen Ensemble vor.Diese Systeme sind von ihrer Umwelt isoliert. Da die elementaren physikali-schen Prozesse invariant unter Zeitumkehr sind, erwartet man dies auch fur einSystem im Gleichgewicht, in dem der Anfangszustand “vergessen” ist. Auchkanonische und makrokanonische Ensembles haben detailliertes Gleichgewicht,da man man diese Systeme gemeinsam mit ihrem Warme- oder Teilchenbadbetrachten kann, und dann liegt wieder ein isoliertes System vor. Offene Sy-steme dagegen wechselwirken bestandig mit ihrer Umwelt und konnen daherpermanent fern von thermischen Gleichgewicht sein. Also kann man in ihnenkein detailliertes Gleichgewicht erwarten.

Mastergleichungen lassen sich meist nicht analytisch losen. Daher simuliertman sie ublicherweise am Computer, oder macht man Naherungen und gehtz.B. zur Fokker-Planck-Gleichung uber (s.u.).

2.3.3 Exkurs: Monte-Carlo-Simulationen zur Bestimmung

statistischer Eigenschaften

Das folgende Beispiel ist im Gegensatz zu den meisten anderen Beispielen diesesKapitels kein eindimensionales System. Die Zustande i und j seien Zustandeeines Systems im kanonischen Ensemble. Wir wollen durch eine Computer-simulation bestimmen, wie haufig welcher Zustand im kanonischen Ensemblevorkommt oder was der Wert einer Observablen im kanonischen Ensemble ist.Hierbei ist das Konzept von Ubergangsraten und Mastergleichungen sehr hilf-reich.

Das statistische Gewicht eines Zustands i im Ensemble ist e−E(i)/kBT /Z,wobei E die Energie des Zustands ist und Z die kanonische Zustandssumme.Die direkte Auswertung der Observablen durch Summation uber alle quanten-mechanischen Zustande des Systems ist normalerweise nicht durchfuhrbar, undman muss sich stattdessen mit einer reprasentativen Stichprobe begnugen. Aller-dings ist es wenig effizient, eine Stichprobe von Zustanden zufallig zu erzeugenund die Observable nur bezuglich dieser Zustande zu ermitteln, da die mei-sten zufallig erzeugten Zustande eine hohe Energie haben und daher zusammennur ein kleines statistisches Gewicht haben. Viel effizienter ist es, einen stocha-stischen Prozess zu definieren, der einen reprasentativen Satz von Zustandenerzeugt. Man nennt solche Computersimulationen Monte-Carlo-Simulationen.

Zu diesem Zweck wahlt man die Ubergangsraten zwischen den Zustandenso, dass detailliertes Gleichgewicht vorliegt, also dass

wij

wji=

peq(i)

peq(j)= e−(E(i)−E(j))/kBT (2.14)

ist. Damit ist gewahrleistet, dass in einer unendlich langen Zeitreihe jeder Zu-stand genau mit der Haufigkeit vorkommt, die seinem statischen Gewicht ent-spricht.

Zwei beliebte Moglichkeiten, die Ubergangsraten zu wahlen, sind die Glauber-Dynamik

wij = aij1

2

[

1 + tanh

(

−E(i) − E(j)

2kBT

)]

(2.15)

12

und die Metropolis-Dynamik

wij = aijmin(1, e−(E(i)−E(j))/kBT ) . (2.16)

Die Parameter aij = aji haben in der Praxis fur viele Zustandspaare den Wert0 und fur die anderen Zustandspaare einen konstanten Wert.

In einer Computersimulation beginnt man mit einem beliebigen Anfangszu-stand i0 und erzeugt dann die weiteren Zustande i1, i2 etc mit Ubergangsratenentsprechend einer der letzten beiden Formeln. Dabei muss man sicherstellen,dass die Zeitreihe so lang wird, dass die in ihr erzeugte Stichprobe auch re-prasentativ ist. Dabei ist es ratsam, den ersten Teil der Zeitreihe zu ignorieren,da der Anfangszustand wenig typisch fur das kanonische Ensemble ist, wenn ernicht mit viel Einsicht ausgewahlt wurde.

Als Beispiel wahlen wir das Ising-Modell auf dem zweidimensionalen Git-ter, E = −J

∑

α,β Nachbarnσασβ , wobei die Indizes Gitterplatze durchzahlen.

Die Summe geht hier uber nachste-Nachbarpaare, und die σα haben die Werte±1. Das Ising-Modell ist ein Modell fur den Magnetismus, da im Grundzustandalle Spins parallel ausgerichtet sind. Ein Zustand des Gesamtsystems ist durchdie Spinkonfiguration σα gegeben. Wir wahlen aij = 1, wenn sich die Kon-figurationen i und j nur um einen Spin unterscheiden, und sonst aij = 0. DieEnergieanderung des Gesamtsystems bei Flippen des Spins α betragt ∆E =2Jσα

∑

β σβ , wobei die Summe uber die 4 Nachbarspins des Spins α lauft.Eine Monte-Carlo-Simulation des Ising-Modells sieht dann so aus:

1. Erzeuge eine Anfangskonfiguration σα.

2. Wahle einen Spin α zufallig aus allen aus.

3. Berechne die Energieanderung ∆E, die mit dem Umkehren dieses Spinsverbunden ist.

4. Berechne das wij dafur, dass der Spin umgedreht wird, unter Verwendungder Glauber- oder Metropolis-Dynamik.

5. Erzeuge eine Zufallszahl zwischen 0 und 1.

6. Wenn die Zufallszahl kleiner ist als wij , dann drehe den Spin um.

7. Gehe zu Schritt 2.

Bei der Erstellung dieser Regeln haben wir es umgangen, Zeitintervalle zuwahlen, die so klein sind, dass die Wahrscheinlichkeit fur einen Ubergang wahrendeines Zeitintervalls klein ist. Stattdessen haben wir eine Regel gewahlt, die si-cherstellt, dass dann wenn der nachste Ubergang passiert (egal wann er pas-siert), er mit den richtigen relativen Wahrscheinlichkeiten in einen der benach-barten Zustande geht. Mit diesem Trick konnen wir die Simulation mit maxi-malem Tempo ablaufen lassen und erhalten trotzdem eine korrekte Abfolge vonZustanden zu unserer Master-Gleichung. Allerdings verzerren wir die Zeitskaladabei.

Eine alternative Methode, die kleinen Ubergangswahrscheinlichkeiten amComputer zu behandeln, ohne dabei die Zeitskala zu verzerren, lernen wir inden Ubungen kennen.

13

Aufgaben

1. Diskretisiere die Zeit in der Mastergleichung und zeige, dass sich dann dieZeitentwicklung einer Markovkette, also eine Zeitentwicklung der Form

~rt+1 = Q~rt

ergibt.

2. Betrachte den in der Zeit kontinuierlichen Random Walk, bei dem derPfad mit einer Rate µl bzw µr einen Schritt nach links bzw einen Schrittnach rechts geht. (a) Schreibe fur den in der Zeit kontinuierlichen RandomWalk die Mastergleichung auf. (b) Wie wurdest Du einen kontinuierlichenRandom Walk am Computer simulieren?

3. In einer Bakterienkolonie aus n Individuen entstehen mit einer Rate r(n−n2/K) neue Individuen, und sie sterben mit einer Rate µ pro Individuum.Interpretieren Sie die Parameter r und K. Wie lautet die entsprechendeMastergleichung?

2.4 Kramers-Moyal-Entwicklung und Fokker-Planck-

Gleichung

Ab jetzt betrachten wir wieder eindimensionale Prozesse. Die Zustande i und jersetzen wir daher durch die Variable x. Wir gehen weiterhin davon aus, dassdie Ubergange uberwiegend zu nahegelegenen Werten von x stattfinden.

Solche eindimensionalen Prozesse konnen genau wie der Random Walk je-weils in der Zeit und im Ort diskret oder kontinuierlich sein. Im Unterschiedzum Random Walk sind bei allgemeinen eindimensionalen Prozessen die Uber-gangswahrscheinlichkeiten (Variante 1), die Ubergangsraten (Variante 2), dieSchrittgroßenverteilung (Variante 3) oder die Ratendichten (Variante 4) vomOrt x abhangig.

Die Mastergleichung (2.11) wurde fur Variante 2 formuliert. Fur Variante 4sieht sie folgendermaßen aus:

∂P (x, t)

∂t=

∫

µ(x − ∆x; ∆x)P (x − ∆x, t)d(∆x) − P (x, t)

∫

µ(x; ∆x)d(∆x) .

(2.17)Da die Ratendichten µ vom Ort abhangen, ist als erstes Argument jeweils dieAusgangsposition des betreffenden Schrittes angegeben.

Wir konnen diese Gleichung auch anders interpretieren, wenn wir viele ein-zelne Schritte zu einem Schritt zusammenfassen. Dies bedeutet, dass wir diex-Achse aus so großer Entfernung betrachten, dass benachbarte Werte der Zu-fallsvariablen sehr dicht liegen, d.h. dass im Intervall dx, das weiter unten in denIntegralen steht, viele Werte der Variablen liegen. Auch auf der Zeitskala fuhrenwir eine solche Naherung durch. Die “infinitesimalen” Zeitintervalle dt seien alsoso groß, dass in ihnen viele Ubergange stattfinden, so dass wir µ(x; ∆x)dt alseine Schrittgroßenverteilung g(x; ∆x) interpretieren konnen.

Das Ziel dieses Abschnittes ist es, die Fokker-Planck-Gleichung aus der Ma-stergleichung herzuleiten. Die Fokker-Planck-Gleichung ist die verallgemeiner-te Variante der Diffusionsgleichung (2.8). Eine wichtige Voraussetzung fur die

14

Gultigkeit der Herleitung ist, dass die Ratendichte µ(x; ∆x) sich stetig mit xandert, und dass die Verteilung p(x, t) sich ebenfalls stetig mit x andert. Wenndiese Voraussetzungen erfullt sind konnen wir in den Integralen in Gleichung(2.17) eine Taylorentwicklung in ∆x machen, und erhalten

∂P (x, t)

∂t=

∞∑

n=1

(−1)n

n!

∂n

∂xn[an(x)P (x, t)] (2.18)

mit

an(x) =

∫ ∞

−∞

(∆x)nµ(x; ∆x)d∆x (2.19)

Dies ist die Kramers-Moyal-Entwicklung der Mastergleichung.Wenn die Anderung von µ und P mit x langsam ist (bezogen auf die typische

Schrittgroße), genugt es, nur die ersten beiden Terme dieser Entwicklung zuberucksichtigen, und wir erhalten die Fokker-Planck-Gleichung

∂P (x, t)

∂t= − ∂

∂x[a1(x)P (x, t)] +

1

2

∂2

∂x2[a2(x)P (x, t)] . (2.20)

Man nennt a1(x) den Drift-Koeffizienten. Er ist die mittlere Anderung der Va-riablen x pro Zeiteinheit. Wenn µ von x unabhangig ist, ist der stochastischeProzess ein Random Walk, und a1 ist identisch mit der Geschwindigkeit v. a2(x)ist der Diffusionskoeffizient. Er ist die mittlere quadratische Anderung von xpro Zeiteinheit. Wenn µ von x unabhangig ist, ist der stochastische Prozess einRandom Walk, und a2 ist identisch mit 2D, also mit zweimal der Diffusionskon-stanten.

Damit die Fokker-Planck-Gleichung eine eindeutige Losung hat, benotigtman die noch Anfangsbedingungen und die Randbedingungen.

Wenn der zugrunde liegende stochastische Prozess einer der anderen Varian-ten entspricht, kann man a1(x) und a2(x) naturlich direkt aus diesem Prozessausrechnen, unter Verwendung der Ubergangswahrscheinlichkeiten bzw -raten.a1dt ist die mittlere Anderung von x in der Zeit dt, und diese ist dieselbe beieiner diskreten oder einer kontinuierlichen Betrachtung des Prozesses. Wichtigist nur, dass auf der Lange a1dt sich der Wert von a1 praktisch nicht andert, sodass die Funktion a1(x) sauber definiert ist, und dies war ja die Voraussetzung,die wir zu Beginn dieses Teilkapitels gemacht haben. a2dt ist die Varianz derAnderung von x in der Zeit dt, und sie ist proportional zu dt unabhangig davon,ob wir auf Zeitskalen sind, auf denen die Diskretheit der x-Werte wahrnehmbarist oder nicht.

Die Vernachlassigung der Terme hoherer Ordnung ist umso mehr gerecht-fertigt, je langsamer P (x, t) und µ(x; ∆x) mit der Zufallsvariablen x variieren.Dann wird nicht nur die Fokker-Planck-Gleichung eine immer bessere Nahe-rung, sondern dann darf man auch immer mehr Einzelschritte beim Ubergangzum Kontinuum zusammenfassen. Nach dem zentralen Grenzwertsatz nahertsich dann die Schrittgroßenverteilung immer mehr einer Gaußverteilung.

Aufgaben

1. Der letzte Satz dieses Teilkapitels lautete “Nach dem zentralen Grenz-wertsatz nahert sich dann die Schrittgroßenverteilung immer mehr einer

15

Gaußverteilung.” Schreibe diese Gaußverteilung explizit hin unter Verwen-dung von a1(x) und a2(x). Wie hangen die an(x) von a1dt und a2dt ab?Was ergibt sich daraus im Grenzfall dt → 0?

2. Berechne a1 und a2 fur die letzte Aufgabe (Bakterienkolonie) des vorigenAufgabenblocks.

3. Berechne a1 und a2 fur den im Abschnitt 2.3.1 beschriebenen radioaktivenZerfall von n Atomen.

4. Leite die Fokker-Planck-Gleichung fur einen stochastischen Prozess derVariante 1 auf direktem Weg gemaß der folgenden Aleitung her. Der Pro-zess sei ein Random Walk mit ortsabhangigen Werten von p und q. DieAusgangsgleichung ist also

p(m, N + 1) = pm−1 p(m − 1, N) + qm+1 p(m + 1, N) . (2.21)

Fuhre auf beiden Seiten der Gleichung eine Taylorentwicklung durch, gehedann zu den Variablen x = m∆x und t = N∆t uber und drucke dieVorfaktoren auf der rechten Seite durch a1(x) und a2(x) aus. Fuhre denGrenzubergang ∆x → 0 und ∆t → 0 aus.

5. Wie verallgemeinert man die Fokker-Planck-Gleichung auf mehrdimensio-nale Zufallsvariablen?

2.4.1 Absorbierende Randbedingungen

Es gibt noch eine andere Variante der Fokker-Planck-Gleichung. Um diese Va-riante zu motivieren, gehen wir wieder zum Random Walk zuruck. In vielenBeispielen ist der Random Walk zuende, wenn ein bestimmter Wert uberschrit-ten wird. Wenn Peter oder Paul sein ganzes Vermogen verloren hat, konnendie beiden nicht weiter spielen, und wenn der Besoffene zuhause angekommenist, braucht er nicht mehr weiterzulaufen. Wir definieren daher die Wahrschein-lichkeit G(xf , x, t) dafur, dass der Random Walk zur Zeit t noch nicht bei xf

angekommen ist, wenn er zur Zeit 0 bei x gestartet ist. Fur einen allgemeineneindimensionalen stochastischen Prozess konnen wir eine entsprechende Wahr-scheinlichkeit definieren. Wir geben im Folgenden eine Master- und eine Fokker-Planck-Gleichung fur die Wahrscheinlichkeit G(xf , x, t) an. Wir benutzen zudiesem Zweck wieder Variante 4. Der entscheidende Trick besteht dabei darin,den ersten Zeitschritt und nicht den letzten zu betrachten, so dass man eineDiffusionsgleichung bzgl der Anfangsposition x erhalt. Die Zeitableitung vonG(xf , x, t) konnen wir dann schreiben als

∂G(xf , x, t)

∂t=

∫

µ(x; ∆x) (G(xf , x + ∆x, t) − G(xf , x, t)) d(∆x) . (2.22)

Die Taylorentwicklung in x und Abbrechen nach dem 2.Term ergibt

∂G(xf , x, t)

∂t≃ a1(x)

∂G(xf , x, t)

∂x+

a2(x)

2

∂2G(xf , x, t)

∂x2. (2.23)

Dies ist die Fokker-Planck-Gleichung bezuglich der Anfangsposition. a1 und a2

sind genauso definiert wie bisher. Man achte darauf, dass sie diesmal vor den

16

Ableitungen stehen und nicht dahinter. Der zweite Unterschied zur normalenFokker-Planck-Gleichung besteht darin, dass der Driftterm ein positives Vorzei-chen hat.

Wir kommen nun wieder zuruck zum Random Walk und setzten a1 = vund a2 = 2D. Die Losung von (2.23) hat diesmal nicht die Form (2.7), dadie Anfangs- und Randbedingungen ganz anders sind. Es gibt jetzt auch keineNormierungsbedingung fur G. Die Randbedingung lautet

G(xf , xf , t) = 0 (2.24)

und die Anfangsbedingung lautet

G(xf , x, 0) = 1 (2.25)

fur x 6= xf .Wir wahlen unsere x-Achse so, dass xf = 0 ist und dass der Startpunkt im In-

tervall x ≥ 0 liegt. Im Folgenden werden wir die Falle v > 0 und v < 0 zunachstseparat behandeln. Im Fall v > 0 besteht eine nichtverschwindende Wahrschein-lichkeit, dass der Random Walk nie beim Ursprung ankommt. Die Wahrschein-lichkeit hierfur, die “Entkommwahrscheinlichkeit” G(0, x,∞) bekommen wir aus(2.23) mit der Bedingung ∂G/∂t = 0. Die Losung zur richtigen Randbedingung(G = 0 fur x = 0 und G = 1 fur x = ∞) ist dann

G(0, x,∞) = 1 − e−vx/D . (2.26)

Im Fall v < 0 ist die mittlere Zeit, bis der Random Walk beim Ursprungankommt, endlich. Wir nennen diese Zeit τ (x). Es gilt

τ (x) = −∫ ∞

0

t∂

∂tG(xf , x, t) dt =

∫ ∞

0

G(xf , x, t) dt . (2.27)

Dies folgt daraus, dass der Random Walk mit Wahrscheinlichkeit

G(xf , x, t) − G(xf , x, t + dt) = −dt ∂G/∂t

im Zeitintervall [t, t+ dt] bei xf ankommt. Integration der Gleichung (2.23) vont = 0 bis t = ∞ ergibt dann

−1 = vdτ

dx+ D

d2τ

dx2. (2.28)

Die allgemeine Losung mit der Eigenschaft τ (0) = 0 lautet

τ(x) =x

|v| + A(e|v|x/D − 1) . (2.29)

mit einem noch zu bestimmenden Parameter A. Im Limes x → ∞ darf dieserAusdruck nicht exponentiell divergieren, sondern muss sich der Form τ (x) =x/|v| nahern, so dass wir A = 0 erhalten.

Die Wahrscheinlichkeitsdichte dafur, dass ein Random Walk, der bei x star-tet, zum Zeitpunkt t zum ersten Mal den Ursprung erreicht, betragt

f(t) =x√

4πDt3exp

[

− (x + vt)2

4Dt

]

. (2.30)

17

Fur v ≤ 0 ist diese Wahrscheinlichkeitsdichte auf 1 normiert, fur v > 0 istdas Integral

∫ ∞

0 f(t)dt durch (2.26) gegeben. Man nennt die Verteilung (2.30)die inverse Gaussverteilung, und sie wurde im Jahr 1915 von Schrodinger undSmoluchowski hergeleitet. Anschaulich kann man diese Formel folgendermaßenbegrunden: Die Wahrscheinlichkeitsdichte dafur, dass ein Random Walk, der beix startet, zur Zeit t am Ursprung ist, ist proportional zum Wert der zu diesemDiffusionsprozess gehorenden Gaußverteilung am Ursprung zur Zeit t, also zu

1√4πDt

exp

[

− (x + vt)2

4Dt

]

.

Je großer t ist, desto wahrscheinlicher ist es allerdings, dass der Random Walkzur Zeit t nicht zum ersten Mal beim Ursprung ist. Ein doppelt so lang dauernderWalk von x nach 0 uberschreitet im Durchschnitt den Ursprung doppelt so oft.Also ist die Wahrscheinlichkeit, dass er zur Zeit t zum ersten Mal bei 0 ist,proportional zu 1/t. Wenn man dies an die Gaußverteilung multipliziert unddiese dann so normiert, dass

∫ ∞

0f(t)dt = 1 ist (fur v < 0), erhalt man die

inverse Gaussverteilung.

Anwendung: Wahrscheinlichkeit fur die Fixierung einer vorteilhaften

Mutation

Wir betrachten eine biologische Population aus N Individuen und nehmen an,dass in einem Individuum eine Mutation stattgefunden hat, die seine Fitness einwenig erhoht. Intuitiv gehen viele Leute davon aus, dass sich solch eine vorteil-hafte Mutation im Laufe der Generationen in der gesamten Population durch-setzen wird. Wir zeigen im Folgenden, dass diese Intuition falsch ist und dassdie Wahrscheinlichkeit, dass die Mutation sich durchsetzt, recht klein ist. Deranschauliche Grund hierfur ist, dass durch zufallige Einflusse das Individuummit der Mutation oder seine Nachkommen evtl. keine Kinder bekommen, undsomit wurde die Mutation wieder aus der Population verschwinden. Die Zahlder Individuen, die die Mutation haben, ist also ein stochastischer Prozess. Wirdefinieren ihn folgendermaßen: Wir gehen davon aus, dass die Generationen alledieselbe Populationsstarke N haben, und dass sie nicht uberlappen. Es werdenalso in jedem Zeitschritt alle Individuen der vorigen Generation durch die Indi-viduen der nachsten Generation ersetzt. Der Einfachheit halber betrachten wirungeschlechtliche Fortpflanzung, so dass jedes Individuum nur einen Elternteil,die Mutter, hat. Die Mutation vererbt sich von Mutter zu Tochter. Die Tochter-generation bestimmen wir nach folgender Regel aus der Elterngeneration: Wirlosen Nmal aus, welches Individuum der Elterngeneration Mutter eines Indivi-duums der Tochtergeneration wird. Dabei haben Trager der Mutation eine umden Faktor (1+s) großere Wahrscheinlichkeit, gewahlt zu werden, als Individuenohne die Mutation. Unsere stochastische Variable x ist der Anteil der Indivi-duen mit der Mutation an der Gesamtpopulation. Wenn die PopulationsstarkeN groß ist, konnen wir x als kontinuierliche Variable interpretieren und eineFokker-Planck-Gleichung aufstellen. Wir benotigen hierzu die Terme a1(x) unda2(x). Der Erwartungswert der stochastischen Variablen in der nachsten Gene-ration ist

〈x(t + 1)〉 =x(t)(1 + s)

x(t)(1 + s) + (1 − x(t))=

x(t) + x(t)s

1 + x(t)s≃ x(t) + sx(t) − sx2(t) .

18

Im letzten Schritt haben wir benutzt, dass s klein ist. Es ergibt sich also a1(x) =x(1 − x)s. a2 nahert sich fur kleine s einem von s unabhangigen Wert. Wirberechnen zunachst die Varianz der Zahl der Trager der Mutation (mit s = 0),und teilen sie am Ende durch N2, um die Varianz des Anteils der Trager derMutation zu bekommen. Die Varianz der Zahl der Trager der Mutation ist Nmal die Varianz fur ein Kind, x ∗ (1 − x)2 + (1 − x) ∗ x2 = x(1 − x), worausa2 = x(1 − x)/N folgt.

Uns interessiert die Wahrscheinlichkeit u(x, t) = (1 − G(1, x, t = ∞)) dafur,dass wenn zur Zeit 0 der Anteil x aller Individuen der Population die Mutationhat, bis zur Zeit t = ∞ dieser Anteil den Wert xf = 1 erreicht hat. G(xf , x, t)genugt der Fokker-Planck-Gleichung bzgl. des Anfangswertes, Gleichung (2.23).Fur t = ∞ setzten wir die linke Seite dieser Gleichung zu Null. Einsetzen vona1 und a2 fuhrt dann zu

∂2u(x,∞)

∂x2= −2sN

∂u(x,∞)

∂x

was mit den Randbedingungen u(0, t) = 0 und u(1, t) = 1 auf die Losung

u(x,∞) =1 − e−2Nsx

1 − e−2Ns

fuhrt. Dies vereinfacht sich zu 2Nsx wenn Ns groß ist und x von der Großen-ordnung 1/N . Wenn es anfangs nur ein Individuum mit der Mutation gibt, wirddie Mutation mit der sehr kleinen Warscheinlichkeit 2s fixiert.

Aufgaben

1. Angenommen, der Besoffene hat einen Heimweg von einem KilometerLange. Wenn er um Mitternacht losgeht und alle 2 Sekunden einen Schrittder Große 1m macht – leider nur mit Wahrscheinlichkeit 1/2 jeweils in dierichtige Richtung–, mit welcher Wahrscheinlichkeit kommt er vor Sonnen-aufgang zuhause an?

2. Leite eine zu (2.27) ahnliche Beziehung fur die hoheren Momente τn derWahrscheinlichkeitsverteilung von τ her.

2.5 Langevin-Gleichungen

Den Ubergang zu kontinuierlicher Zeit und kontinuierlicher Zufallsvariable, denwir bei der Herleitung der Fokker-Planck-Gleichung vollzogen haben, kann manauch direkt in der Formulierung des stochastischen Prozesses machen.

Wir schreiben diesen Prozess in der folgenden Form:

x = v(x(t)) + b(x(t))η(t) (2.31)

In dieser Gleichung ist v(x) ein deterministischer Term, der einfach eine Teil-chengeschwindigkeit ist, wenn wir x als den Ort eines Teilchens interpretieren.Der Einfluss des Zufalls auf die Bewegung steckt in η. Wir wahlen v(x) so, dassder Erwartungswert 〈η(t)〉 des Rauschens verschwindet. In vielen Fallen wirdweißes Rauschen verwendet, und dann gilt

〈η(t)η(t′)〉 = δ(t − t′) . (2.32)

19

Hier fehlt noch eine Angabe uber die Wahrscheinlichkeitsverteilung von η(t).Vorher ist allerdings noch zu klaren, wie uberhaupt (2.31) zusammen mit (2.32)zu verstehen ist. Es macht ja zunachst keinen Sinn, auf der linken Seite x zuschreiben und folglich davon auszugehen, dass sich x kontinuierlich in der Zeitandert, und dann auf der rechten Seite einen Ausdruck zu schreiben, der sichin jedem noch so kleinen Zeitintervall unendlich oft andert. Wir multiplizierenzunachst (2.31) mit dt und erhalten

dx = v(x(t))dt + b(x(t))η(t)dt ≡ v(x(t))dt + b(x(t))dW . (2.33)

Jetzt ist η(t)dt = dW als Anderung von x im Zeitintervall [t, t + dt] durchdas Rauschen zu interpretieren. (Bemerkung: Die Regeln, die wir hier fur dasRechnen mit dW verwenden, sind der sogenannte Ito-Formalismus. Eine alter-native Methode ist der Stratonovich-Formalismus.) Fur dW setzen wir einenWiener-Prozess an, also die Wahrscheinlichkeitsverteilung (2.7) mit v = 0 undD = 1,

p(dW ) =1

√

2π(dt)e−(dW )2/2(dt) . (2.34)

Die Idee dahinter ist naturlich wieder die, dass in der Zeit dt eigentlich sehr vieleeinzelne Schritte des stochastischen Prozesses stattgefunden haben, so dass siesich nach dem zentralen Grenzwertsatz zu einer Gaußverteilung aufaddieren.

Fur nicht infinitesimale Zeiten t ergibt sich daraus

p(W ) =1√2πt

e−W 2/2t . (2.35)

Wenn man nur die Bewegungsgleichung x = η hat und bei x = 0 startet, ist dieVerteilung der Position x nach der Zeit t folglich eine Gaußverteilung der Breitet. Das Rauschen, das wir verwenden, nennt man Gaußsches weißes Rauschen.

Zu diesem stochastischen Prozess gibt es eine Fokker-Planck-Gleichung (2.20).Wir berechnen hierfur nun die Koeffizienten a1(x) und a2(x). Die mittlere Ande-rung der Variable x pro Zeiteinheit dt ist offensichtlich

a1(x) = v(x) .

Die Varianz von ∆x pro Zeit dt ist

a2(x) = (b(x))2〈dW 2〉/dt = (b(x))2 .

Die zu dem durch die Langevin-Gleichung beschriebenen stochastischen Pro-zess gehorende Sprunggroßenverteilung g(x; ∆x) ist

g(x; ∆x) =e− (∆x−vdt)2

2b2(dt)

√

2πb2(dt)(2.36)

Wenn wir die Zeitskala feiner wahlen, wird diese Verteilung immer schmaler,und im Limes dt → 0 verschwinden alle hoheren Momente 〈(∆x)n〉/dt. Diesist ein Vorteil der Kontinuumsdarstellung. Die Vernachlassigung von Termenhoherer Ordnung in der Kramers-Moyal-Entwicklung ist eine Naherung, wenndie Zufallsvariable diskret ist, und sie wird exakt, wenn die Zufallsvariable konti-nuierlich ist. Allerdings sollte man sich bewusst machen, dass die Kontinuums-formulierung ebenfalls eine Idealisierung ist, da ein physikalisch verursachtes

20

Rauschen naturlich immer eine endliche Korrelationszeit hat, so dass der “in-finitesimale” Zeitschritt dt eine nichtverschwindende Große hat und dass derTerm δ(t − t′) in (2.32) daher als δt,t′/dt interpretiert werden darf. Das Rau-schen verkorpert all diejenigen Freiheitsgrade und Einflusse, die nicht explizitin der Rechnung berucksichtigt werden und die daher “zufallig” sind. Um dasRauschen durch eine Delta-Funktion nahern zu durfen, muss die Dynamik dernicht berucksichtigten Einflusse schnell sein auf der betrachteten Zeitskala. Diestrifft zum Beispiel auf thermisches Rauschen zu, das durch die thermische Be-wegung der Umgebung entsteht. Fur ein Brownsches Teilchen kommt das ther-mische Rauschen von der Bewegung der Flussigkeits- oder Gasmolekule, die dasTeilchen umgeben. Nicht jedes Rauschen ist Gaußsches weißes Rauschen. Nurwenn das Rauschen aus der Uberlagerung genugend vieler elementarer Prozessemit wohldefinierter Varianz resultiert, ist es Gaußsch. Denn nur dann gilt derzentrale Grenzwertsatz. Wenn der elementare Prozess eine Verteilung mit diver-gierendem zweiten Moment hat, dann resultiert ein Rauschen, bei dem extremeEreignisse nicht vernachlassigbar sind.

Nachdem wir uns nun uberlegt haben, wie mit der Gleichung (2.31) zu rech-nen ist, betrachten wir noch einige Anwendungen.

2.5.1 Brownsche Bewegung

Wir beginnen mit der Brownschen Bewegung

mv(t) = −mζv(t) +√

λη(t) (2.37)

mit Gaußschen weißen Rauschen. Die formale Losung dieser Gleichung ist

v(t) = v0e−ζt + e−ζt

∫ t

0

dτeζτ

√λη(τ)

m.

Damit ist

〈v(t)v(t′)〉 = v20e−ζ(t+t′)+e−ζ(t+t′)

∫ t

0

dτ

∫ t′

0

dτ ′eζ(τ+τ ′) λδ(τ − τ ′)

m2≃ e−ζ|t−t′| λ

2ζm2.

Hierbei haben wir berucksichtigt, dass 〈η(τ)〉 = 0 ist und dass in dem Doppelin-tegral wegen der delta-Funktion die beiden Integrationsgrenzen t und t′ durchden kleineren dieser beiden Werte ersetzt werden durfen. Im letzten Schritt ha-ben wir außerdem angenommen, dass t, t′ ≫ ζ−1 ist. ζ−1 ist eine Relaxations-zeit, namlich die Zeitskala der Reibung, durch die die Geschwindigkeit verringertwird.

Mit derselben Naherung berechnen wir

〈x(t)2〉 =

∫ t

0

dτ

∫ t

0

dτ ′〈v(τ)v(τ ′)〉 ≃ λt

ξ2m2.

Die Fokker-Planck-Gleichung zur Brownschen Bewegung ist

∂P (v, t)

∂t=

∂

∂v(ζvP ) +

λ

2m2

∂2P

∂v2. (2.38)

Im Gleichgewicht wird die Geschwindigkeitsverteilung stationar, und wir ha-ben dann ∂P/∂t = 0 und

Peq(v) ∝ e−m2ζv2/λ (2.39)

21

Aus der statistischen Physik wissen wir, dass dies die Boltzmannverteilung∝ e−E/kBT sein muss (mit E = mv2/2), und daraus schießen wir, dass dieBeziehung

mζ

λ=

1

2kBT(2.40)

gelten muss. Wie zu erwarten, besteht ein Zusammenhang zwischen der Tem-peratur und der Rauschstarke λ. Solch einen Zusammenhang finden wir immerdann, wenn wir Systeme betrachten, die durch das Rauschen ins thermischeGleichgewicht gebracht werden. Wenn das Rauschen nicht thermischer Naturist, besteht dieser Zusammenhang nicht.

2.5.2 Brownsches Teilchen im Kraftfeld

Nun nehmen wir an, dass es eine ortsabhangige Kraft K(x) = −dV (x)/dx gibt,so dass die Bewegungsgleichung fur das Teilchen die Form

mx = −mζx + K(x) +√

λη(t) (2.41)

annimmt.Wir betrachten zunachst den Grenzfall starker Dampfung, in dem der Term

mx auf der linken Seite vernachlassigt werden kann. Dann vereinfacht sich dieBewegungsgleichung zu

x = −ΓdV

dx+

√

2ΓkBTη(t) (2.42)

mit Γ = 1/mζ. Die Starke des Rauschterms wird durch die Bedingung gegeben,dass die Gleichgewichtsverteilung Peq(x) ∝ e−V (x)/kBT ist. Die Fokker-Planck-Gleichung fur den Grenzfall starker Dampfung ist

∂P (x, t)

∂t= − ∂

∂x(ΓK(x)P ) + ΓkBT

∂2P

∂x2. (2.43)

Man nennt diese Gleichung Smoluchowski-Gleichung.Wenn man nicht den Grenzfall starker Dampfung betrachtet, hangt die

Fokker-Planck-Gleichung von den beiden Variablen x und v ab. Anhand vondiesem Beispiel zeigen wir, wie man eine Fokker-Planck-Gleichung fur einenzweidimensionalen Prozess herleiten kann. Zunachst schreiben wir (2.41) in derForm

v = −ζv +K(x)

m+

√

2ζkBT

mη(t)

x = v

und leiten dann die Fokker-Planck-Gleichung uber eine Taylorentwicklung derGleichung

∂P (x, v, t)

∂t=

∫ ∫

[

P (x − ∆x, v − ∆v, t)µ(x − ∆x, v − ∆v; ∆x, ∆v)

−P (x, v, t)µ(x, v; ∆x, ∆v)]

d(∆x)d(∆v)

in ∆x und ∆v her. Dies ergibt

∂P (x, v, t)

∂t= −v

∂P

∂x− K

m

∂P

∂v+ ζ

∂

∂v(vP ) +

ζkBT

m

∂2P

∂v2. (2.44)

22

Da der Rauschterm nur in der Gleichung fur v steht, enthalt der Diffusionstermnur Ableitungen bzgl. v. Man pruft leicht nach, dass die stationare Losung

Peq(x, v) ∝ e−(V (x)+mv2/2)/kBT (2.45)

ist.

2.5.3 Brownsche Bewegung im Doppelmuldenpotenzial: Kra-

mers Problem

Wir spezifizieren nun, dass das Potenzial V (x) ein Doppelmuldenpotenzial ist,mit einem weniger tiefen Minimum bei xA, einem tieferen Minimum bei xC >xA, und einem Maximum bei xB zwischen den beiden Minima. Wir beginnenmit einem Anfangszustand, in dem das Brownsche Teilchen im linken (hoheren)Minimum sitzt und fragen, was die Ubergangsrate in das tiefere Minimum beixC ist. Das Inverse der Ubergangsrate ist die durchschnittliche Zeit τ , die manwarten muss, bis der Ubergang passiert.

Kramers Anwendungsbeispiel war eine chemische Reaktion. Zustand A ent-spricht den Ausgangsprodukten und C den Endprodukten. Damit die Reaktionstattfindet, muss eine Energiebarriere uberwunden werden, die Aktivierungs-energie. Wir gehen davon aus, dass der Reaktionsverlauf durch eine eindimen-sionale Achse beschrieben werden kann. x ist also eine “Reaktionskoordinate”.Die Achse wird so skaliert, dass die durch das Potenzial verursachten Kraftedurch −V ′(x) gegeben sind. Weiterhin gehen wir davon aus, dass die Reaktionthermisch aktiviert ist, also dass die Reaktionspartner sich in einer thermalisier-ten Umgebung befinden und die Energie zum Uberwinden der Barriere durchthermische Stoße bekommen, die wir durch Gaußsches weißes Rauschen model-lieren konnen. Die Barriere ist so hoch, dass die Zeit τ viel großer ist als dieZeit, die das System braucht, um im Ausgangsminimum bei xA die Gleichge-wichtsverteilung zu erreichen, also V (xB) − V (xA) ≫ kBT .

Wir betrachten den Grenzfall starker Dampfung, so dass unser System durchdie beiden Gleichungen (2.42) und (2.43) beschrieben wird. Starke Dampfungbedeutet, dass die Reibungskraft x/Γ viel starker als die durch das Potenzi-al verursachte Kraft V ′′(x)(x − xA) ≡ mω2

A(x − xA) ist. Ohne Reibung waredie Geschwindigkeit von der Großenordnung aωA, wobei a die Amplitude derSchwingung um das Potenzialminimum ist. Die durch das Potenzial verursachteBeschleunigung ware also von der Großenordnung aω2. Dies ist vernachlassigbargegenuber der Beschleunigung durch Reibungskrafte aω/mΓ, wenn

mΓω ≪ 1 (2.46)

ist.Nun berechnen wir die mittlere Zeitdauer τ , die das Teilchen benotigt, um

von A nach C zu gelangen. Zur Beantwortung dieser Frage benotigen wir dieFokker-Planck-Gleichung bzgl der Anfangsposition, Gl. (2.23), wobei xf durchxC zu ersetzen ist.

Bevor wir die fur unser Problem zutreffenden Ausdrucke fur a1 und a2 ein-setzen, berechnen wir den Ausdruck fur τ zunachst allgemein. Ausgehend vonder Beziehung (2.27) erhalten wir eine zu (2.28) analoge Gleichung

−1 = a1(x)dτ

dx+

1

2a2(x)

d2 τ

dx2(2.47)

23

mit den Randbedingungen

τ (xC) = 0 , limx→−∞

d

dxτ (x) = 0 . (2.48)

Die Losung, die diese Randbedingungen erfullt, ist

τ(x) = 2

∫ xC

x

dx′φ−1(x′)

∫ x′

−∞

dx′′ φ(x′′)

a2(x′′)(2.49)

mit

φ(x) = exp

[∫ x

−∞

dx′ 2a1(x′)

a2(x′)

]

.

Wenn wir nun die Ausdrucke aus (2.43) a1 = −ΓV ′(x) und a2 = 2ΓkBT einset-zen, erhalten wir fur Kramers Problem

τ (x) =1

ΓkBT

∫ xC

x

dx′ exp

[

V (x′) − V (xA)

kBT

] ∫ x′

−∞

dx′′ exp

[

−V (x′′) − V (xA)

kBT

]

.

(2.50)Nun nutzen wir noch die Naherungen aus, die wir oben gemacht haben. DaV (xB) − V (xA) ≫ kBT ist, kommt der Hauptbeitrag zum ersten Integral ausder unmittelbaren Umgebung von xB. Wir entwickeln daher das Potenzial inder ersten Exponentialfunktion bis zur zweiten Ordnung in (x − xB),

V (x′) = V (xB) − 1

2mω2

B(x′ − xB)2 .

Im zweiten Integral kann man folglich die obere Integrationsgrenze gleich xB

setzen. Hier kommt der Hauptbeitrag zum Integral von der Umgebung des lo-kalen Minimums bei xA, so dass wir das Potenzial im zweiten Integral naherndurch

V (x′′) = V (xA) +1

2mω2

A(x′ − xA)2 .

Da die Integrale von der Umgebung der Extrema dominiert werden, konnen wiralle Integrationsgrenzen nach ±∞ schicken. Wir haben also nur noch Gauß-Integrale auszufuhren und erhalten

τ ≃ 2πm

ωAωBe(V (xB)−V (xA))/kBT . (2.51)

Die Ubergangsrate von xA nach xC ist das Inverse dieses Ausdrucks. Fur dieumgekehrte Ubergangsrate von xC nach xA erhalt man den entsprechendenAusdruck mit C statt A.

Dies ist Kramers beruhmtes Ergebnis fur die Reaktionsrate. Die Zeitskalawird von einem thermisch aktivierten Prozess bestimmt, wobei die Energie-barriere die Hohe V (xB) − V (xA) hat. Ubergangsraten, die exponenziell in derHohe der Barriere (geteilt durch kBT ) sind, findet man in vielen Anwendungen.Man nennt dies Arrhenius-Verhalten.

24

Kapitel 3

Nichtlineare Dynamik

Wir betrachten in diesem Kapitel dynamische Systeme, die durch Gleichungender Form

dx(1)

dt= F1(x

(1), x(2), . . . , x(N)) ,

dx(2)

dt= F2(x

(1), x(2), . . . , x(N)) ,

...

dx(N)

dt= FN (x(1), x(2), . . . , x(N)) . (3.1)

beschrieben werden. Wenn man den Anfangszustand ~x(t0) kennt, kann man mitHilfe dieser Gleichungen die zeitliche Entwicklung von ~x berechnen.

Wenn die Fi keine linearen Funktionen in den x(j) sind, liegt nichtlineareDynamik vor, die oft schwer vorherzusagen und sehr komplex sein kann. Einigeder stochastischen Systeme des letzten Kapitels werden zu deterministischennichtlinearen Systemen, wenn wir nur die Dynamik der Mittelwerte betrachten.Ein Beispiel hierfur ist das Wachstum einer Bakterienkolonie (siehe Aufgabe 3 inAbschnitt 2.3), das bei Vernachlassigung des Rauschens und bei Umbenennungder Parameter durch

n = rn(1 − n/K) (3.2)

beschrieben wird. In diesem Kapitel werden wir noch mehr Anwendungsbeispie-len aus der Populationsdynamik begegnen.

Auch Systeme, in denen Krafte und damit Beschleunigungen (also zweiteZeitableitungen) auftreten, lassen sich auf die Form (3.1) bringen. Hierzu wirddie Geschwindigkeit als weitere Variable eingefuhrt. Damit wird aus der Glei-chung d2x/dt2 = F (x)/m das Gleichungssystem

dx

dt= v

dv

dt= F (x)/m .

Wenn die Kraft zeitabhangig ist, F = F (x, t), fuhrt man eine dritte Variableein, um diese explizite Zeitabhangigkeit zu eliminieren und die Form (3.1) zu

25

erhalten:

ds

dt= 1

dx

dt= v

dv

dt= F (x, s)/m .

Der Raum, der von den Variablen des Systems aufgespannt wird, ist der Pha-

senraum. Die Bahn, die von dem Zustand des Systems im Phasenraum zuruck-gelegt wird, nennt man Trajektorie oder Orbit. Wir haben vorhin gesehen, dassdie zeitliche Entwicklung deterministisch und fur gegebenen Anfangszustandeindeutig ist. Daraus folgt unmittelbar, dass sich Trajektorien im Phasenraumnicht schneiden konnen.

Abb. 3.1 zeigt Phasenraumtrajektorien fur das ideale (d.h. reibungfsreie)Pendel, das durch die Gleichungen

dϕ

dt= v

dv

dt= −a sinϕ

(mit a = g/l) beschrieben ist.

-2 0 2ϕ

-3

-2

-1

0

1

2

3

v

Abbildung 3.1: Phasenraumtrajektorien des Pendels fur verschiedene Anfangs-geschwindigkeiten und fur a = 1.

Man unterscheidet zwischen konservativen und dissipativen dynamischen Sy-stemen. In konservativen Systemen kann man die Variablen so wahlen, dass dieGroße eines Phasenraumvolumenelements unter der zeitlichen Entwicklung kon-stant bleibt. Das Volumenelement andert nur seine Form. Hierzu gehoren dieHamiltonschen Systeme, die ublicherweise in der Klassischen Mechanik behan-delt werden. In dissipativen Systemen gibt es Gebiete des Phasenraums, in denenein Phasenraumvolumenelement mit der Zeit immer kleiner wird. Dieses Kapi-tel beschaftigt sich nur mit dissipativen Systemen. Ob ein Phasenraumvolumenunter der Dynamik gleich bleibt oder sich andert, kann man an dem Ausdruck∑

i ∂Fi/∂x(i) ablesen. Um dies zu zeigen, machen wir den Ansatz V = l1l2 . . . lN

26

(mit infinitesimalen li). Es ist l1 = x(1)e − x

(1)a ≃ l1∂F1/∂x(1). Damit erhalten

wir

V = l1l2 . . . lN + l2l1l3 . . . lN + . . .

=∂F1

∂x(1)l1l2 . . . lN +

∂F2

∂x(2)l2l1 . . . lN + . . .

=∑

i

∂Fi

∂x(i)l1l2 . . . lN

= V ~∇ · ~F . (3.3)

Wenn ∇ · ~F = 0 ist, bleibt das Phasenraumvolumenelement in seiner Großekonstant, wenn ∇ · ~F < 0 ist, schrumpft es, ansonsten wachst es.

Dissipative Systeme haben Attraktoren im Phasenraum. Dies sind Teilmen-gen des Phasenraums, denen sich alle Trajektorien annahern, die in einem be-stimmten endlichen Phasenraumvolumenelement starten. Attraktoren konnenPunkte sein, wie z.B. fur das gedampfte Pendel, das durch die Gleichungen

dϕ

dt= v

dv

dt= −a sinϕ− ηv

beschrieben wird. Der Fixpunkt v = 0, ϕ = 0 ist der Attraktor fur samtliche An-fangsbedingungen. Abb. 3.2 zeigt eine Phasenraumtrajektorie fur das gedampftePendel.

0 2 4 6 8ϕ

-2

-1

0

1

2

3

v

Abbildung 3.2: Phasenraumtrajektorie fur das gedampfte Pendel. v(t = 0) = 3,a = 1, η = 0.15.

Die Attraktoren konnen auch Grenzzyklen sein, wie z.B. in Abb. 3.3 fur dasperiodisch getriebene Pendel

dϕ

dt= v

dv

dt= −a sinϕ− ηv + c sin(ωs)

s = t mod 2π/ω .

27

-5 0 5 10 15ϕ

-2

-1

0

1

2

v

Abbildung 3.3: Zwei Phasenraumtrajektorien fur das gedampfte periodisch ge-triebene Pendel, projiziert in die ϕ − v-Ebene. a = 1, η = 0.15, ω = 1, c = 1.

Wenn die Dynamik eines dissipativen Systems chaotisch ist, konnen dieAttraktoren auch kompliziertere Struktur haben. Man nennt sie dann seltsa-

me Attraktoren. Ein Beispiel fur einen seltsamen Attraktor ist in Abbildung3.28 gezeigt. Auch unser getriebenes gedampftes Pendel hat seltsame Attrak-toren. Um komplexe Trajektorien oder Attraktoren besser zu erfassen, ist esoft zweckmaßig, einen Poincare-Schnitt zu machen. Man tragt dann nur dieSchnittpunkte der Trajektorie mit einer bestimmten (Hyper)-Ebene auf. Diedreidimensionale Trajektorie unseres getriebenen Pendels wird dann zu einerzweidimensionalen Ansammlung von Punkten, wie in Abb. 3.4 gezeigt. Die Tra-jektorie wiederholt sich nie genau. Wenn man den Poincare-Schnitt vergroßert,sieht man immer weitere Details des Attraktors (wenn die Rechengenauigkeitgenau genug war und die Laufzeit lang genug....). Mit dem Poincare-Schnitthat man aus einem dreidimensionalen dynamischen System eine zweidimensio-nale diskrete Abbildung gemacht, die jedem Punkt der Ebene einen Folgepunktzuordnet. Diskrete Abbildungen konnen also auch Chaos zeigen, ebenso wiekontinuierliche dynamische Systeme. Auf diskrete Abbildungen werden wir indieser Vorlesung allerdings nicht eingehen.

Bei Veranderung eines oder mehrerer Parameter eines nichtlinearen Systemskonnen Bifurkationen auftreten. Dies sind qualitative Anderungen der Attrak-toren und ihrer Stabilitat.

In diesem Kapitel werden zunachst mit eindimensionalen Systemen begin-nen und uns, zum Teil mit graphischen Methoden, mit Fixpunkten und ihrerStabilitat und den moglichen Bifurkationen befassen. Dann wenden wir uns Bei-spielen in zwei Dimensionen zu. Als wichtigstes neues Phanomen kommen hierGrenzzyklen hinzu. In drei Dimensionen schließlich kann es seltsame Attrakto-ren geben, fur die wir einige Beispiele sehen werden.

28

-2 0 2ϕ

-5

-4

-3

-2

-1

0

v

Abbildung 3.4: Poincare-Schnitt einer chaotische Trajektorie mit der s = 0–Ebene fur das gedampfte periodisch getriebene Pendel. v(t = 0) = 1.2, a = 1,η = 0.15, ω = 0.4, c = 1.5.

3.1 Nichtlineare Dynamik in einer Dimension

3.1.1 Phasenportraits

Wir behandeln zunachst Systeme, deren Dynamik in der Form

x = f(x) (3.4)

geschrieben werden kann. Uber ein solches System konnen wir eine Menge aus-sagen, ohne die Bewegungsgleichung explizit zu losen. Jedem Punkt der x-Achsekonnen wir einen Pfeil zuordnen, der das Vorzeichen von x angibt. Wenn f(x)positiv ist, zeigt der Pfeil nach rechts, wenn f(x) negativ ist, zeigt der Pfeilnach links. Wenn f(x) verschwindet, haben wir einen Fixpunkt. Ein Bild wieAbbildung 3.5, das alle qualitativ verschiedenen Trajektorien zeigt, nennt manein Phasenportrait.

3.1.2 Stabilitat von Fixpunkten

Um die Stabilitat eines Fixpunktes zu ermitteln, untersuchen wir, ob Trajekto-rien, die in unmittelbarer Umgebung des Fixpunktes starten, in ihn hinein odervon ihm weglaufen. Wir bezeichnen den Fixpunkt mit x∗ und den Startpunktder Trajektorie mit x∗ + u mit kleinem u. Taylorentwicklung von f(x) bis zurersten Ordnung gibt

du

dt= f ′(x∗)u . (3.5)

Wenn f ′(x∗) positiv ist, ist der Fixpunkt instabil, wenn f ′(x∗) negativ ist, ist erstabil. Wenn f ′(x∗) = 0 ist, konnen verschiedene Situationen vorliegen, die derLeser sich anhand der folgenden Beispiele veranschaulichen kann: (a) x = −x3,(b) x = x3, (c) x = x2, (d) x = 0.

29

x

f(x)

Abbildung 3.5: Beispiel fur ein Phasenportrait fur ein eindimensionales dynami-sches System. Stabile Fixpunkte werden mit gefullten Kreisen markiert, instabileFixpunkte mit leeren Kreisen.

Aufgaben

1. Welches Langzeitverhalten, abgesehen vom Hineinlaufen in einen Fixpunkt,kann ein eindimensionales dynamisches System noch zeigen?

2. Die logistische Gleichung

N = rN

(

1− N

K

)

(3.6)

wurde 1838 von Verhulst vorgeschlagen, um das Wachstum der Welt-bevolkerung zu beschreiben. Zeichnen Sie das Phasenportrait fur diesesModell. Diese Gleichung passt tatsachlich ganz gut fur einfache Organis-men wie Bakterien oder Hefe, wenn man sie im Labor unter konstantenBedingungen wachsen lasst. Diskutieren Sie, warum das Modell fur Insek-ten nicht gut ist.

3. Der Allee-Effekt: Fur viele Organismen ist die relative WachstumsrateN/N am großten, wenn N nicht zu klein ist. Dies kann zum Beispieldaher kommen, dass es fur kleine N schwierig ist, einen Partner zu finden,oder dass durch Inzucht wenig lebensfahige Individuen enstehen. ZeigenSie, dass dieser Effekt durch Gleichungen der Form

N

N= r − a(N − b)2 (3.7)

beschrieben wird. Welche Bedingungen mussen an a und b gestellt werden?Finden Sie alle Fixpunkte und diskutieren sie ihre Stabilitat. Skizzieren Sieden qualitativen Verlauf von N(t) fur verschiedene Anfangsbedingungen.Vergleichen Sie ihre Ergebnisse mit denen fur die logistische Gleichung.

4. Wissen Sie, wie man die Gleichung x = f(x) am Computer integriert?

5. Autokatalyse: Betrachten Sie die chemische Reaktion

A + X

k+−→k−←−

2X

30

Stellen Sie eine dynamische Gleichung fur die Konzentration x der X-Molekule auf. Bestimmen Sie die Fixpunkte und ihre Stabilitat. SkizzierenSie den zeitlichen Verlauf von x fur verschiedene Anfangsbedingungen.

3.1.3 Sattelknotenbifurkationen

Bei Veranderung der Parameter eines eindimensionalen dynamischen Systemskonnen Fixpunkte entstehen und vergehen, und ihre Stabilitat kann sich andern.Parameterwerte, an denen solche Bifurkationen auftreten, heißen Bifurkations-punkte.

Durch eine Sattelknotenbifurkation entstehen Fixpunkte bzw werden Fix-punkte zerstort. Wenn man einen Parameter verandert, bewegen sich zwei Fix-punkte aufeinander zu, kollidieren miteinander und vernichten einander. Wennman den Parameter in entgegengesetzter Richtung verandert, entstehen plotz-lich aus dem Nichts zwei Fixpunkte, die zunachst sehr nahe beieinander sindund dann auseinanderlaufen.

Das Standardbeispiel fur eine Sattelknotenbifurkation ist das System

x = r + x2 (3.8)

mit dem Parameter r, wie in Abbildung 3.6 gezeigt. Wenn r durch den Wert 0

f(x)

x

f(x)

x

f(x)

x

r=0r<0 r>0

Abbildung 3.6: Die Sattelknotenbifurkation. Bei r = 0 ist der Fixpunkt “halb-stabil”.

geht, geht das Minimum der Parabel uber die x-Achse, und die beiden Fixpunkteverschwinden. Allerdings ist der “Geist” dieser Fixpunkte noch da, wenn r kleinist, da x in der Umgebung des Ursprungs sehr klein wird. Das System verweiltalso lange in der Umgebung der verschwundenen Fixpunkte.

Das Bifurkationsdiagramm fur die Sattelknotenbifurkation ist in Abbildung3.7 gezeigt. Ihren Namen bekommt die Sattelknotenbifurkation daher, dass diebeiden Fixpunkte in zweidimensionalen dynamischen Systemen ein Sattel (alsoein Fixpunkt mit einer instabilen und einer stabilen Richtung) und ein Knoten(ein vollstandig stabiler Fixpunkt) sind.

Als weiteres Beispiel fur eine Sattelknotenbifurkation betrachten wir dasSystem

x = r − x− e−x . (3.9)

Hier kann man die Werte der Fixpunkte nicht sofort angeben. Durch eine gra-phische Methode kommt man am schnellsten zum Ziel. Die Fixpunkte sind dort,wo sich die Kurven y = r − x und y = e−x schneiden. (Der Leser zeichne diesselbst. Ich bin zu faul, dies fur das Skript zu zeichnen....). Die Richtung der Pfei-le im Phasenportrait zeigt nach rechts, wenn die Gerade r − x uber der Kurve

31

r

x

Abbildung 3.7: Bifurkationsdiagramm fur die Sattelknotenbifurkation. Der sta-bile Fixpunkt ist durch eine durchgezogene Linie markiert, der instabile durcheine gestrichelte.

e−x liegt, und sonst nach links. Daraus ermittelt man sofort die Stabilitat derFixpunkte.

Die Bifurkation ist dort, wo die Kurve r− x gerade die Tangente an e−x ist,und das ist offensichtlich bei r = 1 der Fall. Der halbstabile Fixpunkt ist dannbei x = 0. Wenn wir die Gleichung (3.9) fur r = 1 + ǫ (mit kleinem ǫ) und furkleine x jeweils bis zur fuhrenden nichtverschwindenden Ordnung entwickeln,erhalten wir

x ≃ ǫ− x2

2. (3.10)

Bis auf das Vorzeichen von x2 ist dies identisch mit unserem Standardbeispiel(3.8), wenn wir die Variable y = x/2 einfuhren und r = ǫ/2 setzen. In derUmgebung einer Sattelknotenbifurkation erhalt man in der Tat immer ein dy-namisches System der Form r± x2. Der Grund hierfur ist, dass zwei Fixpunktedann entstehen bzw vergehen, wenn sich ein lokales Minimum oder Maximumahnlich wie in Abbildung 3.6 dargestellt uber die x-Achse schiebt. In der Um-gebung dieses Minimums oder Maximums ergibt eine Taylorentwicklung abergenau die genannte quadratische Form. Man nennt sie daher auch die Normal-form. Den Vorfaktor vor dem x2 kann man durch eine Umskalierung von x zumVerschwinden bringen.

3.1.4 Die transkritische Bifurkation

Transkritische Bifurkationen treten in Systemen auf, in denen ein Fixpunktfur alle Parameterwerte existiert und nicht zerstort werden kann. Dies ist zumBeispiel in der Populationsdynamik der Fall, wo N = 0 ein Fixpunkt sein muss.Solch ein Fixpunkt kann aber seine Stabilitat andern, und dies passiert bei einertranskritischen Bifurkation. Die Normalform lautet

x = rx − x2 . (3.11)

Wenn r sich andert, sieht das Phasenportrait so aus, wie in Abbildung 3.8gezeigt. Das r − x-Diagramm ist in Abbildung 3.9 gezeigt. Fur eine biologischePopulation macht naturlich nur der Bereich x ≥ 0 Sinn. Aber es kann trotzdemhilfreich sein, auch den Bereich x < 0 aufzutragen, um besser zu sehen, was ander Bifurkation mit den Fixpunkten passiert.

32

x x x

f(x) f(x) f(x)

r<0 r=0 r>0

Abbildung 3.8: Die transkritische Bifurkation. Bei r = 0 ist der Fixpunkt “halb-stabil”.

r

x

Abbildung 3.9: Bifurkationsdiagramm fur die transkritische Bifurkation. Derstabile Fixpunkt ist durch eine durchgezogene Linie markiert, der instabile durcheine gestrichelte.

Als weiteres Beispiel betrachten wir ein einfaches Lasermodell, das 1983 vonH. Haken eingefuhrt wurde. Die dynamische Variable ist die Zahl n(t) von Pho-tonen in dem Laser. Sie entwickelt sich gemaß

n = GnN − kn .

N ist die Zahl der angeregten Atome im Laser. Der erste Term beschreibt also dieZunahme von Photonen durch Wechselwirkung mit den angeregten Atomen, derzweite Term das Entkommen eines Teils der Photonen am Laserende. Es fehltnoch eine Gleichung, die N durch n ausdruckt. Durch Pumpen des Lasers kanneine maximale Zahl N0 von angeregten Atomen erreicht werden, und ihre Zahlverringert sich durch die stimulierte Emission von Photonen. Also setzen wir an

N(t) = N0 − αn

und erhalten damitn = (GN0 − k)n− (αG)n2 . (3.12)

Daraus erhalten wir die Normalform fur eine transkritische Bifurkation, wennwir x = αGn und r = GN0−k setzen. Der Fixpunkt n∗ = 0 wird instabil, wennN0 die Schwelle k/G uberschreitet. Es entsteht ein stabiler Fixpunkt n∗ > 0,und der Laser strahlt Licht ab.

3.1.5 Die Heugabelbifurkation

Eine Heugabelbifurkation tritt in Systemen auf, in denen f(x) eine Symmetriebesitzt, so dass Fixpunkte in symmetrischen Paaren auftreten und verschwin-

33

den. Es gibt zwei Arten von Heugabelbifurkationen, die superkritische und diesubkritische.

Die superkritische Heugabelbifurkation hat die Normalform

x = rx − x3 . (3.13)

Hier erkennen wir die Symmetrie bezuglich des Vorzeichens von x.Wenn r sich andert, sieht das Phasenportrait so aus, wie in Abbildung 3.10

gezeigt. Das r − x-Diagramm ist in Abbildung 3.11 gezeigt.

x x x

f(x) f(x) f(x)

r<0 r=0 r>0

Abbildung 3.10: Die superkritsche Heugabelbifurkation.

r

x

Abbildung 3.11: Bifurkationsdiagramm fur die superkritsche Heugabelbifurka-tion.

Als weiteres Beispiel betrachten wir das System

x = −x + β tanhx , (3.14)

das in magnetischen Systemen und neuronalen Netzen auftritt. Ahnlich wie beimBeispiel fur die Sattelknotenbifurkation kommen wir hier am schnellsten zumZiel, wenn wir die rechte Seite in zwei Terme zerlegen, deren Form uns vertrautist. Die Fixpunkte sind die Schnittpunkte der Kurven y = x und y = β tanhx.Wenn β < 1 ist, gibt es nur den Schnittpunkte am Ursprung, wenn β > 1 ist,gibt es zwei weitere Schnittpunkte. (Auch bei diesem Beispiel uberlasse ich esdem Leser, die Bilder zu malen.) Durch Entwicklung um β = 1 und x = 0erhalten wir x ≃ ǫx − x3/3, was bis auf die Skalierung von x der Normalformentspricht.

Die subkritische Heugabelbifurkation hat die Normalform

x = rx + x3 . (3.15)

34

r

x

Abbildung 3.12: Bifurkationsdiagramm fur die subkritsche Heugabelbifurkation.

Jetzt sind die Stabilitaten der Fixpunkte umgekehrt, wie man im Bifurkations-diagramm sieht (Abb. 3.12).

Fur r > 0 gibt es in diesem Modell uberhaupt keinen Fixpunkt mehr. DieVariable x lauft nach ±∞. In einem echten System passiert dies naturlich nicht.Wenn wir uns in Erinnerung rufen, dass die Normalform aus einer Taylorent-wicklung in der Umgebung der Bifurkation resultiert, dann ist es naheliegend,den Term der nachsten Ordnung mitzunehmen. Um die Symmetrie zu wahren,nehmen wir den Term x5 und erhalten dann die erweiterte Normalform

x = rx + x3 − x5 . (3.16)

Durch geeignete Skalierung von x und t kann man immer erreichen, dass dieKoeffizienten der beiden nichtlinearen Terme 1 sind. Das Bifurkationsdiagrammsieht jetzt aus wie in Abbildung 3.13 gezeigt. Zusatzlich zur Heugabelbifurkati-

r

x

Abbildung 3.13: Bifurkationsdiagramm fur die subkritsche Heugabelbifurkationmit dem Term funfter Ordnung.

on sieht man auch zwei Sattelknotenbifurkationen, durch die die Fixpunkte zugroßen Werten von ±x entstehen.

Fur einen endlichen Bereich von r-Werten gibt es drei stabile Fixpunkte(zwei mit x ≥ 0). Wenn wir den Parameter r hoch- und runterfahren, kannalso Hysterese auftreten. Wir beginnen mit einem Wert von r, der so klein ist,dass es nur den Fixpunkt x∗ = 0 gibt. Wenn r den Wert 0 uberschreitet, wirddieser Fixpunkt instabil, und das System lauft auf einen von Null verschiedenen

35

Fixpunkt. Wenn wir jetzt wieder r runterfahren, wird irgendwann der Punkterreicht, wo dieser Fixpunkt instabil wird, und das System springt wieder aufx∗ = 0. Zwischen den beiden Sprungpunkten befindet sich das System mal aufdem Fixpunkt x∗ = 0 und mal auf dem Fixpunkt x∗ 6= 0, je nachdem, welcheTrajektorie der Parameter r zuletzt durchlaufen hat.

Wenn ein stabiler Fixpunkt bei Anderung der Parameter instabil wird unddas System auf einen anderen, weit davon entfernten Fixpunkt springt, ist daseine Katastrophe. Bei der superkritischen Heugabelbifurkation gibt es keine Ka-tastrophe, da der neue stabile Fixpunkt zunachst in unmittelbarer Nahe desinstabil gewordenen liegt. Bei der subkritischen Heugabelbifurkation ist diesvollig anders.

3.1.6 Bifurkationen mit mehr als einem Parameter

Wenn das System mehr als einen Parameter hat, treten auch mehr Bifurkatio-nen auf, und es kann leicht Katastrophen geben. Als erstes Beispiel betrachtenwir eine “unvollkommene” Heugabelbifurkation, wie sie auftritt, wenn die Sym-metrie bzgl x nicht perfekt ist. Dann haben wir eine Gleichung der Form

x = h + rx− x3 . (3.17)

Fur h 6= 0 ist die Symmetrie bzgl x gebrochen. Um herauszufinden, wann es wieviele Fixpunkte gibt, zeichnen wir die beiden Kurven y = −h und y = rx − x3

und ermitteln, wie oft sie sich schneiden.Aus Abbildung 3.14 kann man ablesen, dass es genau dann drei Fixpunkte

gibt, wenn r > 0 und |h| <√

8r3/27 ist. Sonst gibt es einen Fixpunkt. Die

x

r>0

x

r<0

rx−x3

−h

Abbildung 3.14: Die unvollkommene Heugabelbifurkation.

Bereiche im Parameterraum, in denen es einen bzw drei Fixpunkte gibt, sindim Stabilitatsdiagramm gezeigt, Abbildung 3.15. Das Bifurkationsdiagramm istin Abbildung 3.16 gezeigt. Man kann auch den Fixpunkt gegen h auftragen, wiein Abbildung 3.17. Aus den beiden letzten Abbildungen kann man erkennen,dass sowohl beim Variieren von h als auch beim Variieren von r Katastrophenauftreten konnen.

Als weiteres Beispiel betrachten wir ein Modell fur den Ausbruch einer In-sektenplage, des “spruce budworm” in Kanada. Die Zeitentwicklung der Insek-tenpopulation wird durch die Gleichung

N = RN

(

1− N

K

)

− BN2

A2 + N2. (3.18)

beschrieben. Der erste Term beschreibt logistisches Wachstum der Population,die die Nadeln des “balsam fir tree” fressen. Der zweite Term beschreibt das

36

r

h

0

0

1 Fixpunkt

3 Fixpunkte

Abbildung 3.15: Stabilitatsdiagramm fur die unvollkommene Heugabelbifurka-tion.

r

x x

r

h=0 h>0

Abbildung 3.16: Bifurkationsdiagramm fur die unvollkommene Heugabelbifur-kation.

Sterben von Insekten. Sie werden von Vogeln gefressen. Erst wenn die Insekten-population großer als A wird, wenden sich die Vogel ernsthaft diesen Insektenzu, und dann wird der zweite Term groß. Der Parameter B ist proportional zurZahl der Vogel.

Wir machen die Gleichung dimensionslos durch die Wahl x = N/A, τ =Bt/A, r = RA/B, k = K/A und erhalten dann

dx

dτ= rx

(

1− x

k

)

− x2

1 + x2. (3.19)

Wir fuhren wieder nur eine graphische Analyse durch. Die Fixpunkte undihre Stabilitat konnen wir ablesen, wenn wir die Kurven y = r(1 − x/k) undy = x/(1+x2) auftragen, wie in Abbildung 3.18. Wenn k klein ist, gibt es fur je-den (positiven) Wert von r nur einen Schnittpunkt der beiden Kurven, und zwarfur recht kleine Werte von x. Die Insektenpopulation ist also unter Kontrolle.Fur großere k gibt es mit wachsendem r zunachst eine Sattelknotenbifurkation,die zwei weitere Fixpunkte erzeugt, und dann eine Sattelknotenbifurkation, beider der ursprungliche stabile Fixpunkt verschwindet und nur der Fixpunkt mitgroßem x∗ bleibt. Dieser zweite Ubergang fuhrt im doppelten Sinn zu einer Ka-tastrophe: die Insektenplage bricht plotzlich aus, da die Population von einemkleinen auf einen großen Wert ansteigt. Die Baume werden kahlgefressen. Da-durch werden die Parameter r und k irgendwann wieder klein, und die Plage

37

xx

h

r>0

h

r<0

Abbildung 3.17: Bifurkationsdiagramm bei festem r.

x

r

k

1+xx

2

Abbildung 3.18: Die beiden Terme der Insektendynamik.

ist zu Ende – bis der Wald wieder so dicht ist, dass die nachste Katastrophepassiert...

Das Stabilitatsdiagramm fur dieses Modell ist in Abbildung 3.19 skizziert.Eine genauere, quantitative Untersuchung wird dem Leser uberlassen.

Aufgaben

1. Es hilft oft der Anschauung, das dynamische System in der Form x =−dV (x)/dx zu schreiben. V (x) kann man als ein “Potenzial” auffassen, indem ein “Teilchen” mit der Position x sich bergab bewegt. Die Minimavon V entsprechen den stabilen Fixpunkten, die Maxima den instabilenFixpunkten. Zeichne solch ein Potenzial fur eine superkritische und einesubkritische Heugabelbifurkation.

2. Zeigen Sie, dass die Zeitdauer, wahrend der nach einer Sattelknotenbifur-kation eine Trajektorie in der Umgebung des “Geists” gefangen bleibt, wie1/√

ǫ skaliert.

3. Ein einfaches Modell fur die Zahl von Fischen in einem Fischteich ist

N = rN(1 − N

K)−H

N

A + N

mit H > 0 und A > 0. Der erste Summand beschreibt das Wachstum der

38

r

k

Plage

bistabil

harmlos

Abbildung 3.19: Stabilitatsdiagramm fur das “spruce budworm” Modell.

Population, wenn nicht gefischt wird. Der zweite Summand beschreibt denEffekt des Fischens auf die Population.

(a) Gib eine biologische Interpretation des Parameters A.

(b) Zeige, dass man die Gleichung auf die Form

dx

dτ= x(1− x) − h

x

a + x

bringen kann.

(c) Zeige, dass das System einen, zwei oder drei Fixpunkte hat, je nachdem Wert von a und h. Untersuche die Stabilitat all dieser Fixpunkte.

(d) Untersuche die Dynamik nahe x = 0 und zeige, dass fur h = a eineBifurkation auftritt. Was ist es fur eine Bifurkation?

(e) Zeige, dass eine weitere Bifurkation fur h = (a+1)2/4 auftritt, wenna < ac ist. Bestimme ac und die Art der Bifurkation.

(f) Zeichne das Stabilitatsdiagramm in der a-h-Ebene. Kann es irgendwoHysterese geben?

4. Ein biochemischer Schalter: Zebrastreifen und Muster auf Schmetterlings-flugeln sind Beispiele fur biologische Musterbildung. Eine wichtige Zutatzu Modellen fur solche Musterbildung ist ein biochemischer Schalter, indem ein Gen G durch ein Signalmolekul S aktiviert wird. Wenn die Kon-zentration von S eine Schwelle ubersteigt, wird das Gen eingeschaltet undein Pigment wird produziert. Die Konzentration g des Pigments andertsich mit der Zeit gemaß

g = k1s0 − k2g +k3g

2

k24 + g2

(3.20)

mit positiven Konstanten ki. Hier ist s0 die Konzentration von S (die sichauf der betrachteten Zeitskala nicht andert), der zweite Term beschreibtdie Degradation des Pigments, der dritte eine Feedbackschleife, bei derdie Produktion des Pigments sich selbst nichtlinear verstarkt, bis zu einerSattigung.

39

(a) Skalieren Sie die Variablen so um, dass sich die Gleichung zu

dx

dτ= s− rx +

x2

1 + x2(3.21)

vereinfacht.

(b) Zeigen Sie, dass es fur s = 0 zwei positive Fixpunkte gibt, wennr < rc ist. Was ist der Wert von rc?.

(c) Der Anfangswert sei g(0) = 0, und s wird nun langsam von 0 hochge-fahren. Was passiert mit g(t)? Geht g auf Null zuruck, wenn s wiederauf Null zuruckgeht?

(d) Welche Bifurkationen konnen auftreten?

(e) Skizzieren Sie das Stabilitatsdiagramm.

5. Synchronisation von Leuchtkafern: In manchen Teilen Sudostasiens ver-sammeln sich tausende mannlicher Leuchtkafer nachts in Baumen undblinken synchron, um die Weibchen anzulocken. Wir modellieren die Re-aktion eines Leuchtkafers auf ein periodisches Leuchtsignal durch die fol-gende Gleichung:

θ = ω + A sin(θ −Θ) . (3.22)

ω ist die Eigenfrequenz des Kafers und Θ = Ωt das periodische Signal. θist die Phase, in der sich der Leuchtzyklus des Kafers gerade befindet.

(a) Wahlen Sie die Phasendifferenz θ−Θ als Variable und schreiben Siefur diese eine dimensionslose Gleichung auf.

(b) Zeichnen Sie das Phasenportrait fur verschiedene Werte des Parame-ters µ = (Ω − ω)/A und ermitteln Sie diejenigen Werte von µ, furdie der Leuchtkafer eine feste Phasenbeziehung mit dem Signal hat.Interpretieren Sie die Ergebnisse anschaulich.

3.2 Nichtlineare Dynamik in zwei Dimensionen

In zwei Dimensionen ist die Vielfalt an Fixpunkten und Bifurkationen großerals in einer Dimension. Außerdem gibt es Grenzzyklen.

3.2.1 Fixpunkte

Wir betrachten das Gleichungssystem d~x/dt = ~F (~x). und untersuchen das Ver-halten einer Trajektorie, die in der Nahe eines Fixpunktes ~x∗ startet. Dannist

~x(t) = ~x∗ + ~u(t) ,

und

dxi/dt = dui/dt = Fi(~x∗ + ~u(t)) ≃

∑

j

(

∂Fi

∂xj

)

~x∗

uj(t) .

Dies ist eine Matrixgleichung der Form

~u = A~u (3.23)

40

mit

Aij =

(

∂Fi

∂xj

)

~x∗

.

Diese Gleichung hat Losungen der Form

~u(t) =

N∑

k=1

~vkeλkt ,

wobei ~vk die Losung der Eigenwertgleichung

A~vk = λk~vk

ist.Da die Matrixelemente von A reell sind, gibt es zu jedem komplexen Eigen-

wert einen komplex konjugierten Partner.In zwei Dimensionen haben wir zwei Eigenwerte

λ1,2 =1

2

(

τ ±√

τ2 − 4∆)

, ∆ = λ1λ2, τ = λ1 + λ2 .

τ ist die Spur von A und ∆ die Determinante. Wir unterscheiden die folgendenFalle (siehe Abb. 3.20):

• ∆ > 0 und τ2 − 4∆ > 0: Beide Eigenwerte sind reell und haben dasselbeVorzeichen. Der Fixpunkt ist ein Knoten. Seine Stabilitat wird durch dasVorzeichen von τ bestimmt. Bei einem stabilen Knoten schmiegen sich dieTrajektorien mit der Zeit immer starker an die “langsame” Eigenrichtung(d.h. die mit dem kleineren |λ|) an, bei einem instabilen Knoten richtetsich die Trajektorie mit der Zeit an der “schnellen” Eigenrichtung aus.

• ∆ < 0: Die Eigenwerte sind reell und haben entgegengesetztes Vorzeichen.Der Fixpunkt ist ein Sattelpunkt. Es gibt eine stabile und eine instabileRichtung.

• ∆ > 0 und τ2 − 4∆ < 0: Die Eigenwerte sind komplex konjugiert. DerFixpunkt ist eine Spirale. Seine Stabilitat wird durch das Vorzeichen vonτ bestimmt. Wir setzen λ1 = λ′ + iλ′′ und λ2 = λ′ − iλ′′ an. Dann ist

~u(t) = ~v1e(λ′+iλ′′)t + ~v2e

(λ′−iλ′′)t .

~u ist reell, und folglich sind ~v1 und ~v2 komplex konjugiert, ~v1 = ~v′ + i~v′′

und ~v2 = ~v′ − i~v′′. Dies gibt

~u(t) = 2eλ′t (~v′ cosλ′′t− ~v′′ sin λ′′t) .

Die relative Orientierung von ~v′ und ~v′′ und das Vorzeichen von λ′′ be-stimmen den Drehsinn der Spirale.

• τ = 0: Der Fixpunkt ist ein Zentrum. Die Trajektorien laufen in geschlos-senen Bahnen um ihn herum. Zentrum und Sattelpunkt sind die einzigmoglichen Fixpunkte in Hamiltonschen Systemen, in denen die Energieerhalten ist. In dissipativen Systemen kann eine kleine Nichtlinearitat dasZentrum zerstoren, und es entsteht eine Spirale.

41

• τ2−4∆ = 0: In diesem Fall sind entweder alle Richtungen Eigenrichtungenzum selben Eigenwert und wir haben einen Stern, oder es gibt nur eineEigenrichtung, und dann ist der Fixpunkt entartet. In beiden Fallen be-findet sich der Fixpunkt an der Grenze zwischen Knoten und Spirale, undeine kleine Parameteranderung oder eine Nichtlinearitat fuhrt zu einemKnoten oder einer Spirale.

• ∆ = 0: In diesem Fall ist ein Eigenwert 0. Aus der entsprechenden Eigen-richtung lauft man weder aus dem Fixpunkt heraus, noch in ihn hinein.Es gibt also eine ganze Linie von Fixpunkten. Eine kleine Storung machtaus der Fixpunktlinie einen Sattel oder Knoten.

stabile Spirale

stabiler Knoten

Sattelpunkt

instabiler Knoten Stern

entarteter FixpunktFixpunktlinie

instabile SpiraleZentrum

Abbildung 3.20: Die erwahnten Arten von Fixpunkten in 2 Dimensionen.

3.2.2 Bifurkationen in zwei Dimensionen

Sattelknotenbifurkation

In zwei Dimensionen gibt es ebenfalls die Sattelknoten- und die Heugabelbi-furkation. Wir nehmen zur Normalform in einer Dimension noch die zweiteDimension hinzu. Fur die Sattelknotenbifurkation erhalten wir dann

x = r + x2

y = −y . (3.24)

42

Das Phasenportrait ist in Abbildung 3.21 gezeigt. Die Verallgemeinerung auf

x x

r<0 r=0 r>0

x

y

Abbildung 3.21: Die Sattelknotenbifurkation in zwei Dimensionen.

weitere Dimensionen ist offensichtlich.

Heugabelbifurkation

Die Normalform der Heugabelbifurkation ist in zwei Dimensionen

x = rx− x3

y = −y . (3.25)

Das Phasenportrait ist in Abbildung 3.22 gezeigt. Naturlich gibt es auch in zwei

x

y

x

r=0

x

r<0 r>0

Abbildung 3.22: Die Heugabelbifurkation in zwei Dimensionen.

Dimensionen die subkritische Variante, bei der die Stabilitat aller Fixpunkteinvertiert ist.

Hopf-Bifurkation

Eine Bifurkation, die es in einer Dimension nicht gibt, ist die Hopf-Bifurkation.Hier wird eine Spirale instabil, wenn ein Parameter variiert wird. Statt desFixpunktes gibt es nach der Bifurkation einen Grenzzyklus. Die Normalform ist(in Polarkoordinaten)

r = µr − r3

ϑ = ω + br2 . (3.26)

Fur µ < 0 gibt es eine stabile Spirale am Ursprung. Das Vorzeichen von ωbestimmt den Drehsinn der Spirale. Fur µ > 0 gibt es eine instabile Spirale

43

am Ursprung und einen stabilen Grenzzyklus bei r =√

µ. (Siehe Abb. 3.23.)Der Realteil der beiden komplex konjugierten Eigenwerte der Stabilitatsmatrix

µ>0µ<0

Abbildung 3.23: Die Hopfbifurkation.

wechselt bei einer Hopfbifurkation das Vorzeichen. In unserem Beispiel ist µdieser Parameter.

Ebenso wie fur die Heugabelbifurkation gibt es auch fur die Hopfbifurkationeine subkritische Variante (die normale Hopfbifurkation wird oft auch die su-perkritische Hopfbifurkation genannt). Hierbei verschmilzt ein instabiler Grenz-zyklus mit einer stabilen Spirale und lasst eine instabile Spirale ubrig. Nachder Bifurkation muss die Trajektorie zu einem entfernten Attraktor springen.Dieser kann ein Fixpunkt, ein Grenzzyklus, oder ein chaotischer Attraktor sein.Ein Beispiel ist durch die Gleichungen

r = µr + r3 − r5

ϑ = ω + br2 . (3.27)

gegeben. (Siehe Abb. 3.24.)

µ>0µ<0

Abbildung 3.24: Beispiel fur eine subkritische Hopfbifurkation.

Wir betrachten ein weiteres Beispiel zur Hopf-Bifurkation. Gegeben ist dasGleichungssystem

x = µx− y + xy2

y = x + µy + y3 .

Wir betrachten die Umgebung des Fixpunktes x = y = 0. Die Stabilitatsmatrix

44

ist dort

A =

(

µ −11 µ

)

,

d.h. τ = 2µ und ∆ = µ2 + 1. Weil τ2 − 4∆ negativ ist, ist der Fixpunkt eineSpirale. Fur µ < 0 ist sie stabil, fur µ > 0 instabil. Um zu entscheiden, ob dieHopf-Bifurkation subkritisch oder superkritisch ist, mussen wir die Radialbewe-gung genauer untersuchen. Es ist

1

2

d

dtr2 =

1

2

d

dt(x2 + y2) = xx + yy = (µ + y2)r2

undr = µr + y2r > µr .

Wenn µ > 0 ist, ist r auf jeden Fall positiv, und es gibt keinen Grenzzyklus inder Umgebung der instabilen Spirale fur kleine µ. Damit ist die Hopf-Bifurkationsubkritisch. Fur µ < 0 muss es also einen instabilen Grenzzyklus geben, der furµ→ 0 mit dem Fixpunkt verschmilzt.

Weitere Bifurkationen von Grenzzyklen

Grenzzyklen konnen noch durch 3 weitere Bifurkationen entstehen oder zerstortwerden. Sie sind schwerer zu finden als eine Hopf-Bifurkation, weil bei ihnennicht nur die Umgebung eines Fixpunktes, sondern ein großerer Bereich desPhasenraums beteiligt ist.

• Sattelknotenbifurkation periodischer Bahnen: Wir betrachten das Glei-chungssystem fur die subkritische Hopf-Bifurkation (3.27). Hier gibt esneben der Hopf-Bifurkation bei µ = 0 noch eine zweite Bifurkation. Furµc = −1/4 hat die radiale Gleichung eine Sattelknotenbifurkation. Bei die-sem Wert entstehen aus dem Nichts ein stabiler und ein instabiler Grenz-zyklus, zusatzlich zu der stabilen Spirale am Ursprung. Das Anfertigender Bilder hierzu wird dem Leser uberlassen.

• “Infinite-period”-Bifurkation (Sattelknotenbifurkation auf einem Grenzzy-klus): Hierbei wird die Bewegung auf dem Grenzzyklus an einem Punktimmer langsamer, bis schließlich dort eine Sattelknotenbifurkation auftritt.Ein Beispiel ist das Gleichungssystem

r = r(1 − r2)

θ = µ− sin θ (3.28)

mit µ > 0. Die Bifurkation geschieht bei µ = 1, wie in Abbildung 3.25gezeigt.

• Homokline Bifurkation: Beispiel fur eine subkritische Hopfbifurkation. Die-se Bifurkation passiert, wenn ein Grenzzyklus in die Nahe eines Sattel-punkts wandert und ihn schließlich beruhrt. Dies passiert bei µc ≃ −0.8645im Gleichungssystem

x = y

y = µy + x− x2 + xy . (3.29)

Qualitativ sieht das Phasenportrait so aus wie in Abbildung 3.26 darge-stellt.

45

µ>1 µ<1

Abbildung 3.25: Sattelknotenbifurkation auf einem Grenzzyklus.

x

y

x

y

x

y

µ=µ µ>µµ<µc c c

Abbildung 3.26: Homokline Bifurkation

3.2.3 Phasenportraits und topologische Uberlegungen

Die quantitative Analyse eines zweidimensionalen nichtlinearen System kommtselten ohne Computer aus, aber dennoch kann man viele qualitative Aussagenuber das System durch einfache Uberlegungen und Zeichnungen machen. Nebender Bestimmung von Fixpunkten und einer linearen Stabilitatsanalyse in ihrerUmgebung ist auch eine Bestimmung der Nullklinen hilfreich. Dies sind Linien,langs derer x = 0 oder y = 0 ist. Sie trennen diejenigen Bereiche voneinander, indenen x oder y verschiedenes Vorzeichen haben. Damit weiß man, wo die Pfeileim Phasenportrait nach rechts oben, nach rechts unten, nach links oben, undnach links unten laufen.

Der Index einer Kurve

Um zu entscheiden, ob in einem bestimmten Bereich des Phasenraums Fixpunk-te oder Grenzzyklen liegen, ist es hilfreich, den Index gschlossener Kurven zubestimmen. Der Index einer geschlossenen Kurve C ist eine ganze Zahl, die misst,wie viele Umdrehungen ein Vektorfeld langs dieser Kurve macht. Sei ~x = ~F (x)ein stetiges Vektorfeld im zweidimensionalen Phasenraum, und sei C eine sichselbst nicht schneidende geschlossene Kurve, die durch keinen Fixpunkte geht.Dann hat das Vektorfeld an jedem Punkt der Kurve einen wohldefinierten Win-kel φ ∈ [0, 2π) mit der positiven x-Achse. Wenn man gegen den Uhrzeigersinn

46

langs der Kurve wandert, andert sich der Winkel φ kontinuierlich (wenn man 0und 2π als identisch betrachtet). Nach einem Umlauf hat man wieder den ur-sprunglichen Winkel. Also ist die uber den Umlauf integrierte Winkelanderungein Vielfaches von 2π. Der Index IC ist diese Winkelanderung, geteilt durch 2π.Beispiele fur IC = 1 und IC = −1 sind in Abbildung 3.27 gezeigt. Der Index

Abbildung 3.27: Beispiele fur IC = 1 (links) und IC = −1 (rechts).

hat einige nutzliche Eigenschaften:

1. Zwei Kurven, die ineinander deformiert werden konnen, ohne dabei einenFixpunkt zu uberqueren, haben denselben Index. Grund: Bei der Defor-mation muss sich der Index kontinuierlich andern. Da er eine ganze Zahlist, geht das nur, wenn er gleich bleibt.

2. Eine Kurve, die keinen Fixpunkt einschließt, hat den Index Null. Grund:Wir konnen sie auf einen winzig kleinen Kreis zusammenziehen, ohne dabeiden Index zu andern. Aber auf einem winzig kleinen Kreis andert sich einstetiges Vektorfeld praktisch nicht.

3. Der Index andert sich nicht, wenn man alle Pfeile des Vektorfelds rum-dreht.

4. Wenn C eine Trajektorie des dynamischen Systems ist, ist IC = 1.

5. Man kann einem Fixpunkt auch einen Index zuordnen. Dies ist der In-dex derjenigen Kurven, die nur diesen Fixpunkt einschließen und keinenanderen.

6. Wenn C mehrere isolierte Fixpunkte einschließt, ist der Index dieser Kur-ve die Summe der Indizes der eingeschlossenen Fixpunkte. Grund: Wirkonnen die Kurve auf lauter kleine Kreise um die Fixpunkte und ganzschmale “Brucken” zwischen den Fixpunkten zusammenschrumpfen las-sen. Die Winkelanderung des Vektorfelds langs dieser deformierten Kur-ve ergibt sich aus den Winkelanderungen auf den Kreisen und auf denBrucken (dort ist sie Null, da die Brucken in beide Richtungen durchlau-fen werden).

7. Jede geschlossene Trajektorie eines dynamischen Systems muss Fixpunkteeinschließen, die in der Summe einen Index 1 haben.

8. Alle isolierten Fixpunkte in Abbildung 3.20 außer dem Sattelpunkt habenden Index I = 1, der Sattelpunkt hat I = −1.

47

Wie man periodische Bahnen ausschließen kann

Wenn man den Verdacht hat, dass es keine geschlossene Trajektorie in einemSystem gibt, kann man mehrere Methoden ausprobieren, um dies zu zeigen.

• Mit Hilfe der Indextheorie: Eine geschlossene Bahn muss Fixpunkte ein-schließen, deren Index insgesamt 1 ist. In einigen Situationen ist es moglichzu zeigen, dass es solche Bahnen nicht gibt (Siehe Aufgabe 1 unten).

• Wenn ein Gradientensystem vorliegt: Wenn das dynamische System inder Form ~x = −∇V geschrieben werden kann, kann es keine geschlossenenBahnen geben. Grund: Die Anderung von V langs der geschlossenen Bahn

ware ∆V =∫ T

0∇V · ~xdt < 0. Aber es muss ∆V = 0 sein, weil V eindeutig

ist.

• Mit Hilfe einer Lyapounov-Funktion: Eine Lyapounov-Funktion L(~x) isteine Funktion, die langs aller Trajektorien abnimmt. Also ist dL/dt < 0uberall, außer an Fixpunkten. In diesem Fall gibt es offensichtlich keinegeschlossenen Trajektorien wegen der Eindeutigkeit von L.

• Dulac’s Kriterium: Wenn es eine stetig differenzierbare Funktion g(~x) gibt,

so dass ∇(g~x) in einem einfach zusammenhangenden PhasenraumgebietR uberall dasselbe Vorzeichen hat, kann es in R keine geschlossenen Tra-jektorien geben. Um dies zu zeigen, integriert man ∇(g~x) uber die von

der Trajektorie eingeschlossene Flache A,∫

A∇(g~x)d2x =

∮

Cg~x · ~ndl. Die

rechte Seite muss Null sein, da ~x auf die Normale ~n sekrecht steht.

Das Poincare-Bendixon-Theorem

Das Poincare-Bendixon-Theorem besagt, dass eine Trajektorie, die in einemabgeschlossenen, beschrankten Bereich des zweidimensionalen Phasenraums ge-fangen ist, in dem es keine Fixpunkte gibt, in eine periodische Trajektorie laufenmuss. Dieses Theorem ist sehr plausibel, aber der Beweis ist schwierig, und daherverzichten wir hier auf ihn.

Das Poincare-Bendixon-Theorem hat wichtige Konsequenzen. Zum einemkonnen wir es benutzen, um die Existenz eines periodischen Orbits zu zeigen.Wir mussen hierzu nur einen ringformigen Bereich finden, der keine Fixpunkteenthalt, und an dessen Randern das Vektorfeld uberall nach innen zeigt. Dannbleiben namlich Trajektorien, die in diesen Bereich laufen, fur immer dort, undwir konnen das Theorem anwenden und folgern, dass es in dem Bereich eineperiodische Trajektorie geben muss.

Außerdem folgt aus dem Poincare-Bendixon-Theorem, dass es in zweidimen-sionalen dynamischen Systemen keine auf ein endliches Gebiet beschranktenchaotischen Attraktoren geben kann. Also ist Chaos erst in drei Dimensionenmoglich.

Aufgaben

1. Ein Modell zweier Spezies, die sich von derselben Nahrung ernahren (“Scha-fe und Kaninchen”) ist zum Beispiel durch die folgenden Gleichungen ge-

48

geben:

x = x(3 − x− 2y)

y = y(2− x− y) .

Die erste Gleichung ist fur die Kaninchen, die zweite fur die Schafe.

(a) Erklare die anschauliche Bedeutung der verschiedenen Terme und derGroßenverhaltnisse der Parameter.

(b) Finde die vier Fixpunkte des Systems und untersuche ihre Stabilitat.

(c) Zeichne die Fixpunkte und die Nullklinen in der x− y-Ebene.

(d) Zeige mit Hilfe der Indextheorie, dass dieses System keine geschlos-senen Trajektorien hat.

(e) Skizziere ein paar Trajektorien in der x-y-Ebene, so dass klar ist,wie sich die Populationen von welchem Anfangspunkt aus entwickelnwerden.

(f) Formuliere in Worten die Schlussfolgerung der Rechnung: Konnenzwei Spezies, die miteinander um dieselbe Nahrung konkurrieren, ko-existieren?

2. Nun verallgemeinere das Modell der vorigen Aufgabe, so dass die Koeffi-zienten beliebige Werte annehmen konnen:

N1 = r1N1(1−N1/K1)− b1N1N2

N2 = r2N2(1−N2/K2)− b2N1N2 .

Alle 6 Parameter seien positiv. Was ist die biologische Interpretation dieserParameter? Wie viele dieser Koeffizienten kann man durch eine geeigneteWahl der Skalierung auf 1 setzen? Zeige, dass es jetzt 4 qualitativ ver-schiedene Phasenportraits gibt. Was sind die Bedingungen dafur, dass diebeiden Spezies koexistieren?

3. Ein einfaches Modell fur die Ausbreitung einer Epidemie ist

x = −kxy

y = kxy − ly .

x und y sind die Zahl der gesunden und kranken Leute.

(a) Was bedeuten die Terme in den beiden Gleichungen anschaulich?

(b) Finde alle Fixpunkte und berechne ihre Stabilitat.

(c) Zeichne die Linien, auf denen x = 0 oder y = 0 ist.

(d) Finde eine Erhaltungsgroße. Bilde zu diesem Zweck eine Differenzi-algleichung fur dy/dx und integriere sie nach Variablenseparation.

(e) Zeichne das Phasenportrait. Was passiert fur t→∞ ?

(f) Sei (x0, y0) die Anfangsbedingung. Man sagt, dass eine Epidemie auf-tritt, wenn y(t) zunachst ansteigt. Unter welchen Bedingungen tritteine Epidemie auf?

49

4. Ein einfaches Modell fur oszillierende chemische Reaktionen ist der Brus-selator

x = 1− (b + 1)x + ax2y

y = bx− ax2y .

x und y sind die dimensionslos gemachten Konzentrationen zweier Sub-stanzen. a und b sind positive Parameter.

(a) Was bedeuten die Terme in den beiden Gleichungen anschaulich?

(b) Finde alle Fixpunkte und berechne ihre Stabilitat.

(c) Zeichne die Nullklinen. In welcher Region mussen sich die Trajekto-rien nach genugend langer Zeit befinden?

(d) Zeige, dass fur einen Parameterwert b = bc eine Hopf-Bifurkationauftritt. Finde bc.

(e) Gibt es den Grenzzyklus fur b > bc oder b < bc? Benutze zur Beant-wortung der Frage das Poincare-Bendixon-Theorem.

(f) Finde naherungsweise die Periode des Grenzzyklus nahe bei bc.

5. Ein stetiges Vektorfeld im zweidimensionalen Phasenraum habe genauzwei geschlossene Trajektorien, von denen eine innerhalb der anderen liegtund im entgegengesetzten Sinn umlauft. Ist die folgende Aussage falschoder richtig: Es muss mindestens einen Fixpunkt zwischen den beiden pe-riodischen Orbits geben. Wenn die Aussage richtig ist, beweise sie. Wennsie falsch ist, bringe ein Gegenbeispiel.

3.3 Seltsame Attraktoren

Dynamische Systeme mit mehr als zwei Dimensionen konnen Chaos zeigen. Esgibt keine allgemein akzeptierte Definition von Chaos, aber die folgenden For-mulierung erfasst das Wesentliche am Chaos: “Chaos ist aperiodisches Langzeit-verhalten in einem deterministischen System, das sehr empfindlich von den An-fangsbedingungen abhangt.” Normalerweise ist mit der empfindlichen Abhangig-keit ein exponenziell in der Zeit anwachsender Abstand infinitesimal benachbar-ter Trajektorien verbunden. Chaotische Attraktoren nennt man auch “SeltsameAttraktoren”. einen seltsamen Attraktor, namlich den fur das periodisch getrie-bene Pendel haben wir schon kennengelernt. Im Folgenden betrachten wir nochein historisch bedeutendes Modell. Wer sich genauer uber die Methoden zurAnalyse seltsamer Attraktoren informieren mochte, kann sich eines der vielenBucher uber Chaos besorgen oder mal nachschauen, was in der Gruppe Bennergeforscht wird.

Der Meteorologie-Professor Edward Lorenz vom MIT war in den 60er Jahreneiner der ersten, die einen seltsamen Attraktor in einem System aus drei Varia-blen fanden und die empfindliche Abhangigkeit des Systems von den Anfangsbe-dingungen zeigten. Er wollte den Grunden fur die schwere Vorhersagbarkeit desWetters nachgehen. Er begann mit dem Rayleigh-Benard-System als Modell furKonvektionsrollen in der Athmosphare. Im Rayleigh-Benard-Experiment wirdeine Flussigkeit zwischen zwei horizontalen Platten eingesperrt und von untengeheizt. Wenn der Temperaturgradient einen Schwellwert uberschreitet, bilden

50

sich in der Flussigkeit Konvektionsrollen aus. Wenn der Temperaturgradientnoch großer wird, findet man Mehrfachrollen, und schließlich sogar Chaos. Diehydrodynamischen Gleichungen, die dieses System beschreiben, kann man imFourierraum als ein dynamisches System aus unendlich vielen Variablen (denFourierkomponenten) ausdrucken. Lorenz betrachtete nur die drei fuhrendenFourierkomponenten und erhielt das folgende Gleichungssystem:

x = s(−x + y)

y = rx − y − xz

z = −bz + xy . (3.30)

Die zeitliche Anderung des Phasenraumvolumens berechnet sich mit Gl. (3.3)zu

dV

dt= −(s + 1 + b)V < 0 ,

d.h. das Phasenraumvolumen nimmt exponenziell in der Zeit ab. Es gibt indiesem Modell Fixpunkte und Grenzzyklen. Außerdem findet man fur weiteParameterbereiche seltsame Attraktoren, wie den in Abb. 3.28 gezeigten. Der

-20

20 -30

300

50

Abbildung 3.28: Der Lorenz-Attraktor fur r = 26, s = 9, b = 2.6.

Attraktor hat die Form zweier Spiralen. Die Trajektorie lauft zunachst mehrfachin der einen Spirale um, wobei sie immer weiter nach außen kommt, dann wech-selt sie auf die andere Seite und macht dort mehrere Umlaufe, etc. Die Zahl derUmlaufe auf den beiden Seiten andert sich dabei auf scheinbar erratische Weise.Wenn man zwei Trajektorien vergleicht, die anfangs nahe beieinander liegen,wird nach einiger Zeit die Zahl der Umlaufe auf den beiden Seiten verschieden.

3.4 Die fraktale Dimension

Eine wichtige Eigenschaft eines Attraktors oder eines Satzes von Punkten istseine Dimension. Ein Punktattraktor hat die Dimension 0, und ein Grenzzy-klus hat die Dimension 1. Seltsame Attraktoren haben auch auf beliebig kleinen

51

Skalen eine Unterstruktur, und deshalb ist es eine nicht triviale Aufgabe, ihneneine Dimension zuzuordnen. In diesem Abschnitt fuhren wir das Konzept frak-taler, d.h. nicht ganzzahliger Dimensionen ein und stellen Methoden zu ihrerBerechnung vor. Der Name “Fraktale” wurde von dem Mathematiker BenoitMandelbrot eingefuhrt.

3.4.1 Kastchenzahlmethode

Betrachten wir eine Menge, die in einem N -dimensionalen Raum liegt. Wiruberdecken den Raum mit einem Gitter aus N -dimensionalen “Wurfeln” derKantenlange ǫ. (In einem eindimensionalen Raum sind die “Wurfel” Intervalle,in zwei Dimensionen sind sie Kastchen.) Wir zahlen dann die Zahl von WurfelnN(ǫ), die benotigt werden, um die Menge zu uberdecken. Wir wiederholen diesfur immer kleinere Werte von ǫ. Die Dimension der Menge ist definiert als

D0 = limǫ→0

ln N(ǫ)

ln(1/ǫ). (3.31)

Als erstes Beispiel betrachten wir einige einfache Mengen in einem zweidimen-sionalen Raum (siehe Abb.3.29). Fur die beiden Punkte benotigt man N(ǫ) = 2Kastchen fur alle Werte von ǫ, die kleiner sind als der Abstand der beiden Punk-te. Also ist D0 = 0. Fur das Kurvensegment benotigt man im Grenzfall kleinerǫ eine Zahl N(ǫ) ∼ l/ǫ an Kastchen, wobei l die Lange des Kurvensegments ist.Also ist D0 = 1. Fur die Flache benotigt man fur kleine ǫ ein Zahl N(ǫ) ≃ A/ǫ2

an Kastchen, wobei A die Flache ist. Also ist D0 = 2. Fur nichtfraktale Mengenergibt sich also die richtige Dimension.

Nun betrachten wir ein Beispiel fur eine fraktale Menge, die Cantor-Menge.Man erzeugt sie wie folgt: Man beginnt mit dem Intervall [0, 1] und entfernt dasmittlere Drittel (1/3, 2/3). Nun hat man zwei Intervalle der Lange 1/3 von denenman jeweils wieder das mittlere Intervall entfernt. Dann hat man 4 Intervalleder Lange 1/9. Wenn man diesen Prozess unendlich oft fortsetzt, erhalt man dieCantor-Menge. (Siehe Abb. 3.30). Um ihre fraktale Dimension zu bestimmen,wahlen wir die Folge ǫn = 3−n von ǫ-Werten. Da sich die Zahl von Intervallenbei jedem Schritt verdoppelt, haben wir N(ǫn) = 2n und

D0 = ln 2/ ln 3 = 0.63... .

3.4.2 Das Spektrum der fraktalen Dimensionen

Um die Notwendigkeit fur weitere Definitionen der Dimension eines seltsamenAttraktors zu erkennen, bedecken wir den Attraktor mit einem Gitter aus klei-nen Wurfeln. Wir bestimmen dann, welchen Anteil µi der Zeit eine Trajektoriein dem Wurfel i verbringt, wenn die Zeit gegen unendlich geht. Wenn µi fur al-le Anfangspunkte (bis auf einen verschwindenen Anteil) im Einzugsbereich desAttraktors gleich ist, nennen wir die µi das naturliche Maß der Wurfel. Im All-gemeinen ist es so, dass der großte Anteil des Maßes nur von einem kleinen Anteilder Wurfel kommt, die den chaotischen Attraktor bedecken. Sie werden von derTrajektorie viel haufiger besucht als die anderen Wurfel. Die Kastchenzahldi-mension gibt uns aber keine Auskunft hieruber, da sie nur die Zahl der Wurfelbestimmt, die den Attraktor bedecken, und nicht berucksichtigt, dass verschiede-ne Wurfel verschiedenen wichtig sind. Um zu berucksichtigen, dass verschiedene

52

(a) (b)

(c)

Abbildung 3.29: Illustration von N(ǫ) fur Mengen aus (a) zwei Punkten, (b)einem Kurvensegment, (c) der Flache innerhalb einer geschlossenen Kurve. DerWert von N(ǫ) in diesem Beispielen ist 2, 11, 15.

Wurfel ein verschiedenes Maß haben, haben Grassberger (1983) und Henschelund Procaccia (1983) eine allgemeinere Dimension Dq eingefuhrt, die von einemkontinuierlichen Parameter q abhangt:

Dq =1

1− qlimǫ→0

ln I(q, ǫ)

ln 1/ǫ, (3.32)

wobei

I(q, ǫ) =

N(ǫ)∑

i=1

(µi)q

ist, und die Summe ist uber alle N(ǫ) Wurfel der Kantenlange ǫ zu nehmen, dieden Attraktor bedecken. Fur q > 0 bekommen diejenigen Wurfel mehr Gewicht,die ein großeres Maß haben. Fur q = 0 erhalten wir I(0, ǫ) = N(ǫ), und damitist D0 identisch mit der Kastchenzahldimension. Wenn alle µi gleich sind, istµi = 1/N(ǫ) und ln I(q, ǫ) = (1− q) ln N(ǫ). In diesem Fall ist Dq = D0 fur alleq. Wir definieren D1 uber den Limes D1 = limq→1 Dq und erhalten (z.B. unterVerwendung der Regel von l’Hopital)

D1 = limǫ→0

∑N(ǫ)i=1 µi ln µi

ln ǫ≡ lim

ǫ→0

S(ǫ)

ln(1/ǫ). (3.33)

53

1/3

1/3

2/3

2/31/9 2/9 7/9 8/90

0

0

1

1

1

Abbildung 3.30: Die ersten drei Schritte der Konstruktion einer Cantor-Menge.

Dies ist die sogenannte Informationsdimension. Die Bezeichnung kommt daher,

dass S(ǫ) = −∑N(ǫ)

i=1 µi ln µi die Information ist, die wir gewinnen, wenn wirdie µi kennen und erfahren, dass die Trajektorie im Wurfel Nummer j ist. DieInformationsdimension sagt uns also, wie dieser Informationsgewinn wachst,wenn ǫ gegen Null geht.

Aus der statistischen Physik kennen wir S(ǫ) als Entropie. Je großer dieEntropie ist, desto mehr Mikrozustande gibt es zu einem gegebenen Makrozu-stand. Wenn wir den Makrozustand kennen, ist die Entropie ein Maß fur unsereUnwissenheit uber den Mikrozustand. Wenn wir aber erfahren, in welchem Mi-krozustand sich das System befindet, ist unser Informationsgewinn um so großer,je großer vorher die Ungewissheit (also die Entropie) war. So kommt es, dassder Ausdruck fur den Informationsgewinn der gleiche wie der fur die Entropieist.

Aufgaben

1. Bestimme die fraktale Dimension derjenigen Cantormenge, die entsteht,wenn man jeweils von jedem Intervall einen Abschnitt der halben Inter-valllange (statt wie vorher eine Drittel Intervalllange) aus der Mitte ent-fernt.

2. Bestimme die fraktale Dimension der Kochkurve, die dadurch entsteht,dass man jeweils uber dem mittleren Drittel jedes Intervalls die Kanteneines gleichseitigen Dreiecks einzeichnet.

3. Konstruiere ein fraktales Objekt nach der folgenden Regel: Beginne miteinem Quadrat. Teile es in neun gleiche Quadrate. Nimm alle außer denvier Eckquadraten weg. Wiederhole mit jedem Eckquadrat diese Prozedur,etc. Was ist die fraktale Dimension dieses Objekts?

54

Kapitel 4

Phasenubergange und

kritische Phanomene

In diesem und dem nachsten Kapitel behandeln wir raumlich ausgedehnte Sy-steme. Bevor wir im nachsten Kapitel die Ausbildung von raumlichen und raum-zeitlichen Mustern im Nichtgleichgewicht behandeln, befassen wir uns zunachstmit Strukturen, die im Gleichgewicht auftreten. Wir konzentrieren uns hierbeiauf Phasenubergange zwischen verschiedenen Strukturen, wie sie bei Anderungeines Kontrollparameters, wie z.B. der Temperatur, des Drucks oder der Kon-zentrationsverhaltnisse passieren. Unser Hauptinteresse gilt hierbei den konti-

nuierlichen Phasenubergangen, die mit kritischem Verhalten verbunden sind.Phasenubergange zeigen gewisse Parallelen zu den Bifurkationen, die wir im

vorigen Kapitel behandelt haben. Nach einer allgemeinen Einfuhrung werden wirmit der Landau-Theorie beginnen, die formal der Normalform einer Heugabel-bifurkation gleicht. Diese fuhrt uns auf die Definition der kritischen Exponenten.Dann werden wir explizit die raumliche Struktur berucksichtigen und Fluktua-tionen zunachst in Gaußscher Naherung behandeln. Da kritisches Verhalten miteiner Selbstahnlichkeit des Systems einhergeht, bietet sich die Methode der Re-

normierung zur quantitativen Analyse an, die wir in ihrer einfachsten Variantevorstellen werden. Hier begegnen wir wieder dem Konzept von Trajektorien undFixpunkten, wenn auch mit einer ganz anderen Bedeutung als bei dynamischenSystemen. Am Ende dieses Kapitels behandeln wir einen speziellen Phasenuber-gang, den Perkolationsubergang, der rein mathematisch definiert werden kannund der in mancherlei Hinsicht der denkbar einfachste Phasenubergang ist.

4.1 Einfuhrung in Phasenubergange

Phasenubergange finden statt, wenn die freie Energie oder eine ihrer Ableitun-gen eine Singularitat hat. Die Singularitat kann ein Sprung sein, ein Knick,oder eine Divergenz ins Unendliche. Man beobachtet dort meist einen scharfenUbergang in den Eigenschaften einer Substanz. Ein Beispiel ist der flussig-Gas-Ubergang, der in Abbildung 4.1 im Phasendiagramm des Wassers zu sehen ist.Langs der Grenzlinien findet ein Phasenubergang erster Ordnung statt, d.h. einPhasenubergang, bei dem eine erste Ableitung der freien Energie unstetig ist.

55

Abbildung 4.1: Phasendiagramm des Wassers, schematisch.

In unserem Beispiel hat die Entropie

S = −(

∂H

∂T

)

P

einen Sprung, und folglich wird bei der Kondensation oder beim Frieren dielatente Warme ∆Q = T∆S frei. Das Hauptinteresse dieses Kapitels gilt demkritischen Punkt, an dem der flussig-Gas-Phasenubergang verschwindet. An die-sem Punkt findet ein kontinuierlicher Phasenubergang (d.h. ein Phasenubergangzweiter oder hoherer Ordnung) statt, und jenseits dieses Punktes gibt es nurnoch eine Phase, die man als sehr dichtes Gas oder eine sehr energiehaltigeFlussigkeit verstehen kann.

Ein weiteres bekanntes Beispiel fur einen Phasenubergang ist der Ferroma-gnet. Sein Phasendiagramm ist in Abbildung 4.2 gezeigt. Die beiden Phasen,zwischen denen der Phasenubergang erster Ordnung stattfindet, sind diejeni-gen mit positiver und mit negativer Magnetisierung. Oberhalb der kritischenTemperatur gibt es diesen Phasenubergang nicht mehr.

Abbildung 4.2: Phasendiagramm eines einfachen Ferromagneten.

Weitere Beispiele fur Phasenubergange sind der Ubergang zwischen Normal-

56

CuZn

Abbildung 4.3: Links: Die beiden moglichen Orientierungen des NH+4 -Tetraeders

von NH4Cl in der CsCl-Struktur. Unterhalb der kritischen Temperatur bevorzu-gen die Tetraeder eine der beiden Orientierungen. Rechts: β-Messing unterhalb(oben) und oberhalb (unten) von Tc = 738 K.

leiter und Supraleiter, der Ubergang zum Antiferromagneten, die Entmischungvon binaren Flussigkeiten, Ordnungs-Unordnungsubergange in Kristallen, diePhasenubergange von Flussigkristallen und von Lipidmolekulen in Wasser unddie Mikrophasenseparation von Diblock-Copolymeren. Die Abbildungen 4.3 bis4.6 illustrieren einige dieser Beispiele.

Im Allgemeinen wachst der Grad der Ordnung unterhalb des kritischenPunktes. Die Ursache hierfur ist die Konkurrenz zwischen der Entropie undder Energie. Meist bedeutet mehr Ordnung niedrigere Entropie und niedrigereEnergie. Aus der Beziehung F = E − TS erkennt man, dass bei tieferen Tem-peraturen die freie Energie weniger stark von der Entropie abhangt, so dass diegeordnete Phase mit der niedrigeren Energie auch die niedrigere freie Energiehat. Der Ubergang zu hoherer Ordnung bedeutet das Brechen einer Symme-trie. Beim Ferromagneten ist dies die Symmetrie zwischen den verschiedenenOrientierungen der Magnetisierung.

Abbildung 4.4: Schematische Anordnung der Flussigkristallmolekule in den ver-schiedenen Phasen eines Flussigkristalls.

Phasenubergange haben gewohnlich einen Ordnungsparameter, der unter-halb der kritischen Temperatur von Null verschieden ist und oberhalb ver-schwindet. Fur den flussig-Gas-Ubergang ist dies die Dichtedifferenz zwischenden beiden Phasen, und fur den Ferromagneten ist dies die Magnetisierung. Oft

57

Abbildung 4.5: Schematische Darstellung der Struktur von in Wasser gelostenLipidmolekulen. Von links oben nach rechts unten wachst die Molekulkonzen-tration.

gibt es ein “Feld”, das einen Ubergang zwischen den verschiedenen moglichenTieftemperaturphasen veranlasst. Im Ferromagneten ist dies das Magnetfeld, imflussig-Gas-System der Druck. Oberhalb der kritischen Temperatur ruft das Feldebenfalls eine Symmetriebrechung hervor (sofern der Phasenubergang mit einerSymmetriebrechung verbunden ist). In vielen Systemen kann man das “Feld”theoretisch definieren, aber oft lasst es sich nicht auf einfache Weise experimen-tell realisieren. Man denke z.B. an den Ordnungsubergang in β-Messing, wo das“Feld” die Zinkatome auf das eine Untergitter und die Kupferatome auf dasandere Untergitter ziehen musste.

4.2 Mean-Field-Theorie

Eine Mean-Field-Theorie, auf Deutsch auch oft Molekularfeldtheorie genannt,ist i.A. die einfachste Theorie, die man zu Phasenubergangen machen kann. Diegrundlegende Annahme hierbei ist, dass jedes Teilchen nur einen gemitteltenEffekt der anderen Teilchen spurt, so dass raumliche Fluktuationen ignoriertwerden konnen. Dies ist eine Naherung, die meist nicht zu exakten Ergebnissenfuhrt, aber oft die Phanomenologie des Phasenubergangs richtig erfasst. Be-kannte Beispiele sind die van-der-Waals-Theorie fur Gase und die Weiss’scheTheorie fur Magnetismus, die wir beide in den Ubungen wiederholen werden. Indiesem Teilkapitel werden wir die Landau-Theorie fur Phasenubergange behan-deln, die die “Normalform” fur einen Phasenubergang mit symmetrischem Ord-nungsparameter darstellt. Bis ca 1960 waren Mean-Field-Theorien (MFT) dievorherrschenden Theorien fur Phasenubergange. Da die experimentell bestimm-ten Werte der kritischen Exponenten von denen der MFT abwichen, wurden imfolgenden Jahrzehnt weitergehende Theorien etwickelt und auch Computersimu-lationen durchgefuhrt. Diese Theorien beruhten zum einen auf storungstheore-tischen Berechnungen, bei denen in einer kleinen Große systematisch entwickelt

58

Abbildung 4.6: Das Phasendiagramm fur Diblock-Copolymere. Auf der x-achseist das Mischungsverhaltnis aufgetragen, auf der y-Achse die Starke der fur dieEntmischung verantwortlichen Wechselwirkung.

wurde, zum anderen wurden Skalentheorien eingefuhrt, die die Idee, dass ein kri-tisches System selbstahnlich ist, implementierten. Diese Entwicklungen gipfeltenin der Renormierungsgruppentheorie, die 1971 von K.G. Wilson vorgestellt wur-de, und fur die er 1982 den Nobelpreis erhielt.

4.2.1 Landau-Theorie fur Phasenubergange

Eine einfache und elegante phanomenologische Theorie fur die Umgebung deskritischen Punktes wurde 1937 von Landau vorgestellt. Er nahm an, dass diefreie Energie sich in der Nahe des Phasenubergangs allein als Funktion des Ord-nungsparameters m schreiben lasst, wobei die Koeffizienten von den Kontrollpa-rametern des Systems abhangen. Wir betrachten den Fall, dass der Ordnungs-parameter eine skalare Große ist, und dass eine Symmetrie bzgl. des Vorzeichensdes Ordnungsparameters vorliegt. Eine solche Situation ist z.B. fur einen Ferro-magneten mit einer Vorzugsachse gegeben, oder fur den flussig-Gas-Ubergang.In der Nahe des kritischen Punktes ist der Ordnungsparameter klein, und eineEntwicklung in m ergibt als erste Terme

F

kBTN= f0 + rm2 + um4 . (4.1)

Die ungeraden Potenzen treten nicht auf, da die freie Energie nicht vom Vorzei-chen von m abhangt. Damit die freie Energie nach unten beschrankt ist, mussu > 0 sein. Die freie Energie hat also die in Abbildung 4.7 gezeigte Form.

In einem System mit festem Volumen und fester Teilchenzahl, das auf kon-stanter Temperatur gehalten wird, nimmt die Freie Energie im thermischen

59

m

F(m)

r=0r<0r>0

Abbildung 4.7: Die freie Energie als Funktion des Ordnungsparameters in derLandau-Theorie fur verschiedene r. Die Kurven fur r < 0 haben zwei Minima.

Gleichgewicht ihr Minimum an.Man erkennt sofort, dass fur r > 0 die freie Energie ihr Minimum bei m = 0

hat, so dass wir dort in der Hochtemperaturphase sind. Es ist also r = a(T −Tc)mit einem konstanten Koeffizienten a und der kritischen Temperatur Tc. Furr < 0 gibt es die beiden Losungen

m = ±√

−r

2u.

Wir sehen also, dass der Ordnungsparameter in der Nahe des kritischen Punktesmit der Temperatur uber ein Potenzgesetz der Form

m ∝ |T − Tc|β (4.2)

verbunden ist, wobei der Wert des Exponenten β = 1/2 ist.Als nachstes fuhren wir ein Feld h ein und addieren den den Term −mh zu

(4.1),F

kBTN= f0 + rm2 + um4 − mh . (4.3)

Minimieren der Freien Energie ergibt jetzt

2rm + 4um3 − h = 0 .

Fur T = Tc (also r = 0) finden wir

h ∝ mδ (4.4)

mit δ = 3. Wir haben also ein weiteres Potenzgesetz gefunden. Fur die Suszep-tibilitat erhalten wir

(

∂m

∂h

)

T

=1

2r + 12um2=

12r

fur r > 01

4|r| fur r < 0

Es gilt alsoχ ∝ |T − Tc|−γ (4.5)

60

mit γ = 1.Schließlich berechnen wir noch die Warmekapazitat, die ja ebenfalls als zwei-

te Ableitung der freien Energie geschrieben werden kann. Wir erhalten sie, wennwir in (4.1) die Losung fur den Ordnungsparameter einsetzen und −T∂2F/∂T 2

berechnen. Zum Einen gibt es den Beitrag, der aus f0 resultiert, und den wirim Folgenden ignorieren, da er bei Tc regular ist. Der zweite Beitrag resultiertaus den vom Ordnungsparameter abhangigen Termen. Fur T > Tc ist er 0, undfur T < Tc ist er

a2NkBT

2u[T + 2(T − Tc)].

Damit liegt ein Sprung in der Warmekapazitat vor. Wenn wir einen Exponentenfur das Verhalten von CH bei Annaherung an Tc von links oder rechts einfuhrenwollten, so hatte dieser den Wert α = 0.

4.2.2 Dynamik

Wenn der Ordnungsparameter nicht seinen Gleichgewichtswert hat, wird er sichmit der Zeit andern, bis er ihn erreicht hat. Die einfachste Form einer solchenDynamik ist durch die Gleichung

m = −λdF

dm(4.6)

gegeben mit einer Relaxationskonstanten λ. Der Ordungsparameter bewegt sichbergab in der Freie-Energie-Landschaft, bis er ein lokales Minimum erreicht.Durch Einsetzen der Gleichung (4.1) erhalten wir

m = −λNkBT (2rm + 4um3) , (4.7)

was formal einem dynamischen System gleicht, das bei r = 0 eine Heugabelbi-furkation macht. Wenn wir das Feld hinzunehmen, ist das analoge dynamischeSystem eine unvollkommene Heugabelbifurkation.

4.2.3 Kritische Exponenten

Wir haben am Beispiel der Landau-Theorie gesehen, dass sich in der Nahe deskritischen Punktes eine Reihe von kritischen Exponenten definieren lassen. Diesgilt nicht nur in Mean-Field-Theorie, sondern allgemein, wobei der Werte derExponenten i.A. allerdings von denen der Mean-Field-Theorie abweichen. DerGrund fur all diese Potenzgesetze ist, dass es am kritischen Punkt große, langle-bige Fluktuationen des Ordnungsparameters gibt. Eine Flussigkeit am kritischenPunkt sieht milchig aus, weil die Dichtefluktuationen sogar auf Langenskalen sogroß wie die Lichtwellenlange auftreten. Diese Lange betragt ungefahr 4000 A,also das 4000fache des Teilchenabstands! Da in diesen Potenzgesetzen oft derAbstand von der kritischen Temperatur als Parameter auftritt, definiert mandie dimensionslose Große

t =T − Tc

Tc

. (4.8)

Die folgende Tabelle gibt einige Potenzgesetze fur den kritischen Punkt desflussig-Gas-Ubergangs und des einfachen Ferromagneten.

61

Große Gas Magnet

Ordungsparameter ρfl − ρgas ∼ (−t)β M ∼ (−t)β

Warmekapazitat CV ∼ |t|−α CH ∼ |t|−α

Suszeptibilitat κT ∼ |t|−γ χT ∼ |t|−γ

Feld |p − pc| ∼ |ρfl − ρgas|δ |H | ∼ |M |δKorrelationslange ξ ∼ |t|−ν ξ ∼ |t|−ν

Korrelationsfunktion 〈ρ(~x)ρ(~x1)〉 − 〈ρ(~x)〉2 〈m(~x)m(~x1)〉 − 〈m(~x)〉2G(~x − ~x1) ∼ |~x − ~x1|−d+2−η ∼ |~x − ~x1|−d+2−η

Die Korrelationslange ist ein Maß fur die Ausdehnung der großten Fluktua-tionen. Auf genugend großen Langen fallt die Korrelationsfunktion exponentiellmit der Entfernung ab,

G(~x − ~x1) ≃ Ce−|~x−~x1|/ξ .

Das Zeichen ∼ liest man “skaliert wie”. Ein Ausdruck der Form

F (t) ∼ |t|λ (4.9)

bedeutet

λ = limt→0

ln |F (t)|ln |t| .

Die Beziehung (4.9) steht also nicht nur fur einfache Potenzgesetze F = A|t|λ,sondern auch fur kompliziertere Ausdrucke der Form

F (t) = A|t|λ(1 + btλ1 + . . . )

mit einem positiven λ1. Genugend nah an t = 0 dominiert der fuhrende Term,aber weiter entfernt gibt es i.A. Korrekturen zu dem einfachen Potenzgesetz.

In Experimenten zeigte sich, dass total verschiedene Systeme dieselben kri-tischen Exponenten haben. So haben z.B. der Ferromagnet YFeO3 (hat eineVorzugsachse), die Flussigkeiten CO2 und Xe, der Antiferromagnet FeF2, derOrdnungs-Unordnungsubergang in β-Messing und in NH4Cl und das dreidi-mensionale Isingmodell dieselben Werte der Exponenten, α ≃ 0.11, β ≃ 0.33,γ ≃ 1.24, δ ≃ 4.8, ν ≃ 0.63, η ≃ 0.03. Eine weitere Universalitatsklasse beinhal-tet das superfluide Helium, magnetische Systeme mit einer Vorzugsebene undden nematisch-smektischA-Ubergang in Flussigkristallen. Dort haben die Expo-nenten die Werte α ≃ −0.008, β ≃ 0.35, γ ≃ 1.32, δ ≃ 4.8, ν ≃ 0.67, η ≃ 0.03.

Experimente zeigen auch, dass die kritischen Exponenten bei Annaherungvon unten und oben an Tc denselben Wert haben. Die Werte der kritischenExponenten sind nicht alle unabhangig, sondern es gelten die Skalenrelationenγ = −β(1− δ), α = 2− 2β − γ, γ = ν(2− η) und 2− α = dν. Diese und andereEigenschaften kritischer Systeme werden wir spater aus der Renormierungsana-lyse herleiten. Eine entscheidende Annahme hierbei ist, dass fur das kritischeVerhalten mikroskopische Details des Systems unwichtig sind, und dass nur die-jenigen Eigenschaften des Systems relevant sind, die sich auch noch auf großenLangenskalen auswirken. Anders lasst sich auch die Universalitat der kritischenExponenten nicht erklaren. Ihre Werte hangen von folgenden Eigenschaften desSystems ab:

• Die raumliche Dimension

• Die Zahl der Komponenten und die Symmetrie des Ordnungsparameters

62

• Die An- oder Abwesenheit von langreichweitigen Wechselwirkungen

• Die Erhaltungsgroßen, an die der Ordungsparameter koppelt

Der letzte Punkt ist fur die dynamischen kritischen Eigenschaften relevant, dieaber hier nicht behandelt werden. Der wichtigste Effekt der raumlichen Dimen-sion ist, dass in hoheren Dimensionen Fluktuationen um die Mittelwerte wenigerwichtig sind. Wir werden das spater explizit zeigen. Daraus folgt, dass es i.A.eine untere kritische Dimension gibt, unterhalb derer der Phasenubergang uber-haupt nicht mehr auftreten kann, weil die Tieftemperaturphase selbst bei sehrniedrigen Temperaturen durch Fluktuationen zerstort wird. Fur das Isingmo-dell ist dies die Dimension 1. Außerdem gibt es eine obere kritische Dimension,oberhalb derer die Fluktuationen nur kleine Storungen sind, so dass die in derMean-Field-Theorie berechneten Exponenten gultig sind.

Ubungsaufgaben

1. Die freie Energie ist gegeben durch

F = −kBT ln∑

i

e−Ei/kBT ,

wobei die Summe uber alle mikroskopischen Zustande i des Systemsgeht. Dies ist offensichtlich eine analytische Funktion, da die Energie ana-lytisch von den Systemparametern Volumen, Magnetfeld, etc, abhangt.Wie kann es sein, dass F an einem Phasenubergang Singularitaten hat?

2. Was ist der Ordnungsparameter fur die in der Vorlesung erwahnten Ord-nungs-Unordnungsubergange? Fur einen Antiferromagneten? Fur Flussig-kristalle? Fur Supraleiter?

3. Gehe aus vom Phasendiagramm fur Wasser und betrachte den Ubergangzwischen der flussigen und gasformigen Phase. Zeichne in der ρ−P -Ebene(ρ ist die Dichte, die ja in der Flussigkeit hoher ist als im Gas) den quali-tativen Verlauf von Linien konstanter Temperatur (also Isothermen) ein.Wie unterscheiden sie sich ober- und unterhalb von Tc? (Hinweis: Erinneredich an die Behandlung des van-der-Waals-Gases in der Vorlesung uberStatistische Physik.)

4. Gehen vom Phasendiagramm eines Ferromagneten aus. Zeichne qualitativden Verlauf von Isothermen in der M -H-Ebene. Zeichne qualitativ M(T )fur verschiedene Werte von H .

5. Kritische Exponenten in der van-der-Waals-Theorie: Die van-der-Waals-Gleichung lautet

(

p + aN2

V 2

)

(V − Nb) = NkBT . (4.10)

(a) Schreibe diese Gleichung so um, dass du eine kubische Gleichung inV erhaltst und bestimme den kritischen Punkt(pc, Vc, Tc), an demdiese Gleichung sich auf die Form

0 = (V − Vc)3 = V 3 − 3V 2Vc + 3V V 2

c − V 3c (4.11)

reduziert.

63

(b) Fuhre die reduzierten Großen

π =p

pc

, v =V

Vc

, t =T

Tc

(4.12)

ein, und schreibe damit die van der Waals-Gleichung (4.10) in einerForm, die keine Materialkonstanten mehr enthalt und deshalb fur alleSubstanzen Gultigkeit haben sollte.

(c) Fuhre die Variablen ω = v − 1, ǫ = t − 1 und q = π − 1 ein, lose dievan-der-Waals-Gleichung nach q(ǫ, ω) auf und leite daraus schließlichdie folgenden Beziehungen ab

κT ∼ |T − Tc|−γ mit γ = 1 (4.13)

|ρ − ρc| ∼ (Tc − T )β mit β = 1/2 (4.14)

p − pc ∼ |ρ − ρc|δ mit δ = 3 (4.15)

6. Die Weiss’sche Theorie fur den Magnetismus: Die Grundannahme dieserTheorie ist, dass jeder Spin dasselbe lokale (mittlere) Feld sieht, das wi-derum von den Spins selbst erzeugt wird und folglich proportional zurmittleren Spinorientierung ist. Der Mittelwert eines Spins ist also

〈σ〉 =ea〈σ〉 − e−a〈σ〉

ea〈σ〉 + e−a〈σ〉

wobei a proportional zu 1/kBT ist. Zeige, dass diese Theorie nahe bei Tc

mit der Landau-Theorie identisch ist.

7. Zeichne die freie Energie der Landau-Theorie in Gegenwart eines Feldes.Was passiert mit den Minima der freien Energie, wenn bei T < Tc das Feldvon −h auf +h kontinuierlich hochgefahren wird? Erklare damit qualitativdie Hystereseschleife eines Magneten.

64

4.3 Ginzburg-Landau-Funktional und Fluktua-

tionen in gaußscher Naherung

In der MFT haben wir angenommen, dass der Ordnungsparameter nicht vomOrt abhangt, sondern uberall seinen Mittelwert annimmt. Jetzt wollen wir raum-liche Fluktuationen im Ordnungsparameter berucksichtigen. Wir schreiben alsovon nun an m = m(~x), wobei m(~x) z.B. fur einen Magneten die Magnetisie-rung, also das magnetische Moment pro Einheitsvolumen ist. Wir gehen auch indiesem Abschnitt davon aus, dass der Ordnungsparameter ein Skalar ist. Damitwir zum Kontinuumslimes ubergehen konnen, mussen wir das Einheitsvolumengroß wahlen verglichen mit dem Volumen pro Spin.

Wir gehen davon aus, dass sich die Wahrscheinlichkeit fur das Auftreten einerKonfiguration [m(~x)] allein als Funktion des Ordnungsparameters ausdruckenlasst,

P [m(~x)] ∝ e−H[m(~x)]/kBT (4.16)

mit

H[m(~x)]

kBT=

∫

ddx[

rm2(~x) + um4(~x) + c (∇m(~x))2 − hm(~x)

]

. (4.17)

d ist die Dimension des Systems.Wie bei der Landau-Theorie haben wir bei Gleichung (4.17) angenommen,

dass m klein ist, und unter Berucksichtigung der Symmetrie bezuglich des Vor-zeichens von m im Ordnungsparameter entwickelt. Außerdem haben wir an-genommen, dass auch der Gradient von m klein ist und auch hier nur denfuhrenden Term mitgenommen. Der Ausdruck (4.17) wird das Ginzburg-Landau-

Funktional genannt und spielt die Rolle der Energie. Statt es phanomenologischherzuleiten, wie wir es taten, kann man es auch explizit aus dem Isingmodellableiten, indem man uber Zellen, die mehrere Spins enthalten, mittelt, zumKontinuumslimes ubergeht und die fuhrenden Terme behalt.

Ebenso wie in der Landau-Theorie konnen wir auch hier eine dynamischeGleichung hinschreiben. Eine einfache Relaxationsdynamik ist gegeben durch

dm(~x, t)

dt= −λ

δHδm(~x, t)

+ bη(~x, t)

= −λkBT[

2rm(~x, t) + 4um3(~x, t) − 2c∆m(~x, t) − h]

+ bη(~x, t) .

(4.18)

Da der Ordnungsparameter aufgrund der Temperatur lokalen Fluktuationenunterworfen ist, haben wir weißes Rauschen addiert. Wir werden uns in denfolgenden Rechnungen allerdings auf zeitunabhangige Großen beschranken, diemit Hilfe der kanonischen Zustandssumme ausgewertet werden konnen.

Die Zustandsumme ist

Z =

∫

D[m(~x)]e−H[m(~x)]/kBT , (4.19)

wobei∫

D[m(~x)] =

∫

dm(~x1) . . . dm(~xN )

65

mit einer genugend feinen Diskretisierung ist. Einen von der Diskretisierungabhangigen Normierungsfaktor brauchen wir nicht einzufugen, da er sich beider Berechnung physikalischer Großen immer wegkurzt.

Die freie Energie ist F = −kBT ln Z. Sie ist identisch mit Landaus frei-er Energie, wenn der Ordnungsparameter raumlich homogen ist, d.h. wenn ∇mverschwindet. Dann ist namlich H/kBT = N(rm2+um4−hm), und im thermo-dynamischen Limes wird die Zustandssumme von demjenigen m dominiert, dasdie Funktion rm2 + um4 − hm minimiert. Wenn die raumlichen Fluktuationendes Ordnungsparameters um den wahrscheinlichsten Wert klein sind, konnenwir erwarten, dass die MFT gilt. Wir werden sehen, dass in Dimensionen d < 4die Fluktuationen so stark sind, dass die Ergebnisse der MFT ungultig sind.

Im folgenden betrachten wir Fluktuationen um die wahrscheinlichste Konfi-guration fur den Fall h = 0. Wir schreiben

m(~x) = m0 + m1(~x) ,

wobei m0 der aus der Landau-Theorie resultierende Ordnungsparameter ist. Wirdiskutieren zunachst den Fall r > 0, fur den m0 = 0 ist. Wenn die Fluktuationenklein sind, konnen wir den Term m4

1 vernachlassigen und erhalten das Ginzburg-Landau-Funktional

H0/kBT =

∫

ddx(

rm21 + c(∇m1)

2)

. (4.20)

Den Index 0 verwenden wir immer dann, wenn wir ein Ginzburg-Landau-Funk-tional nur bis zur zweiten Ordnung in m1 betrachten. Bei der Berechnung derZustandssumme und von Erwartungswerten treten in diesem Fall nur gaußscheIntegrale auf, die ausintegriert werden konnen.

Um die Zustandssumme auszuwerten, gehen wir zunachst in den Fourier-raum. Hier erweist sich mal die diskrete und mal die kontinuierliche Varianteals zweckmaßiger. Wir schreiben

m1(~x) = L−d∑

~k

m~kei~k~x,

wobei L die Ausdehnung des Systems in jeder Richtung ist und Ld das Volumen.Wenn wir fur das System periodische Randbedinungen annehmen, haben dieKomponenten des Wellenvektors nur Werte, die ein Vielfaches von 2π/L sind.Außerdem gibt es einen großtmoglichen Wert fur die Komponenten des Wel-lenvektors, der von der mikroskopischen Gitterkonstante des Systems abhangt.Wir vereinfachen die Rechnungen, indem wir fur den Betrag des Wellenvektors(statt fur jede der Komponenten getrennt) einen Cutoff Λ einfuhren, da wir

meist in Kugelkoordinaten integrieren. Damit konnen wir die Summe uber ~kdurch ein Integral nahern:

1

Ld

∑

~k

≃ 1

(2π)d

∫

ddk ≃ Kd

∫ Λ

0

dk kd−1.

Der letzte Schritt wird gemacht, wenn der Integrand nur vom Betrag k desWellenvektors abhangt.

Kd ≡∫

dΩd

(2π)d=

(

2d−1πd/2Γ

(

d

2

))−1

66

ist die Oberflache der d-dimensionalen Einheitskugel, dividiert durch (2π)d.Wir erhalten also

H0/kBT =

∫

ddx

L2d

∑

~k,~k1

ei(~k+~k1)~x(

r − c(~k~k1))

m~km~k1

= L−d∑

~k,~k1

δ~k,−~k1

(

r − c(~k~k1))

m~km~k1

= L−d∑

~k

(r + ck2)m~km−~k

= L−d∑

~k

(r + ck2)|m~k|2 . (4.21)

Damit berechnet sich die Zustandssumme in gaußscher Naherung zu

Z0 = C∏

~k

[∫

dm~ke−L−d(r+ck2)|m~k

|2]

= C∏

~k

(

π

r + ck2

)1/2

.

Wir haben Potenzen von Ld in den konstanten Vorfaktor C absorbiert. Mansieht, dass fur r = 0 Singularitaten bei ~k = 0 auftreten. Wir wollen zunachst

die Warmekapazitat CH = −T ∂2F0

∂T 2 auswerten. Mit

F0 = −kBT ln Z0 = −kBT ln C +∑

~k

[

−kBT

2ln π +

kBT

2ln(r + ck2)

]

erhalten wir (mit r = a(T − Tc))

CH = −∑

~k

kBTa

r + ck2+∑

~k

kBT 2a2

2(r + ck2)2.

Im thermodynamischen Limes L → ∞ nahert sich der Betrag des kleinstenWellenvektors, 2π/L, dem Wert Null, und die spezifische Warme divergiert beider kritischen Temperatur. Der letzte Term enthalt die fuhrende Divergenz, undwir werten im folgenden nur ihn aus:

CsingH =

kBT 2a2Ld

2

∫

ddk

(2π)d

1

(r + ck2)2

=kBT 2a2KdL

d

2r2

∫ Λ

0

dkkd−1

(1 + crk2)2

=kBT 2a2KdL

d

2c2

( c

r

)2−d/2∫ Λ

√c/r

0

dqqd−1

(1 + q2)2

≡ kBT 2a2KdLd

2c2ξ4−d

∫ Λξ

0

dqqd−1

(1 + q2)2

≡ kBT 2a2KdLd

2c2ξ4−dI .

Die Große√

c/r ist eine Lange, die wir ξ genannt haben. Es wird sich nachherzeigen, dass es die Korrelationslange ist. Damit ist der Exponent ν berechnetund betragt ν = 1/2.

67

Fur d < 4 konvergiert das Integral I, und wir haben

CH ∼ (T − Tc)−(4−d)/2

und α = (4 − d)/2. Die spezifische Warme hat also eine Divergenz bei Tc.Die Fluktuationen modifizieren somit das kritische Verhalten des Systems ind < 4. In 3 Dimensionen ist α = 1/2. Dies ist großer als der tatsachliche Wert0.11 von α. Die Fluktuationen in gaußscher Naherung uberschatzen also diesenExponenten. Mit Hilfe einer Storungstheorie, die uber die gaußsche Naherunghinausgeht, erhalt man einen realistischeren Wert fur α.

In d = 4 hat I eine logarithmische Divergenz, und es ist

CH ∼ ln(T − Tc) .

In d > 4 ist

I =

∫ ξΛ

0

qd−5dq −∫ ξΛ

0

dqqd−5 + 2qd−3

(1 + q2)2=

(ξΛ)d−4

d − 4−∫ ξΛ

0

dqqd−5 + 2qd−3

(1 + q2)2.

Der zweite Summand ist ein konvergentes Integral (solange d < 6 ist), und wirerhalten

CH = A − B(T − Tc)(d−4)/2 .

Die Warmekapazitat hat bei T = Tc eine Spitze, oft auch “Cusp” genannt.Fur r < 0 ist m0 =

√

|r|/2u, und wir erhalten

H0/kBT =

∫

ddx(

2|r|m21 + c(∇m1)

2)

. (4.22)

Der einzige Unterschied zu Gleichung (4.20) besteht also in der Ersetzung r →2|r|, und die Formeln fur die Korrelationslange und die Warmekapazitat konnenmit einer entsprechenden Modifikation der Vorfaktoren von der Rechnung furr > 0 ubernommen werden.

Als Nachstes berechnen wir die Korrelationsfunktion

G(~x − ~x1) ≡ 〈m(~x)m(~x1)〉 − 〈m(~x)〉〈m(~x1)〉= 〈(m(~x) − m0)(m(~x1) − m0)〉= 〈m1(~x)m1(~x1)〉 . (4.23)

Hierbei haben wir ausgenutzt, dass das System translationsinvariant ist, so dassG nur von der Differenz ~x − ~x1 abhangt und nicht von beiden Orten separat,und dass 〈m(~x)〉 = 〈m(~x1)〉 = m0 ist.

Wir haben bei der Berechnung der Warmekapazitat gesehen, dass die Rech-nungen einfacher im Fourierraum durchzufuhren sind. Wir berechnen also dieGroße

〈m~km~k1

〉 =

∫

ddxddx1e−i~k~x−i~k1~x1〈m(~x)m(~x1)〉

=

∫

ddxddx1e−i~k(~x−~x1)−i~x1(~k+~k1)〈m(~x)m(~x1)〉

=

∫

ddx2ddx1e

−i~k~x2−i~x1(~k+~k1)〈m(~x2 + ~x1)m(~x1)〉

= Ld

∫

ddx2e−i~k~x2G(~x2)δ~k,−~k1

= δ~k,−~k1LdG~k

. (4.24)

68

Im vorletzten Schritt haben wir ausgenutzt, dass wegen der Translationsinvari-anz 〈m(~x2 + ~x1)m(~x1)〉 unabhangig von ~x1 ist. Wenn wir dieselbe Umformungmit kontinuierlichem Impuls (also unendlichem Volumen) wiederholen, erhaltenwir

〈m~km~k1

〉 = (2π)dδd(~k + ~k1)G~k.

In gaußscher Naherung haben wir also den Mittelwert

〈m~km~k1

〉0 ≡ Z−10

∏

~k2

[∫

dm~k2

]

m~km~k1

e−H0[m~k2

]/kBT

zu berechnen. Dies geht am einfachsten, wenn wir ein Hilfsfeld h~keinfuhren und

Hh/kBT = H0/kBT −∑

~k

h~km−~k

definieren. Dann ist namlich (fur r > 0)

〈m~km−~k

〉0 =

(

∂2 ln Zh

∂h~k∂h−~k

)

[h~k]=0

, (4.25)

wobei Zh die mit Hh berechnete Zustandssumme ist,

Zh = C∏

~k

[∫

dm~kexp

−(r + ck2)L−dm~km−~k

+ m~kh−~k

]

= C∏

~k

[

∫

dm~kexp

−(r + ck2)L−d

(

m~km−~k

−Ldm~k

h−~k

r + ck2

)]

= C∏

~k

[

∫

dm~kexp

−(r + ck2)L−d

(

m~k− Ldh~k

2(r + ck2)

)(

m−~k−

Ldh−~k

2(r + ck2)

)

exp

h~kh−~k

Ld

4(r + ck2)

]

= Z0 exp

∑

~k

h~kh−~k

Ld

4(r + ck2)

.

Also ist (mit (4.25))

〈m~km~k1

〉0 =Ldδ~k,−~k1

2(r + ck2)≃ (2π)dδd(~k + ~k1)

2(r + ck2). (4.26)

Unter Verwendung von (4.23) und (4.24) erhalten wir die Korrelationsfunktion

69

im Ortsraum in gaußscher Naherung

G0(~x) ≡ 〈m1(~x)m1(0)〉0=

1

L2d

∑

~k

〈m~km−~k

〉0ei~k~x

=1

2Ld

∑

~k

ei~k~x

(r + ck2)

≃ 1

2c

∫

ddk

(2π)d

ei~k~x

ξ−2 + k2. (4.27)

Fur T < Tc ergibt sich dasselbe Ergebnis, allerdings mit ξ =√

c/2|r|. In dreiDimensionen ergibt sich

G0(~x) =π2

(2π)3c

e−x/ξ

x. (4.28)

Dies ist die Ornstein-Zernicke-Korrelationsfunktion. Bei T = Tc ist ξ = ∞, unddie Korrelationsfunktion zerfallt wie 1/x mit der Entfernung x. Der Exponent ηhat also den Wert 0 in gaußscher Naherung in drei Dimensionen. In der Nahe vonTc sieht man auf Abstanden großer als die Korrelationslange im Wesentlicheneinen exponentiellen Zerfall der Korrelationen.

Zum Abschluss schatzen wir ab, unter welchen Bedingungen die in die-sem Abschnitt verwendete gaußsche Naherung eine gute Naherung ist. Diesist nur dann der Fall, wenn die betrachteten Fluktuationen kleine Korrekturenzur Landau-Theorie machen. Wie wir am Beispiel der Warmekapazitat gesehenhaben, ist dies unterhalb von 4 Dimensionen nicht mehr der Fall. Die Landau-Theorie gab einen Sprung ∆c in der spezifischen Warme, und der Beitrag derFluktuationen zur spezifischen Warme ist nur dann eine kleine Korrektur, wenner viel kleiner als dieser Sprung ist, also wenn

∆c ≫ kBT 2c LdIKda

d/2

2cd/2(T − Tc)

(d−4)/2 .

Durch Umformung erhalt man daraus das Ginzburg-Levanyuk-Kriterium

( |Tc − T |Tc

)(4−d)/2

≫ Nξd0∆c

, (4.29)

wobei N ein numerischer Faktor ist und ξ0 =√

c/aTc. Genugend nahe an Tc

wird diese Ungleichung in weniger als 4 Dimensionen immer verletzt. Aber etwasweiter entfernt von Tc gibt es oft einen Bereich, in dem die gaußsche Naherunggut ist.

Ubungsaufgaben

1. Konnen die Parameter λ und b in (4.18) unabhangig voneinander sein?wenn Nein, welche Beziehung muss zwischen ihnen gelten?

2. Leite Gleichung (4.22) her.

70

3. Fur d < 4 istCH ≃ A+(T − Tc)

−(4−d)/2

fur T > Tc undCH ≃ A−(T − Tc)

−(4−d)/2

fur T < Tc. Berechne das Amplitudenverhaltnis A−/A+.

4. Werte das Integral in (4.27) in 3 Dimensionen aus und zeige damit dieRichtigkeit von (4.28).

4.4 Ortsraumrenormierung und Skalengesetze

4.4.1 Ortsraumrenormierung

Eine sehr leistungsfahige Methode in der Analyse kritischer Systeme ist dieRenormierung. Sie beruht auf der Annahme, dass kritische Systeme skalenin-

variant sind, dass sie also auf verschiedenen Langenskalen gleich aussehen. Furdiese Annahme gibt es eine Reihe von Hinweisen:

• In der Umgebung des kritischen Punktes gelten Potenzgesetze, die erstjenseits der Korrelationslange zusammenbrechen. Potenzgesetze bedeutenaber Skaleninvarianz.

• Die Werte der kritischen Exponenten werden durch makroskopische Ei-genschaften des Systems bestimmt und hangen nicht von mikroskopischenDetails ab. Das bedeutet, dass auf verschiedenen Langenskalen dieselbenMechanismen fur die kritischen Eigenschaften des Systems verantwortlichsind.

Bei der Renormierungsmethode, die in diesem Abschnitt vorgestellt wird,der Ortsraumrenormierung, wird die Langenskala des Systems geandert. Ska-leninvarianz bedeutet, dass das daraus resultierende System im Wesentlichenidentisch mit dem ursprunglichen ist. Diese Bedingung liefert Fixpunkte derRenormierungstransformation, die kritische Punkte sind oder ganz homogeneSysteme. Wir werden dies nun am Beispiel des zweidimensionalen Isingmodellsvorfuhren.

Wir beginnen mit dem in Abbildung 4.8 gezeigten quadratischen Gitter.Die Zustandssumme ist

Z ≡∑

σi

e−H/kBT =∑

σi

e(J/kBT )P

<ij> σiσj ≡ eKP

<ij> σiσj ,

wobei die Summe in der Exponentialfunktion uber nachste-Nachbar-Paare <ij > genommen wird. Wir andern die Langenskala des Systems, indem wir dieSpur uber die Spins auf den hell gezeichneten Gitterplatzen nehmen. Es ent-steht eine Zustandssumme mit einer neuen Hamiltonfunktion, die nur noch vonden verbleibenden Spins abhangt. Das verbleibende Gitter hat eine Gitterkon-stante, die um den Faktor

√2 großer ist. Skaleninvarianz bedeutet, dass nach

wiederholter Anwendung dieser Transformation (nachdem die mikroskopischenDetails eliminiert sind) die im Exponenten von Z stehende Hamiltonfunktionsich nicht weiter andert.

71

σ

σ01

2

3

4

σ

σ σ

Abbildung 4.8: Dezimierung auf dem Quadratgitter. Uber die Spins auf den hellgezeichneten Platzen wird die Spur gebildet. Das neue, gestrichelte Gitter hateine Gitterkonstante, die um den Faktor

√2 großer ist als die des alten Gitters.

Wir greifen einen Spin, σ0, heraus und untersuchen explizit, welche Termedurch seine Elimination entstehen: Der Beitrag zu Z, der σ0 enthalt, ist derFaktor

eKσ0(σ1+σ2+σ3+σ4) .

Spurbildung uber die beiden Einstellungen dieses Spins liefert

A ≡ eK(σ1+σ2+σ3+σ4) + e−K(σ1+σ2+σ3+σ4) = 2 cosh(K(σ1 + σ2 + σ3 + σ4)) .

Dieser Ausdruck kann nur drei verschiedene Werte annehmen:

(i) A1 ≡ e4K + e−4K = 2 cosh 4K ,

wenn alle vier Spins dasselbe Vorzeichen haben;

(ii) A2 ≡ e2K + e−2K = 2 cosh2K ,

wenn drei Spins dasselbe Vorzeichen haben und der vierte Spin das andere;

(iii) A3 ≡ 2 ,

wenn zwei Spins das eine Vorzeichen haben und zwei das andere.Wir wollen A nun so umschreiben, dass man die Beitrage zur neuen Hamil-

tonfunktion ablesen kann, also dass es die Form

B = ea+b(σ1σ2+σ1σ3+σ1σ4+σ2σ3+σ2σ4+σ3σ4)+cσ1σ2σ3σ4

hat. Die im Exponenten auftretenden Beitrage sind unter allen moglichen dieje-nigen, die mit den bestehenden Symmetrien in A vereinbar sind. A ist invariant,

72

wenn alle Spins geflippt werden oder wenn zwei Spins vertauscht werden, unddasselbe muss folglich fur B gelten. Die drei Gleichungen zur Bestimmung vona, b und c sind die folgenden:

A1 = ea+6b+c ,

A2 = ea−c ,

A3 = ea−2b+c .

Die Division A1/A3 liefert

b =1

8ln cosh 4K .

A1A33A

−42 liefert

c =1

8(ln cosh 4K − 4 ln cosh 2K) .

A1A33A

42 liefert

a =1

8(ln cosh 4K + 4 ln cosh 2K + ln 128) .

Elimination der Spins auf den hell gezeichneten Gitterplatzen erzeugt also nichtnur Kopplungen zwischen nachsten Nachbarn (wie σ1 und σ2) auf dem neuenGitter, sondern auch zwischen ubernachsten Nachbarn (wie σ1 und σ3), undVierspinprodukte. Man kann sich vorstellen, dass beim nachsten Eliminations-schritt noch mehr neue Kopplungen und noch hohere Spinprodukte erzeugtwerden. In einer exakten Rechnung mussten all diese Kopplungen berucksich-tigt werden. An dieser Stelle wird nun eine Naherung gemacht. Im Folgen-den beschranken wir uns auf die Zweispinkopplungen zwischen nachsten undubernachsten Nachbarn und vernachlassigen alle hoheren Kopplungen. In unse-rem Beispiel kann man nachrechnen, dass der erzeugte Beitrag c zur Vierspin-kopplung am kritischen Punkt viel kleiner ist als die Zweispinkopplung b, so dassihre Vernachlassigung gerechtfertigt werden kann. Durch diese Naherung erhal-ten wir auch nur genaherte Werte fur den kritischen Punkt und die kritischenExponenten. Wir schreiben also

H/kBT = −K∑

NN

σiσj − L∑

uNN

σiσj

fur die Hamiltonfunktion vor dem Eliminierungsschritt und

H′/kBT = −K ′∑

NN

σiσj − L′∑

uNN

σiσj

fur die Hamiltonfunktion nach dem Eliminierungsschritt. NN steht fur nachste-Nachbar-Paare auf dem betrachteten Gitter, und uNN fur ubernachste-Nachbar-Paare. Im Prinzip mussen wir nun die obige Rechnung unter Berucksichtigungder Kopplung L wiederholen und nach der Elimination samtlicher Spins auf denhellen Gitterplatzen den Wert von K ′ und L′ ablesen. Aber auch diese Rech-nung ist noch recht kompliziert, und wir machen zwei weitere Vereinfachungen:Erstens entwickeln wir in K, so dass in fuhrender Ordnung b = K2 ist. Zweitenslesen wir direkt aus Abbildung 4.8 ab, welches die wichtigsten Beitrage zu K ′

und L′ sind: Jede Kopplung K ′ bekommt zweimal einen Beitrag b, da ja auf

73

beiden Seiten ein Nachbar eliminiert wird. Außerdem gibt es noch einen BeitragL, da die ubernachste-Nachbar-Kopplung nach der Elimination zur nachsten-Nachbar-Kopplung geworden ist. Also ist

K ′ ≃ 2K2 + L

L′ ≃ K2 . (4.30)

Diese Transformation hat drei Fixpunkte (K∗, L∗):

1. K∗ = L∗ = 0: Dieser Fixpunkt beschreibt ein System aus ungekoppeltenSpins. Solch ein System liegt bei unendlicher Temperatur vor.

2. K∗ = L∗ = ∞: Dies bedeutet T = 0. Die Kopplungen sind so stark, dassalle Spins parallel sind.

3. K∗ ≡ K∗c = 1/3, L∗ ≡ L∗

c = 1/9: Dies ist der kritische Punkt.

In Abbildung 4.9 ist das aus (4.30) resultierende Flussdiagramm dargestellt.Ein Pfeil verbindet einen Punkt (K, L) mit dem sich nach zwei Iterationen erge-

Abbildung 4.9: Das aus (4.30) resultierende Flussdiagramm.

benden Punkt. (Punkte nach jeder Iteration zu zeichnen, wurde das Diagrammunubersichtlich machen, da aufgrund eines negativen Eigenwertes (s.u.) Oszilla-tionen auftreten.) Man kann in diesem Diagramm alle drei Fixpunkte erkennen.Man sieht auch, dass die beiden trivialen Fixpunkte attraktiv sind, d.h. dassdie Iteration die Kopplungen immer naher zu diesem Fixpunkten bringt. Derkritische Punkt ist instabil, denn selbst wenn die Kopplungen anfangs in seinerNahe sind, laufen sie irgendwann weg von ihm. Dieses Ergebnis macht Sinn:Je mehr Iterationen wir machen, auf desto großeren Langenskalen betrachtenwir das System. Wenn das System nicht genau am kritischen Punkt ist, hat

74

es eine endliche Korrelationslange, und auf genugend großen Skalen ist das Sy-stem nicht mehr kritisch. Befindet man sich in der Hochtemperaturphase, sindauf genugend großen Abstanden die Spinorientierungen unabhangig voneinan-der, und man ist folglich am Hochtemperaturfixpunkt K = L = 0. Befindetman sich in der Tieftemperaturphase, sind auch auf großen Abstanden die Spi-norientierungen gekoppelt. Dass die Kopplungsstarke unter der Iteration (4.30)gegen Unendlich geht bedeutet, dass Abweichungen von der Vorzugsspinorien-tierung nur auf kurzen Skalen (in Form von kleinen “Tropfchen” der falschenOrientierung) auftreten und man sie folglich auf genugend großen Skalen nichtmehr sieht.

Um die kritischen Exponenten zu erhalten, linearisieren wir die Rekursi-onsrelationen (4.30) in der Umgebung des kritischen Punktes, K = K∗

c + δK,L = L∗

c + δL, und erhalten

(

δK ′

δL′

)

=

(

4K∗c 1

2K∗c 0

)(

δK ′

δL′

)

=

(

4/3 12/3 0

)(

δK ′

δL′

)

Die Eigenwerte der Transformationsmatrix sind

λ1 =1

3(2 +

√10) ≃ 1.7208

und

λ2 =1

3(2 −

√10) ≃ −0.3874 ,

und die Eigenvektoren sind ~v1 = (1, (√

10 − 2)/3) und ~v2 = (1, (−√

10 − 2)/3).Da |λ2| < 1 ist, wird nach einigen Iterationen der Vektor (δK, δL) parallel zu~v1. Wir haben dann

(δK ′, δL′) = λ1(δK, δL) .

Den Exponenten ν erhalten wir, wenn wir uns bewusst machen, dass eineIteration die Korrelationslange ξ auf ξ′ = ξ/

√2 verkurzt, und dass δK ≃

−J(T − Tc)/kBT 2c ist. Also ist (δK ′)−ν = (δK)−ν/

√2, und

ν = ln 2/2 lnλ1 ≃ 0.638 .

Zum Vergleich: der exakte Wert von ν ist 1.Betrachten wir noch den Effekt eines endlichen Magnetfeldes auf die Orts-

raumrenormierung, und addieren einen Term −h∑

i σi zu H/kBT . Eine genaher-te Rekursionsrelation kann wieder intuitiv aufgestellt werden. Auf die verblei-benden Spins wirkt direkt das Feld h und durch die orientierende Wirkung desFeldes auf die eliminierten Nachbarspins ein Zusatzfeld Kh, also insgesamt

h′ = h + Kh .

Die Fixpunktwerte dieser Rekursionsrelation sind h∗ = 0 und und h∗ = ∞. Inder Umgebung des kritischen Punktes ergibt sich

h′ = (1 + K∗c )h =

4

3h .

Wir haben also einen weiteren Eigenwert λh = 4/3 der linearisierten Rekur-sionsrelationen. Da λh > 1 ist, wachst ein anfanglich kleines Feld an unter

75

Renormierung, ebenso wie das mit |T − Tc| passiert. Auch wenn man auf klei-nen Langenskalen ein schwaches Magnetfeld kaum wahrnimmt, sieht man aufgroßeren Skalen, dass die Spins eine Vorzugsrichtung haben. Im nachsten Ab-schnitt werden wir sehen, wie man aus λ1 und λh die kritischen Exponentenund die Skalenrelationen zwischen ihnen bestimmen kann.

Die in diesem Abschnitt vorgerechnete Ortsraumrenormierung ist nur einevon vielen moglichen. Andere Ortsraumrenormierungen fuhren Blockspins ein,d.h. sie fassen eine Gruppe von Spins zusammen und ersetzen sie durch einenneuen Spin, der die Orientierung der Mehrzahl der Spins in der Gruppe hat. Dieneuen Kopplungen zwischen den neuen Spins werden naherungsweise berechnet.Wieder andere Renormierungen vermeiden das Erzeugen der hoheren Kopplun-gen, indem sie zuerst Verbindungen zwischen Spins parallel verschieben so dassjeder zu eliminierende Spin nur noch zwei Nachbarn hat, und dann die Spuruber die zu eliminierenden Spins bilden. Dieses Verfahren ist die sogenannteMigdal-Kadanoff-Naherung und ist in Abbildung 4.10 illustriert. Alle diese Re-

Abbildung 4.10: Die Schritte in der Migdal-Kadanoff-Renormierungstrans-formation fur ein Quadratgitter. (a) Das ursprungliche Gitter mit der Kopp-lungskonstanten J . (b) Die Halfte der Verbindungen werden enfernt, um dieLangenskala des Gitters um den Faktor 2 zu vergroßern. Zum Ausgleich fur dieentfernten Kopplungen wird jede der verbleibenden Kopplungen verdoppelt. (c)Die mit einem Kreuz markierten Platze werden durch Spurbildung uber die aufihnen sitzenden Spins entfernt, und es resultiert eine neue KopplungskonstanteJ ′.

normierungsverfahren sind um so besser, je kleiner die Terme sind, die von ihnenvernachlassigt werden. Da aber in vielen Fallen die Große der vernachlassigtenTerme schwer abzuschatzen ist, ist die Ortsraumrenormierung oft unkontrol-liert. Eine alternative Methode, die Impulsraumrenormierung, wird daher meistbevorzugt.

Ubungsaufgaben

1. Bestimme den aus (4.30) resultierenden Wert fur J/kBTc fur das zweidi-mensionale Isingmodell. Zum Vergleich: der exakte Wert von J/kBTc furdas zweidimensionale Isingmodell ist 0.4406.

2. Was passiert, wenn man in der in der Vorlesung vorgefuhrten Ortsraumre-normierung nicht nur die Vierspinkopplung, sondern auch die ubernachste-Nachbarkopplung vernachlassigt?

3. Fuhre eine Ortsraumrenormierung fur das eindimensionale Isingmodelldurch, indem du in der Zustandssumme uber jeden zweiten Spin die Spur

76

bildest. Welche Fixpunkte ergeben sich? Was ist ihre physikalische Bedeu-tung? Diese Ortsraumrenormierung ist exakt!

4.4.2 Skalengesetze

Obwohl die eben dargestellte Dezimierungsprozedur mit nur wenigen Parame-tern keine quantitativ befriedigenden Resultate liefert und auch fur die Berech-nung von Korrelationsfunktionen ungeeignet ist, zeigt sie doch die allgemeineStruktur von Renormierungsgruppentransformationen, die wir jetzt benutzen,um die Skalengesetze abzuleiten.

Eine allgemeine RenormierungsgruppentransformationR bildet die ursprung-liche Hamiltonfunktion H auf eine neue ab

H′ = RH .

Wir haben am Beispiel der Ortsraumrenormierung gesehen, dass die neue Ha-miltonfunktion neue Kopplungen und Felder hat. Außerdem beinhaltet dieseTransformation auch die Reskalierung der Langen um einen Faktor b, der inunserem Beispiel b =

√2 war. Die Zahl der Freiheitsgrade (z.B. Gitterplatze)

reduziert sich dabei um einen Faktor b−d. In unserem Beispiel war dieser Faktor1/2, da jeder zweite Spin eliminiert wurde. Die Fixpunkt-Hamiltonfunktion istdurch

RH∗ = H∗

bestimmt, und in ihrer Umgebung gelten linearisierte Rekursionsrelationen

LδH = δH′ .

Die Eigenoperatoren dieser linearen Transformation sind durch die Eigenwert-gleichung

LδHi = λiδHi

bestimmt. In der Nahe des Fixpunktes konnen wir schreiben

H = H∗ + tδHt + hδHh +∑

i≥3

ciδHi .

Die beiden Storungen δHt und δHh haben Eigenwerte, deren Betrag großer als1 ist, |λt| ≡ byt > 1 und |λh| ≡ byh > 1, und sie sind mit der Temperatur t =(T −Tc)/Tc und dem Feld h assoziiert. Man nennt sie relevante Storungen. Alleanderen Storungen haben negative Eigenwerte und sind irrelevant. Wir habendies explizit am Beispiel der Ortsraumrenormierung des Isingmodells gesehen,aber es scheint allgemein fur kritische Phanomene in der Gleichgewichtsphysikzu gelten.

Fur die freie Energie pro Teilchen gilt dann in der Umgebung des kritischenPunktes nach n Iterationen

FN (t, h, c3, . . . )/N ≡ f(t, h, c3, . . . ) = b−dnf(tbytn, hbyhn, c3by3n, . . . ) . (4.31)

Hier haben wir konstante additive Terme, die bei jeder Iteration entstehen,weggelassen, da sie die folgende Herleitung der Skalengesetze nicht beeinflussen.Beziehung (4.31) gilt fur jeden Wert von t, h und n, solange man nahe genug

77

am kritischen Punkt ist. Wir konnen nun n so wahlen, dass |t|bytn = 1 ist. Dannergibt sich

f(t, h, c3, . . . ) = |t|d/yt f±(h|t|−yh/yt , . . . ). (4.32)

Nahe bei Tc kann die Abhangigkeit von c3 etc. vernachlassigt werden. Die Punktein (4.32) stehen also fur irrelevante Beitrage. Wir konnen nun die kritischenExponenten als Funktion von yt und yh bestimmen: Die spezifischen Warmeberechnet sich durch zweimaliges Ableiten von (4.32) nach t bei h = 0, woraus

α = 2 − d/yt

resultiert. Der spontane Ordnungsparameter ergibt sich durch einmaliges Ab-leiten nach h und der Ersetzung h = 0, woraus

β =d − yh

yt

folgt. Zweimaliges Ableiten nach h gibt die Suszeptibilitat und

γ =2yh − d

yt

.

Den Exponenten δ erhalten wir am einfachsten, wenn wir in (4.31) die Ersetzungbn = h−1/yh machen, was auf

f(t, h, c3, . . . ) = hd/yh f(th−yt/yh , . . . )

fuhrt. Ableitung nach h bei t = 0 gibt

m ∼ h(d−yh)/yh

undδ =

yh

d − yh

.

Wir hatten schon im Zusammenhang mit der Ortsraumrenormierung argu-mentiert, dass

ν = 1/yt

ist. Eine Beziehung fur den Exponenten η kann man in Verbindung mit derWilson-Renormierung ableiten.

Aus der Skalenform (4.32) fur die freie Energie erhalten wir durch Ableitennach dem Feld h eine Skalenform fur den Ordnungsparameter:

m(t, h) = |t|(d−yh)/yt f±(h|t|−yh/yt) = |t|β f±(h|t|−βδ) . (4.33)

Wenn man den Ordnungsparameter m als Funktion des Feldes h in der Nahedes kritischen Punktes misst und m/|t|β gegen h/|t|βδ auftragt, fallen alle Datenfur t > 0 auf eine Kurve und alle Daten fur t < 0 auf eine andere Kurve. Kurvenm(h), die man fur verschiedene Temperaturen gemessen hat, konnen also mitdieser skalierten Auftragung zum Zusammenfallen gebracht werden. DerartigeSkalenkollapse sind sehr nutzlich zur Bestimmung der kritischen Exponenten,denn die Kurven fallen naturlich nur dann zusammen, wenn man die richtigenWerte fur die Exponenten verwendet.

78

Skalenformen wie (4.32) und (4.33) spielen eine zentrale Rolle in der Unter-suchung kritischer Phanomene. Wir haben in dieser Vorlesung die Skalenformenzwar aus der Renormierung hergeleitet, aber sie waren schon vor der Entwick-lung der Renormierungsgruppentheorie bekannt und waren phanomenologisch,mit guter Intuition, aufgestellt worden. Solche phanomenologischen Skalenuber-legungen sind insbesondere dann wichtig, wenn man physikalische Situationenuntersucht, fur die es noch keine renormierungsgruppentheoretische Behandlunggibt.

Ubungsaufgaben

1. Ausgehend von (4.33) beantworte folgende Fragen:(a) Wie hangt m von t ab, wenn h = 0 ist?(b) Wie hangt m von h ab, wenn t = 0 ist?

(c) Wie ist also das asymptotische Verhalten von f±(x) fur große x undfur kleine x?

2. Leite aus der freien Energie der Landau-Theorie f = rm2+um4−hm eineSkalenform fur m ab und lies daran die Exponenten β und δ ab. Wenndu mit dem Computer umgehen kannst, dann trage in einem Diagrammm(h) fur verschiedene (kleine) t auf, und in einem anderen Diagramm denSkalenkollaps. Dies muss fur t > 0 und t < 0 getrennt gemacht werden.

4.5 Perkolation

4.5.1 Definition und Beispiele

Perkolation ist ein kritsches Phanomen, das rein geometrischer Natur ist, aberviele Anwendungen in physikalischen Systemen hat. Sie lasst sich am einfachstenmithilfe eines Modells auf einem Gitter veranschaulichen, wie in Abbildung 4.11dargestellt.

Jeder Gitterplatz ist mit einer Wahrscheinlichkeit p besetzt und mit Wahr-scheinlichkeit 1 − p leer. p ist der einzige Parameter des Systems. Die Perkola-tionstheorie befasst sich mit Clustern besetzter Platze.

Alle besetzten nachsten Nachbarn eines besetzten Platzes gehoren zu demsel-ben Cluster wie dieser Platz. In Abbildung 4.11 sind ein Cluster aus drei Platzenund ein Cluster aus 18 Platzen gekennzeichnet. Unter der Große s eines Clu-sters versteht man die Zahl der besetzten Platze, die zu diesem Cluster gehoren.Solange p kleiner ist als ein kritischer Wert pc, der auf dem Quadratgitter denWert pc ≃ 0.5927 hat, gibt es selbst auf einem unendlich großen Gitter nurCluster endlicher Große. Jenseits von pc gibt es einen unendlichen Cluster, dervon einem Ende des Systems zum anderen Ende reicht. Den kritischen Wert pc

nennt man die Perkolationsschwelle. Wie bei anderen kritischen Phanomenen,so ist auch der Perkolationsubergang nur im thermodynamischen Limes unend-licher Gittergroße ein scharfer Ubergang. Kleinere Gitter haben manchmal auchunterhalb der Perkolationsschwelle schon einen Cluster, der von einem Endezum anderen reicht, und sie haben manchmal oberhalb der Perkolationsschwellekeinen solchen Cluster.

Neben der eben geschilderten Besetzungsperkolation, bei der Gitterplatzebesetzt oder leer sind, gibt es auch die Verbindungsperkolation, bei der Ver-

79

bindungen zwischen Nachbarplatzen entweder geschlossen oder geoffnet sind.Wenn der Prozentsatz geoffneter Verbindungen genugend groß ist, gibt es einenPfad durch das ganze System, der uber geoffnete Verbindungen lauft. Beide Ar-ten Perkolation haben dasselbe kritische Verhalten, und wir werden deshalb imfolgenden nur die Besetzungsperkolation beschreiben.

Fur Perkolation gibt es zahlreiche Beispiele in der Natur. Man denke anein Feuer, das sich von einem brennenden Baum zu Baumen, die genugend na-he sind, ausbreiten kann. Wenn die Baumdichte gering ist, wird das Feuer baldverloschen. Wenn sie groß ist, kann es sich durch den gesamten Wald ausbreiten.Ahnliches gilt fur die Ausbreitung von Krankheiten, wenn man die Baume durchIndividuen ersetzt, die nicht gegen die Krankheit immun sind. Ein weiteres Bei-spiel fur Perkolation ist das Eindringen einer Flussigkeit in ein poroses Medium.Wenn es genugend Verbindungen zwischen den Poren gibt, durch die die Flussig-keit fließen kann, kann die Flussigkeit tief in das Medium eindringen. Andernfallsbleibt sie nahe der Oberflache. Auch in Materialien, die sich aus elektrisch lei-tenden und nichtleitenden Elementen zusammensetzten, tritt Perkolation auf.Wenn der Prozentsatz der leitenden Elemente oberhalb eines Schwellwertes ist,kann ein Strom durch das Material fließen. Ein Beispiel fur Verbindungsper-kolation ist der Sol-Gel-Ubergang. Solange die Molekule eines solchen Systemsnur wenige Verbindungen untereinander haben, konnen sich sich gegeneinanderbewegen. Wenn die Dichte der Verbindungen so hoch geworden ist, dass sich dasVerbindungsnetzwerk uber das gesamte System erstreckt, ensteht ein Gel. Diespassiert z.B. beim Kochen eines Eis oder eines Puddings. Allerdings lassen sichnicht alle diese Beispiele durch das einfache Modell unkorrelierter Perkolationdarstellen, da oft die Wahrscheinlichkeiten zweier Nachbarplatze oder Nachbar-verbindungen, besetzt zu sein, nicht unabhangig voneinander sind.

4.5.2 Kritische Exponenten und Skalenverhalten

Der Ordnungsparameter fur den Perkolationsubergang ist der Prozentsatz Pder Platze, die zum unendlichen Cluster gehoren. Er ist Null unterhalb pc undwachst oberhalb der Perkolationsschwelle wie

P ∼ (p − pc)β . (4.34)

Der Wert von β auf dem zweidimensionalen Gitter ist 5/36.Wir definieren die Clustergroßenverteilung np(s) als die Zahl der Cluster

der Große s, geteilt durch die Gesamtzahl der Gitterplatze. Den unendlichenCluster, der oberhalb pc vorhanden ist, nehmen wir aus dieser Verteilung heraus.Ein allgemeiner Skalenansatz fur die Clustergroßenverteilung nahe pc ist

np(s) = s−τC±(s|p − pc|1/σ) . (4.35)

Bei p = pc bedeutet dies, dass die Clustergroßenverteilung ein Potenzgesetzist. Dies ist auch tatsachlich der Fall, wenn man die ganz kleinen Cluster nichtberucksichtigt. Auf kleinen Langenskalen spurt man noch die mikroskopschenDetails des Gitters, aber auf großeren Langenskalen liegt ein Potenzgesetz vor.Wenn p nicht genau dem kritischen Wert pc entspricht, gibt es einen Cutoffzu diesem Potenzgesetz. Die Funktion C±(x) ist konstant fur kleine x und falltfur große x exponentiell ab wie e−ax. Die Clustergroße, bei der dieser Cutoffeinsetzt, ist

smax ∼ |p − pc|−1/σ .

80

Auf dem zweidimensionalen Gitter ist τ = 187/91 und σ = 36/91. Alle dieseWerte sind aus analytischen Berechnungen bekannt.

Die Besetzungswahrscheinlichkeit p lasst sich durch die Clustergroßenvertei-lung ausdrucken und ist

p =∑

s

snp(s) + P .

Der Ordnungsparameter P lasst sich durch die Clustergroßenverteilung aus-drucken, wenn wir davon ausgehen, dass alle Cluster, die jenseits des Cutoffssmax zu liegen kamen, Teil des unendlichen Clusters sind:

P =∑

s

s(

s−τC+(0) − s−τC+(s|p − pc|1/σ))

≃∫ ∞

1

s1−τ(

C+(0) − C+(s|p − pc|1/σ))

ds

≃ [(p − pc)1/σ]τ−2

∫ ∞

0

x1−τ (C+(0) − C+(x)) dx

∼ (p − pc)(τ−2)/σ .

Damit haben wir die Skalenrelation

β = (τ − 2)/σ .

Die mittlere Große des Clusters, in dem sich ein willkurlich ausgewahlterbesetzter Platz befindet, ist (wenn man den unendlichen Cluster ausschließt)

S =

∑

s s2np(s)∑

s snp(s)=

1

p − P

∑

s

s2np(s)

≃∫

s2s−τC±(s|p − pc|1/σ)ds

∼ |p − pc|(τ−3)/σ ≡ |p − pc|−γ .

Der Exponent γ hangt also in der Perkolationstheorie mit der mittleren Großeendlicher Cluster zusammen.

Den Radius Rs eines Clusters der Große s definieren wir durch die Beziehung

R2s =

s∑

i=1

|~xi − ~x0|2/s

mit der Schwerpunktskoordinate

~x0 =s∑

i=1

~xi/s .

Der mittlere quadratische Abstand zweier Platze in demselben Cluster ist

∑

ij

|~xi − ~xj |2/s2 = 2R2s .

Dieses Ergebnis erhalt man am einfachsten, wenn man den Schwerpunkt in denUrsprung des Koordinatensystems setzt, ~x0 = 0.

81

Schließlich definieren wir noch die Korrelationslange ξ als die Wurzel desmittleren quadratischen Abstands zweier Platze, die zu demselben Cluster gehoren,

ξ2 =2∑

s R2ss

2np(s)∑

s s2np(s)∼ R2

smax.

Die fraktale Dimension eines (großen) Clusters ist durch die Beziehung

Rs ∼ s1/D

definiert, und wir erhalten somit

ξ ∼ |p − pc|−ν

mitν = 1/σD .

Fur das zweidimensionale Gitter ist ν = 4/3 und D = 91/48.Ein großer fraktaler Cluster der Verbindungsperkolation ist in Abbildung

4.12 gezeigt.

4.5.3 Ortsraumrenormierung

Mit einer Ortrsraumrenormierung kann man genaherte Resultate fur die Perko-lationsschwelle und den Exponenten ν bekommen. Wir zeigen dies am Beispieleines Dreiecksgitters, Abbildung 4.13. Je drei Gitterplatze werden bei einemTransformationsschritt durch einen ersetzt. Wir mussen nun ein Kriterium fest-legen, unter welchen Umstanden der neue Platz leer oder besetzt ist. Eine sinn-volle Festlegung ist, dass der neue Platz besetzt ist, wenn es einen Pfad besetzterPlatze durch das Dreieck gibt, also wenn zwei oder drei der alten Platze besetztsind, und dass er ansonsten leer ist. Damit erhalten wir

p′ = p3 + 3p2(1 − p) .

Die Fixpunkte dieser Rekursionsrelation sind p∗ = 1 (vollig besetztes Gitter),p∗ = 0 (vollig leeres Gitter), und p∗c = 1/2 (Perkolationsschwelle). Der Wert furdie Perkolationsschwelle ist exakt, was aber Zufall ist, da die Rechnung keineexakte Berechnung der Perkolationsschwelle ist. Unter der Renormierung lauftder Fluss nach p∗ = 0, wenn der Anfangswert unterhalb p∗c liegt, und nachp∗ = 1, wenn der Anfangswert oberhalb p∗c liegt. Die lineariserte Rekursionlautet

δp′ =3

2δp .

Bei einem Renormierungsschritt andert sich die Gitterkonstante um den Faktorb =

√3. Damit ist

ν =ln√

3

ln 3/2≃ 1.355 .

Dieser Wert ist recht nah am exakten Wert 4/3.

82

Ubungsaufgaben

1. Perkolation in einer Dimension: Leite die folgenden Ergebnisse her:(a) pc = 1.(b) np(s) = ps(1 − p)2.(c) ν = 1.(d) γ = 1.Anmerkung: Die Skalenrelation β = (τ −2)/σ ist in einer Dimension nichterfullt, da np(s) einen Faktor (pc − p)2 enthalt, wahrend in hoheren Di-mensionen np(s) fur kleinere s (d.h. dort wo die Skalenfunktion C konstantist) nicht von p abhangt.

2. Perkolation auf dem Bethe-Gitter: Das Bethe-Gitter, auch Cayley-Baumgenannt, ist ein Gitter, das sich immer weiter verzweigt und keine geschlos-senen Wege hat. Man konstruiert es folgendermaßen: Vom Gitterplatz amUrsprung gehen z Aste aus, and deren Enden wieder Gitterplatze liegen,von denen wieder jeweils z − 1 neue Zweige ausgehen, usw., wie in Abbil-dung 4.14 dargestellt. Leite die folgenden Ergebnisse her:(a) pc = 1/(z − 1).(b) γ = 1. (Hinweis: Berechne eine Beziehung fur die mittlere Clustergroßein einem Zweig.)(c) β = 1. (Hinweis: Berechne die Wahrscheinlichkeit, dass ein gegebenerPlatz durch einen gegebenen Zweig nicht mit nach Unendlich verbundenist. Diese Rechnung geht am einfachsten fur z=3.)(d) σ = 1/2. (Hinweis: Die Rechnung geht am einfachsten fur z = 3. Zeigezuerst, dass ein Cluster der Große s genau s(z−2)+2 leere Nachbarplatzehat. Dann berechne fur z = 3 einen Ausdruck fur np(s)/npc

(s) im Grenz-fall kleiner p − pc.)(e) τ = 5/2. (Hinweis: Berechne τ aus γ und σ.)

83

Abbildung 4.11: Perkolation auf dem Quadratgitter. Bei einer Besetzungswahr-scheinlichkeit p = 0.25 gibt es nur kleine Cluster (Oben). Naher an der Per-kolationsschwelle treten schon großere Cluster auf (Mitte), und oberhalb derPerkolationsschwelle sind Platze uber das gesamte Gitter hinweg verbunden(Unten).

84

Abbildung 4.12: Ein großer Perkolationscluster.

Abbildung 4.13: Ortsraumrenormierung auf dem Dreiecksgitter. Die neuen Git-terplatze liegen im Zentrum der schraffierten Dreiecke.

Abbildung 4.14: Bethe-Gitter mit Koordinationszahl z = 3.

85

Kapitel 5

Musterbildung in raumlich

ausgedehnten Systemen

In diesem Kapitel verbinden wir Elemente der vorigen drei Kapitel. Bisher habenwir nichtlineare Dynamik nur in “einfachen” Systemen betrachtet, die wenigeVariablen enthalten und in denen die raumliche Ausdehnung keine Rolle spielt.Im letzten Kapitel haben wir zum ersten Mal raumlich ausgedehnte Systemebetrachtet, aber nur im thermischen Gleichgewicht. Das Interessanteste, was insolchen Systemen passieren kann, sind Phasenubergange. Wenn wir nun nichtli-neare Dynamik mit raumlicher Ausdehnung koppeln, werden neue, interessantePhanomene moglich. Zu ihnen gehoren Wellen und raumzeitliche Musterbil-dung. Formal manifestiert sich die raumliche Ausdehnung dadurch, dass dieVariablen nun ortsabhangig sind, und dass ihre zeitliche Anderung auch einenDiffusionsterm beinhaltet. Da Diffusion ein Hauptthema im zweiten Kapitel war,verbindet dieses Kapitel somit alle bisherigen Themen.

5.1 Wellen

In vielen Systemen entfernt vom Gleichgewicht, insbesondere auch in biologi-schen Systemen, treten Wellen auf. Beispiele sind die Ausbreitung einer Popula-tion in ein neues Gebiet, die Ausbreitung einer Epidemie und biochemische Wel-len in der Embryonalentwicklung. Wir betrachen im Folgenden solche Wellen,die ihre Form wahrend der Ausbreitung nicht andern und die daher geschriebenwerden konnen als

u(x, t) = U(x− ct) ≡ U(z) . (5.1)

Die eben genannten Beispiele konnen als “Reaktions-Diffusions-Systeme”verstanden werden. Es gibt einen Ausbreitungsmechanismus (Diffusion), der al-leine noch keine Wellen machen kann, sondern einfach eine immer breiter und fla-cher werdende Konzentration bewirkt. Außerdem gibt es einen “Reaktions”mecha-nismus, durch den die Konzentrationen u(x, t) sich lokal andern konnen. Wennwir eine eindimensionale Darstellung wahlen, haben wir dann Gleichungen derForm

∂u

∂t= f(u) +D

∂2u

∂x2. (5.2)

86

Eine eindimensionale Darstellung ist dann angemessen, wenn die Welle sich imWesentlichen in einer Richtung ausbreitet.

5.1.1 Fisher-Kolmogoroff-Gleichung

Die einfachste Gleichung der Form (5.2), die auf sich ausbreitende Wellen fuhrt,ist die Fisher-Kolmogoroff-Gleichung

∂u

∂t= ru(1 − u) +D

∂2u

∂x2. (5.3)

Der “Reaktions”term beschreibt logistisches Wachstum (3.6), d.h. die Konzen-tration (z.B. einer Population) wachst an, bis sie den Sattigungswert 1 erreicht.Der Diffusionsterm beschreibt das Abwandern von Mitgliedern der Populationin Nachbargebiete.

Durch die Variablentransformation

t→ t/r , x→ x√

D/r (5.4)

wird die Fisher-Kolmogoroff-Gleichung dimensionslos,

∂u

∂t= u(1 − u) +

∂2u

∂x2. (5.5)

Zwei offensichtliche Losungen der Fisher-Kolmogoroff-Gleichung sind diekonstanten Losungen u = 0 und u = 1. Dies sind die beiden Fixpunkte der logi-stischen Gleichung. Wir suchen im Folgenden nach Losungen, die die Form ei-ner in positiver x-Richtung wandernden Wellenfront haben. Anschaulich konnenwir uns vorstellen, dass sich eine Population in ein bisher von ihr unbesiedeltesGebiet ausbreitet. Dort, wo die Wellenfront noch nicht war, ist also die Konzen-tration u = 0, dort wo die Wellenfront schon lange vorbeigekommen ist, hat dieKonzentration den Sattigungswert 1. An der Front selbst fallt die Konzentrationmit wachsendem x von 1 auf 0 ab. Um die Form der Front zu bestimmen, setzenwir z = x− ct und u(x, t) = U(z). eingesetzt in (5.5) ergibt das

U ′′ + cU ′ + U(1 − U) = 0 . (5.6)

Nun haben wir nur noch Ableitungen bezuglich einer Variablen (z) und konnendie Gleichung (5.6) formal wie ein dynamisches System behandeln. Wir fuhrendaher V = U ′ ein und erhalten somit ein Gleichungssystem mit Ableitungenerster Ordnung:

U ′ = V

V ′ = −cV − U(1 − U) . (5.7)

Dieses System beschreibt, wie sich die Konzentration U und ihre Steigung mitwachsendem z andern. Es hat zwei Fixpunkte. Der eine Fixpunkt ist (U, V ) =(0, 0). Eine lineare Stabilitatsanalyse gibt die Eigenwerte

λ =1

2

(

−c±√

c2 − 4)

, (5.8)

und folglich ist dieser Fixpunkt fur c2 > 4 ein stabiler Knoten und fur c2 < 4eine stabile Spirale. Da nur positive U -Werte Sinn machen (denn U ist eine

87

Konzentration), mussen wir die Spirale ausschließen, denn sie wurde bedeuten,dass das vordere Ende der Front (wo die Konzentration noch in der Nahe von 0ist) um die u-Achse oszilliert, mit abnehmender Amplitude bei wachsendem z.Wir folgern daher, dass c ≥ 2 sein muss, und dass der erste Fixpunkt ein stabilerKnoten ist. (Naturlich gibt es auch die Losung c ≤ −2, und sie entspricht einernach links laufenden Wellenfront. Sie ergibt sich aus der nach rechts laufendenLosung durch Spiegelung z → −z, sie bringt also nichts Neues.)

Der zweite Fixpunkt von (5.7) ist (U, V ) = (1, 0). Er hat die Eigenwerte

λ =1

2

(

−c±√

c2 + 4)

, (5.9)

d.h. er ist ein Sattelpunkt. Aus der Gleichung U ′ = V folgern wir, dass Tra-jektorien in der (U, V )-Ebene oberhalb der U -Achse nach rechts und unterhalbnach links laufen. Damit wissen wir, wie das Phasenportrait qualitativ ausse-hen muss (Abb. 5.1). Diejenige Trajektorie, die die beiden Fixpunkte verbindet,

0 U

V

1

Abbildung 5.1: Qualitatives Phasenportrait fur das System (5.7).

beschreibt die gesuchte Wellenfront. Sie fallt mit wachsendem z kontinuierlichvon U = 1 auf U = 0 ab. Die genaue Form von U(z) muss numerisch bestimmtwerden, und sie hangt von c ab. Welche Geschwindigkeit die Wellenfront hat,hangt von der Anfangsbedingung u(x, 0) fur x → ±∞ ab. Um dies zu sehen,betrachten wir die Anfangsbedingung

u(x, 0) ∼ Ae−ax (5.10)

fur x → ∞. Also ist dort u(x, t) = Ae−a(x−ct) fur Zeiten t > 0. Da u in diesemBereich von x-Werten sehr klein ist, konnen wir den quadratischen Term in (5.5)vernachlassigen und erhalten sofort einen Zusammenhang zwischen c und a:

c = a+1

a. (5.11)

Der kleinstmogliche Wert c = 2 wird fur a = 1 angenommen. Dies ist die Ge-schwindigkeit, die sich nach einiger Zeit einstellt, wenn am Anfang u = 0 istfur alle x > x0 fur ein beliebiges x0. Der Beweis hierzu ist aufwendig. Aberman kann sich diesen Sachverhalt folgendermaßen plausibel machen: Am Uber-gang bei x0 zu u = 0 hat man a = ∞. Dieser Teil der Front bewegt sich also

88

sehr schnell nach rechts und sorgt auf diese Weise dafur, dass der Knick weni-ger scharf wird. Nach einer Weile hat der vorderste Teil der Front dann einenkleineren Wert von a. Auch dieser Teil wandert schnell nach rechts, und dernun dominante Teil hat ein noch kleineres a, ...., bis der Wert a = 1 und dielangsamste Geschwindigkeit erreicht sind.

5.1.2 Storungstheorie fur die Form der Front

Auch wenn wir die Form der Front nicht analytisch berechnen konnen, konnenwir sie zumindest naherungsweise bekommen. Es gibt namlich eine kleine Große,in der wir entwickeln konnen. Um dies zu sehen, machen wir die Ersetzungenz = cξ und U(z) = g(ξ) in (5.6) und erhalten

ǫd2g

dξ2+dg

dξ+ g(1 − g) = 0 (5.12)

mit ǫ = 1/c2 ≤ 0.25. Die Losung ist eindeutig, wenn wir die Bedingungeng(−∞) = 1, g(∞) = 0 und g(0) = 1/2 auferlegen. Die letzte Bedingung legt diePosition der Front langs der ξ-Achse fest.

Da ǫ maximal den Wert 0.25 hat, behandeln wir es als eine kleine Große, inder wir entwickeln. Wir machen deshalb den Ansatz

g(ξ, ǫ) = g0(ξ) + ǫg1(ξ) + . . . . (5.13)

Solch ein Ansatz funktioniert nur, wenn die Gleichung mit ǫ = 0, die nun eineGleichung erster Ordnung ist, die vorgegebenen Randbedingungen erfullen kann.Dies funktioniert in der Tat, wie man durch explizites Berechnen von g0 findet:In Ordnung ǫ0 haben wir die Gleichung

dg0dξ

= −g0(1 − g0) . (5.14)

Durch Separation der Variablen erhalt man

g0(ξ) =A

1 +Aeξ, (5.15)

und aus der Bedingung g0(0) = 1/2 folgt A = 1. In Ordnung ǫ erhalten wir

dg1dξ

+ (1 − 2g0)g1 = −d2g0dξ2

, (5.16)

was wir umformen konnen zu

g′1 −g′′0g′0g1 = −g′′0 . (5.17)

Wenn man diese Gleichung genugend lange anschaut, erkennt man, dass derAnsatz g1 = g′0f mit einer zu bestimmenden Funktion f weiterhelfen konnte.Eingesetzt ergibt dies

f ′ = −g′′0/g′0woraus wir f = A − ln(−g′0) erhalten mit einer Konstanten A. Das Vorzeichenhinter dem Logarithmus ergibt sich aus der Einsicht, dass g′0 negativ ist. Die

89

Randbedingungen sind g1(−∞) = g1(0) = g1(∞) = 0, woraus A = − ln 4 folgt.Also haben wir schließlich

g1 = −g′0 ln(−4g′0) =eξ

(1 + eξ)2ln

[

4eξ

(1 + eξ)2

]

. (5.18)

Ausgedruckt durch die ursprunglichen Variablen erhalten wir dann in ersterOrdnung in ǫ den folgenden Ausdruck fur die Front:

U(z, ǫ) =1

1 + ez/c+

ez/c

c2(1 + ez/c)2ln

[

4ez/c

(1 + ez/c)2

]

+ O(c−4) . (5.19)

Selbst fur c = 2 ist die Naherung nullter Ordnung (also der erste Term) nurwenige Prozent von der exakten Losung entfernt, und die Naherung bis zurersten Ordnung ist somit erst recht erstaunlich gut.

5.1.3 Stabilitat der Front

Um die Stabilitat einer Front zu bestimmen, gehen wir in das mitbewegte Be-zugssystem und schreiben u(x, t) = u(z, t) mit z = x − ct. Dann ist die Bewe-gungsgleichung fur die Front

ut = u(1 − u) + cuz + uzz . (5.20)

Nun betrachten wir eine kleine Storung um die Front zur Geschwindigkeit c,

u(z, t) = uc(z) + ωv(z, t) (5.21)

mit einem kleinen Parameter ω. Wir gehen davon aus, dass die Storung v nurauf einem endlichen Intervall [−L,L] von Null verschieden ist. Fur v erhaltenwir durch Einsetzen die Gleichung

vt = [1 − 2uc(z)]v + cvz + vzz . (5.22)

Wenn v fur t → ∞ gegen 0 oder gegen duc(z)/dz geht, ist die Front stabilgegenuber der Storung. (Der zweite Fall bedeutet namlich nur eine Verschiebungder Front um ω, nicht eine Anderung ihrer Form.) Wir losen die Gleichung furv mit denselben Methoden wie die Schrodingergleichung. Wir machen zunachstden Ansatz

v(z, t) = g(z)e−λt (5.23)

und eliminieren die erste Ableitung mit dem Ansatz

g(z) = h(z)e−cz/2 . (5.24)

Dann bleibt fur h die Differenzialgleichung

−h′′ +[

2uc(z) +c2

4− 1

]

h = λh , (5.25)

die die Form einer zeitunabhangigen Schrodingergleichung hat. Das ”Potenzial”dieser ”Schrodingergleichung” ist uberall ≥ 0. Es gibt daher keine Zustandemit negativer Energie. Also sind alle Eigenwerte λ positiv, und aus (5.23) folgtsomit, dass v mit der Zeit auf 0 abklingt.

90

Aufgabe

Gegeben ist das Gleichungssystem

∂U

∂t= AU

(

1 − U

K

)

−BUV

∂V

∂t= CUV −DV +D2∇2V , (5.26)

das die Populationen einer Nahrungspflanze (U) und eines Pflanzenfressers (V )beschreibt. Alle Konstanten A,B,C,D,D2 seien positiv.

1. Interpretieren Sie die Terme und bringen Sie durch Umskalieren das Glei-chungssystem auf eine Form mit nur zwei Parametern.

2. Finden Sie die konstanten Losungen und interpretieren Sie diese.

3. Im Folgenden suchen wir Wellenfronten. Fuhren Sie die Variable z = x−ctein und bestimmen Sie die resultierenden Differenzialgleichungen. BringenSie diese auf die Form eines dynamischen Systems (also nur Gleichungenmit ersten Ableitungen.)

4. Uberlegen Sie anschaulich, welche Art von Wellen es geben konnte (alsowelche konstanten Werte von U und V bei z = ±∞). Wahlen Sie diejenige,die an keinem der beiden Enden die weniger interessante Losung U = V =0 hat.

5. Dann suchen Sie die Fixpunkte des dynamischen Systems und machen Sieeine lineare Stabilitatsanalyse in der Umgebung der beiden gewunschtenFixpunkte. Was ergibt sich daraus fur die minimale Ausbreitungsgeschwin-digkeit der Front? In der Umgebung des einen Fixpunktes bekommen Sieein Polynom dritten Grades p(λ) als Bestimmungsgleichung fur die Ei-genwerte. Betrachten Sie dieses Polynom zunachst fur den Fall, dass derkonstante Term verschwindet und skizzieren Sie den Verlauf des Polynomsqualitativ. Aufgrund dieser Skizze konnen Sie das Polynom auch fur allge-meine Werte des konstanten Termes zeichnen. Leiten Sie daraus ab, dasses Fronten geben kann, die monoton in den Wert bei z = −∞ einmunden,als auch solche, die oszillierend in ihn einmunden. Skizzieren Sie quali-tativ, wie diese Fronten aussehen sollten. Warum darf die Losung in derUmgebung des einen Fixpunktes oszillieren, aber nicht in der Umgebungdes anderen Fixpunktes?

6. Bestimmen Sie numerisch am Computer Beispiele fur diese beiden Losungs-typen, indem Sie jeweils geeignete Parameterwerte und die zugehorige mi-nimale Frontgeschwindigkeit c wahlen.

5.2 Turing-Muster

5.2.1 Allgemeiner Formalismus

In der Natur gibt es Prozesse, durch die ein raumliches Muster einmal entstehtund dann lange so bleibt. Beispiele sind Tierfellmuster oder Schmetterlingsflugel.

91

Auch wenn man bis heute die genauen Prozesse, durch die diese Muster entste-hen, nicht aufgeklart hat, gibt es uberzeugende Modelle, die derartige Musterproduzieren konnen. Diese Modelle sind Reaktions-Diffusions-Systeme mit zwei(oder mehr) chemischen Substanzen. Der Mechanismus fur die Musterbildung indiesen Modellen ist eine diffusionsgetriebene Instabilitat. Um dies zu verstehen,betrachten wir das Gleichungssystem

u = γf(u, v) + ∇2u

v = γg(u, v) + d∇2v . (5.27)

u und v sind die Konzentrationen der beiden Substanzen. Durch eine geeig-nete Umskalierung kann man Reaktions-Diffusions-Gleichungen immer auf dieangegebene Form bringen, mit nur zwei Parametern d und γ (außer den in fund g auftretenden Parametern). d ist das Verhaltnis der beiden Diffusions-koeffizienten. γ ist die relative Starke der Reaktionsterme verglichen mit denDiffusionstermen. Wir werden sehen, dass die Funktionen f und g bestimmteBedingungen erfullen mussen, damit sich Muster bilden konnen. Wir gehen da-von aus, dass es genau einen homogenen stationaren Zustand (u0, v0) gibt. Es istalso f(u0, v0) = g(u0, v0) = 0. Außerdem verlangen wir, dass dieser stationareZustand stabil ist gegenuber homogenen Storungen. Die zum Fixpunkt (u0, v0)gehorige Stabilitatsmatrix muss also eine negative Spur und eine positive De-terminante haben, d.h.

τ0γ

=∂f

∂u

∣

∣

∣

∣

u0,v0

+∂g

∂v

∣

∣

∣

∣

u0,v0

≡ fu + gv < 0 (5.28)

und∆0

γ2=∂f

∂u

∂g

∂v

∣

∣

∣

∣

u0,v0

− ∂f

∂v

∂g

∂u

∣

∣

∣

∣

u0,v0

≡ fugv − fvgu > 0 (5.29)

Nun pragen wir diesem homogenen Zustand eine kleine inhomogene Storungauf und untersuchen, ob eine derartige Storung wieder abklingt. Wir betrachteneine Storung, die nur bezl einer Raumrichtung eine Inhomogenitat hat. Dies seidie x-Richtung. Solange solche Storungen klein sind, uberlagern sie sich linear.Wir konnen sie also in ihre Fourierkomponenten zerlegen und diese Komponen-ten getrennt untersuchen. Es genugt daher, eine raumlich periodische Storung

(

δu(t = 0)δv(t = 0)

)

=

(

A cos(kx− ϕ)B cos(kx− ϕ)

)

(5.30)

um den homogenen Fixpunkt zu betrachten. Die Phase ϕ hangt von den gewahl-ten Randbedingungen ab. Wenn am Rand kein Molekul das System verlassenkann, ist dort (~n · ∇)u = (~n · ∇)v = 0, wobei ~n senkrecht auf den Rand steht.Die erlaubten k-Werte hangen von der Ausdehnung des Systems ab.

Die Bewegungsgleichung fur solch eine periodische Storung ist

(

˙δu

δv

)

=

(

γf(u0 + δu, v0 + δv) − k2δuγg(u0 + δu, v0 + δv) − dk2δv

)

=

(

γfu − k2 γfv

γgu γgv − dk2

)(

δuδv

)

(5.31)Die Spur der Stabilitatsmatrix ist jetzt also

τ = τ0 − (1 + d)k2 (5.32)

92

und ist negativ, weil wir τ0 < 0 gefordert haben. Die Determinante der Stabi-litatsmatrix ist

∆ = ∆0 − γk2(gv + dfu) + dk4 . (5.33)

Wenn diese Determinante negativ ist, ist der homogene Fixpunkt instabil ge-genuber einer Storung der Wellenzahl k. Weil ∆0 positiv ist, ist eine notwendigeBedingung fur negatives ∆

gv + dfu > 0 . (5.34)

Zusammen mit der oben hergeleiteten Bedingung fu + gv < 0 folgt daraus, dassgv und fu das entgegengesetzte Vorzeichen haben mussen. Wir wahlen fu > 0und gv < 0. Dann muss d > 1 sein. Die zweite chemische Substanz muss alsoschneller diffundieren als die erste. Eine weitere notwendige Bedingung fur einenegatives ∆ ist

(dfu + gv)2

4d>

∆0

γ2. (5.35)

Diese Bedingung ist fur ein genugend großes Verhaltnis d > dc der Diffusions-koeffizienten erfullt. Der großere Eigenwert

λ =τ

2+

√

τ2

4− ∆ (5.36)

der Stabilitatsmatrix ist in Abbildung 5.2 skizziert. wenn d > dc ist, gibt es ein

d<dc

k

d=dd>d

λ

0

c

c

2

Abbildung 5.2: Der großere Eigenwert der Stabilitatsmatrix als Funktion derWellenzahl k2 der Storung fur verschiedene d.

Intervall von k-Werten, fur das die Storung anwachst. Sie wachst exponentiell inder Zeit, also wie eλt. Wenn wir keine kosinusformige Storung anlegen, sonderneine beliebige Storung, dann wachst diejenige Fourierkomponente am schnell-sten an, die zum großten λ gehort. Nach einiger Zeit dominiert also diejenigeWellenzahl k, die zum Maximum der Kurve λ(k2) in Abbildung 5.2 gehort.

Naturlich setzt sich diese exponenzielle Wachstum nicht unbegrenzt fort. Wirhaben ja eine lineare Storungstheorie um den homogenen Fixpunkt (u0, v0) ge-macht. Sobald die Storung nicht mehr klein ist, werden die nichtlinearen Termeder dynamischen Gleichungen wichtig, und das Muster (u(x, t), v(x, t)) erreichteine asymptotische Form. Wichtig hierfur ist naturlich, dass die Funktionenf und g so beschaffen sind, dass kein unbegrenztes Anwachsen von u oder vmoglich ist.

93

5.2.2 Zwei Klassen von Systemen

Wir hatten festgelegt, dass fu > 0 und gv < 0 ist und somit d > 1. Wenn u großerist als u0, wachst es also an, d.h. u beschreibt die Konzentration eines Aktivators,der sich selbst verstarkt. Entsprechend ist v die Konzentration eines Inhibitors,der sich selbst abbaut, und der viel schneller diffundiert als der Aktivator.

Wir uberlegen nun, welche Vorzeichen die anderen beiden Eintrage der Sta-bilitatsmatrix des homogenen Fixpunktes, fv und gu haben konnen. Die Be-dingung ∆0 > 0 fuhrt zu fvgu < 0. Wir betrachten der Reihe nach die beidendaraus resultierenden moglichen Falle.

gu > 0, fv < 0

In diesem Fall aktiviert u nicht nur sich selbst, sondern auch v. Der Inhibitorv inhibiert nicht nur sich selbst, sondern auch u. Dies ist in Abbildung 5.3illustriert.

+

+

−

−u v v

u

Abbildung 5.3: Darstellung des Situation gu > 0, fv < 0. Rechts ist der damitverbundene zeitliche Verlauf der beiden Konzentrationen skizziert.

Der Verlauf der Nullklinen in der Umgebung des homogenen Fixpunktes liegtdamit auch fest. Sie genugen den Gleichungen

fu(u − u0) + fv(v − v0) = 0

gu(u − u0) + gv(v − v0) = 0 .

Beide Nullklinen haben folglich positive Steigung. Aus ∆0 > 0 folgt fu/|fv| <gu/|gv|, d.h. die Kurve f = 0 ist flacher als die Kurve g = 0, so wie in Abbildung5.4 dargestellt.

u

vg=0

f=0

u

vg=0

f=0

Abbildung 5.4: Verlauf der Nullklinen fur den Fall gu > 0, fv < 0 (links) undfur den Fall gu < 0, fv > 0 (rechts).

94

gu < 0, fv > 0

In diesem Fall inhibiert u die Variable v, wahrend v die Produktion von u ak-tiviert. Dies ist in Abbildung 5.5 illustriert. Die beiden Nullklinen haben jetzt

+ −

u v

u−

+ v

Abbildung 5.5: Darstellung des Situation gu < 0, fv > 0. Rechts ist der damitverbundene zeitliche Verlauf der beiden Konzentrationen skizziert.

negative Steigung, wobei g = 0 wieder steiler verlauft als f = 0, wie in Abbil-dung 5.4 rechts gezeigt.

Wir werden im Folgenden konkrete Beispiele fur diese beiden Situationenbetrachten.

5.2.3 Das Schnakenberg-Modell

Das Schnakenberg-Modell (1979) ist das einfachste Modell fur chemische Reak-tionen zweier Spezies, das Oszillationen zeigen kann. Damit solche Oszillationenauftreten, muss es einen Reaktionsterm mit drei Molekulen geben. Dass dreiMolekule gleichzeitig am selben Ort zusammentreffen und direkt miteinanderreagieren, ist sehr unwahrscheinlich. Aber Reaktionsterme mit drei Molekulenergeben sich in Systemen mit mehr Partnern, wenn man sie auf die Beschrei-bung zweier Molekulsorten reduziert. Man denke zum Beispiel an Reaktionen,bei denen Enzyme als Katalysatoren fungieren. Das Schnakenberg-Modell istdurch die folgenden Gleichungen beschrieben:

A = k1 − k2A+ k3A2B +DA∇2A

B = k4 − k3A2B +DB∇2B . (5.37)

Beide Spezies werden mit konstanter Rate erzeugt. Die erste Spezies wird in ei-ner weiteren Reaktion abgebaut, und dann gibt es noch die Dreimolekulreaktion2A + B → 3A. Uns interessiert hier nicht der Parameterbereich, in dem Oszil-lationen auftreten, sondern derjenige, in dem es eine stabile homogene Losunggibt.

Wir bringen die Gleichungen nun auf die Form (5.27). Zunachst multipli-zieren wir beide Gleichungen mit

√

k3/k2 und defineren u =√

k3/k2A und

v =√

k3/k2B. Dann erhalten wir

u = k2

(

k1

k2

√

k3

k2− u+ u2v

)

+DA∇2u

v = k2

(

k4

k2

√

k3

k2− u2v

)

+DB∇2v .

95

Jetzt mussen wir noch alle Großen dimensionslos machen. Die Reaktionskon-stanten ki haben die Dimension einer inversen Zeit, und die Diffusionskonstan-ten haben die Dimension (Lange)2/Zeit. Wir messen alle Langen ~x in Einheiteneiner Lange L, ersetzen also ~x durch L~x, so dass ~x jetzt eine dimensionsloseLange ist. Spater wird es sich als zweckmaßig erweisen, fur L die raumlicheAusdehnung des Systems zu wahlen. Entsprechend messen wir die Zeit in Ein-heiten von L2/DA, ersetzen also t durch tL2/DA. Wenn wir die dimensionslosen

Parameter γ = k2L2/DA, d = DB/DA, a = k1

k2

√

k3

k2und b = k4

k2

√

k3

k2einfuhren,

erhalten wir schließlich

u = γ(

a− u+ u2v)

+ ∇2u

v = γ(

b− u2v)

+ d∇2v . (5.38)

Der homogene Fixpunkt ist u0 = a + b, v0 = b/(a + b)2. Damit er positiv ist,mussen b und a+ b positiv sein. An diesem Fixpunkt ist

fu =b− a

b+ a, fv = (a+ b)2 > 0, gu =

−2b

a+ b< 0 ,

gv = −(a+ b)2 < 0 , ∆0/γ2 = fugv − fvgu = (a+ b)2 > 0 . (5.39)

Damit fu > 0 ist, muss b > a sein. Des Weiteren fuhrt die Bedingung τ0 < 0auf

0 < b− a < (a+ b)3 (5.40)

Damit die diffusionsgetriebene Instabilitat auftreten kann, muss dfu + gv > 0und (dfu + gv)

2/4d > ∆0/γ2 sein, was auf die beiden Bedingungen

d(b− a) > (a+ b)3 ,

[d(b− a) − (a+ b)3]2 > 4d(a+ b)4 (5.41)

fuhrt. Damit ist ein Bereich im Raum der Parameter a, b, d festgelegt, in demsich raumliche Muster ausbilden konnen. Diese Analyse zeigt auch, dass diezweite im vorigen Teilkapitel geschilderte Situation vorliegt.

Nun betrachten wir eine Storung in x-Richtung der Form (5.30). Die Ausdeh-nung des Systems in x-Richtung sei von x = 0 bis x = p. Wenn wir verlangen,dass am Rand kein Molekul das System verlassen kann, haben wir ϕ = 0 undk = nπ/p mit positiven ganzzahligen n. Die Bedingung, dass die Determinan-te ∆ aus Gleichung (5.33) positiv sein soll, fuhrt dazu, dass k2 zwischen zweiWerten k2

1 und k22 liegen muss, die gegeben sind durch die Ausdrucke

k21,2 = γ

[d(b− a) − (a+ b)3] ∓ [d(b − a) − (a+ b)3]2 − 4d(a+ b)41/2

2d(a+ b).

(5.42)All diejenigen n, fur die k2 = (nπ/p)2 in diesem Intervall liegt, entsprecheninstabilen Storungen, die exponentiell wie eλ(n)t anwachsen, wobei λ(n) durch(5.36) gegeben ist. Das explizite Ausschreiben dieses Ausdrucks sparen wir uns...

Die Storung bildet also ein Streifenmuster, in dem sich hohe u-Konzentrationenmit niedrigen u-Konzentrationen langs der x-Achse abwechseln. Die Zahl derNulldurchgange von δu ist n, und die Zahl der Streifen ist somit n + 1, wobeisich “helle” und “dunkle” Streifen (also solche mit δu > 0 und δu < 0) abwech-seln. Wenn wir die Ansdehnung p des Systems verdoppeln, verdoppeln wir die

96

Werte n, fur die sich das Muster ausbildet. Großere Systeme bilden also mehrStreifen aus. Außerdem gibt es in großeren Systemen naturlich mehr Werte vonk, fur die Musterbildung moglich ist, da die k-Werte der Fourierreihe dichterliegen.

Statt p zu verdoppeln, konnen wir auch die fundamentale Langeneinheit Lmit der Ausdehnung des Systems gleichsetzen. Dann ist immer p = 1, und eineVerdopplung von L fuhrt zu einer Erhohung von γ um den Faktor 4, da γ jadurch γ = k2L

2/DA definiert war. Dies fuhrt wegen der Beziehung (5.42) zueiner Verdopplung der Grenzwerte des k-Intervalls, in dem Musterbildung auf-tritt, also wieder zu doppelt soviel Streifen. Indem wir am Parameter γ drehen,konnen wir also die Ausdehnung des Systems andern.

Eine analoge Betrachtung konnen wir fur Storungen in y-Richtung durchfuhren.Wenn die Ausdehnung des Systems in y-Richtung q betragt, haben wir alsoky = mπ/q. Allgemeine Storungen in x- und y-Richtung konnen in einem recht-eckigen System dann in den Moden

cos(nπx/p) cos(mπy/q)

entwickelt werden. Dann ist

k2 = π2

(

n2

p2+m2

q2

)

,

und dieses k2 muss in den durch (5.42) gegebenen Grenzen liegen, wenn dieStorung anwachsen soll. Die Storungen entsprechen also vertikalen und horizon-talen Streifen und “Schachbrett”mustern mit rechteckigen Feldern. Welches dermoglichen Muster bei einer allgemeinen Storung am Ende gewinnt, kann nichtdurch die hier durchgefuhrte lineare Stabilitatsanalyse bestimmt werden.

5.2.4 Das Gierer-Meinhardt-Modell

Gierer und Meinhardt fuhrten 1972 das folgende Aktivator-Inhibitor-Modell ein,das wir gleich in dimensionsloser Form anschreiben:

u = γ

(

a− bu+u2

v(1 + ku2)

)

+ ∇2u

v = γ(

u2 − v)

+ d∇2v . (5.43)

Der Aktivator vermehrt sich autokatalytisch, solange die Konzentration des In-hibitors klein ist. Gleichzeitig fuhrt eine hohe Aktivatorkonzentration zur Zu-nahme des Inhibitors. Also haben wir hier ein Beispiel fur die erste Klasse vonSystemen. Wir setzen im Folgenden k = 0. Der homogene Fixpunkt ist dannu0 = (a + 1)/b, v0 = u2

0. Also mussen a + 1 und b positiv sein. An diesemFixpunkt ist

fu = b1 − a

1 + a, fv = − b2

(1 + a)2< 0, gu =

2(a+ 1)

b> 0 ,

gv = −1 , ∆0/γ2 = fugv − fvgu = b > 0 . (5.44)

Damit fu > 0 ist, muss a < 1 sein. Des Weiteren fuhrt die Bedingung τ0 < 0auf

b <1 + a

1 − a(5.45)

97

Damit die diffusionsgetriebene Instabilitat auftreten kann, muss dfu + gv > 0und (dfu + gv)

2/4d > ∆0/γ2 sein, was auf die beiden Bedingungen

b1 − a

1 + a>

1

d,

db <1 + a

(1 − a)2

(

3 + a− 2√

2(1 + a))

oder db >1 + a

(1 − a)2

(

3 + a+ 2√

2(1 + a))

(5.46)

fuhrt. Damit ist ein Bereich im Raum der Parameter a, b, d festgelegt, in demsich raumliche Muster ausbilden konnen.

Die Bestimmung des Wellenzahl-Intervalls, in dem die Musterbildung ge-schieht, bleibt dem Leser uberlassen.

Ein weiteres, haufig verwendetes Modell fur diffusionsgetriebene Instabilitatist das Thomas-Modell, das ebenfalls zur ersten Klasse von Systemen gehort. Eswurde als Modell fur ein experimentell untersuchtes Substrat-Inhibitor-System(Sauerstoff und Uricinase) eingefuhrt und ist spezifiziert durch

f(u, v) = a− u− ρuv

1 + u+Ku2

g(u, v) = α(b − v) − ρuv

1 + u+Ku2. (5.47)

5.2.5 Tierfellmuster

Wenn man eines der Aktivator-Inhibitormodelle am Computer implementiert,kann man erforschen, was jenseits der linearen Storungstheorie passiert. Manerhalt ausgehend von einer kleinen Storung des homogenen Zustands nach ei-niger Zeit ein stationares Muster. Je nach Wert der Parameter konnen dabeiFlecken oder Streifen entstehen, wie sie fur Tierfellmuster charakteristisch sind.Man geht davon aus, dass die Grundlage fur diese Muster relativ fruh in derEmbryonalentwicklung gelegt werden. Diese Muster wachsen dann mit dem Tiermit. Wir haben gesehen, dass eine kleine raumliche Ausdehnung nur eine oderwenige Perioden des Musters zulasst. Je kleiner der Embryo wahrend der Mu-sterbildung ist, desto “grober” wird also am Ende das Muster sein. Interessantist auch, dass im Schwanz von Raubkatzen oft ein Ubergang von Flecken zuStreifen stattfindet. Streifen quer zum Schwanz bedeuten, dass in der kurzenRichtung (die des Schwanzumfangs) uberhaupt keine Musterbildung moglichwar, sondern nur in der langen Richtung.

Ubungsaufgaben

1. Diskutieren Sie, warum im thermodynamischen Gleichgewicht keine raum-lichen Muster vorliegen konnen.

2. Langs der geradlinigen Kuste eines Landes sei der Fischfang innerhalbvon H Kilometern streng reguliert. Jenseits dieser Schutzzone wird aberso intensiv gefischt, dass die Fischdichte dort praktisch Null ist. Die Fi-sche vermehren sich durch logistisches Wachstum und breiten sich durcheinen Diffusionsprozess aus. Begrunden Sie, dass das folgende Modell diebeschriebene Situation beschreibt:

u = ru(

1 − u

K

)

− Eu+Duxx (5.48)

98

mit den Randbedingungen u(H) = 0 und ux(0) = 0. Zeigen Sie, dass die

Schutzzone großer als π2

√

Dr−E

Kilometer sein muss, damit der Fischbe-

stand in Kustennahe nicht auf Null absinkt.

3. Betrachte das folgende Aktivator-Inhibitor-System:

u =u2

v− bu+ ∇2u

v = u2 − v + d∇2v (5.49)

mit positiven Konstanten b und d. Welches ist der Aktivator und welches der In-hibitor? Finde die nichttriviale stationare Losung des homogenen Systems undzeige, dass sie fur b < 1 stabil ist. Ermittle die Bedingungen dafur, dass dieserhomogene Zustand durch Diffusion instabil wird. Zeige, dass der Parameterbe-reich fur die diffusionsgetriebene Instabilitat durch die Ungleichungen 0 < b < 1und db > 3 + 2

√2 gegeben ist. Markiere diesen Bereich in der (b, d)-Ebene. Am

Rand dieses Bereichs findet die Bifurkation zur Instabilitat statt. Dort gibt esgenau eine Wellenzahl, die instabil wird. Zeige, dass diese kritische Wellenzahlkc die Bedingung k2

c = (1 +√

2)/d erfullt.

5.3 Die komplexe Ginzburg-Landau-Gleichung

Im vorigen Teilkapitel haben wir ein System betrachtet, in dem die homogeneLosung instabil wird, wenn die Parameter gewisse kritische Werte uberschreiten.Mit der linearen Storungstheorie konnten wir nur fur kurze Zeiten beschreiben,wie die Storung mit der Zeit anwachst. Da die lineare Storungstheorie ein rei-nes exponentielles Anwachsen beschreibt, kann sie nicht die fur langere Zeiteneintretende Sattigung der Storung beschreiben.

Die komplexe Ginzburg-Landau-Gleichung ist die einfachste denkbare Glei-chung, die uber die lineare Theorie hinausgeht. Wir gehen davon aus, dass dieParameterwerte in der Nahe der Bifurkation liegen, so dass die Amplitude derStorung immer klein bleibt. Dann konnen wir eine Taylorentwicklung in Po-tenzen der Amplitude der Storung durchfuhren und nach den ersten Termenabbrechen.

Dieses Vorgehen haben wir sowohl im Zusammenhang mit der nichtlinearenDynamik, als auch bei der Theorie kritischer Phanomene gewahlt. Im Rahmender nichtlinearen Dynamik haben wir die “Normalform” verschiedener Bifur-kationen eingefuhrt, die sich aus einer Taylorentwicklung der Variablen ergibt.Im Rahmen der Theorie kritischer Phanomene haben wir die Landau-Theorieund die Ginzburg-Landau-Theorie betrachtet. Beide basieren auf einer Entwick-lung in Potenzen des Ordnungsparameters, wobei die Ginzburg-Landau-Theorieauch raumliche Fluktuationen zulasst. Es ist kein Zufall, dass die komplexeGinzburg-Landau-Gleichung ebenfalls nach Ginzburg und Landau benannt ist.Sie ist namlich eine Verallgemeinerung auf komplexe Koeffizienten und Ord-nungsparameter von derjenigen Gleichung, die im Rahmen der Theorie kritischerPhanomene die zeitliche Entwicklung des Ordnungsparameters beschreibt. Da-her konzentrieren wir uns zunachst auf die Ordnungsparameterdynamik in derNahe eines kritischen Punktes.

99

5.3.1 Ordnungsparameterdynamik in der Nahe eines kri-

tischen Punktes

In Kapitel 4 haben wir keine Zeitabhangigkeiten bei kritischen Phanomenen be-trachtet, da die Dynamik im thermischen Gleichgewicht sich ja auf wenig inter-essante Fluktuationen beschrankt. Nun wollen wir den Fall betrachten, dass einSystem in der Nahe des kritischen Punktes außerhalb des thermischen Gleichge-wichts initiiert wird und dann sich selbst uberlassen wird. Dies erreicht man z.B.dadurch, dass man die Temperatur von einem Wert oberhalb von Tc sprunghaftauf einen Wert unterhalb von Tc absenkt. Nach diesem Temperatursprung wirdder Ordnungsparameter mit der Zeit anwachsen, bis er seinen Gleichgewichts-wert erreicht. Wenn der Ordnungsparameter keine Erhaltungsgroße ist, setztman hierfur eine Relaxationsdynamik an,

∂m(~x, t)

∂t= −λδH[m(~x, t)]

δm(~x, t), (5.50)

wobei H das in Gleichung (4.17) gegebene Ginzburg-Landau Funktional ist,

H[m(~x, t)]

kBT=

∫

ddx[

rm2(~x, t) + um4(~x, t) + c (∇m(~x, t))2]

.

Anschaulich bedeutet dies, dass der Ordnungsparameter sich uberall so andert,dass dadurch die Energie erniedrigt wird. Je großer die Energieerniedrigungist, desto schneller ist die Anderung. Man kann sich das so vorstellen, dassder Ordnungsparameter in der Energielandschaft bergab wandert, bis er einEnergieminimum erreicht. Dieses Minimum ist der durch die Landau-Theoriegegebene homogene Ordnungsparameter m = −

√

r/2u, wenn T < Tc ist. Wennwir den Ausdruck fur H in die Gleichung fur den Ordnungsparameter einsetzen,erhalten wir formal eine Reaktions-Diffusionsgleichung:

m(~x, t) = λkBT(

−2rm− 4um3 + 2c∇2m)

. (5.51)

Den “Diffusionsterm” erhalt man am einfachsten, indem man das Integral uberx diskretisiert und dann erst die Funktionalableitung ausrechnet.

5.3.2 Anwachsen einer diffusionsgetriebenen Instabilitat

Als nachstes betrachten wir die diffusionsgetriebenen Instabilitaten des vorigenTeilkapitels. Wir beschranken und auf den Fall, dass die Instabilitat eine ein-deutige und einkomponentige Wellenzahl hat. Das bedeutet, dass die homogeneLosung nur bezuglich einer Raumrichtung instabil wird, was wir dadurch errei-chen konnen, dass die Ausdehnung des Systems in der anderen Richtung sehrklein ist, oder dass die Diffusion anisotrop ist, oder dass der Bereich instabiler~k-Werte so schmal ist und so genau positioniert ist, dass dort ky = 0 ist.

Dann gilt fur δu (und analog fur δv):

δu(~x, t) = A(~x, t)eiqcx +A∗(~x, t)e−iqcx , (5.52)

wobei wir nun die Wellenzahl der instabilen Storung mit qc notiert haben. ImRahmen der linearen Stabilitatsanalyse war A(t) = A(0)eλ(qc)t. Dadurch, dasswir A komplex wahlen, konnen wir die Phase der Storung berucksichtigen. Wir

100

wahlen im Folgenden die Zeitskala so, dass λ(qc) = 1 ist. Dann istA(t) = A(0)et,bzw. A = A. Wenn wir nun davon ausgehen, dass die Korrekturen zu dieserlinearen Theorie aus einer Taylorentwicklung resultiert, haben wir die folgendenfuhrenden Terme in A und seiner Ortsableitung:

A = A− |A|2A+ ∇2A . (5.53)

Durch geeignete Skalierung von A und der Langenskala konnen wir alle Ko-effizienten zu 1 machen. Es treten nur ungerade Potenzen von A auf, da dasVorzeichen von A keine Rolle spielen darf. Auch eine Umkehr der Raumachsendarf diese Gleichung nicht andern, und deshalb tritt der ∇-Operator in der zwei-ten Potenz auf. Der Sattigungswert von A darf nicht von seiner Phase abhangen,deshalb stehen im zweiten Term die Betragsstriche. Genau wie Gleichung (5.51)lasst sich diese Gleichung durch die Ableitung eines Energiefunktionals aus-drucken,

∂A(~x, t)

∂t= −δH[A(~x, t)]

δA∗(~x, t), (5.54)

mit

H[A(~x, t)] =

∫

ddx

[

−|A|2(~x, t) +1

2|A|4(~x, t) + |∇A(~x, t)|2

]

.

Die Amplitude der Storung wachst uberall an, bis sie ihren Sattigungswert 1erreicht hat.

5.3.3 Die komplexe Ginzburg-Landau-Gleichung

Die beiden bisher betrachteten Situationen beschreiben Spezialfalle der allge-meineren komplexen Ginzburg-Landau-Gleichung. Sie wurde von mehreren Ar-beitsgruppen im Jahr 1971 eingefuhrt und enthalt zwei komplexe Koeffizienten.Damit kann man nicht nur Storungen betrachten, die inhomogen im Raum sind,sondern auch Storungen, die zeitlich oszillieren, also zusatzlich zu einem qc auchdurch ein ωc charakterisiert sind. Sie lautet

A = A− (1 + ic)|A|2A+ (1 + ib)∇2A . (5.55)

Eine Variable u, deren homogener Fixpunkt unter einer Storung mit qc und ωc

instabil wird, andert sich nun gemaß

δu(~x, t) = A(~x, t)eiqcx−iωct +A∗(~x, t)e−iqcx+iωct , (5.56)

Hier ist auch der Spezialfall qc = 0, ωc 6= 0 enthalten, der eine raumlich homoge-ne und zeitlich oszillierende Instabilitat beschreibt, wie sie bei manchen chemi-schen Reaktionen auftritt. Der Fall qc 6= 0, ωc 6= 0 tritt zum Beispiel beim Ein-setzen der Rayleigh-Benard-Konvektion auf. Die komplexe Ginzburg-Landau-Gleichung tritt in einer Vielfalt von Systemen auf. Hierzu gehoren auch Supra-leitung und Superfluiditat, Stringtheorie, Flussigkristalle und Bose-Einsteinkon-densation.

Jetzt lasst sich die rechte Seite der Bewegungsgleichung fur A nicht mehr alsAbleitung eines Energiefunktionals schreiben, und es kann eine große Vielfaltvon Phanomenen auftreten. Hierzu gehoren ebene und spiralformige Wellen mitQuellen und Senken, verschiedene Arten von raumzeitlichem Chaos, und glas-artiges Verhalten (also ungeordnet im Raum und extrem langsam in der Zeit).

101

Es gibt lange Ubersichtsartikel und Bucher zur komplexen Ginzburg-Landau-Gleichung. Wir wollen im Folgenden nur zwei Themen aus der großen Vielfaltanschneiden, namlich ebene Wellen und eindimensionale koharente Strukturen.

Ebene Wellen

Eine ebene WelleA =

√

1 −Q2ei[ ~Q~x−ωp(Q)t+φ] (5.57)

lost die Gleichung (5.55), wenn man Q2 < 1 und ωp(Q) = c(1−Q2)+bQ2 wahlt.Hierbei ist φ eine willkurliche Phase. Wir fuhren die Abkurzung F 2 = 1−Q2 ein.Solche ebenen Wellen konnen naturlich nur dann beobachtet werden, wenn siestabil gegenuber kleinen Storungen sind. Um dies zu uberprufen, fuhrt man einelineare Stabilitatsanalyse der ebenen Welle durch. Man macht also den Ansatz

A =[

F + δa+eλt+i~k~r + δa−e

λt−i~k~r]

ei[ ~Q~x−ωp(Q)t+φ] , (5.58)

mit kleinen Amplituden δa±. Eingesetzt in die komplexe Ginzburg-Landau-Glei-chung erhalt man eine Beziehung fur λ. Wenn der Realteil von λ fur alle k nega-tiv ist, ist die betrachtete Welle stabil. Die entsprechende Rechnung fuhren wirhier nicht durch. Man kann sie z.B. in dem Ubersichtsartikel von Aranson undKramer (Review of Modern Physics 74) aus dem Jahr 2002 finden. Insbesonderekann man durch eine Entwicklung in Potenzen von k2 die Stabilitat der Wellegegenuber langwelligen Storungen untersuchen.

Eindimensionale koharente Strukturen

Eindimensionale koharente Strukturen haben die Form

A(x, t) = e−iωkt+iφ(x−vt)a(x− vt) . (5.59)

Wenn man diesen Ansatz in die komplexe Ginzburg-Landau-Gleichung (5.55)einsetzt, ergeben sich folgende Gleichungen fur die Funktionen a und ψ = dφ/dx:

ax = s ;

sx = ψ2a− γ−1[(1 + bωk)a+ v(s+ bψa) − (1 + bc)a3] ;

ψx = −2sψ/a+ γ−1[b− ωk + v(sb/a− ψ) − (b− c)a2] . (5.60)

Wir haben hierbei die Abkurzung γ = 1+ b2 verwendet. Diese drei Gleichungenkonnen wir als dynamisches System mit drei Variablen interpretieren. Die Fix-punkte dieses dynamischen Systems sind der triviale Zustand a = 0 und ebeneWellen A(x, t) mit konstantem a und ψ = dφ/dx.

Trajektorien, die in diesen Fixpunkten starten und enden, beschreiben lo-kalisierte koharente Strukturen. Sie konnen homokline oder heterokline Orbitssein. Homokline Orbits enden in demselben Fixpunkt, in dem sie gestartet sind,heterokline Orbits verbinden zwei verschiedene Fixpunkte. Wenn der Fixpunktein Sattelpunkt in zwei Dimensionen ist, lauft das homo- bzw heterokline Or-bit langs einer instabilen Eigenrichtung aus dem Fixpunkt heraus und langseiner stabilen Eigenrichtung hinein. Wenn der Ausgangsfixpunkt N allgemeinnN instabile Eigenrichtungen hat, hat man nN − 1 Parameter frei, um die An-fangsrichtung der Trajektorie festzulegen. Außerdem hat man noch ωk und v als

102

freie Parameter, zusammen also nN +1 freie Parameter. Wenn der EndfixpunktL allgemein nL instabile Eigenrichtungen hat, gibt es nL Bedingungen, die dieeinlaufende Trajektorie so einschranken, dass sie zu all diesen Richtungen ortho-gonal ist. Die Zahl der Freiheitsgrade, die unsere Trajektorie charakterisieren,ist also n = nN − nL + 1. Fur n ≥ 1 gibt dies eine ganze Schar von koharentenStrukturen, fur n = 0 gibt es nur einen diskreten Satz solcher Strukturen, undfur n < 0 gibt es sie gar nicht.

Man klassifiziert sie folgendermaßen: Pulse entsprechen dem homoklinen Or-bit zum trivialen Fixpunkt a = 0. Fronten sind heterokline Orbits zwischeneiner ebenen Welle und dem trivialen Fixpunkt. Domanengrenzen entsprechenheteroklinen Orbits zwischen zwei verschiedenen ebenen Wellen. Homoklonen

entsprechen homoklinen Orbits zu einer ebenen Welle. Sie konnen periodischoder auch nichtperiodisch angeordnet sein.

Literaturhinweise

Zu dem Thema jedes Kapitels gibt es eine Reihe von Lehrbuchern und Uber-sichtsartikeln. Die folgenden Bucher und Artikel wurden fur diese Vorlesungkonsultiert und sind in der Fachbereichs-Bibliothek bzw Online verfugbar:

• Kapitel 2: L. Reichl, A modern course in statistical physics; Paul undBaschnagel, Stochastic processes from physics to finance.

• Kapitel 3: S. Strogatz, Nonlinear dynamics and Chaos. J.D. Murray, Ma-thematical Biology, Bd.1.

• Kapitel 4: F. Schwabl, Statistische Physik. M. Kardar, Statistical physicsof fields. D. Stauffer und A. Aharony, percolation theory.

• Kapitel 5: J.D. Murray, Mathematical Biology, Bd.2. I.S. Aranson, L. Kra-mer, “The world of the complex ginzburg-Landau equation”, Reviews ofModern Physics 74 (2002), ab S.99.

103