Inhaltsverzeichnis

Kapitel 1. Einführung 5 1.1. Gliederung 7

Kapitel 2. Oberflächendynamik in partiellen Differentialgleichungen 9 2.1. Langevingleichungen 11 2.1.1. Einfluß des Rauschens 11 2.1.2. Skalenverhalten 12 2.1.3. Ein Beispiel: Die lineare Gleichung 13 2.1.4. Galilei–Invarianz 14 2.2. Strukturbildung und -vergröberung 15 2.3. Chaotische Dynamik 18

Kapitel 3. Die Kuramoto-Sivashinsky–Gleichung 21 3.1. Grundlegende Bemerkungen 21 3.1.1. Nichtdissipative Dynamik 21 3.1.2. Linearisierung und Strukturbildung 22 3.1.3. Verhalten auf großen Skalen 22 3.2. Ein Teilchenmodell 25 3.2.1. Definition des Modells 25 3.2.2. Statistik des Teilchenmodells auf großen Skalen 29 3.2.3. Zusammenfassende Bewertung 33 3.3. Die anisotrope KS–Gleichung 34 3.3.1. Einführung 34 3.3.2. Anisotropie im linearen Term 35 3.3.3. Anisotropie im nichtlinearen Term 37 3.3.4. Vorzeichenwechsel von β 38 3.4. Zusammenfassung 40

Kapitel 4. Molekularstrahlepitaxie 41 4.1. Homoepitaxie 43 4.2. Stufenbarrieren 45 4.3. Wachstum mit Stufenbarrieren 47 4.3.1. Oberflächen mit schwacher Neigung 47 4.3.1.1. Pyramidenwachstum 47 4.3.1.2. Lagenwachstum 49 4.3.2. Oberflächen mit starker Neigung 50

Kapitel 5. Kontinuumsmodelle für Oberflächen 53 5.1. Oberflächen im lokalen Gleichgewicht 54 5.2. Bemerkungen zum Nichtgleichgewicht 57 5.3. Chemisches Potential für Adatome 57 5.4. Dynamik der Oberfläche 59 5.4.1. Oberflächendiffusion 59 5.4.2. Evaporation und Kondensation 61

Kapitel 6. Stufenbarrieren im Kontinuumsbild 63 6.1. Nichtgleichgewichtsstrom 63 6.1.1. Thermodynamische Grundüberlegung 63 6.1.2. Neigungsabhängigkeit des Oberflächenstroms 64 6.2. Vergleich zur Kinetik der Phasenordnung 68 6.3. Lineare Stabilitätsanalyse 71 6.4. Strukturvergröberung auf singulären Oberflächen 74 6.4.1. Skalenrelationen 74 6.4.2. Sphärische Näherung 76 6.4.3. Weitere Näherungsverfahren 80 6.5. Instabilität auf vizinalen Oberflächen 82 6.5.1. Skalenabschätzung 83 6.5.2. Numerische Ergebnisse 85 6.5.2.1. Das Solid-on-Solid–Modell 86 6.5.2.2. Längen- und Zeitskalen in den Simulationen 87 6.5.2.3. Die Übergangsphase 88 6.6. Instabiles Wachstum, eine Schlußbetrachtung 89

Kapitel 7. Das Kristallgitter im Kontinuumsmodell 93 7.1. Die erhaltene Sinus-Gordon–Gleichung 93 7.2. Renormierung nach Nozières und Gallet 95 7.2.1. Die Methode 95 7.2.2. Anwendung auf die erhaltene Sinus-Gordon–Gleichung 96 7.3. Dämpfung der Wachstumsoszillationen 99 7.4. Die erhaltene KPZ–Nichtlinearität 102

Kapitel 8. Zusammenfassung 105

Anhang A. Strukturfaktor in sphärischer Näherung 107

Anhang B. Allgemeine Erläuterungen 109

INHALTSVERZEICHNIS 3

B.1. Verwendete Begriffe 109 B.2. Numerik 111

Literaturverzeichnis 113

Ein verwandtes Gefallen fand er an Eisblumen [ …]. Ich möchte sagen: alles wäre gut gewesen und man hätte darüber zur Tagesordnung übergehen können, wenn die Erzeugnisse sich, wie es ihnen zukam, im Symmetrisch-Figürlichen, streng mathematischen und Regelmäßigen gehalten hätten. Aber daß sie mit einer gewissen gaukelnden Unverschämtheit Pflanzliches nachahmten, aufs wunderhübschesten Farrenwedel, Gräser, die Becher und Sterne von Blüten vortäuschten, daß sie mit ihren eisigen Mitteln im Organischen dilettierten, das war es, worüber Jonathan nicht hinwegkam, und worüber seines gewissermaßen mißbilligenden, aber auch bewunderungsvollen Kopfschüttelns kein Ende war. Bildeten, so lautete seine Frage, diese Phantasmagorien die Formen des Vegetativen vor oder bildeten sie sie nach? Keines von beidem, erwiderte er wohl sich selbst; es waren Parallelbildungen. Doch die schöpferisch träumende Natur träumte hier und da dasselbe, und durfte von Nachahmung die Rede sein, so gewiß nur von wechselseitiger. Sollte man die wirklichen Kinder der Flur als die Vorbilder hinstellen, weil sie organische Tiefenwirklichkeit besaßen, die Eisblumen aber bloße Erscheinungen waren? Aber ihre Erscheinung war das Ergebnis keiner geringeren Kompliziertheit stofflichen Zusammenspiels als diejenige der Pflanzen.

Thomas Mann [1]

KAPITEL 1 Einführung

Wie ist es möglich, daß in der Natur wie von außen bewußt geordnete Strukturen entstehen, obwohl ihre Gesetze im Kleinen wirken und für die größere Ordnung der Bestandteile blind zu sein scheinen? In der unbelebten Natur erscheint dieser Gegensatz besonders unverständlich. Wir waren und sind deshalb immer wieder geneigt, auch ihr Eigenschaften des Lebens zuzusprechen.

Auch an Oberflächen zwischen verschiedenen Medien beobachtet man geordnete Strukturen: Eine Wasseroberfläche bildet ein Streifenmuster, wenn Wind über sie hinwegbläst, Sand im seichten Wasser am Strand wird von den ankommenden Wellen zu einem Rippenprofil zusammengeschoben, Schnee auf einer ebenen Fläche schmilzt nicht gleichmäßig ab, sondern bildet getrennt von scharfen Graten viele kleine Senken gleicher Größe.

Physikalische Experimente zeigen Strukturbildung in allen Größen. Mit bloßem Auge sichtbar sind z.B. Konvektionszellen oder die Finger zwischen zwei Flüssigkeiten unterschiedlicher Viskosität. In dieser Arbeit werden Strukturen auf wachsenden Kristalloberflächen untersucht, die sich nur über zehn bis einige hundert Atomradien erstrecken, wie das in Bild 1.1 gezeigte auf einer Platinoberfläche gewachsene Muster aus Rippen und Buckeln.

Wie in ihrer Bildung, so zeigt sich auch bei der weiteren Entwicklung der Strukturen eine Ordnung im Makroskopischen, die von rein mikroskopischen Regeln bedingt wird. Starke Ausprägung einer solchen Ordnung ist die Skaleninvarianz, die in vielen Experimenten und Simulationen beobachtet wird, und die als Hypothese vielen theoretischen Auswertungen und Interpretationen zugrunde liegt. In allgemeinster Formulierung und im Hinblick auf die in dieser Arbeit untersuchten Systeme fordert die Skalenhypothese, daß die Konfigurationen zu zwei verschiedenen Zeitpunkten “gleich aussehen” – vor geeigneten Längenmaßstäben und in geeignet angepaßtem Zeitfluß.

Die genannten Beispiele zeigen eine Gemeinsamkeit der zu untersuchenden Systeme: Sie befinden sich nicht im Gleichgewicht, sondern sind von Abbildung 1.1. Diese rastertunnelmikroskopische Aufnahme einer Pt(111)–Oberfläche nach Homoepitaxie von fünf Monolagen bei T = 400 K zeigt exemplarisch Strukturen, deren Bildung und Wachstum in dieser Arbeit untersucht werden. Auf einer Makrostufe (geneigtes Stück einer Kristalloberfläche mit geraden Stufen gleichen Abstands) haben sich Rippen und Buckel gebildet, auf einem Stück einer Hochsymmetriefläche (oben links) gleichförmige Hügel (mit freundlicher Genehmigung von T. Michely, [2]).

außen getrieben. An den in dieser Arbeit behandelten wachsenden Kristallen ist das offensichtlich. Im Formalismus ihrer Modellierung gleichen sie einer anderen Ausprägung der Nichtgleichgewichtsdynamik, nämlich Spinsystemen mit sich entmischenden Phasen.

Die Vorgehensweise ist phänomenologisch. Die Modelle sind Kontinuumsgleichungen, die das Geschehen auf einer mesoskopischen Skala erfassen. Ihre Konstruktion beruht auf Plausibilitätsüberlegungen, etwa ob sie Symmetrien und Erhaltungsgrößen der zu beschreibenden Systeme berücksichtigen. Ihre Rechtfertigung erhalten sie durch den Vergleich mit Experimenten oder “realistischeren” Simulationen, in die mehr der mikroskopischen physikalischen Eigenschaften eingehen. Die Ergebnisse haben daher zweierlei Bedeutung: Sie sind einerseits Vorhersagen für durchzuführende detailliertere Untersuchungen. Andererseits reduziert die phänomenologische Modellbildung die Gesamtheit des Geschehens auf wenige vereinfachte

1.1. GLIEDERUNG 7

Prozesse – wie jede physikalische Modellierung – und trägt so zum Verständnis der Vorgänge bei, wenn die Phänomene des vereinfachten Modells mit der Wirklichkeit übereinstimmen.

1.1. Gliederung

Die Arbeit gliedert sich im weiteren Verlauf wie folgt: In einem ersten Teil wird der theoretische Rahmen abgesteckt, in dem Oberflächendynamik untersucht wird (Kapitel 2 und 3). Ihm gegenüber steht ein zweiter Teil, der mehr den experimentellen Hintergrund des Kristallwachstums berücksichtigt.

Verschiedene Typen von Modellen zur Strukturbildung und Oberflächendynamik werden zusammen mit gängigen Methoden zur Untersuchung und Beschreibung in Kapitel 2 vorgestellt. Im folgenden Kapitel 3 werden sie auf zwei Modelle angewandt, um deren theoretisch interessante Aspekte herauszuarbeiten.

Daran schließt sich der zweite Teil an, der Oberflächenmodelle auf theoretische Fragestellungen im Zusammenhang mit Kristallwachstum unter Molekularstrahlepitaxie untersucht. Deren experimenteller Hintergrund und Grundkonzepte theoretischer Modellierung der herausragenden Phänomene werden in Kapitel 4 vorgestellt. Im Detail werden Kontinuumsmodelle, partielle Differentialgleichungen und alle ihre Elemente zur Beschreibung der Oberflächendynamik, in Kapitel 5 konstruiert. Charakteristische Nichtgleichgewichtseigenschaften und ihr Beitrag zur Strukturbildung und -entwicklung auf Oberflächen wachsender Kristalle sind Gegenstand von Kapitel 6. Kapitel 7 schließt den zweiten Teil mit einem stochastischen Modell ab, das den Übergang zwischen schichtweisem und rauhem Wachstum beschreibt. Der zweite, eher experimentell motivierte Teil, kann mit gelegentlicher Bezugnahme auf die vorhergehenden rein theoretischen Kapitel für sich allein gelesen werden.

Am Ende steht eine Zusammenfassung, gefolgt von einem Anhang mit der Darstellung einiger Rechnungen, einigen Definitionen und Erläuterungen zu an vielen Stellen der Arbeit verwendeten Begriffen, sowie einer Skizze des verwendeten Algorithmus zur numerischen Integration nichtlinearer partieller parabolischer Differentialgleichungen.

8 1. EINFÜHRUNG

KAPITEL 2 Oberflächendynamik in partiellen Differentialgleichungen

In diesem Kapitel werden allgemeine Konzepte der theoretischen Beschreibung der Oberflächendynamik eingeführt, wie sie im zweiten Teil der Arbeit zur Modellierung physikalischer Systeme benötigt werden. Der physikalische Hintergrund soll nur am Rande in erläuternden Bemerkungen berücksichtigt werden.

Die erste Wahl einer geeigneten Parametrisierung einer Fläche im Raum

r = r(x, t), (1)

bezeichnet mit x∈ℝ^{d} eine innere Koordinate der Fläche, und mit r ∈ ℝ^{d+1} die Koordinate eines Oberflächenpunktes im umgebenden Raum. Physikalische Anwendungen bietet die Einbettung einer d=2–dimensionalen Fläche in den ℝ^{3}, aber auch der Fall d = 1, etwa zur Beschreibung einer Kristallstufe auf einer Oberfläche.

Die Dynamik ist in allen physikalischen Anwendungen dieser Arbeit überdämpft. Effekte wie die einer Masse oder Trägheit der Oberfläche sollen vernachlässigt werden. Stattdessen wird die Dynamik hinreichend beschrieben, wenn man zu jedem Punkt r(x, t) eine Geschwindigkeit v kennt, mit der er durch den umgebenden Raum verschoben wird. v = v_{k}+v_{e} spaltet sich auf in einen konfigurationsabhängigen und einen externen Anteil. Ersterer berücksichtigt die Wirkung am Punkt r, die von der Konfiguration der Oberfläche an anderen Punkten r′ ausgeht, etwa die Tendenz zur Relaxation nach Verformungen, der zweite den Einfluß des umgebenden Mediums auf die Oberfläche, meist eine zufällige Kraft, die einerseits thermisch bedingt sein kann, andererseits auf Unordnung im umgebenden Medium zurückgeht.

∂_{t}r(x, t) = v({r(x, t)}, x, t) (2)

wäre die allgemeinste Form dieser Dynamik.

Eine günstige Vereinfachung ergibt sich aus der Monge–Darstellung der Oberfläche: Hat sie eine Vorzugsrichtung, bezüglich der es keine Überhänge gibt, kann man jeden Oberflächenpunkt auf einen “Grundflächenpunkt” x

10 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

projizieren und ihm eine eindeutige Höhe h(x, t) zuordnen. In diesem Fall ist die Oberfläche

r(x, t) = (x, h(x, t)) (3)

der Graph der Höhenfunktion h. Die Dynamik wird angegeben durch die jeweilige Verschiebung in h–Richtung über einem Punkt x:

∂_{t}h(x, t) = v_{k}{h} + v_{e}(h(x, t), x, t). (4)

Prinzipiell kann v_{k} über einem Fußpunkt x von der gesamten Konfiguration h im Raum abhängen, oft jedoch nur von h in einer kleinen Umgebung um x. Dann wird die Abhängigkeit v_{k} sinnvoll in räumliche Ableitungen von h entwickelt und (4) erhält die Form einer parabolischen partiellen Differentialgleichung.

Die Struktur der Gleichungen muß die Symmetrien des zugrundeliegenden physikalischen Problems widerspiegeln. Allen Formen von Oberflächendynamik gemein ist die Translationsinvarianz, sei es unter einer diskreten Gruppe von Translationen, wenn man die Gitterstruktur eines Kristalls berücksichtigt, oder andernfalls auch für beliebige Verschiebungen. Alle Gleichungen vom Typ von (4) sind daher invariant unter der Translation h → h + c mit c ∈ ℝ bzw. mit a einer Gitterkonstanten (c/a) ∈ ℤ. Die unten angeführte Gleichung (5) fällt nicht in diese Klasse, beschreibt auch keine Oberflächendynamik, weshalb die Konfigurationen hier mit ϕ statt h bezeichnet werden.

In den folgenden drei Abschnitten sind partielle Differentialgleichungen für Oberflächendynamik nach verschiedenen Typen klassifiziert. Beispiele dafür sind:

• Relaxierende deterministische Gleichungen wie die Diffusionsgleichung

∂_{t}h(x, t) = Δh(x, t)

– Δ bezeichnet in dieser Arbeit immer den Laplace–Operator ∇^{2} – oder die Ginzburg-Landau–Gleichung zur Relaxation des ϕ^{4}–Modells bei Temperatur T = 0 (Modell A) [3, 4], ein Standardbeispiel der Strukturbildung und des Domänenwachstums,

∂_{t}ϕ(x, t) = Δϕ(x, t) + ϕ(x, t) − ϕ(x, t)^{3}. (5)

Das Interesse liegt an der von einem ungeordneten Anfangszustand ausgehenden Entwicklung solcher Gleichungen.

2.1. LANGEVINGLEICHUNGEN 11

• Gleichungen mit einer zufälligen externen Kraft (Langevingleichungen), wie [5]

∂_{t}h(x, t) = ν Δh(x, t) + η(x, t) (6)

und die Kardar-Parisi-Zhang–Gleichung (KPZ) [6]

∂_{t}h(x, t) = ν Δh(x, t) + \frac{λ}{2}(∇h(x, t))^{2} + η(x, t). (7)

Anders als beim letzten Typ von Gleichungen untersucht man hier die Dynamik ausgehend von einer Anfangskonfiguration h ≡ 0, also die Aufrauhung einer ursprünglich glatten Fläche.

• Chaotische deterministische Gleichungen wie [7, 8]

∂_{t}h(x, t) = −ν Δh(x, t) − K Δ^{2}h(x, t) + \frac{λ}{2}(∇h(x, t))^{2}, (8)

bei denen die sich aus einer Anfangskonfiguration entwickelnde Struktur in ständiger Auflösung und Erneuerung begriffen ist. Sie verhalten sich auf großen Skalen in vielem ähnlich wie die verrauschten Gleichungen [9, 10, 11] und stellen ein Bindeglied zwischen den beiden erstgenannten Typen von Gleichungen dar.

Gleichungen vom ersten Typ werden in den Kapiteln 5 und 6 behandelt, eine Langevingleichung in Kapitel 7. Verallgemeinerungen von Gleichung (8) folgen im nächsten Kapitel.

Zu ihrer Untersuchung nützliche Konzepte, wie die Symmetrie unter bestimmten Transformationen, Skaleninvarianz ihrer statistischen Eigenschaften, und die Stabilität oder Instabilität einer stationären Konfiguration und ihre Bedeutung für die Strukturbildung, werden in den folgenden Abschnitten dieses Kapitels vorgestellt.

2.1. Langevingleichungen

2.1.1. Einfluß des Rauschens. Der Rauschterm in Langevingleichungen modelliert eine zufällige auf die Oberfläche wirkende Kraft. Ausgehend von der glatten Anfangskonfiguration h(x, t = 0) ≡ 0 rauht sich die Oberfläche unter dem Einfluß des Rauschens auf. Da in die zum Zeitpunkt t auf den Punkt (h(x, t), x) einwirkende Kraft viele mikroskopische Details eingehen, nimmt man allgemein gaußisches weißes Rauschen an, zur Vereinfachung meist mit Mittelwert ⟨η(x, t)⟩ = 0. In dieser Arbeit sollen nur kurzreichweitige Korrelationen in Raum und Zeit untersucht werden. Die Kovarianz des Rauschens wird angenommen als

⟨η(x, t) η(x′, t′)⟩ = 2 D Δ^{n}δ(x−x′) δ(t−t′). (9)

12 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

Im Falle von n = 1 bleibt beim Rauschen das Volumen unter der Oberfläche ∫ d^{d}x h erhalten, man spricht von “erhaltendem Rauschen”. Mit dieser Kovarianz ist das Rauschen η =∇·j die Divergenz eines zufälligen Oberflächenstroms. In der Anwendung in Kapitel 7 ist das Rauschen dagegen nicht volumenerhaltend, n=0, was aus der beschriebenen experimentellen Situation einsichtig ist: Der auf die Kristalloberfläche auftreffende Teilchenstrahl ist nicht homogen sondern besteht aus zufälligen Einzelereignissen, nämlich dem Auftreffen einzelner Teilchen. Der Rauschterm berücksichtigt deshalb die Schwankungen der insgesamt aufgedampften Masse um den jeweiligen Mittelwert.

Das Rauschen steht im Wechselspiel mit der glättenden Tendenz des deterministischen (konfigurationsabhängigen) Anteils der Dynamik. In den Gleichungen (6) und (7) kann man den Term Δh als von einer Oberflächenspannung getriebenes Bestreben interpretieren, den Flächeninhalt zu verringern (siehe dazu auch Kapitel 5). Ohne Rauschen würde die Oberfläche in eine stationäre Konfiguration – etwa eine Form minimalen Flächeninhalts – relaxieren.

Wie dieses Wechselspiel die Form der Oberfläche zeitlich verändert, sieht man z.B. in der mittleren Breite

w(t) := \sqrt{⟨h^{2}⟩ − ⟨h⟩^{2}}, (10)

besser aber in der Verteilung der Höhendifferenzen über zwei verschiedenen Fußpunkten x und y. Ihr zweites Moment

G_{2}(x, y, t) = ⟨(h(x, t) − h(y, t))^{2}⟩ (11)

wird dabei oft zur Auswertung verwendet. Mit der Wahl einer glatten Anfangsbedingung und gleicher Rauschstärke an jedem Raumpunkt (s.o.) bleibt ⟨h(x, t) − h(y, t)⟩ ≡ 0, so daß die Definition (11) sinnvoll ist. Die Breite nimmt von w(0)=0 ausgehend zu.

Invarianz der Dynamik unter bestimmten Tranformationen findet man wieder in den Korrelationsfunktionen. Sind Anfangsbedingungen und Rauschen invariant gegen Translation des Fußpunkts x, so ist die Höhendifferenzkorrelationsfunktion

G_{2}(x, y, t) = G_{2}(x−y, t)

nur von der Differenz x−y abhängig.

2.1.2. Skalenverhalten. Über einen weiten Zeitbereich ist der Zustand des Systems nicht stationär, z.B. wächst w(t) an. Dennoch besteht

2.1. LANGEVINGLEICHUNGEN 13

in vielen Fällen eine Ähnlichkeit zwischen den Zuständen zu zwei aufeinanderfolgenden Zeiten, die man als statistische Skaleninvarianz bezeichnet. Sie besagt, daß die Oberfläche h(x, t) und ihr reskaliertes Gegenstück

h̃(x, t) = b^{−ζ} h(bx, b^{z}t) (12)

dieselben statistischen Eigenschaften besitzen. Dies gilt über einen weiten Bereich von Längenskalen zwischen der mikroskopischen Länge der Oberfläche (Gitterkonstante eines Kristalls, Persistenzlänge einer Membran) a und ihrer makroskopischen Ausdehnung L ≫ a. Über die Skaleneigenschaften sind auch der Rauhigkeitsexponent ζ und der dynamische Exponent z definiert. Es folgt inbesondere, daß G_{2} eine Skalenform besitzt, denn

G_{2}(x, t) = b^{−2ζ}G_{2}(bx, b^{z}t) = t^{2ζ/z}g_{2}(x/t^{1/z}) (13)

mit der speziellen Wahl b = t^{−1/z} und der Definition g_{2}(x) := G_{2}(x, 1). In das Argument der Skalenform geht die räumliche Abhängigkeit der Korrelationen im Verhältnis zu einer zeitlich anwachsenden Korrelationslänge ξ(t) ∼ t^{1/z} ein.

Bei Anwachsen von G_{2} und der Form von g_{2} unterscheidet man konventionelles und anomales Skalenverhalten. Bei konventionellem Skalenverhalten nimmt zur Zeit t die Verteilung der Höhendifferenzen auf Abständen |x−y| < ξ(t) einen stationären Zustand an. Bei anomalem Skalenverhalten saturiert die Verteilung nicht, solange ξ(t) anwächst (bis schließlich ξ die Systemgröße erreicht). In dieser Arbeit treten nur Beispiele konventionellen Skalenverhaltens auf, für eine allgemeinere Übersicht sei auf [12] verwiesen.

Aus der Zeitunabhängigkeit auf kleinen Abständen bei konventionellem Skalenverhalten folgt für g_{2}(x)

g_{2}(x) ∼ x^{2ζ}. für x ≪ 1.

Auf großen Abständen sind die Höhen h(x, t) und h(y, t) unkorreliert, und die Korrelationsfunktion G_{2}(x−y, t)=2 w(t)^{2} hängt nicht weiter vom Abstand der Fußpunkte x und y ab:

g_{2}(x) → const für x → ∞.

2.1.3. Ein Beispiel: Die lineare Gleichung. Verdeutlichen wir das Skalenverhalten am Beispiel der linearen Gleichung (6). Man betrachtet sie nach Fouriertransformation

∂_{t}h(q, t) = −ν q^{2}h(q, t) + η(q, t) (14)

14 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

mit dem transformierten Höhenprofil

h(q, t) = ∫ d^{d}x e^{−iq·x}h(x, t).

Die Transformierte des Rauschterms hat die Kovarianz

⟨η(q, t) η(q′, t′)⟩ = 2 D (q^{2})^{n} δ(q+q′) δ(t−t′).

(Für Gleichung (6) ist natürlich n = 0.) Nach der Anfangsbedingung h ≡ 0 ergibt sich die Höhenkonfigration zu späteren Zeiten aus der Geschichte des Rauschens

h(q, t) = ∫_{0}^{t}ds e^{−νq^{2}(t−s)}η(q, s), (15)

und daraus die Korrelationsfunktion

⟨h(q, t)h(q′, t′)⟩ = \frac{D}{νq^{2}} (1 − e^{−2νq^{2}t}) δ(q+q′). (16)

Die Korrelationsfunktion der Höhendifferenzen im Ortsraum ergibt die entsprechende Rücktransformation

(17) G_{2}(x, t) = \frac{2}{(2π)^{2d}} ∫ d^{d}q (1 − cos q·x)\frac{D}{νq^{2}} (1 − e^{−2νq^{2}t}) .

Die Skalenfunktion in G_{2}(x, t) = ξ(t)^{2−d}g_{2}(x/ξ(t)) ist

g_{2}(x/ξ) = \frac{2}{(2π)^{2d}} \frac{D}{ν} ∫ d^{d}κ (1 − cos \frac{κ·x}{ξ})\frac{1 − e^{−κ^{2}}}{κ^{2}} (18)

mit der Korrelationslänge ξ(t) = \sqrt{2νt}. Die im letzten Abschnitt definierten Skalenexponenten liest man in diesem Beispiel aus (17) ab:

z = 2 und ζ = \frac{2 − d}{2}. (19)

2.1.4. Galilei–Invarianz. Offensichtliche Invarianzen der Gleichungen (6) und (7) betreffen Transformationen jeweils im Raum (Translation, Rotation, Inversion) oder der Zeit (Translation). Gleichung (7) ist darüber hinaus invariant unter einer Transformation, die Raum und Zeit verknüpft. Für das Vektorfeld u(x, t)=−∇h(x, t)/λ wird Gleichung (7) zur Burgers– Gleichung mit Rauschen

∂_{t}u + (u · ∇)u = ν Δu − \frac{1}{λ} ∇η. (20)

Die linke Seite bezeichnet die Zeitableitung des Feldes an mit der Geschwindigkeit u(x, t) mitfließenden Punkten. Gleichung (20) ist daher invariant gegen Galilei–Transformationen

u(x, t) → u(x−tu_{0}, t) + u_{0}, (21)

2.2. STRUKTURBILDUNG UND -VERGRÖBERUNG 15

außer wenn die Korrelationen von η(x, t) diese Symmetrie brechen. Bei Korrelationen wie in (9) sind alle Raum- und Zeitpunkte gleichwertig und unkorreliert, und die Invarianz erstreckt sich auch auf das Rauschen.

Das Gradientenfeld u wie in (21) in einem bewegten System zu betrachten ist äquivalent zu einer Verkippung (genauer: Scherung) des Höhenprofils

h(x, t) → h(x−tu_{0}, t) − \frac{x · u_{0}}{λ} + \frac{1}{2λ} u_{0}^{2} t, (22)

mit zusätzlicher gleichmäßiger Erhöhung der mittleren Wachstumsgeschwindigkeit. Das transformierte Feld in (22) erfüllt ebenso wie h Gleichung (7).

Galilei-Invarianz in Verbindung mit Skaleninvarianz hat eine wichtige Konsequenz für die Exponenten z und ζ. Die Transformation (22) muß auch für ein reskaliertes Höhenfeld gelten, so daß unter Reskalierung (12) der Koeffizient λ konstant bleibt. Der Vektor u_{0} hat zwei Bedeutungen. Einerseits bezeichnet er eine Neigung, wird also unter (12) reskaliert wie

u_{0} → b^{ζ−1}u_{0},

andererseits hat er die Bedeutung einer Geschwindigkeit, was

u_{0} → b^{1−z}u_{0}

erfordert. Beide Skalenfaktoren müssen identisch sein, weshalb die Relation

z + ζ = 2 (23)

gilt [6, 13]. Daraus ergibt sich (außer für Dimension d=2), daß Gleichung (7) andere Skalenexponenten besitzt als die der linearen Gleichung (19).

Die “Verbindung” von Raum und Zeit in der Galileitransformation zeigt sich im Vergleich der Dynamik auf mit u_{0}/λ verkippten Oberflächenstücken gegenüber waagrechten Teilen. Strukturen entwickeln sich gleich, nur auf der verkippten Oberfläche laufen sie mit der zusätzlichen Geschwindigkeit u_{0} bergab. Gleichung (8)besitzt dieselbe Symmetrie, was im folgenden Kapitel eine wichtige Rolle spielen wird.

2.2. Strukturbildung und -vergröberung

Strukturbildung, wie sie an Kontinuumsgleichungen in dieser Arbeit untersucht wird, setzt die Existenz eines instabilen Fixpunkts im Kofigurationsraum voraus, bei Gleichung (5) ist dies ϕ(x, t)≡0. Man verfolgt die Dynamik ausgehend von einer “kleinen” Störung dieses Fixpunkts, einer Anfangsbedinung ϕ_{0}(x) überall “nahe bei” 0, etwa mit einer Gaußverteilung mit ⟨ϕ_{0}(x)⟩ = 0 und ⟨ϕ_{0}(x)ϕ_{0}(y)⟩ = D δ(x−y) mit D ≪ 1. Die Allen-Cahn–Gleichung (5) besitzt zwei stabile Fixpunkte ϕ(x, t)≡±1, von denen

16 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

ϕ an jedem Punkt x zunächst einem zustrebt. Es entstehen durch Wände getrennte Domänen, in denen das Profil jeweils nahe bei einem der stabilen Werte liegt.

Die Anfänge der Strukturbildung werden von einer Linearisierung um den instabilen Fixpunkt erfaßt. Die Fouriertransformierte gehorcht

∂_{t}ϕ(q, t) = (1 − q^{2})ϕ(q, t). (24)

Die Anteile in ϕ_{0} mit Wellenvektoren q^{2} < q_{c}^{2} = 1 wachsen (zunächst exponentiell) an, während kurzwelligere Anteile weggedämpft werden. Die Domänen der anfangs entstehenden Strukturen haben daher eine Ausdehnung L>2π/q_{c}.

Die weitere Entwicklung besitzt ein Lyapunov-Funktional, die Freie Energie einer Konfiguration ϕ,

F\{ϕ\} = ∫ d^{d}x [\frac{1}{2} (∇ϕ)^{2} + V (ϕ)] ,

wobei V (ϕ) = −\frac{1}{2}ϕ^{2}+ \frac{1}{4}ϕ^{4}. Gleichung (5) ist glechbedeutend mit

∂_{t}ϕ = − \frac{δF\{ϕ\}}{δϕ}, (25)

bezeichnet also die von F getriebene Relaxation des Systems in einen Zustand minimaler Freier Energie.

Selbst für sehr große Domänen, die man bei fortschreitender Relaxation beobachtet, bleibt die Dicke der Wände endlich. Eine spezielle stationäre Lösung von (25) macht dies deutlich, nämlich eine “Stufe”, ϕ(x, t) = 0 auf der Hyperebene x_{1} = 0 mit ϕ(x, t) → ±1 für x_{1} → ±∞. Fernab der Wand nähert sich das stationäre Profil wie exp (−x_{1}/V′′(±1)^{1/2}) den Minimalwerten an. Die Dicke der Wand ist daher von der Größenordnung 1. Für eine detailliertere Darstellung, auch der nun folgenden Bemerkungen, verweisen wir auf [3].

Bei einer solchen Stufenkonfiguration ist der wesentliche Beitrag zu F{ϕ} in der Nähe der Wand lokalisiert. Das System relaxiert vor allem dadurch, daß es den Flächeninhalt der Domänenwände verringert. Wie für ein anderes Problem in Kapitel 5 im Zusammenhang mit Gleichung (60) hergeleitet wird, geschieht dies dadurch, daß sich die Wände überall entlang ihrer Normalen mit einer Geschwindigkeit proportional zur Krümmung verschieben. Bild 2.1 zeigt die Relaxation von (5) im Anfangsstadium, und zu drei späteren Zeiten, wo sie hauptsächlich durch Verschiebung der Domänenwände vorangetrieben wird. Dargestellt ist jeweils die “Höhenlinie” ϕ ≡ 0. Nimmt man nun an, daß in der Konfiguration der Wände nur eine Abbildung 2.1. Relaxation nach Gleichung (5) dargestellt durch die “Höhenlinie” ϕ = 0 zu verschiedenen Zeiten: Oben links die anfangs entstandene Struktur, gut beschrieben durch Gl. (24). Sind die Domänen angewachsen und deutlich größer als die trennenden Wände (unten links und rechts), ist die Skalenannahme gut erfüllt.

18 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

Längenskala ξ(t) bestimmend ist, was Simulationen gut bestätigen [3], so ist ein Richtwert für die Krümmung der Wände 1/ξ, und

ξ̇ ∼ \frac{1}{ξ} bzw. ξ ∼ t^{1/2}

Zu späten Zeiten ist diese Längenskala viel größer als die Wanddicke, ξ ≫1. Sie ist die einzige Länge, die in die Verteilung der Oberflächenkonfigurationen eingeht. Dies legt die Annahme einer Skalenhypothese nahe, die durch Experimente und Simulationen gerechtfertigt ist und für die weitere analytische Behandlung hilfreich. Anschaulich bedeutet dies, daß die drei letzten Konfigurationen in Abbildung 2.1 jeweils durch Vergrößerung um einen festen Faktor ineinander übergehen. Wie aus den Folgerungen des Ansatzes (12) ergeben sich daraus für alle Korrelationsfunktionen Skalenformen, so zum Beispiel für die gleichzeitige Zweipunktkorrelation

C(x, t) = ⟨ϕ(x + y, t) ϕ(y, t)⟩ = f(|x|/ξ).

Am Beispiel der Allen-Cahn–Gleichung haben wir die Punkte skizziert, die wir im Kapitel 6 bei der Modellierung von Kristalloberflächen bei epitaktischem Wachstum wiederfinden werden: Linearisierung um einen homogenen, instabilen Anfangszustand liefert eine Vorhersage über die zunächst entstehenden Strukturen. Die Dynamik ist dissipativ insofern, als zu jeder Konfiguration eine Freie Energie definiert ist, die ein Lyapunov-Funktional der Dynamik ist. Im Limes großer Zeiten t ist die Verteilung der Oberflächenkonfigurationen offenbar nur noch von einer Länge ξ(t) abhängig, hier die Ausdehnung der Domänen konstanter Magnetisierung, bei wachsenden Kristallen die Größe der Hügel auf der Oberfläche. Das Skalenverhalten ist (außer in Ausnahmefällen) nicht rigoros aus den Gesetzen der Dynamik ableitbar, wenn auch empirisch beobachtet und wertvolle Grundlage für weitere theoretische Vorhersagen.

2.3. Chaotische Dynamik

Grundsätzlich anders verhalten sich Systeme mit nichtdissipativer Dynamik. Eine analog zum letzten Abschnitt definierte Freie Energie nimmt hier nicht stetig ab, vielmehr gibt es ein Wechselspiel zwischen Energie erzeugenden und dissipierenden Mechanismen. Von chaotischer Dynamik spricht man, wenn das System folgende Eigenschaften aufweist. Zwei Trajektorien im Konfigurationsraum, die zu einem Zeitpunkt t sehr nahe beieinander liegen, gehen im weiteren Verlauf exponentiell auseinander. Den

2.3. CHAOTISCHE DYNAMIK 19

(endlich vielen) Lyapunovexponenten in einem dynamischen System endlicher Dimension entspricht in räumlich ausgedehnten Systemen die Spektraldichte der Eigenwerte des Operators der linearen Antwort (linear response function) auf eine externe Störung [14]. Wenn es so etwas wie “typische” Konfigurationen zu späten Zeiten gibt, liegen sie auf einem chaotischen Attraktor im Konfigurationsraum und nicht auf einem stabilen Grenzzyklus oder Fixpunkt wie im letzten Abschnitt.

Mit diesen allegemeinen Eigenschaften bieten Oberflächen mit chaotischer Dynamik eine Verbindung zwischen den beiden Gruppen der letzten Abschnitte. Die Forderung, daß sich benachbarte Trajektorien voneinander entfernen sollen, bedeutet daß jeder Fixpunkt der Dynamik instabil ist. Eine glatte, unter der Dynamik invariante Oberfläche sollte daher unter einer kleinen Störung komplizierteren Strukturen weichen, gleich dem im letzten Abschnitt beschriebenen Szenario zu frühen Zeiten.

Im Spätstadium gleicht die Dynamik oft der von Langevingleichungen. Dem Rauschen entspricht die “unvorhersehbare” Entwicklung einer einzelnen Trajektorie im Konfigurationsraum, ihre Bindung an einen Attraktor bewirkt eine globale Stabilisierung und entspricht den glättenden Termen der stochastischen Modelle. Man versucht daher zu einer gegebenen chaotischen Gleichung eine effektive stochastische Beschreibung, eine Langevingleichung mit denselben makroskopischen Eigenschaften zu finden.

Die hier angesprochenen Analogien können in dieser Allgemeinheit nicht anders als vage sein, allein weil die Bezeichnung “chaotische Oberflächendynamik” ein weites Feld umfaßt. Präzisere Fragestellungen und Antworten ergeben sich am Einzelbeispiel, was im nächsten Kapitel am Beispiel von Gleichung (8), der sogenannten Kuramoto-Sivashinsky–Gleichung (KS), unter den hier angesprochenen Aspekten – Fehlen eines Lyapunovfunktionals, Strukturbildung im Frühstadium und effektive Dynamik im Spätstadium – illustriert werden soll.

20 2. OBERFLÄCHENDYNAMIK IN PARTIELLEN DIFFERENTIALGLEICHUNGEN

KAPITEL 3 Die Kuramoto-Sivashinsky–Gleichung

3.1. Grundlegende Bemerkungen

Die Kuramoto-Sivashinsky–Gleichung (8) ist ein einfaches Modell, das räumlich ausgedehntes chaotisches Verhalten zeigt. In seinen ursprünglichen Herleitungen beschreibt es z.B. oszillierende chemische Reaktionen [7], Flammenfronten in vorgemischten reagierenden Gasen [8] oder das Profil eines dünnen viskosen Flüssigkeitsfilms, wenn er eine Wand hinabfließt [15]. Im Falle einer Flammenfront ist die Modellierung anschaulich: Das Höhenprofil h(x, t), wie in Abbildung 3.1 dargestellt, beschreibt die Höhe der Flammenfront zur Zeit t über einem Fußpunkt x am Boden des die reagierenden Gase einschließenden Behälters.

3.1.1. Nichtdissipative Dynamik. Ersetzen wir die Höhe, die laterale Länge und die Zeit durch die dimensionslosen Raum- und Zeitvariablen \sqrt{ν/K}x, (λ/ν)h und (ν^{2}/K)t, so wird Gleichung (8) zu

∂_{t}h = −Δ^{2}h − Δh + \frac{1}{2} (∇h)^{2} . (26)

Wie beim Übergang von (7) zu (20) gilt für das Gradientenfeld u=−∇h

∂_{t}u + (u · ∇)u = −Δ^{2}u − Δu,

eine Gleichung, die unter der Galileitransformation

u(x, t) → u(x−tu_{0}, t) + u_{0}

invariant ist, wobei man wiederum u(x, t) als Geschwindigkeitsfeld interpretiert. Diese Interpretation legt die Wahl einer “kinetischen Energie”

E(t) = \frac{1}{2} ∫ d^{d}x u(x, t)^{2} = \frac{1}{2} ∫ d^{d}x (∇h(x, t))^{2}

nahe. Gleichung (8) in Fourierdarstellung zeigt, in welchem Bereich Energie erzeugt und wo sie dissipiert wird,

(27) ∂_{t}h(k, t) = (k^{2} − k^{4}) h(k, t) − ∫ \frac{d^{d}q}{(2π)^{d}} q·(k−q) h(q, t) h(k−q, t).

22 3. DIE KURAMOTO-SIVASHINSKY–GLEICHUNG

Großskalige Fluktuationen mit |k| < 1 werden verstärkt, und so dem System Energie zugeführt. Der nichtlineare Term (∇h)^{2} koppelt über das Faltungsintegral die einzelnen Moden h(k, t) untereinander und transportiert die Anregungen in den Bereich |k|>1, wo Energie dissipiert wird.

3.1.2. Linearisierung und Strukturbildung. Läßt man den Faltungsterm in (27) weg, erhält man eine Linearisierung um den instabilen Fixpunkt h≡const. Anders als bei der linearisierten Allen-Cahn–Gleichung in (24) wächst hier jede Mode mit der Rate q^{2}−q^{4} an, am stärksten solche mit Wellenvektoren |q| = q_{m} = 1/\sqrt{2}. Beginnt die Dynamik in der Nähe des instabilen Fixpunkts, so ist die lineare Näherung von (8) gut, bis die maximalen Gradienten von der Größenordnung (∇h)^{2} ≃1 werden. In dieser Frühzeit der Dynamik sind vor allem die am meisten instabilen Moden mit |q|=1/\sqrt{2} angewachsen und bestimmen die Struktur. Die erste Abbildung in 3.1 zeigt eine solche Konfiguration, die aus unregelmäßig angeordneten Buckeln der Größe ∼ 2π/q_{m} ≃9 besteht.

Die Moden der Wellenlänge q_{m} bestimmen auch weiter die Struktur der Oberfläche, die zu einem späteren Zeitpunkt im zweiten Teilbild von 3.1 gezeigt ist. Qualitativ verläuft die Dynamik der Oberfläche auf mikroskopischer Skala wie folgt: Der nichtlineare Term in (8) läßt die Flanken der Buckel stärker wachsen, die Buckel werden also breiter. Sie verbreitern sich soweit, bis ihre Kappen eine Größe erreichen, die im Bereich der instabilen Wellenlängen liegt, und brechen dann auf.

3.1.3. Verhalten auf großen Skalen. Unregelmäßigkeiten im Anwachsen und Aufbrechen der Buckel führen zu Fluktuationen auf großen Skalen, die zudem vom linearen Anteil der Dynamik verstärkt werden. Dem tritt folgender glättender Effekt entgegen: Von einer großskaligen Wölbung laufen, wie im letzten Abschnitt bei (22) beschrieben wurde, an den Seiten die Buckel hinab. Dadurch werden die Buckel oben auf der Wölbung gedehnt und in den Senken daneben zusammengedrückt [10]. In der Zone gedehnter Buckel ist ⟨(∇h)^{2}⟩ und damit die mittlere Wachstumsgeschwindigkeit kleiner, gestauchte Buckel erhöhen dagegen die Wachstumsgeschwindigkeit in einem Gebiet. Die zurückgebliebenen Senken können daher die vorangegangenen Wölbungen wieder einholen. In 1+1 Dimensionen ist die Stabilität der Lösung von (8) auf einem endlichen Intervall mit periodischen Randbedingungen bewiesen: Abhängig von der Systemgröße existiert eine Schranke für ⟨(∇h)^{2}⟩ [16].

Abbildung 3.1. Zwei “Schnappschüsse” des Höhenprofils h(x, t) der KS–Gleichung (26). Ausgehend von einer nahezu glatten Oberfläche (siehe Text) wachsen zunächst Moden mit q_{m} an, der nichtlineare Term bleibt noch ohne Einfluß (oben). Typische Konfigurationen im späten Stadium (unten) sind wesentlich von der Nichtlinearität geprägt: Die Symmetrie h → −h besteht nicht mehr. Zellen in Senken werden zusammengedrückt, auf Erhebungen verbreitern sie sich und brechen wieder auf (vgl. auch Abbildung 3.2).

24 3. DIE KURAMOTO-SIVASHINSKY–GLEICHUNG

Zu den oben beschriebenen Merkmalen von Langevingleichungen besteht eine qualitative Verwandtschaft, wenn man (8) geglättet, d.h. räumlich gemittelt über mehrere Zellen betrachtet. Das schnelle Aufbrechen und Entstehen neuer Zellen ähnelt dann zufälligem Rauschen η(x, t), der glättende Effekt der Wechselwirkung der Zellen entspricht in führender Ordnung der Relaxation durch eine effektive Oberflächenspannung, ν Δh.

Ein quantitativer Test für den Vergleich einer chaotischen Oberfläche mit einer “effektiven” stochastischen Beschreibung durch eine Langevingleichung ist das Skalenverhalten im infraroten Limes. Besteht eine Skaleninvarianz wie in Gleichung (12), und, wenn ja, bei welcher Langevingleichung tritt sie mit denselben Exponenten auf? Mit der Bezeichnung aus der Theorie der kritischen Phänomene sucht man nach der Universalitätsklasse der Oberfläche.

Ein Indiz für eine Antwort liefert bei Gleichung (8) die Galilei-Invarianz. Sie bedingt die Skalenrelation z + ζ = 2 in Gleichung (7) [6, 13]. Insofern sind außer in Dimension d=2 die Skalenexponenten (19) der linearen Gleichung (6) verboten. Verschiedene Arbeiten postulieren, daß die Kuramoto-Sivashinsky–Gleichung zur Universalitätsklasse von (7) gehört [17, 18]. In 1+1 Dimensionen ist dies durch detaillierte Untersuchung der Dynamik auf großen Skalen [9] sowie der kurzskaligen Wechselwirkung der Zellen [10] bestätigt.

Die Zahlreichen Ansätze auf der Grundlage einer Renormierungsgruppenanalyse analog zu [6] scheinen zur theoretischen Lösung dieses Problems ungeeignet [17, 19, 20]. Zwei grundsätzliche Schwierigkeiten bestehen: Die Ansätze sind perturbativ, auch wenn sie im physikalisch interessanten Fall von einem instabilen linearen Fixpunkt ausgehen müssen. Daraus ergibt sich insbesondere für eine Impulsschalenrenormierung das Problem, wie jene kritische Wellenzahl zu behandeln ist, die den linear stabilen vom instabilen Bereich trennt. Die zweite Schwierigkeit besteht darin, daß zur Mittelung über kurzwellige Fluktuationen ein zusätzlicher Rauschterm angesetzt wird, obwohl unklar ist, ob er die charakteristischen Eigenschaften der Dynamik erhält, und obwohl eigentlich über ihr inhärentes raumzeitliches “Rauschen” (der chaotischen Zellendynamik) gemittelt werden müßte. Großskalige Stabilität und Unordnung der KS–Oberfläche sind wesentlich nichtlineare Effekte, die im Renormierungsgruppensinn von einem Fixpunkt starker Kopplung bestimmt sind, der mit oben angeführten Methoden nicht zugänglich ist.

3.2. EIN TEILCHENMODELL 25

Im folgenden Abschnitt wird auf die Frage der makroskopischen Dynamik im Rahmen eines Modells für die Wechselwirkung der Zellen eingegangen. In 2+1 Dimensionen ist diese Frage mehrfach untersucht und kontrovers aber nicht befriedigend beantwortet [21, 22]. Ein interessanter Ansatzpunkt zur Lösung dieses Gegensatzes bietet sich in der Verallgemeinerung auf eine anisotrope Gleichung, die im übernächsten Abschnitt 3.3 behandelt wird.

3.2. Ein Teilchenmodell

3.2.1. Definition des Modells. Die Dynamik der Kuramoto-Siva-shinsky–Gleichung kann, wie im letzten Abschnitt skizziert, über die Wechselwirkung der Zellen verstanden werden [10]. Ihre Anzahl in einem endlichen System entspricht der effektiven Anzahl von Freiheitsgraden [14]. Es liegt daher nahe, ein dynamisches System zu konstruieren, in dem jede Zelle durch eine (zeitlich veränderliche) Variable repräsentiert wird. Hier soll dies durch ein System von wechselwirkenden klassischen Teilchen geschehen.

Bild 3.2 zeigt die Entwicklung von Gleichung (8) in 1+1 Dimensionen und die Trajektorien der Minima zwischen den Zellen, die wir als “Teilchen” auffassen wollen. Die Wechselwirkung benachbarter Zellen zeigt sich auch in den Trajektorien, sie beeinflussen sich gegenseitig. Wenn sich eine Zelle spaltet, entsteht ein neues Teilchen zwischen zwei bereits bestehenden, wenn eine Zelle verschwindet, verschmelzen zwei Teilchen zu einem.

Ein wichtiger Anhaltspunkt zur Modellierung einer Teilchendynamik ist die Galilei-Invarianz von (8). Neigt man die Oberfläche gleichmäßig, verändert sich die Dynamik nur insofern, als zur Bewegung aller Zellen dieselbe Geschwindigkeit addiert wird, nicht aber ihre Relativbewegungen. Daher muß die Teilchendynamik galileiinvariant formuliert werden: Die Teilchen besitzen eine träge Masse, die aufeinander ausgeübten Kräfte bewirken Beschleunigungen. Bezeichnen wir die Position der Teilchen (Minima zwischen Zellen) mit x_{i}(t) und ihre Geschwindigkeit (bzw. Impuls, wenn man ihnen die Masse m=1 gibt) mit p_{i}(t), so ist die Dynamik in allgemeinster Form

ẋ_{i} = p_{i}
(28)  ṗ_{i} = F (x_{i+1}−x_{i}, x_{i}−x_{i−1}).

Eine Vereinfachung besteht in der Annahme lokaler Wechselwirkungen. Die Translationsinvarianz erfordert dann, daß die auf ein Teilchen wirkende Kraft nur von den Abständen zu den nächsten Nachbarn abhängt, wie in (28) angegeben. Sich teilende und verschwindende Zellen werden wie folgt modelliert: Stoßen zwei benachbarte Teilchen zusammen, werden sie durch Abbildung 3.2. Zeitliche Entwicklung der KS–Gleichung (26) in 1 + 1 Dimensionen. Höhenprofile h(x, t) zu Zeiten t=1, 2, …, 200 in gemeinsamer Darstellung zeigen die Zellendynamik (oben): Sie teilen sich auf Erhebungen, verschinden in Senken, und laufen an Flanken bergab. Trajektorien der Minima zwischen Zellen (unten) ähneln den Trajektorien im Teilchenmodell (Bild 3.3).

ein einzelnes am Ort x = x_{i} = x_{i−1} ersetzt, sind sie “zu weit” voneinander entfernt, x_{i}−x_{i−1} = 1, wird ein neues am Ort x = (x_{i}+x_{i−1})/2 erzeugt. Mit dieser Definition wird die minimale instabile Wellenlänge auf 1 festgelegt.

3.2. EIN TEILCHENMODELL 27

Die den neuen Teilchen zugeordneten Impulse p müssen sinnvoll gewählt werden. Wegen der Galileiinvarianz muß p ein gewichteter Mittelwert aus den Impulsen der kollidierenden bzw. erzeugenden Teilchen sein, die einzig mögliche Wahl werden wir nach einigen Überlegungen treffen.

Aus den Geschwindigkeiten p_{i} der Teilchen definiert man ein makroskopisches Geschwindigkeitsfeld

u(x, t) mit u(x_{i}(t), t) = p_{i}(t),

dessen großskalige Eigenschaften mit (8) verglichen werden. Eine natürliche Wahl besteht in einer stückweise konstanten (Stufen-) Funktion

u(x, t) := p_{i} für x ∈ [\frac{x_{i−1}+x_{i}}{2} , \frac{x_{i}+x_{i+1}}{2} ). (29)

Das Geschwindigkeitsfeld ist der Gradient eines zugehörigen Höhenprofils ∂_{x}h(x, t) = −u(x, t), das daher aus einem Zug geneigter Geradenstücke mit Neigungen −p_{i} und Eckpunkten bei (x_{i}+x_{i+1})/2 besteht.

Die Eigenschaft von h als Höhenprofil bedeutet eine wichtige Einschränkung für die Teilchendynamik, für die Kraft F und für die Regeln nach denen neue Teilchen erzeugt und kollidierende Paare durch ein einzelnes Teilchen ersetzt werden. Die Dynamik soll invariant unter einer Verschiebung h → h+const sein, es ist also

∂_{t}h = Φ(∂_{x}h, ∂_{x}^{2}h, …).

Weil dies äquivalent ist zu

∂_{t}u + ∂_{x}Φ(u, ∂_{x}u, …) = 0,

ist das Geschwindigkeitsfeld eine erhaltene Größe. Auf jedem Intervall hängt die Änderung von ∫_{a}^{b} u(x, t) = h(a, t) − h(b, t) nur lokal von der Konfiguration bei a und b ab, in einem endlichen System mit periodischer Randbedingung ist

∫_{0}^{L}dx u(x, t) = \frac{1}{2} ∑_{i}(x_{i} − x_{i−1})(p_{i} + p_{i−1}) (30)

eine Konstante der Dynamik.

Dies ist gewährleistet, wenn die neuen Geschwindigkeiten bei einer Kollision von den Abständen zu den Nachbarn abhängt. Mit der Bezeichnung ξ_{i} := x_{i}−x_{i−1} ist dies

p = \frac{ξ_{i+1}p_{i} + ξ_{i−1}p_{i−1}}{ξ_{i+1} + ξ_{i−1}} .

Bei einer Erzeugung eines neuen Teilchens genügt

p = \frac{p_{i} + p_{i−1}}{2}

obiger Bedingung. Auch für die Kraft Φ(ξ_{i+1}, ξ_{i}) in (28) besteht wegen der Erhaltung von u eine Einschränkung, nämlich

Φ(ξ_{i+1}, ξ_{i}) = \frac{φ(ξ_{i+1}) − φ(ξ_{i})}{ξ_{i+1} + ξ_{i}} (31)

mit einer beliebigen “Potentialfunktion” φ. Man prüft leicht nach, daß (30) damit erhalten bleibt.

Je nach Form von φ(ξ) ist eine periodische Folge von Teilchen mit Abständen a stabil oder instabil gegen Störungen. Linearisierung der Bewegungsgleichungen (28) ergibt die Wachstumsrate ω(k) einer Störung mit Wellenzahl k,

ω(k)^{2} = \frac{φ′(a)}{a} (1 − cos ak).

Eine periodische Abfolge von Zellen ist stabil, wenn φ′(a) < 0. In diesem Fall reagiert das System mit Schallwellen auf eine Störung. Der analoge Fall tritt auch bei der Kuramoto-Sivashinsky–Gleichung auf. Bei bestimmten Zellengrößen reagieren die Muster elastisch auf Verzerrungen [23].

Im folgenden wird ein System untersucht, das keinen stabilen Bereich besitzt, um die chaotischen Eigenschaften der Dynamik deutlicher hervorzuheben. Die Wahl φ(ξ^{2}) = ξ^{2} führt zu besonders einfachen linearen Gleichungen der Teilchendynamik,

ẋ_{i} = p_{i}
(32)  ṗ_{i} = 2x_{i} − x_{i+1} − x_{i−1},

mit denen die im folgenden präsentierten Simulationen ausgeführt wurden. In Abbildung 3.3 sind die Trajektorien bei linearer Kraft zwischen benachbarten Teilchen dargestellt. Sie ähneln den Trajektorien der Minima zwischen Zellen in Bild 3.2.

In den Simulationen wurde als Anfangskonfiguration für alle Teilchen p_{i} ≡ 0 mit zufälligen Abständen ξ_{i} ∈ (0, 1) gewählt, letztere unabhängig und in gleicher Weise mit konstanter Verteilung. Die Integration der linearen Dynamik zwischen den Stößen und Paarerzeugungen wurde sowohl in erster als auch in zweiter Ordnung im Zeitschritt τ genähert. Symbolisch geschrieben wurde bei der Dynamik Ẋ = AX der Propagator exp τA ≃ 1 + τA(+(1/2)(τA)^{2}) angenähert, wobei sich zwischen linearer und quadratischer Näherung keine nennenswerten Abweichungen ergaben. Empfindlicher ist die Dynamik gegen die Genauigkeit, mit der die Zeitpunkte vor allem der Kollisionen erfaßt werden. Es erwies sich als angemessen, Abbildung 3.3. Trajektorien der Teilchen mit linearen Kräften (32).

den Zeitschritt nicht größer als τ = 0.01 zu wählen, damit die makroskopischen Größen, die im nächsten Abschnitt betrachtet werden, verläßlich bestimmt werden konnten.

3.2.2. Statistik des Teilchenmodells auf großen Skalen. Die Annahme statistischer Skaleninvarianz der Teilchendynamik bedeutet für den über (29) definierten Höhenzug zunächst, daß mit Reskalierung (12)

⟨(h̃(x, t) − h̃(0, t))^{2}⟩ = b^{−2ζ}⟨(h(bx, b^{z}t) − h(0, b^{z}t))^{2}⟩.

Andererseits geht in die Höhe ein Produkt aus lateraler Länge und Geschwindigkeit ein, also h̃(x, t) = b^{z−2}h(bx, b^{z}t), weshalb die Korrelationsfunktion auch skaliert wie

⟨(h̃(x, t) − h̃(0, t))^{2}⟩ = b^{2(z−2)}⟨(h(bx, b^{z}t) − h(0, b^{z}t))^{2}⟩.

Direkter Vergleich ergibt wie schon bei Gleichung (7) die Skalenrelation z + ζ = 2.

Statistische Skaleninvarianz zeigt sich bei der Kuramoto-Sivashinsky– Gleichung (8) erst ab einer Längenskala von ungefähr dreihundert Zellen [9, 10], im Teilchenmodell schon ab einer lateralen Länge x ∼ 50, was knapp hundert mittleren Teilchenabständen entpricht. Abbildung 3.4 zeigt die Korrelationsfunktion der Höhendifferenzen G_{2} (11) zu Zeiten 10^{1,...,5} in Abbildung 3.4. Korrelationsfunktion der Höhendifferenzen im System der Größe L = 1000. Ab t = 1000 zeigt sich KPZ–Verhalten: G_{2}(r, t) = (D/ν)|r| für |r| ≪ ξ und das Plateau jenseits von ξ nimmt wie t^{2/3} zu.

einem System der Länge L = 1024. Die Daten wurden über fünfzig Läufe gemittelt. Für kleine Abstände x < 10 und frühe Zeiten t < 10 ist das System eindeutig nicht in einem skaleninvarianten Bereich. Zu späten Zeiten wächst G_{2} mit dem Abstand wie G_{2}(x, t) ∼ x^{2ζ}, wobei ζ nahe bei 1/2 liegt, dem Wert der KPZ–Universalitätsklasse [6]. Auf Abständen jenseits der lateralen Korrelationslänge ξ(t) erreicht G_{2} ein Plateau, das zu späten Zeiten t≥10^{3} wie t^{2ζ/z} =t^{2/3} anwächst.

Bestimmen wir nun die effektiven Parameter ν, λ und D, mit denen Gleichung (7) die Teilchendynamik auf großen Skalen wiedergibt. Aufgrund der Konstruktion der Höhenlinie aus einer galileiinvarianten Dynamik ist λ = 1. Das Verhältnis D/ν bestimmt in d = 1 den stationären Zustand in einem endlichen System: Für das statistische Gewicht einer Konfiguration gilt P{h} ∼ exp[−(ν/2D) ∫ dx (∂_{x}h)^{2}/2]. Der Strukturfaktor des stationären Zustands ergibt sich daraus zu

⟨|h(k, t)|^{2}⟩ = \frac{D}{Lk^{2}ν} . (33)

3.2. EIN TEILCHENMODELL 31

Im großskaligen Limes des Teilchenmodells ist diese Größe ungefähr 1.0 × 10^{−2}/(Lk^{2}), wie in Bild 3.5 für Systeme verschiedener Gößen L = 64, 128, 256, 512 dargestellt ist.

Abbildung 3.5. Die stationären Fluktuationen der Oberfläche des Teilchenmodells für verschiedene Systemgrößen. Im skaleninvarianten Bereich kleiner k ist Lk^{2}⟨|h(k, t)|^{2}⟩ ≡ D/ν konstant.

Die Struktur auf kurzen Skalen (für große k) zeigt einen Abfall, der von der willkürlichen mikroskopischen Definition der Höhenlinie herrührt. Der deutlich ausgeprägte Buckel ist vergleichbar mit der Beobachtung der typischen Zellengröße der Kuramoto-Sivashinsky–Gleichung (8) im Strukturfaktor [18]. Da die Teilchen in jeweils benachbarten Paaren aufeinander zulaufen (siehe Bild 3.3) entsteht bei der verwendeten Konstruktion des Höhenprofils eine Zickzackline mit der Periode des doppelten mittleren Zellenabstands 2ξ̄.

Der Quotient D/ν ist auch im Anstieg der Höhendifferenz-Korrelations-funktion enthalten:

G_{2}(x, t) = \frac{D}{ν} |x|, für ξ̄ ≪ |x| ≪ ξ(t),

woraus sich mit den Daten in Bild 3.4 derselbe Wert wie aus dem Strukturfaktor des stationären Zustands ergibt.

Mit der Kenntnis von λ und D/ν läßt sich das Anwachsen der mittleren Breite im transienten Zustand ξ̄ ≪ ξ(t) ≪ L vorhersagen,

w(t)^{2} = c_{2} (\frac{D^{2}}{ν^{2}} λt)^{2/3} , (34)

wobei c_{2} ≃ 0.40 eine universelle Konstante ist [24]. Die Entwicklung der Breite der Höhenlinie des Teilchenmodells folgt genau dieser Vorhersage (Bild 3.6). Auffällig ist, daß sich kein deutliches Übergangsverhalten von einem Regime, in dem die Statistik von der linearen Gleichung (6) beschrieben wird, zum asymptotischen Bereich mit dem hier beschriebenen KPZ–Verhalten zeigt. Anders die Kuramoto-Sivashinsky–Gleichung in 1+1 Dimensionen: Hier wächst w(t)^{2} über etwa zwei Dekaden in der Zeit wie t^{1/2} an, um dann erst zum asymptotischen Anwachsen mit t^{2/3} überzugehen.

Abbildung 3.6. Zunahme der Breite des im Teilchenmodell definierten Höhenprofils. Die Aufrauhung zeigt KPZ– Verhalten w^{2} = 0.4((D^{2}λ/ν^{2})t)^{2/3}.

Beenden wir die statistische Auswertung des Teilchenmodells mit einer Abschätzung der “nackten” Kopplungskonstanten

g = \frac{λ^{2}D}{ν^{3}}, (35)

die in eine effektive Modellierung mit Gleichung (7) eingehen müßte. Die Kopplung g hat in d=1 die Dimension einer inversen Länge, die wir mit der Wellenzahl k_{max} =: 2π/ℓ des Buckels im Spektrum (Bild 3.5) vergleichen.

Aus dem Verhältnis der lateralen Korrelationslänge der KPZ–Gleichung (7)

ξ(t) = \sqrt{Dν}(λt)^{2/3}

zur intrinsischen Längenskala ℓ ergibt sich eine intrinsische Zeitskala

τ := t(ℓ/ξ(t))^{3/2} = \frac{ℓ^{3/2}}{λ\sqrt{D/ν}} .

Mit den oben bestimmten effektiven Parametern hat das Teilchenmodell τ ≈ 36.3. Da KPZ–Skalenverhalten deutlich ab t = 1000 oder t/τ = 27.5 beobachtet wird, schätzt man die effektive KPZ–Kopplungskonstante auf mikroskopischer Skala zu ℓg ≈ 3.5. Das Teilchenmodell liegt auf seiner mikroskopischen Skala “näher” am Fixpunkt starker Kopplung der KPZ– Gleichung als die KS–Gleichung (8), für die man ℓg ≈0.53 findet [10, 9].

3.2.3. Zusammenfassende Bewertung. Chaotisches Verhalten in einem räumlich ausgedehnten System wurde am Beispiel eines Gases wechselwirkender Teilchen vorgestellt. Das Modell wurde mit den Symmetrien der KS–Gleichung, gegen Verschiebung in h–Richtung und Galilei-Invarianz, und mit elementaren Eigenschaften der Dynamik der Zellen einer chaotischen Lösung von (8) konstruiert. Man findet, daß diese beiden Eigenschaften genügen auf großen Skalen stochastisches Verhalten der Burgers-KPZ–Universalitätsklasse zu erzeugen.

Ein Ansatz dieser Art dient dazu, die chaotische Dynamik ausgedehnter Systeme konzeptuell besser zu verstehen. Interessant ist der Gegensatz zwischen der mikroskopisch instabilen Dynamik und einer Stabilität auf großen Skalen, die sich in einer Effektiven Oberflächenspannung ν äußert. Am Beispiel der Kuramoto-Sivashinsky–Gleichung wurde dieser Vorgang im Abschnitt 3.1.3 des letzten Kapitels qualitativ betrachtet. In diesem Zusammenhang wurde der Begriff der durch Defekte vermittelten Stabilität [10] geprägt. Im Teilchenmodell dieses Kapitels sind die “Defekte” der chaotischen Oberflächendynamik unmittelbar greifbar, nämlich die Teilchen selbst. Ihre ungeordnete Bewegung bewirkt eine effektive Dissipation. Die Analogie zu Boltzmanns Herleitung der thermodynamischen Gesetze des idealen Gases, insbesondere zum Stoßzahlansatz fällt auf: Die Irreversibilität der Dynamik aus einem Nichtgleichgewichtszustand heraus, hydrodynamisch gesprochen die Viskosität, wird erzeugt von der Unordnung. Die inelastischen Stöße der Teilchen erzeugen dagegen keine Viskosität bzw. Oberflächenspannung. Man leitet leicht her, daß das Geschwindigkeitsfeld des Teilchenmodells ohne beschleunigende Kräfte, also F ≡ 0 in (28), der

34 3. DIE KURAMOTO-SIVASHINSKY–GLEICHUNG

Burgers–Gleichung ohne Viskosität folgt [11]. Abbildung 3.7 zeigt eine Simulation dieses Modells mit den typischen Schockfronten. Dies gilt auch in höheren Dimensionen, was in d = 3 zur Modellierung des Zusammenklumpens von Gas- und Staubwolken im Universum verwendet wurde. Deren Verdichtungen entsprechen den Schockfronten, und ihre Größenverteilung läßt eine Statistik, wie sie in der Burgers-Dynamik entsteht, plausibel erscheinen [25].

Abbildung 3.7. Geschwindigkeitsprofil u(x, t) nach (29) des Teilchenmodells ohne Kräfte zwischen den Teilchen, d.h. F ≡0 in (28), das Burgers’ Gleichung ohne Viskosität genügt [11].

Zusätzliche Terme, die die Symmetrie von Gleichung (8) brechen, zerstören oft die chaotischen Eigenschaften der Dynamik. Leicht einzusehen ist dies für einen linearen Dämpfugnsterm ṗ_{i} = −κp_{i}. Inwieweit Resultate für (8) verwandte Gleichungen in einem Teilchenmodell reproduziert werden können, wäre zu untersuchen [26, 27].

3.3. Die anisotrope KS–Gleichung

3.3.1. Einführung. Die Stabilisierung einer linear instabilen Oberfläche durch von nichtlinearen Effekten vermittelte Wechselwirkung wurde im Abschnitt 3.1.3 und in letzten Kapitel vorgestellt. In diesem Kapitel soll eine Veränderung von Gleichung (8) untersucht werden, die einen Übergang vom chaotischen zu global instabilem Verhalten vermittelt, bei dem zudem ein der Phasenordnung verwandtes Phänomen beobachtet wird (vgl. Abschnitt 2.2 im letzten Kapitel).

Eine naheliegende Verallgemeinerung von Gleichung (8) in 2+1 Dimensionen ist die Einführung einer Anisotropie

∂_{t}h = −[ν_{x}∂_{x}^{2} + ν_{y}∂_{y}^{2} + KΔ^{2}] h + \frac{1}{2}[λ_{x}(∂_{x}h)^{2} + λ_{y}(∂_{y}h)^{2}]. (36)

Diese Gleichung wurde zuerst hergleitet zur Beschreibung der Erosion von Oberflächen unter Ionenbeschuß [28]. Man erhält sie ebenfalls für die Modellierung von epitaktischem Wachstum auf vizinalen Oberflächen in der Entwicklung kleiner Gradienten der Fluktuationen auf der Oberfläche, was detailliert im Kapitel 6 dargestellt ist. In diesem Zusammenhang gehen anisotrope partielle Differentialgleichungen zurück auf [29] und [30]. Weitere Anisotropien in den linearen Termen, etwa anstelle des Terms KΔ^{2}h (siehe Kapitel 5) oder ein Term der Form ∂_{x}^{3}h [31] werden hier vernachlässigt.

Durch Reskalierung wie bei (26) wird in (36) nur je einer der ν- und der λ-Koeffizienten zu eins, der andere erhält den Wert, ν_{y}/ν_{x} bzw. λ_{y}/λ_{x}. Existiert mindestens eine linear instabile Richtung, z.B. ist ν_{x} > 0, und verschwindet der nichtlineare Term dieser Richtung nicht, also λ_{x} ≠0, erhält man durch entsprechende Reskalierung

∂_{t}h = −[∂_{x}^{2} + α∂_{y}^{2} + Δ^{2}] h + \frac{1}{2}[(∂_{x}h)^{2} + β(∂_{y}h)^{2}]. (37)

Anisotropie besteht, wenn immer α≠1 oder β ≠1, besonders ausgeprägt wenn einer der Koeffizienten negativ ist. In den nächsten Abschnitten wird dies sowohl am linearen (α) als auch am nichtlinearen Term (β) untersucht.

3.3.2. Anisotropie im linearen Term. Die Anisotropie im linearen Term beeinflußt zunächst die anfängliche Instabilität der glatten Oberfläche. Folgt man den Überlegungen zu (24) und (27) in Kapitel 2, so gehorchen kleine Fluktuationen um einen homogenen Zustand

∂_{t}h(q, t) = (q_{x}^{2} + αq_{y}^{2} − |q|^{4}) h(q, t) =: σ(q) h(q, t). (38)

Alle Moden mit σ(q)>0 – dieser Bereich ist in Abbildung 3.8 für mehrere α < 1 skizziert – wachsen exponentiell an. Bestimmend für die zuerst entstehende Struktur sind die q, bei denen σ(q) maximal wird. Bei α=1 sind dies alle q^{2} = q_{max}^{2} = 1/2, es entsteht ein Muster aus Zellen ohne Vorzugsrichtung (siehe Kapitel 2). Für α < 1 ist dies das Paar q_{max} = (±1, 0), für α>1 das Paar q_{max} = (0, ±\sqrt{α/2}). Im ersten Fall entstehen Rippen in y-, im zweiten in x-Richtung, jeweils mit Abstand 2π/|q_{max}|. Das Verhältnis ν_{x}/ν_{y} bestimmt also die Orientierung der durch die lineare Instabilität entstehenden Rippen. Bei Erosion durch Ionenbeschuß (auch Sputtern) kehrt sich diese Relation mit der Änderung des Einfallwinkels des Ionenstrahls um. Man beobachtet einen Wechsel in der Orientierung der Rippen [28, 32].

Abbildung 3.8. Parameterbereich von Gl. (37). Der instabile Bereich β < min(0, α) ist grau unterlegt. Das eingefügte Bild zeigt das Stabilitätskriterium in einem Quadranten der Ebene der Wellenvektoren. Durchgezogene Kurven begrenzen die linear instabilen Bereiche q_{x}^{2}+αq_{y}^{2}−|q|^{4} > 0 für verschiedene α. Die Strichlinie zeigt die “freien Moden” für ein β <0. Instabile freie Moden existieren, wenn sie den schattierten Bereich schneidet, hier für α=0, 0.2, und 1.

Die für α < 1 in y-Richtung verlaufenden Rippen können auch dann noch die Oberfläche bestimmen, wenn die obige lineare Näherung nicht mehr gültig ist. Betrachten wir als Beispiel den Fall starker Anisotropie im linearen Term mit einem isotropen nichtlinearen Term, also α → −∞ und β =1. Unter Reskalierung der y–Richtung um einen Faktor 1/\sqrt{α} mit unveränderter x-Richtung erhält man

∂_{t}h = −∂_{x}^{2}h − ∂_{x}^{4}h + \frac{1}{2}(∂_{x}h)^{2} + ∂_{y}^{2}h,

wenn man alle Terme der Ordnung α^{−1} vernachlässigt. Dies ist eine entlang der y-Richtung elastisch gekoppelte Kette von eindimensionalen KS-Systemen. Sind die Anfangsbedingungen an zwei verschiedenen Schnitten parallel zur x-Achse, h(x, y_{1}) und h(x, y_{2}), leicht verschieden, so werden sie von der chaotischen eindimensionalen Dynamik weiter auseinandergetrieben. Andererseits wirkt der Term ∂_{y}^{2}h entlang y glättend bis zu Abständen von der Größenordnung 1, oder in der originalen Längenskala von der Größenordnung \sqrt{α}. Zu späten Zeiten schnüren sich die Rippen auf der Oberfläche auf Abständen von etwa \sqrt{α} ab.

3.3.3. Anisotropie im nichtlinearen Term. Unter bestimmten Bedingungen entwickelt die Lösung von Gleichung (37) auch Rippen mit einer anderen nichttrivialen Orientierung. Ist nämlich β < 0, so können sich die beiden Teile des nichtlinearen Terms gegenseitig aufheben. Höhenprofile der Form h(x, y, t) = f(x−uy, t), die entlang der Linien x = s+uy konstant sind, ergeben in (37) eingesetzt

∂_{t}f = −(1 + αu^{2})∂_{s}^{2}f − (1 + u^{2})^{2}∂_{s}^{4}f + \frac{1}{2} (1 + βu^{2})(∂_{s}f)^{2}. (39)

Die Nichtlinearität veschwindet, wenn u = u_{c} = ±\sqrt{−1/β}. Das Oberflächeprofil kann exponentiell anwachsen, außer wenn 1+αu_{c}^{2} = 1−α/β < 0. Mit anderen Worten wächst solch ein Profil nur dann an, wenn die Moden verschwindender Nichtlinearität (im folgenden als freie Moden bezeichnet) auch einen Anteil im instabilen Bereich σ(q) > 0 haben. Der instabile Bereich ist gegeben durch β < min(0, α) und in Abbildung 3.8 graphisch dargestellt.

Die freie Mode mit der höchsten Wachstumsrate τ^{−1} = \frac{1}{4}(α−β)^{2}/(1−β)^{2} bildet ein Rippenmuster mit Wellenlänge ℓ = 2π\sqrt{2(1−β)/(α−β)}, die in einem Winkel ± arctan \sqrt{−β} in bezug zur x-Achse liegen.

Im generischen Fall, wenn Gleichung (37) zu einer zufälligen Anfangsbedingung kleiner Amplitude gelöst wird, sind verschiedene freie Moden angeregt. Numerische Simulationen zeigen folgendes Szenario: Existieren instabile Freie Moden (β < min(0, α)), so bilden sich aus einem chaotischen Hintergrund nach einiger Zeit Domänen mit parallel liegenden Rippen in der Orientierung der Freien Moden, ± arctan \sqrt{−β}. Sie dehnen sich annähernd linear in der Zeit aus, bis Bereiche der beiden Orientierungen das gesamte System bedecken. In Konkurrenz zueinander vergrößern sich die Bereiche, ähnlich dem Phasenordnungsprozeß von Gleichung (5), was in Abbildung 3.9 gezeigt ist. Bisher unternommene Versuche einer numerischen Auswertung haben folgende Gestalt: Man definiert einen “Ordnungsparameter”, der an jeder Stelle die Rippenorientierung angibt. Dazu legt man um jeden Gitterpunkt ein Fenster, etwa der Größe 16×16, und betrachtet die Fouriertransformierte des Höhenprofils h_{i,j}^{F} in diesem Fenster. Die relative Stärke des Spektrums in (1,1)–Richtung, s := (s_{/}−s_{\})/(s_{/}+s_{\}) mit s_{/} = ∑_{i=1}^{8}(h_{i,i}^{F})^{2} und s_{\} aus einer Summe über (i, −i) dient als Ordnungsparameter, der an jedem Punkt den Werten ±1 zustrebt. Allerdings ist es hier numerisch wesentlich aufwendiger, das Domänenwachstum quantitativ zu erfassen, so daß die genaue Art dieses Wachstums ungeklärt ist.

Abbildung 3.9. Entwicklung von Gleichung (37) mit α=1 und β =−1 in einem System der Größe 500×500 (Teilansicht). Auf dem chaotischen Hintergrund (a) nukleieren Domänen zusammenhängender Rippen (b), wachsen (c) und bedecken schließlich das gesamte System (d). An Maximum und Minimum der Oberfläche angepaßte Grauskalendarstellung.

3.3.4. Vorzeichenwechsel von β. Eine genaue Analyse der anisotropen Kuramoto-Sivashinsky–Gleichung könnte auch zur Klärung der Frage der effektiven stochastischen Theorie beitragen, deren widersprüchliche Beantwortung in Abschnitt 3.1.3 erwähnt wurde. Die anisotrope Verallgemeinerung der KPZ–Gleichung (7) in 2+1 Dimensionen zeigt einen Übergang beim Vorzeichenwechsel von β = λ_{y}/λ_{x} [30]. Für β < 0 ist die nichtlineare Kopplung irrelevant, und das Skalenverhalten der Oberfläche entspricht

3.3. DIE ANISOTROPE KS–GLEICHUNG 39

dem der linearen Gleichung (6) mit den Exponenten (19). Der logarithmischen Rauhigkeit (ζ = 0) bei entgegengesetztem Vorzeichen steht für β > 0 algebraische Rauhigkeit (ζ >0) wie in der isotropen KPZ–Gleichung gegenüber.

Eine global stabile Lösung mit anisotroper Nichtlinearität hat Gleichung (36) bzw. (37) bei α<β <0, im weißen Bereich links unten in Figur 3.8. Es wäre sinnvoll, diesen Parameterbereich mit Punkten oberhalb der Achse β = 0 zu vergleichen, um zu sehen, ob fundamentale Unterschiede bestehen, wie in der anisotropen Langevingleichung (7). Mikroskopische Details der Dynamik unterstützen eine solche Vermutung. β > 0 bewirkt, daß sich Fluktuationen auf der Oberfläche, etwa durch die lineare Instabilität entstehende Buckel, in beide Richtungen lateral ausdehnen (siehe Abschnitt 3.1.3). Ist β < 0, so breiten sie sich in x–Richtung zwar aus, ziehen sich in y–Richtung aber zusammen.

Im kollektiven Zusammenspiel der Oberflächenstrukturen ausgehend vom zunächst entstandenen Rippenmuster mit Defekten und Abschnürungen beobachtet man zwei Szenarien, die in Abbildung 3.10 dargestellt sind:

Abbildung 3.10. Vorzeichenwechsel von β. Teilansichten eines Systems der Größe 250 × 250, jeweils zu Zeiten t = 84, 88, 92, 94, mit β = −1 (links) und +1 (rechts). α = 50 in beiden Bilderfolgen. Deutliche Unterschiede im Ablauf der Dynamik fallen auf, links verschieben sich bevorzugt die Endpunkte der Rippen, rechts verbreitern sich die Rippen und bewirken so stärkere Aufrauhung.

40 3. DIE KURAMOTO-SIVASHINSKY–GLEICHUNG

Im einen Fall hat jede Zelle die Tendenz, sich in beide Richtungen zu verbreitern. Große Zellen wachsen auf Kosten ihrer Nachbarinnen und brechen auf, wenn sie eine kritische Größe überschritten haben. Die Dynamik ähnelt insofern der der isotropen Gleichung (8).

Im anderen Fall ziehen sich die Zellen in y–Richtung zusammen. Hat sich als Defekt im Rippenmuster eine Rippe zwischen zwei anderen gebildet, so zieht sich ihr Endpunkt aus der Spalte zwischen den beiden anderen zurück. Folgende Anschauung macht dies deutlich: Auf einem Schnitt parallel zur y–Achse durch den Endpunkt der Rippe ist die Höhe im wesentlichen konstant, nur auf der Flanke direkt am Endpunkt ist (∇h)^{2} ≈ (∂_{y}h)^{2} ≠ 0, weshalb dort die Höhe mit β(∂_{y}h)^{2} <0 abnimmt, sich das Ende der “eingeklemmten” Rippe also mit konstanter Geschwindigkeit v ≈ β|∂_{y}h| zurückzieht.

Beide Szenarien wirken grundsätzlich veschieden, und deuten darauf hin, daß auch bei der Kuramoto-Sivashinsky–Gleichung ein Vorzeichenwechsel im Quotienten der Koeffizienten des nichtlinearen Terms einen Einfluß auf die Rauhigkeit und das statistische Verhalten auf großen Skalen hat. Die genauere Untersuchung des Übergangs bei β = 0 wäre daher ein möglicher Ansatz zur Klärung der Frage anch dem effektiven stochastischen Modell zur Beschreibung von Gleichung (8) auf großen Skalen.

3.4. Zusammenfassung

In diesem Kapitel wurden Spezialfälle räumlich ausgedehnter Systeme mit komplexer Dynamik vorgestellt. Sie sind einerseits im Zusammenhang mit der Klassifizierung im Kapitel 2 zu sehen: Sie berühren die gemeinsamen Eigenschaften von chaotischen und stochastischen Modellen. Bestimmte Änderungen zerstören das chaotische Verhalten und prägen den Systemen eine Dynamik sich ordnender Strukturen auf (siehe Abschnitt 2.2).

Andererseits dient dieses Kapitel wegen der vorgestellten Methoden des Umgangs mit Oberflächendynamik als Vorbereitung auf den nun folgenden zweiten Teil der Arbeit. Die Modellierung mit partiellen Differentialgleichungen, die Auswertung mit statistischen Methoden werden hier zur theoretischen Beschreibung von Kristallwachstum unter Molekularstrahlepitaxie angewandt. Dieser physikalische Aspekt war auch – an verschiedenen Zwischenbemerkungen erkennbar, aber ohne daß näher darauf eingegangen wurde – Grundlage der zuvor behandelten Modelle.

KAPITEL 4 Molekularstrahlepitaxie

Molekularstrahlepitaxie (abgekürzt MBE vom englischen molecular beam epitaxy) bezeichnet die Kristallisierung durch Kondensation eines Dampfes im Hochvakuum. Seit ihrer Entwicklung um die Mitte der siebziger Jahre [33] wurde sie zu einer weit verbreiteten Technik, die es erlaubt Kristallwachstum präzise zu kontrollieren und damit z.B. komplizierte Halbleiterelemente zu entwickeln.

Bild 4.1 zeigt eine MBE–Anlage: In der Mitte einer Ultrahochvakuumkammer (UHV–Kammer) befindet sich ein drehbarer Träger, auf dem mehrere Substrate (wafer) angebracht werden können. Eine Heizung kontrolliert ihre Temperatur, und sie können durch den Drehmechanismus abwechselnd dem aus den Evaporationsquellen austretenden Molekularstrahl ausgesetzt werden.

In den Quellen befindet sich die aufzudampfende Materie in der Gasphase, etwa wie im Inneren eines Ofens. Öffnet man die Verschlüsse (shutter), tritt ein Strahl aus mit entsprechend der Temperatur in den Quellen (nach Maxwell) verteilten Geschwindigkeiten. Die Teilchen treffen ballistisch auf das Substrat auf, kollidieren und reagieren nicht während des Flugs. Teilchen, die nicht sofort auf das Substrat treffen, werden durch die Vakuumpumpen abgesaugt. Die Menge der austretenden Materie läßt sich genau durch die Öffnungszeit des Verschlusses steuern. So lassen sich Schichten von der Dicke eines Atoms, Multilagen mit atomar scharfen Schichtgrenzen herstellen, was bei Wachstum aus einem diffusen oder viskosen Fluß unmöglich wäre.

Eine Ionenpumpe und verschiedene Kryopumpen sorgen für ein Ultrahochvakuum von weniger als 5 × 10^{−11} Torr. Das Hintergrundgas besteht größtenteils aus Wasserstoff, einem zur Erhaltung der Reinheit von Oberflächen “gutartigen” Stoff. Die Partialdrücke von Stoffen, die zu Verunreinigungen im Kristall führen, etwa verschiedener Kohlenstoffverbindungen, liegen unterhalb der Nachweisgrenze eines zur Kontrolle in die Anlage eingebauten Quadrupolmassenspektrometers, also unter 10^{−14} Torr.

Abbildung 4.1. Eine von einem kommerziellen Anbieter (Riber) erhältliche MBE–Anlage [34]: (A) UHV–Kammer, (B) Verschluß (shutter) vor einer Evaporationsquelle (C), (D) drehbare Plattenstation für mehrere Wafer, (E)–(H) Vorrichtungen zur Präparation und zum Transport der Substrate, (I) und (L) Kryopumpen, (J) Ionenpumpen.

Mittels Elektronenstreuung (RHEED, reflection high energy electron diffraction) wird die Menge der aufgedampften Materie überwacht. Hauptsächlich an den Stufen werden Elektronen getreut, die Abschwächung des spekulären Strahls ist daher ein Maß für die Stufendichte auf der Oberfläche. Im Fall schichtweisen Wachstums beobachtet man typische Oszillationen im RHEED–Signal (siehe Bild 7.1) mit Maxima bei kompletter Füllung einer Schicht (kaum Stufen auf der Oberfläche) und Minima bei halber Füllung (maximaler Stufendichte). So lassen sich die aufgedampften Schichten während des Wachstums “abzählen”.

Die hohe Reinheit des wachsenden Kristalls und die genaue Kontrolle der Wachstumsbedingungen erlauben die Untersuchung intrinsischer Effekte des Kristallwachstums. Auftretende Instabilitäten sind inhärent und nicht auf Verunreinigungen und dadurch bedingte Störungen der Kristallstruktur zurückzuführen. Wachstumsexperimente unter MBE stellen daher die Grundlage und bieten die Überprüfungsmöglichkeit von Vorhersagen der in dieser Arbeit untersuchten Wachstumsmodelle.

4.1. HOMOEPITAXIE 43

4.1. Homoepitaxie

Anders als die genaue Kontrolle der Wachstumsstruktur in vertikaler Richtung beim schichtweisen Aufdampfen, stellt die Erzeugung lateraler Strukturen eine Herausforderung dar. Wachstumsinstabilitäten erzeugen lateral variierende Strukturen und bieten möglicherweise einen Lösungsweg. Eine solches Beispiel sind elastische Verspannungen bei Heteroepitaxie, die bereits nach wenigen aufgedampften Monolagen zur Bildung ausgeprägter Inseln mit einer scharfen Größenverteilung führt, wenn sich die Gitterkonstante des aufgedampften Materials von der des Substrats unterscheidet [35].

Aber auch unter Homoepitaxie, wenn das aufgedampfte gleich dem Substratmaterial ist, treten generisch Instabilitäten auf, wie im folgenden Abschnitt dargestellt wird. Der theoretischen Beschreibung ist Homoepitaxie leichter zugänglich, weshalb sich Modelle der Homoepitaxie zur Charakterisierung und zum grundlegenden Verständnis der Phänomene des instabilen Wachstums anbieten.

Bei Homoepitaxie bilden sich in der ersten auf einer atomar glatten Fläche wachsenden Lage Inseln, wenn diffundierende Adatome aufeinandertreffen und einen stabilen Keim für eine Insel bilden. Sie spielen eine wesentliche Rolle bei instabilem Wachstum, wie im nächsten Abschnitt dargestellt wird. Man beobachtet diese Keime, und damit auch die Inseln der ersten Lage in einem typischen Abstand l_{D} voneinander, der als Diffusionslänge bezeichnet wird. Dies ist die Entfernung, die ein einzelnes Adatom in seiner Diffusionsbewegung auf der Kristalloberfläche im Mittel zurücklegt, bevor es auf ein anderes trifft. Ist die gesamte Substratfläche mit Inseln bedeckt, die im Mittel den Abstand l_{D} voneinander haben, können alle Adatome zu den Rändern der bereits bestehenden Inseln gelangen ohne auf andere Einzelatome zu treffen, und keine weiteren Inseln nukleieren. Es entsteht bei einer relativen Füllung der ersten Lage von 0.2 bis 0.5 das typische Bild von etwa gleich großen Inseln mit einheitlichen Abständen zueinander (Abbildung 4.2).

Eine detaillierte Darstellung der Dichte der Nukleationspunkte von Inseln und ihrer Korrelationen über mehrere Kristallschichten findet sich in [36]. Eine einfache Abschätzung der Abhängigkeit der Diffusionslänge l_{D} von der Intensität des einfallenden Strahls und der Diffusionsrate D auf der Oberfläche sei hier kurz erwähnt [37, 38]. Dazu leitet man zunächst die Wahrscheinlichkeit P (l) her, daß auf einem Oberflächenstück S_{l} der Abbildung 4.2. Die Diffusionslänge l_{D} entspricht dem mittleren Abstand der auf einer singulären Oberfläche nukleierten Inseln. Links zwei Konfigurationen einer (kinetischen Monte-Carlo–) Simulation bei einer Bedeckung von 0.2 Monolagen, oben mit niedriger, unten mit hoher Diffusionsrate im Verhältnis zur Aufdampfrate [39]. Rechts eine STM–Aufnahme eines (4000Å)^{2} großen Stücks einer Pt(111)– Oberfläche. Niedrige Diffusionsrate bei niedriger Temperatur (425 K, oben) führt mit l_{D} = 540±4Å zu einer niedrigeren Diffusionslänge als unten (628 K, l_{D} = 1700±250Å) [40].

linearen Ausdehnung l (alle Längen in dieser Betrachtung sind in Einheiten der Gitterkonstanten a angegeben) pro Zeiteinheit eine Insel nukleiert. Der typische Abstand der Inseln l_{D} ist dann gegeben durch die Beziehung P (l_{D})/F = 1, den im entsprechenden Flächenstück findet während des Wachstums einer Monolage im Mittel eine Nukleation statt.

P (l) ist gegeben durch die Rate, mit der Atome auftreffen, F l^{d} mal der Wahrscheinlichkeit p, daß ein Atom ein anderes in S_{l} trifft. p ist ungefähr gleich der Anzahl der in S_{l} während einer Diffusionsbewegung besuchten Gitterplätze – bis auf logarithmische Korrekturen ld̂ mit d̂= min(2, d) – mal

4.2. STUFENBARRIEREN 45

der mittleren Adatomdichte, also p = ld̂⟨ρ⟩. Die Inseln, die die diffundierenden Adatome absorbieren, liegen im Abstand l_{d}, also ist die Adatomdichte zwischen den Inseln ⟨ρ⟩ = F l_{D}^{2}/D. Daraus folgt

P (l) = \frac{F^{2}}{D} l^{d+d̂} l_{D}^{2}

und schließlich

l_{D} ≃ (\frac{D}{F})^{\frac{1}{2+d+d̂}} . (40)

Der Abstand der Nukleationspunkte hängt also ab vom Verhältnis D/F. Je größer F, desto näher liegen die Nukleationspunkte beisammen.

4.2. Stufenbarrieren

Ehrlich und Hudda untersuchten als erste die Selbstdiffusion einzelner Wolframatome auf verschiedenen Oberflächen eines Wolframkristalls im Feldionenmikroskop [41, 42]. Am Emitter wurde ein Kristall mit atomar flachen Oberflächen niedriger Miller–Indizes präpariert, auf denen sich von einem geheizten Wolframdraht ausgetretene Atome anlagern. Das zunächst regelmäßige Muster der Feldionenmikroskopaufnahme weicht teilweise ungeordneten Punkten, die von den Adatomen auf den Facetten des Emitters herrühren (Bild 4.3).

Die Aufnahmen werden bei einer Temperatur gemacht (20 K), die weit unter der Aktivationsenergie für Adatomsprünge auf benachbarte Gitterplätze der Oberfläche liegt. Merkliche Diffusion berichten die Autoren erst bei 288 K auf der 211– und bei 337 K auf den 110– und 321–Oberflächen. Sie nutzen dies zur Beobachtung einzelner Sprünge der Adatomdiffusion: Zwischen den einzelnen Aufnahmen einer Serie wird jeweils die Probe für einige Minuten über die Aktivationsschwelle der Diffusion geheizt. Ist die Wartezeit nicht zu lang, lassen sich beim Vergleich aufeinanderfolgender Bilder die Bewegungen einzelner Atome identifizieren. In der Mitte der Facetten, weit genug entfernt von den Rändern, führen die Adatome Irrfahrten (random walks) aus, und man kann die Diffusionskonstante ihrer Bewegung bestimmen.

Ein wichtiges Ergebnis der Messung betrifft die Ränder der Terrassen: Die Autoren messen die Rate, mit der die Adatomdichte auf den Facetten abnimmt und vergleichen sie mit der Vorhersage, die sich aus der Lösung der Diffusionsgleichung mit absorbierenden Rändern ergibt. Absorbierende Stufen für Adatome sowohl auf der oberen als auch auf der unteren Terrasse Abbildung 4.3. Feldionenmikroskopische Aufnahme einer Wolframspitze. Auf den atomar glatten Facetten (oben) sind nach einiger Zeit von einem geheizten Draht ausgetretene Wolframatome aufgetroffen, die als zusätzliche helle Punkte im unteren Bild erscheinen [41].

wurden von Burton, Cabrera und Frank [43] zur Modellierung von Kristallwachstum angenommen. Die von Ehrlich und Hudda gemessene Rate liegt deutlich unter der Vorhersage. Die Stufen müssen also (zumindest teilweise) Adatome auf der oberen Terrasse reflektieren und somit dort “gefangen” halten, oder anders ausgedrückt: Abwärtssprünge über eine Stufekante benötigen eine höhere Aktivierungsenergie als Sprünge zwischen zwei Plätzen auf der Terrasse. Eine qualitative Erklärung liefert ein einfaches Kugelmodell des Kristalls: Rollt man eine (ein Adatom repräsentierende) Kugel über eine Stufenkante hinab, ist sie in einer Zwischenstellung von ihren Nachbarinnen im Kristall weiter entfernt als beim Übergang zwischen zwei

4.3. WACHSTUM MIT STUFENBARRIEREN 47

Plätzen auf der Terrasse. Die Konsequenzen dieses Effekts für Kristallwachstum liegen auf der Hand und wurden bereits von Ehrlich und Hudda erkannt:

… at low temperatures atoms do not roll over the edge of a plane – instead they are reflected. A growing step can therefore only acquire atoms from one side. [ …] it is in this regime that growth would be most affected.

Die charakteristischen Merkmale von Kristallwachstum mit Stufenbarrieren sind Gegenstand des nächsten Abschnitts.

4.3. Wachstum mit Stufenbarrieren

4.3.1. Oberflächen mit schwacher Neigung. Ehrlich und Hudda erkannten die Stufenbarrieren daran, daß sie diffundierende Atome auf einer Insel in der obersten Kristallage gefangen halten. Inseln können bei epitaktischem Wachstum entstehen, wenn die Terrassenstücke zu Beginn hinreichend Platz bieten also bei genügend kleiner mittlerer Neigung, |∇h| ≪ a/l_{D}. Andernfalls ist der ursprüngliche Stufenabstand kleiner als l_{D} und Adatome werden von den bereits vorhandenen Stufen absorbiert, bevor Inseln nukleieren können.

4.3.1.1. Pyramidenwachstum. Unter Wachstumsbedingungen, wenn ein ständiger Teilchenstrom auf die Oberfläche trifft, werden bevorzugt auf bereits bestehenden Inseln neue kleinere Inseln nukleieren. Dort ist die Adatomdichte größer als auf tieferen Terrassen, von denen aufwärtsgerichtete Stufen die Adatome absorbieren. Eine ursprünglich glatte Oberfläche ist also instabil gegen die Bildung von pyramidenartigen Strukturen aus mehreren übereinanderliegenden Inseln, oft als wedding cake instability bezeichnet.

Der Limes unendlich starker Stufenbarrieren, in dem keine Materie Stufen nach unten überquert, kann in relativ einfachen Modellen beschrieben werden. Hier ein Beispiel:

Die Stufenbarrieren halten die Materie in jeder Lage fest, unter idealen MBE–Bedingungen (ohne Evaporation) kann die Bedeckung der n–ten Lage θ_{n} also nicht abnehmen. Sie nimmt zu in dem Maße wie die darunterliegende Terrasse (Lage n − 1) dem einfallenden Teilchenstrahl – seine Intensität sei F – ausgesetzt ist,

\frac{∂}{∂t}θ_{n} = F (θ_{n−1} − θ_{n}). (41)

48 4. MOLEKULARSTRAHLEPITAXIE

Das Wachstum soll auf einer perfekten Kristalloberfläche beginnen, bei t=0 seien alle Lagen unterhalb der ersten gefüllt und oberhalb leer, also θ_{n} = 1 (0) für n≤0 (n>0).

Die Wahrscheinlichkeit, daß der Kristall über einem beliebigen Fußpunkt zur Höhe n (in Einheiten der Gitterkonstanten a) gewachsen ist, wird bezeichnet mit p_{n} = θ_{n} − θ_{n+1}. Mit obiger Anfangsbedingung, p_{n} = δ_{n 0} bei t=0 ist p_{n} zu späteren Zeiten poissonverteilt

p_{n}(t) = \frac{(F t)^{n}}{n!} e^{−F t n}. (42)

Die mittlere gewachsene Höhe nimmt linear mit der Zeit zu, wie zu erwarten ist bei konstantem äußerem Teilchenstrom:

h(t) = ∑_{n=0}^{∞} n p_{n}(t) = F t.

Die unendlichen Stufenbarrieren bewirken jedoch eine andere Entwicklung der mittleren Breite der Oberfläche, als wir sie aus Kapitel 2 kennen. Da die Höhen poissonverteilt sind, ist

w(t)^{2} = ∑_{n=0}^{∞} n^{2} p_{n}(t) − h(t)^{2} = h(t) = F t.

Die mittlere Breite der Oberfläche nimmt also wie die Wurzel der Zeit zu, w(t) ∼ t^{1/2}, und damit schneller als in den in Kapitel 2 vorgestellten Modellen. Die Stufenbarrieren beschleunigen die Aufrauhung.

Abbildung 4.4. Simulation des im Text beschriebenen 1+1dimensionalen Gittermodells mit unendlich hohen Stufenbarrieren. Die Breite der Säulen in der n-ten Kristallebene ist durch (42) gegeben. Mit freundlicher Genehmigung von M. Schimschak [44].

Dasselbe Aufrauhungsverhalten zeigt folgendes stark vereinfachte Wachstumsmodell: Der (d + 1–dimensionale) Kristall bestehe aus Würfeln, die auf Gitterplätze über der (d–dimensionalen) Substrat(hyper)ebene in

4.3. WACHSTUM MIT STUFENBARRIEREN 49

einem Poissonprozeß mit Rate F “hinabregnen” und dann dort liegenbleiben, egal ob sie auf eine Spitze oder in eine Mulde treffen. Jede Säule wächst unabhängig, die mittlere Breite der Oberfläche folgt aus der Varianz der zufälligen Bewegungen der Säulenhöhen, nimmt also auch zeitlich zu wie w^{2} ∼ h ∼ t.

Gleich den Säulen sind bei epitaktischem Wachstum mit vollständig reflektierenden Stufenbarrieren die in der ersten Lage entstandenen Inseln unabhängig voneinander, insofern als sie keine Materie mehr austauschen können. Die Aufrauhung der Oberfläche wird auch hier bestimmt von den Schwankungen im Wachstum der einzelnen Pyramiden, die sich auf den Inseln der ersten Lage gebildet haben.

Die Füllung der einzelnen Lagen ist im oben beschriebenen Modell (41,42) leicht zu erfassen, nicht so jedoch die Morphologie der entstehenden Pyramiden. In höheren Dimensionen, d+1 mit d>1 können die Stufen zwischen den Lagen vielfältige Formen einnehmen, die davon abhängen, wie Adatome zu den Stufen und an ihnen entlang diffundieren können. Insbesondere können also die Höhenlinien der Pyramiden oder die Ränder ihrer Kristallschichten fraktal oder kompakt sein. Einzig in d = 1 tragen die Stufen selbst keine Freiheitsgrade sondern sind punktförmig, bzw. die Pyramiden haben keine laterale Struktur. Die Form der Pyramiden hängt direkt mit dem Füllungszustand der Lagen zusammen [45]: Die Flanken der Pyramiden bestehen aus Terrassenstücken, deren Längen gerade die p_{n} aus Gleichung (42) sind. Zu späten Zeiten, wenn F t≫1, nähern sich die p_{n} einer Gaußverteilung an, was den Flanken die Form der zugehörigen Verteilungsfunktion, der Fehlerfunktion gibt. Diese Form besitzt eine Symmetrie, die Spitzen der Pyramiden haben dieselbe Form wie die sie trennenden Gräben, insbesondere werden beide unendlich dünn (Siehe Bild 4.4).

4.3.1.2. Lagenwachstum. Anders als in den zuletzt besprochenen Modellen sind im Experiment die Stufenbarrieren nicht unendlich hoch. So können singulär dünne, hohe Spitzen nicht entstehen. Tatsächlich ist es sogar möglich, die Oszillationen im RHEED–Signal, also das schichtweise Wachstum der Oberfläche (auf englisch layer by layer growth) über viele Lagen aufrechtzuerhalten. Ausgehend von einer zunächst glatten Oberfläche bedeutet dies, daß auf den in der ersten Lage darüber gebildeten Inseln (fast) keine erneuten Inseln nukleieren, bevor sich die erste Lage vollständig geschlossen hat.

Durch einen Trick in der Steuerung des einfallenden Teilchenstroms läßt sich Nukleation auf Inseln zusätzlich unterdrücken [46]: Auf die glatte Terrasse treffen die Teilchen mit hoher Intensität (großem F). Dadurch liegen

50 4. MOLEKULARSTRAHLEPITAXIE

die Inseln in der ersten Lage nahe beieinander (l_{D} ist klein). Nachdem sie überall auf dem Substrat nukleiert sind, drosselt man den Teilchenstrom (kleines F), so daß die Diffusionslänge dadurch deutlich größer wird als der Abstand und der Durchmesser der Inseln. Wenn sich die erste Lage in größeren Bereichen geschlossen hat erhöht man wiederum den Teilchenstrom und erhält in der zweiten Lage erneut Inseln mit geringem Abstand.

Auch ohne Stufenbarrieren sind Nukleationen in der zweiten Lage vor Vollendung der ersten möglich und treten unter generischen Bedingungen so häufig auf, daß sich die wachsende Oberfläche nach einer bestimmten Zeit aufrauht und über mehrere Kristallschichten erstreckt. Im RHEED– Signal verschwinden die Oszillationen und die Intensität des spekulären Strahls nimmt monoton ab. Wie lange schichtweises Wachstum bestehen kann, hängt ab von der Größe der Probe, ihrer Temperatur und damit der Diffusionsrate der Adatome, sowie von der Einfallrate der Teilchen [47], siehe auch Kapitel 7 dieser Arbeit.

4.3.2. Oberflächen mit starker Neigung. Auf vizinalen Oberflächen mit Neigung |∇h| ≫ a/l_{D} ist der mittlere Stufenabstand kleiner als l_{D}, und auf den Terrassenstücken finden keine Nukleationen statt. Die Stufenbarrieren bewirken, daß Adatome bevorzugt von unten zu den Stufen diffundieren. Indem sie Adatome absorbieren, bewegen sich die Stufen gleichmäßig in Richtung der tieferliegenden Terrasse. Ist es möglich, daß bei diesem Wachstum mit fließenden Stufen, step flow growth, eine Konfiguration mit geraden Stufen gleichen Abstands stabil ist?

Die Antwort ist unterschiedlich für Fluktuationen der Stufenabstände und Fluktuationen ihrer Form. Vernachlässigen wir zunächst Fluktuationen der Stufenform, wie etwa Mäanderungen, und betrachten nur die mittlere Position der Stufen auf einer Koordinatenachse in Richtung der Oberflächenneigung. Je größer die unter einer Stufe liegende Terrasse ist, umso mehr Materie wird sich an ihr anlagern und umso schneller wird sie sich beim Wachsen vorwärtsbewegen. Bleibt durch eine Fluktuation eine Stufe zurück, vergrößert sich ihr Einzugsgebiet und sie holt ihren Rückstand auf. Umgekehrt fällt eine weiter vorgedrungene Stufe zurück. Diese Stabilisierung der Abstände zwischen den Stufen wurde zuerst von Schwoebel und Shipsey [48] und später in einem größeren Zusammenhang von Villain [29] beschrieben. Es ist jedoch leicht einzusehen, daß eine Folge von Stufen gleicher Abstände nur eine metastabile Konfiguration ist: Entsteht durch eine Fluktuation ein Terrassenstück der Größe l_{D}, so wird dort eine Insel nukleieren. Die Insel zieht Adatome an und bremst so die Bewegung

4.3. WACHSTUM MIT STUFENBARRIEREN 51

der dahinterliegenden Stufe. Dieser Mechanismus führt damit letztlich zur Bildung von Pyramiden wie auf einer Oberfläche schwacher Neigung [44].

Gegenüber transversalen Fluktuationen ist eine Folge gerader, paralleler Stufen instabil [49]. An jede Stufe lagern sich Adatome aus einer diffusiven Bewegung heraus an. Auf ein durch eine Fluktuation vorgeschobenes Stufenstück, eine “Halbinsel”, werden daher mehr Teilchen treffen als auf ein gerades Stück, auf ein zurückgesetztes Stufenstück, eine “Bucht”, entsprechend weniger. Diese Instabilität hängt wesentlich von den Stufenbarrieren ab, da sich bei beidseitiger Anlagerung der Überschuß an Adatomen von der einen Seite und der Mangel von der Anderen ausgleichen würden. Diffusion entlang der Stufen glättet diese Fluktuationen aus, allerdings nur unterhalb einer bestimmten Wellenlänge. Ausgleichender Massentransport entlang der Stufe ist über zu große Entfernungen nicht wirksam. Aus den zunächst geraden Stufen wird so eine Folge von gewellten Stufen, das Oberflächenprofil gleicht einem Wellblechdach mit Rinnen in Richtung der Neigung, wie in Abbildung 1.1 zu sehen ist.

Fassen wir noch einmal die hier beschriebenen Phänomene homoepitaktischen Wachstums mit Stufenbarrieren zusammen: Die “stabilen” Wachstumsformen, also Lagenwachstum auf singulären und gleichmäßig fließende Stufen auf vizinalen Oberflächen sind nur transiente Erscheinungen. Auf hinreichend großen Terrassenstücken und insbesondere auf singulären Oberflächen entstehen Pyramiden, parallele Stufen einer ebenen geneigten Oberfläche verformen sich zu Wellen. Beides ist in Bild 1.1 zu sehen.

Der Rest dieser Arbeit behandelt nun die weitere Entwicklung der Wachstumsinstabilitäten für lange Zeiten. Man beobachtet eine Vergröberung der Oberflächenstrukturen zu weit größeren Skalen als l_{D}. Aus den Pyramiden werden “Berge”, für die wir die Bezeichnung Hügel gebrauchen werden. Zu ihrer Beschreibung sind kontinuierliche Modelle angemessen, also partielle Differentialgleichungen, die im nächsten Kapitel eingeführt werden. Die Stufenbarrieren geben der Diffusion der Adatome eine Vorzugsrichtung, was in diesem Rahmen einem Nichtgleichgewichtsstrom auf der Oberfläche entspricht. Auf der Grundlage des Kontinuumsansatzes werden in Kapitel 6 quantitative Vorhersagen über die Vergröberung der Strukturen, sowie in Kapitel 7 eine theoretische Vorhersage des Übergangs von schichtweisem zu rauhem Wachstum gemacht.

52 4. MOLEKULARSTRAHLEPITAXIE

KAPITEL 5 Kontinuumsmodelle für Oberflächen

Kontinuumsmodelle als Lösungsansatz für Vielteilchenprobleme haben eine lange und erfolgreiche Geschichte. Ein herausragendes Beispiel ist die Hydrodynamik: Bei der Herleitung hydrodynamischer Gleichungen, z.B. der Navier–Stokes–Gleichungen, werden einige fundamentale Annahmen gemacht, die auch zur Beschreibung von Oberflächen sinnvoll sind und hier einführend dargestellt werden sollen.

Jedem Kontinuumsmodell zugrunde liegt das Konzept der Skalentrennung. In die Kontinuumsgleichungen gehen gemittelte Größen ein, wie etwa Dichte, Temperatur, Druck, Strom, oder anderes. Gemittelt werden sie aber nur über einen mesoskopischen Bereich, über Volumen- oder Flächenelemente der linearen Ausdehnung l_{m}, die viel kleiner als das gesamte (makroskopische) System, aber viel größer als die mikroskopischen Volumina oder Flächen sind, die z.B. von der Größenordnung der mittleren freien Weglänge eines Fluids oder der Gitterkonstante eines Kristalls sind.

Ferner ist wichtig, daß die gemittelten Größen eine vollständige Beschreibung der Vorgänge auf mesoskopischer Skala liefern, wie etwa die Zustandsgrößen in der Thermodynamik, die alle “wesentlichen Freiheitsgrade” eines Systems erfassen. Dazu dient die Hypothese des lokalen Gleichgewichts, nach der die gemittelten Größen lokal durch thermodynamische Beziehungen verknüpft sind. In Abschnitt 5.1 werden wir mit dieser Voraussetzung Parameter herleiten, die in die Kontinuumsgleichung eingehen. Das makroskopische System als gesamtes ist hingegen nicht im Gleichgewicht. Die mesoskopisch gemittelten Größen variieren daher räumlich und zeitlich, z.B. die Dichte ρ ≡ ρ(x, t). Das makroskopische System entwickelt sich zeitlich, weil die mesoskopischen Volumenelemente aneinander gekoppelt sind. Dem Konzept der räumlichen Skalentrennung entspricht die Annahme der Existenz einer mesoskopischen Zeitskala τ_{m}, innerhalb derer sich die Einflüsse durch Kopplung nur auf benachbarte Volumenelemente erstrecken. Die Dynamik erhält damit die Form von partiellen Differentialgleichungen, in denen l_{m} als infinitesimale Größe dx und τ_{m} als dt aufgefaßt werden.

54 5. KONTINUUMSMODELLE FÜR OBERFLÄCHEN

Für Kristalloberflächen, den Hauptgegenstand dieser Arbeit, verfährt man bei der mesoskopischen Mittelung wie folgt: Nehmen wir an, der Kristall sei mit einer seiner Hauptsymmetrieachsen entlang der z-Achse eines kartesischen Koordinatensystems orientiert. Die übrigen 2 (im allgemeinen d) Koordinatenachsen spannen eine Hyperebene auf, die als Substrat bezeichnet wird. Vernachlässigt man Leer- und Fehlstellen im Kristall, so ist die Oberfläche zur Zeit t durch eine Höhe h_{x}(t) über jedem Gitterpunkt x in der Substratebene S ⊂ ℝ^{d} gegeben. Glättet man die diskrete Höhenkonfiguration durch einen Faltungskern mit typischer Breite l_{m}, so ergibt sich daraus ein gemitteltes Höhenprofil

h(x, t) := ∑_{x′∈S} χ^{l_{m}}(x − x′)h_{x′}(t). (43)

Als Argument von h(x, t) wird x ∈ ℝ^{d} nunmehr als kontinuierliche Variable aufgefaßt. Das Höhenfeld h ist variiert nur auf Skalen, die größer sind als l_{m}.

Die diskrete Darstellung h_{x} unter Vernachlässigung von Fehlstellen und Überhängen spielt ebenfalls eine wichtige Rolle bei der Modellierung von Kristallstrukturen an Oberflächen. Man stellt sich dabei vor, daß die Einheitszellen des Kristalls ohne Leerräume “aufeinandergestapelt” sind, daher die Bezeichnung Solid-on-solid–Modell oder SOS–Modell. Zu den hier einzuführenden Kontinuumsmodellen verhalten sich SOS–Modelle wie ein Ising–Modell zur Ginzburg-Landau–ϕ^{4}–Theorie: Trotz unterschiedlicher mikroskopischer Natur erwartet und findet man bei geeigneter Wahl der Parameter dasselbe universelle Verhalten auf großen Skalen.

5.1. Oberflächen im lokalen Gleichgewicht

Lokales Gleichgewicht bedeutet hier, daß Oberflächenstücke auf der oben definierten mesoskopischen Längenskala l_{m} ^{1} im Gleichgewicht sind. Unter dieser Voraussetzung definiert man die Oberflächenspannung in ihrer Abhängigkeit von der lokalen Neigung der Oberfläche wie folgt [50]: Für ein einfach kubisches Kristallgitter (andere Gitterformen lassen sich analog konstruieren) betrachtet man ein SOS–Modell, dessen Höhenkonfigurationen gegeben sind durch h_{x} ∈ ℤ über dem Substrat S_{m} := [−l_{m}, l_{m}]^{d} ∋ x. Die mittlere Neigung u wird über helikale Randbedingungen vorgegeben: Die Höhen über zwei einander gegenüberliegenden Fußpunkten y := (y_{1}, …, +l_{m}, …, y_{d}) und y′ := (y_{1}, …, −l_{m}, …, y_{d}) sind verknüpft durch die Beziehung h_{y} − h_{y′} ≡ 2l_{m} u · e_{i}, wobei e_{i} ∈ ℝ^{d} der Einheitsvektor der

1Im folgenden wird l_{m} als ganze Zahl in Einheiten der Gitterkonstanten a aufgefaßt.

5.1. OBERFLÄCHEN IM LOKALEN GLEICHGEWICHT 55

Komponenten mit y_{i} = l_{m} und y′_{i} = −l_{m} ist. Diese Definition verlangt, daß in allen Komponenten 2l_{m}u_{i} eine ganze Zahl ist. Das Konzept der Skalentrennung beinhaltet jedoch l_{m} ≫ 1, und die Wahl von u ist daher praktisch nicht eingeschränkt. Helikale Randbedingungen für die Höhenkonfiguration bedeuten periodische Randbedingungen in den lokalen Höhendifferenzen.

Der Hamiltonian H führt über die zugehörige Gibbsverteilung zu einer Zustandssumme über die Konfigurationen dieses endlichen Systems Z_{m}(u) := ∑_{{h}} exp −βH[h], und damit zur Freien Energie F_{m}(u) := −β^{−1} log Z_{m}(u). Dividiert durch den Flächeninhalt der geglätteten Oberfläche h_{g} := u · y, erhält man die Oberfächenspannung

γ(u) := \frac{1}{\sqrt{1 + u^{2}}} \frac{1}{(2l_{m}−1)^{d}} F_{m}(u), mit l_{m} → ∞ (44)

Die Darstellung der Oberfläche als Höhenkonfiguration im SOS- oder Kontinuumsmodell bietet die Definition der Oberflächenspannung als Freie Energie pro Substratsfläche an, in der Literatur oft als γ̂ bezeichnet. Wir werden im Folgenden jedoch Definition (44) verwenden.

In der Dynamik der Oberflächenmodelle unterscheidet man zwei verschiedene Fällen. Im einen ist die gesamte Masse des Systems M_{m}[h] := ∑_{x∈S_{m}} h_{x} erhalten, im anderen kann sie frei fluktuieren, wird aber durch die Randbedingungen im Mittelwert festgehalten. Im ersten Fall wird Materie entlang der Oberfläche durch Oberflächendiffusion transportiert, im zweiten durch Evaporation und Kondensation im Austausch zwischen Oberfläche und dem darüberliegenden Raum.

Die Lineare Antwort des Systems auf eine infinitesimal kleine treibende Kraft definiert die zur Kraft gehörige Mobilität des Systems. Um diese dynamische Größe zu fassen, führt man zunächst Übergangsraten w_{h→h′} zwischen den Konfigurationen des Systems ein. Sie stellen ein detailliertes Gleichgewicht zwischen den Konfigurationen h bezüglich der Gibbsverteilung her,

w_{h→h′} exp −βH[h] ≡ w_{h′→h} exp −βH[h′]. (45)

Ihre genaue Form, etwa als Metropolisraten, und die Auswahl der erlaubten Übergänge, etwa nur solche, bei denen ein einziges Teilchen bewegt wird, wäre noch zu definieren, soll aber in die allgemeine Betrachtung an dieser Stelle nicht eingehen. Die Bedingung des detaillierten Gleichgewichts führt zu einer natürlichen Konstruktion der Übergangsraten mit treibender Kraft.

Eine treibende Kraft unter Evaporations-Kondensations-Dynamik ist eine Erniedrigung oder Erhöhung des chemischen Potentials für Atome im

56 5. KONTINUUMSMODELLE FÜR OBERFLÄCHEN

Kristall δμ_{0} im Vergleich zum chemischen Potential in der Dampfphase über dem Festkörper. Betrachten wir den Fall der Erniedrigung, δμ_{0} < 0, die den Kristall wachsen läßt: Übergänge, bei denen die Gesamtmasse M_{m} erhöht wird, erhalten eine höhere, solche mit Verminderung der Masse eine niedrigere Rate. Dies geschieht derart, daß die Raten das detaillierte Gleichgewicht der Energieänderungen berücksichtigen, nämlich

w_{h→h′}^{δμ_{0}} = w_{h→h′} exp β(M_{m}[h]−M_{m}[h′])δμ_{0}. (46)

Bei Oberflächendiffusion interessiert eine treibende Kraft, die den Sprüngen der Adatome eine Vorzugsrichtung aufprägt. Experimentell intensiv untersucht [51] ist das Phänomen der Elektromigration: Durch den Kristall fließender elektrischer Strom bewirkt über Streuung der Elektronen an den Adatomen eine Vorzugsrichtung in deren den Diffusionssprüngen. Phänomenologisch ordnet man den Adatomen eine effektive Ladung zu, entsprechend ihrer Reaktion auf das den Strom erzeugende elektrische Feld. Diese Vorstellung entspricht auch der theoretischen Modellierung des Problems [52], bei der der Sprung eines Adatoms um δx neben der Änderung der inneren Wechselwirkungsenergie eine Energieänderung um −q δx · E bewirkt. Bezeichnen h (h′) die Konfigurationen vor (nach) dem Sprung, so ändern sich analog zu (46) die Raten

w_{h→h′}^{E} = w_{h→h′} exp βq δx · E. (47)

Um die lineare Antwort auf kleine Felder zu fassen, genügt die einfachste Wahl, nämlich die Projektionen von δx und E auf die Substratebene.

Ändert man zum Zeitpunkt t = 0 die Raten w_{h→h′} gemäß (46), so nimmt zur Zeit t>0 die mittlere Höhe mit der Geschwindigkeit v(δμ_{0}, t) := (d/dt)⟨M_{m}⟩/(2l_{m} +1)^{d} zu. Analog findet man für die erhaltene Dynamik einen mittleren Strom j(E, t) entlang der Oberfläche. Die jeweiligen linearen Anteile definieren die

• Adatommobilität

σ := ∇_{E} lim_{t→∞} j(E, t)|_{E=0}, (48)

eine wegen der Anisotropie einer vizinalen Oberfläche im allgemeinen tensorielle Größe. Auf Oberflächen mit vier- und sechszähliger Symmetrie (in d = 2 Dimensionen) wird σ skalar. Zur nichterhaltenen Dynamik gehört die

• Oberflächenmobilität

σ_{⊥} := lim_{δμ_{0}→0} lim_{t→∞} \frac{v(δμ_{0}, t)}{δμ_{0}}. (49)

5.3. CHEMISCHES POTENTIAL FÜR ADATOME 57

Die die Dynamik auf mesoskopischer Skala bestimmenden Größen, Oberflächenspannung und Mobilitäten, lassen sich also auf natürliche Weise aus den zugrundeliegenden mikroskopischen Modellen definieren, allerdings nur für Oberflächen im lokalen Gleichgewicht.

5.2. Bemerkungen zum Nichtgleichgewicht

Kristallwachstum unter für Molekularstrahlepitaxie typischen Bedingungen findet fern vom Gleichgewicht statt. Zu den beweglichen Adatomen kommen ständig neue aus dem äußeren Teilchenstrom hinzu, was ausgeglichen wird durch “Verlust” bei Nukleation von Inseln, wenn zwei oder mehr diffundierende Teilchen zusammentreffen, oder bei Anlagerung an bereits vorhandene Kristallstufen. Fern vom Gleichgewicht ist das System, weil dieser Durchsatz mit einer so großen Rate erfolgt, daß der Zustand der Adatome und der Oberfläche bei weitem nicht der Gleichgewichtsverteilung entspricht, wie sie im letzten Abschnitt definiert wurde. Dies äußert sich auch quantitativ: Die mittlere Adatomdichte n_{0} ist viel gößer als der entsprechende Gleichgewichtswert n_{e}q.

Ein thermodynamischer Formalismus für Systeme außerhalb oder gar fern vom Gleichgewicht ist nicht annähernd so weit entwickelt wie für den Gleichgewichtsfall. Ein Kontinuumsmodell, wie es sich für letzteren mikroskopisch herleiten ließ, kann in diesem Fall nur ein phänomenologischer Ansatz sein. Seine Berechtigung erhält es durch Vergleich mit Experimenten oder Simulationen von elementaren Modellen des wachsenden Kristalls.

Bereits ein stationärer Zustand, das Analogon zum Gleichgewichtszustand, ist beim wachsenden Kristall nicht notwendig gegeben. Ein Beispiel dafür, das schichtweise Wachstum (layer by layer growth), werden wir noch eingehend untersuchen. Direkt übertragbar wäre die Herleitung der Nichtgleichgewichts-Oberflächendynamik mit verallgemeinerter Oberflächenspannung oder Freier Energie des Systems, wenn ein stationärer Zustand eine Gibbsverteilung im Konfigurationsraum der Oberfläche wäre. Dies ist keinesfalls gewährleistet, auch wenn Gibbsverteilungen in getriebenen Systemen mit stationärem Zustand eine besondere Rolle spielen [53].

5.3. Chemisches Potential für Adatome

Die mesoskopische Definition der neigungsabhängigen Oberflächenspannung führt zur funktionalen Abhängigkeit der Freien Energie der Oberfläche von einer gegebenen Oberflächenkonfiguration

F\{h\} = ∫_{V} d^{d}x γ(∇h(x, t))\sqrt{1 + ∇h(x, t)^{2}}, (50)

indem man über die Beiträge der einzelnen Oberflächenelemente integriert. Das chemische Potential ist die Energieänderung des Systems unter Zugabe eines Teilchens. Im Kontinuumsbild erscheint dies wie folgt: Variiert man die Höhenfunktion h(x, t) → h(x, t)+ε(x, t), so wird über dem Fußpunkt x die Menge ρ ε(x, t)d^{d}x an Materie zugefügt, wobei die inverse Dichte ρ^{−1} das atomare Volumen des Stoffes bezeichnet. Die Zugabe eines einzelnen Atoms entspricht einer infinitesimalen Änderung, so daß das chemische Potential durch die Funktionalableitung

μ(x, t) = (ρd^{d}x)^{−1} \frac{δF\{h\}}{δh}(x, t) (51)

gegeben ist. Mit der Gestalt von Gleichung (50) folgt

\frac{δF\{h\}}{δh} = −∇ · [\frac{γ(∇h)}{\sqrt{1 + (∇h)^{2}}} ∇h + \sqrt{1 + (∇h)^{2}} ∇_{u}γ|_{u=∇h}] , (52)

wobei ∇_{u}γ nur die Ableitung nach dem Argument von γ (an dieser Stelle als differenzierbar angenommen) bezeichnet. Oft betrachtet man Oberflächen mit kleinen lokalen Steigungen (∇h)^{2} ≪1, und nähert daher linear in ∇h

\frac{δF\{h\}}{δh} ≈ −∂_{i} [δ_{ij}γ(0) + ∂_{ij}^{2}γ(0)] ∂_{j}h. (53)

Bei Oberflächen mit vier- oder sechszähliger Symmetrie ist die Hessesche der Oberflächenspannung skalar, ∂_{ij}^{2}γ(0) = δ_{ij}γ′′(0), und (53) vereinfacht sich zu

\frac{δF\{h\}}{δh} ≈ −γ̃ Δh, (54)

mit der Definition γ̃ ≡ γ(0)+γ′′(0) für die Steifigkeit (englisch stiffness) der Oberfläche. Die letzte Gleichung macht deutlich, daß diese Herleitung auf eine Gibbs-Thomson-Formel für Kristalloberflächen führt: An einer Grenzfläche hängt das lokale chemische Potential von ihrer Kümmung ab.

Der eigenen Erwähnung verdient der Fall, daß γ für bestimmte Orientierungen eine Singularität, etwa eine “Ecke” besitzt. In diesen Orientierungen divergiert ∂_{ij}^{2}γ und daher auch die Steifigkeit in (54). Solch ein Fall tritt auf bei Temperaturen unterhalb des Rauhigkeitsübergangs in diesen Oberflächenorientierungen: Der Singularität in γ entspricht nach der Wulff– Konstruktion eine makroskopische Facette in der Gleichgewichtsform des Kristalls [54, 50]. In der Sprache der Gleichungen (52) bis (54) bedeutet dies, daß auf einer Facette ein “unendlich hohes” chemisches Potential für

5.4. DYNAMIK DER OBERFLÄCHE 59

Adatome herrscht. Tatsächlich wird sich ein auf die Facette aufgebrachtes Adatom an eine Stufe an ihrem Rand anlagern. Bei Wachstum des Kristalls bestehen jedoch keine Rauhigkeitsübergänge [55, 56] (siehe auch Kapitel 7). Das phänomenologische Äquivalent zu γ fern vom Gleichgewicht ist also glatt und die Steifigkeit endlich.

5.4. Dynamik der Oberfläche

Die Dynamik wird getrieben durch das Bestreben, die Freie Energie der Oberfläche zu verringern. Dem “Weg” von h im Raum der Konfigurationen hin zu niedrigeren Freien Energien entspricht der Transport von Materie entlang dem negativen Gradienten des chemischen Potentials, bzw. nach (54) der mittleren Krümmung der Oberfläche. Dabei können unterschiedliche Relaxationsmechanismen, d.h. unterschiedliche Wege des Materietransports, auf unterschiedlichen Zeitskalen ablaufen. In der Relaxationsdynamik sind dann bestimmte Freiheitsgrade “eingefroren”, was verdeutlicht wird beim Vergleich von Evaporations–Kondensations–Dynamik mit Oberflächendiffusionsdynamik.

5.4.1. Oberflächendiffusion. Ein Beispiel eingeschränkter Relaxationsdynamik haben wir bereits kennengelernt: Unter für Molekularstrahlepitaxie typischen Bedingungen bleibt jedes Teilchen auf der Oberfläche haften und wird in den wachsenden Kristall eingebaut. Zunächst wenden wir uns einer vergleichbaren Situation unter lokalem Gleichgewicht zu. Die Unterdrückung der Evaporation hat eine wichtige Konsequenz für die Relaxation. Ist in der Gasphase über der Oberfläche das chemische Potential μ_{0} niedriger als μ(x, t) an einer Stelle auf der Oberfläche, so kann dennoch nicht Materie von diesem Ort sublimieren, und umgekehrt kann der umgebende Raum nicht als Teilchenreservoir zur Rekondensation dienen. Der entsprechende Freiheitsgrad ist unter diesen Bedingungen der Relaxation nicht zugänglich. Dynamik mit erlaubter Evaporation und Rekondensation wird im nächsten Abschnitt untersucht.

Vernachlässigen wir ferner Leerstellen im Innern des Kristalls, entspricht dem Massenerhalt ein Volumenerhalt:

∫_{S} d^{d}x h(x, t) = ρ^{−1}M = const. (55)

Die Dynamik wird daher durch eine Kontinuitätsgleichung beschrieben. In Normalenrichtung verschiebt sich die Oberfläche entsprechend der lokal sich akkumulierenden oder entfernenden Masse, also um die Divergenz eines (tangentialen) Oberflächenstroms. In lokalen Koordinaten ist dies

ρv_{n} = n̂ \frac{1}{\sqrt{|g|}} ∂_{i}(\sqrt{|g|}J^{i}) (56)

Der Strom selbst ist die lineare Antwort auf einen Gradienten im chemischen Potential auf der Oberfläche,

J^{i} = σ^{i}_{j} g^{jk} ∂_{k}μ. (57)

Stellt man die Oberfläche als Graphen einer Funktion h(x, t) dar, so ist der metrische Tensor gegeben durch g^{ij} = δ^{ij} + ∂^{i}h ∂^{j}h und |g| = 1 + ∇h^{2}. Die Mobilität σ^{i}_{j} ist definiert in (48). Die Wachstumsgeschwindigkeit in h–Richtung hängt mit der Normalengeschwindigkeit zusammen über ∂_{t}h = \sqrt{1+(∇h)^{2}} |v_{n}|. Zusammen mit (51) und (52) erhält man daraus eine partielle Differentialgleichung für h.

(58) ∂_{t}h = −∇ · [\sqrt{1 + ∇h^{2}} σ(∇h) g(∇h)

∇ ∇ · [\frac{γ(∇h)}{\sqrt{1 + (∇h)^{2}}} ∇h + \sqrt{1 + (∇h)^{2}} ∇_{u}γ|_{u=∇h}]]

Sie läßt sich einfach verstehen in einer linearen Näherung für kleine Neigungen (∇h)^{2} ≪ 1, die auf Mullins zurückgeht [57]. Die zweite Zeile in (58) wird angenähert nach (54). Für eine lineare Gleichung vernachlässigen wir in den Vorfaktoren (in der ersten Zeile) alle Terme von O(|∇h|) und erhalten mit σ(∇h) ≃ σ(0) = σ𝕀, \sqrt{1 + (∇h)^{2}} ≃ 1 und g(∇h) ≃ g(0) = 𝕀 die Gleichung

∂_{t}h = −K Δ^{2}h (59)

mit K := σγ̃.

Weisen γ und σ eine Singularität bei ∇h=0 auf, also für Kristalloberflächen bei Temperaturen unterhalb des Rauhigkeitsübergangs, so kann die Dynamik nicht in der Form (58) dargestellt werden. Insbesondere würde K in (59) divergieren, oder – wenn σ(∇h=0) = 0 – könnte K formal endlich sein, Gleichung (59) aber dennoch nicht die korrekte Beschreibung der Dynamik liefern. Spohn [58] modelliert die Dynamik in 1 + 1 Dimensionen, indem er umgekehrt u = ∂_{x}h als Funktion von γ̃ auffaßt und eine Zeitentwicklung von ∂_{t}γ̃ herleitet. Nicht für jede Anfangsbedingung beschreibt diese Methode die Relaxation korrekt [59], insbesondere ist die Wahl von σ(∇h = 0) = 0 nicht durch die Herleitung festgelegt. Die Relaxation in facettierte Gleichgewichtszustände ist in der Kontinuumsbeschreibung ein offenes Problem.

5.4. DYNAMIK DER OBERFLÄCHE 61

Wie anschaulich zu erwarten dämpft Oberflächendiffusion kurzskalige Variationen besonders schnell aus. Nach Gleichung (59) nimmt ein Profil mit typischer Wellenslänge k wie exp(−Ktk^{4}) ab. Die Dämpfung langwelliger Variationen des Profils wird dadurch behindert, daß Masse über entsprechend große Entfernungen transportiert werden muß. Von der vierfachen Ableitung Δ^{2} in (59) rühren zwei vom chemischen Potential μ her (siehe (54)), und zwei daher, daß die Höhenänderung die Divergenz eines vom Gradienten von μ getriebenen Stromes ist. Von Relaxation mit Evaporation und Kondensation erwarten wir daher in der zu (59) analogen Gleichung zwei Ableitungen weniger, was Gegenstand des nächsten Abschnitts ist.

5.4.2. Evaporation und Kondensation. Der wesentliche Relaxationsmechanismus läßt sich am einfachsten an einer Oberfläche darstellen, die im Mittel weder wächst noch sublimiert. Das chemische Potential auf einer ebenen Oberfläche μ_{0} ist also gleich dem im umgebenden Raum. Einzige treibende Kraft der Relaxation sind also die Variationen im chemischen Potential gemäß (52), die Materie evaporieren ode kondensieren lassen. Die Oberfläche wird also in ihre Normalenrichtung verschoben mit einer Geschwindigkeit v_{n} = −σ(μ − μ_{0}). Die zugehörige Mobilität σ ist in (49) definiert. Die Dynamik des Höhenprofils h(x, t) folgt damit der partiellen Differentialgleichung

(60) ∂_{t}h = −\sqrt{1 + ∇h^{2}} σ(∇h) ×

∇ · [\frac{γ(∇h)}{\sqrt{1 + (∇h)^{2}}} ∇h + \sqrt{1 + (∇h)^{2}} ∇_{u}γ|_{u=∇h}] .

Wie nach den Überlegungen des vorangegangenen Unterabschnitts zu erwarten besitzt die analog zu (59) linearisierte Form von (60)

∂_{t}h = ν Δh (61)

nur zwei Ortsableitungen mit dem Vorfaktor ν:=σγ̃. Die letzte Gleichung beruht auf demselben Prinzip wie die Grenzflächendynamik bei spinodaler Entmischung ohne Erhaltung der mittleren Magnetisierung (Modell A [4]). Die Grenzflächen verschieben sich lokal mit einer Geschwindigkeit proportional zum Vektor der mittleren Krümmung [3]. Massenerhaltung auf der Oberfläche führt zu einer Modell B [4] verwandten Relaxation. Die Verschiebung der Fläche ist bestimmt durch die Divergenz des Gradienten

62 5. KONTINUUMSMODELLE FÜR OBERFLÄCHEN

des chemischen Potentials [3]. Diese grundlegende Analogie wird im nächsten Kapitel zur Modellierung der Vergröberung der Oberflächenstrukturen (mounds) benutzt.

KAPITEL 6 Stufenbarrieren im Kontinuumsbild

6.1. Nichtgleichgewichtsstrom

6.1.1. Thermodynamische Grundüberlegung. Eine wachsende Oberfläche ist ein System außerhalb des Gleichgewichts, bei Wachstum unter MBE sogar sehr weit vom Gleichgewicht entfernt. Folgen wir der Onsagerschen Theorie der Nichtgleichgewichtsthermodynamik [60], wie sie in Kapitel 14 von [61] dargestellt ist. In der Sprache dieser Theorie bewirkt eine Affinität, d.h. die Differenz einer intensiven Größe an verschiedenen Stellen im System – in einem kontinuierlichen System der Gradient des Feldes einer intensiven Größe –, die Nichtgleichgewichtseigenschaft des Systems. Einfache Beispiele sind ein räumlich variierendes Temperaturfeld auf der Oberfläche, oder ein variierendes chemisches Potential, wie wir es im letzten Kapitel kennengelernt haben. Jede Affinität bewirkt einen Strom in den extensiven Größen. Der Zusammenhang kann für kleine Affinitäten linear angenähert werden

J_{k} = ∑_{j} L_{jk}∇F_{j}. (62)

Entscheidend ist nun, daß die Matrix der kinetischen (oder Onsager-) Koeffizienten im allgemeinen nicht diagonal ist. Eine Affinität in einer intensiven Größe kann also durchaus einen Strom einer anderen als der zugehörigen extensiven Größe induzieren. Für die wachsende Oberfläche bedeutet dies insbesondere, daß eine Differenz zwischen den chemischen Potentialen auf der Oberfläche und im umgebenden Raum (die Affinität) nicht nur einen Massenstrom auf den wachsenden Kristall (Strom der zugehörigen extensiven Größe) sondern auch entlang seiner Oberfläche (Strom einer anderen extensiven Größe) bewirken kann. Bei Wachstum mit Stufenbarrieren läßt sich der induzierte Oberflächenstrom elementar herleiten. Im Folgenden werden wir sehen, daß auf einer gleichmäßig geneigten wachsenden Oberfläche (parallele Stufen gleichen Abstands) ein (Nichtgleichgewichts– )Oberflächenstrom J_{NG} induziert wird. Aber auch in anderen Wachstumsmodellen ohne erkennbare Asymmetrie auf elementarem Niveau kann ein

64 6. STUFENBARRIEREN IM KONTINUUMSBILD

Nichtgleichgewichtsstrom auf geneigten Oberflächen auftreten und wurde tatsächlich beim Wolf-Villain–Modell [62] beobachtet [63].

Umgekehrt folgt aus (62), daß im Gleichgewicht, wenn alle Affinitäten verschwinden, kein Strom einer extensiven Größe herrschen kann. Auf einer gleichmäßig geneigten Oberfläche, in Kontinuumsdarstellung h(x, t) = m · x, fließt im Gleichgewicht kein Oberflächenstrom. J_{NG} ist ein reiner Nichtgleichgewichtseffekt.

6.1.2. Neigungsabhängigkeit des Oberflächenstroms. Wenden wir uns nach diesen allgemeinen Überlegungen der elementaren Modellierung des Oberflächenstroms zu. Die Symmetrie der Oberfläche wird gebrochen von den Stufen. Deren bevorzugte Orientierung hängt ab von der mittleren Neigung, so daß eine Neigungsabhängigkeit des Oberflächenstroms zu erwarten ist.

Burton, Cabrera und Frank [43] beschreiben die diffundierenden Adatome auf der Kristallobefläche durch das Feld ihrer lokalen Dichte ρ(x, t). Liegen die Stufen dicht genug beisammen, findet auf den Terrassenstücken der Oberfläche keine Nukleation statt. Der natürliche Ansatz bei konstanter Aufdampfrate F, Sprungrate D und Gitterkonstante a ist damit

∂_{t}ρ = a^{2}D Δρ + F. (63)

Verluste in der Adatomdichte durch Nukleation würden (in d=1) beschrieben durch einen Term ∼ −ρ^{2}, der bei kleinen ρ gegen die übrigen Terme in (63) vernachlässigt werden kann. Adatome werden also nur an den Stufen aus dem System genommen, während sie auf den Terrassenstücken mit Rate F je Flächeneinheit hinzugefügt werden. Burton et al. [43] argumentieren nun, daß die Stufenbewegung durch Massenanlagerung auf einer wesentlich langsameren Zeitskala erfolgt, als die Diffusion auf der Terrasse. Zu gegebener Stufenkonfiguration kann man also ρ(x, t) quasistatisch berechnen.

Die Absorptionseigenschaften der Stufen bestimmen dabei die Randbedingungen. Ideal reflektierende Ränder führen zu van-Neumann– Randbedingungen, n̂ · ∇ρ = 0, ideal absorbierende Ränder zu Dirichlet– Randbedingungen ρ = 0. Die gemischten Randbedingungen bei endlichen Stufenbarrieren lassen sich mit folgender einfachen Überlegung herleiten.

Im Kontinuumsansatz (63) wird die Gitterstruktur der Kristalloberfläche vernachlässigt. Insbesondere geht dabei die Struktur entlang der Kristallstufen verloren. Zur Herleitung der Randbedingung an Stufen genügt daher eine Betrachtung in 1+1 Dimensionen. Wir weisen jedoch darauf hin, daß die Berücksichtigung der Gitteranisotropie zu anderen Randbedingungen an Stufen als den hier hergeleteten führt.

6.1. NICHTGLEICHGEWICHTSSTROM 65

Der “Strom” oder besser die Durchsatzrate zwischen zwei benachbarten Gitterplätzen ist das Produkt ihrer Besetzungsdichte und der Sprungrate. Auf der Terrasse ist die Sprungrate D, über die Stufenkante nach unten ist sie erniedrigt um einen Arrheniusfaktor, in den die Erhöhung der Energiebarrieren an den Stufen gegenüber denen zwischen benachbarten Gitterplätzen auf der Terrasse eingeht, r_{−} := exp −E_{B}/k_{b}T =: a/l_{−}. Er wird ausgedrückt als das Verhältnis der Gitterkonstanten a zu einer Länge l_{−}, der sogenannten Schwoebellänge. Sie kann interpretiert werden als “Fluchtweite”, nämlich die Entfernung von der Stufe, die ein Teilchen maximal in seiner Diffusionsbewegung in einer unendlich gedachten Halbebene erreicht, nachdem es einmal von der Stufenkante entkommen ist und bevor es dort endgültig absorbiert wird.

Der Strom über die Stufe ist einseitig, hängt also nur von ρ_{o}, der Besetzung des Platzes oberhalb der Stufe ab: j = r_{−}Dρ_{o}. Damit definiert man eine effektive Besetzungsdichte ρ_{eff} des Gitterplatzes unterhalb der Stufe, die ohne Stufenbarriere denselben Strom bewirkte,

D(ρ_{o} − ρ_{eff}) = j = r_{−}Dρ_{o}.

Interpoliert man das diskrete Dichteprofil durch eine glatte Dichtefunktion ρ(x), entspricht der Ableitung an der Stufe ρ′_{eff} = (ρ_{o} − ρ_{eff})/a = ρ_{o}/l_{−}. Eine Länge l_{+} an der aufwärtsgerichteten Stufe konstruiert man analog.

Für das Dichtefeld in (63) gelten daher die Randbedingungen

∇ρ = \frac{ρ}{l_{±}} n̂. (64)

Der Normalenvektor n̂ zeigt wie der Dichtegradient nach innen. Bei ideal reflektierenden Stufen l_{−} = ∞ wird daraus eine van-Neumann–Randbe-dingung ∇ρ = 0 und bei ideal absorbierenden Stufen eine Dirichlet–Randbedingung ρ = 0.

Im folgenden soll der Nichtgleichgewichtsstrom für eine geneigte Oberfläche mit geraden, äquidistanten Stufen berechnet werden. Die mittlere Neigung sei m, der Stufenabstand also l = a/|m|. Transversale Fluktuationen werden vernachlässigt und die Rechnung wird effektiv eindimensional. Das Dichteprofil besteht aus Parabelstücken, die an den Stufen gemäß (64) aneinanderstoßen. Nehmen wir an den Aufwärtsstufen ideale Absorption an, ρ_{+} =0, so ergibt sich an der Stufe gegenüber

ρ_{−} = \frac{F l^{2}}{2a^{2}D} \frac{1}{1 + l/l_{−}} . (65)

Der lokale Diffusionsstrom ist j = −a^{2}D∇ρ, der makroskopische Nichtgleichgewichtsstrom sein räumlicher Mittelwert

J_{NG} = \frac{1}{l}∫_{0}^{l}dx a^{2}Dρ′(x) n̂ = \frac{F}{2} \frac{l}{1 + l/l_{−}} n̂ (66)

Er fließt bergauf, in Richtung von n̂ = ∇h/|∇h|, und sein Betrag ist proportional zu l = 1/|m|, dem Kehrwert der mittleren Neigung. Im Gleichgewicht, bei F = 0 verschwindet J_{NG}. Ein weiterer wesentlicher Unterschied besteht zwischen dem hier definierten J_{NG} und dem Oberflächenstrom (57) in Kapitel 5. Dort war der Strom ein Tangentialvektor an der durch h(x, t) definierten Fläche, in (66) steht dagegen ein Gradient, also ein Vektor in der Substratebene. Für Neigungen |∇h|≪1 unterscheidet sich (57) jedoch kaum von seiner Projektion in die Substratebene. In diesem Fall vereinfacht sich (57) und man kann den Relaxations- und den Nichtgleichgewichtsstrom addieren zu

J = −σ∇μ + J_{NG} (67)

Die Gleichungen (66) und (67) bilden einen Ausgangspunkt für eine mesoskopische Modellierung der Oberflächendynamik mit Nichtgleichgewichtsstrom, nämlich

∂_{t}h = −∇ · [−σ∇μ + J_{NG}] + F. (68)

In diesem Kapitel werden wir Gleichungen vom Typ (68) untersuchen. Fluktuationen in F, sogenanntes Schrotrauschen, und im Oberflächenstrom, das Diffusionsrauschen, sind irrelevant für das Verhalten auf großen Skalen, wie im nächsten Abschnitt argumentiert wird, und werden daher vernachlässigt.

Allerdings ist die Form (66) noch nicht vollständig. Dieser Ansatz vernachlässigt viele Phänomene auf realen Oberflächen, von denen einige hier kurz erwähnt werden sollen.

(a) Der Strom J_{NG} wird genährt von der Fläche vor einer Stufe, aus der heraus Adatome ungehindert zur Stufe diffundieren können. Bei großen Neigungen ist sie proportional zum Stufenabstand (66), bei kleinen Neigungen kann sie nicht größer als l_{D} werden. Dann trägt also nur ein Bruchteil des auf die Oberfläche treffenden Massenstroms zu J_{NG} bei, mithin |J_{NG}| ∼ l_{D}/l ∼ |m|. Für kleine Neigungen ist der Nichtgleichgewichtsstrom eine lineare Funktion der Neigung m.

(b) Neben der Asymmetrie der Absorption der Stufen beobachten man den sogenannten Knock-out–Effekt. Trifft ein Teilchen aus dem äußeren Strahl direkt an einer Stufe auf, so schiebt es ein Atom in der ersten Reihe

6.1. NICHTGLEICHGEWICHTSSTROM 67

an der Stufe in die untere Terrasse hinaus und setzt sich selbst an seinen Platz. Dies bewirkt einen Massentransport nach unten, also entgegengesetzt zum oben beschriebenen Diffusionsstrom. Für sehr große Neigungen, wenn der Stufenabstand nur noch wenige Gitterabstände beträgt, wird Gleichung (66) den Oberflächenstrom nicht mehr korrekt widergeben, sondern wird dominiert von einem Term proportional zu −m.

(c) Der Kontinuumsansatz (63) für die Adatomdichte vernachlässigt die Anisotropie der Diffusionsbewegung und der Stufenorientierung. Zwischen dicht liegenden Stufen muß die diskrete und damit anisotrope Natur der Diffusionsbewegung berücksichtigt werden. Ferner bestehen je nach Orientierung der Stufe unterschiedliche Wege, auf denen Adatome “hinabrollen” können, und damit eine Orientierungsabhängigkeit der Stufenbarrieren, die in die Randbedingungen (64) eingeht.

Die Herleitung eines mesoskopischen J_{NG} = J_{NG}(∇h, Δh, …) aus allen mikroskopischen Eigenschaften ist also nicht möglich. Zum weiteren Vorgehen setzt man den Nichtgleichgewichtsstrom phänomenologisch an, wobei folgende plausible Eigenschaften seiner Abhängigkeit von der lokalen Oberflächenkonfiguration in den Ansatz einbezogen werden:

• J_{NG} = J_{NG}(∇h) ist im wesentlichen eine stetige Funktion der Neigung. Abhängigkeit von der Krümmung und anderen Termen wird vernachlässigt.

Die Auswirkung der Stufenbarrieren auf das Wachstum tritt hervor bei einer Form, die die oben hergeleitete Neigungsabhängigkeit berücksichtigt.

• Für kleine Neigungen |∇h| ≪ a/l_{D} ist J_{NG} ∼ ∇h und insbesondere bergauf gerichtet.

• In einem großen Bereich |∇h| > a/l_{D} bis hin zu |∇h| ≃ O(1) ist |J_{NG}| ∼ 1/|∇h|. Auch in diesem Bereich ist der Strom bergauf gerichtet.

Beide Eigenschaften treten auf in einer interpolierenden Form von J_{NG}, die in den nächsten Abschnitten eingehend behandelt wird:

J_{NG} = \frac{F l_{D}^{2}/a}{1 + (l_{D}/a)^{2}(∇h)^{2}} ∇h. (69)

Bei sehr hoher Stufendichte kommen weitere Effekte zum tragen, wie z.B. der Knock-out–Effekt.

• Für sehr große Neigungen fließt der Nichtgleichgewichtsstrom wieder bergab. Insbesondere hat J_{NG} einen Nulldurchgang mit zunehmendem |∇h|.

68 6. STUFENBARRIEREN IM KONTINUUMSBILD

In einer linearen Stabilitätsanalyse, im übernächsten Abschnitt, wird die Bedeutung von Nullstellen in J_{NG}(∇h) deutlich. Näher untersucht wird weiter unten die Form

J_{NG} = F l_{D}^{2}/a (1 − (l_{D}/a)^{2}(∇h)^{2}) ∇h, (70)

bei der der Effekt der Nullstelle deutlich zutage tritt.

Es ist sinnvoll, den lateralen und vertikalen Längenmaßstab und die Zeit so zu reskalieren, daß die sich ergebende Gleichung vom Typ (68) dimensionslos wird. Die wichtigsten theoretischen Resultate der letzten beiden Abschnitte dieses Kapitels beziehen sich auf isotrope Modelle mit Relaxation wie in (59) und Nichtgleichgewichtströmen wie (69) und (70), die einfach in eine dimensionslose Form gebracht werden können: Definiere dazu die dimensionslose Höhenfunktion ĥ durch die Beziehung

h(x, t) − F t =: A ĥ(x/l, t/τ), (71)

mit den Längenskalen A := aK^{1/2}/(l_{D}^{2}F^{1/2}) und l = K^{1/2}/(l_{D}F^{1/2}) sowie der Zeitskala τ = K/(l_{D}^{4}F^{2}). Außerdem wurde bei ĥ die mittlere Höhe F t abgezogen, die Oberfläche wird also in einem mitbewegten Koordinatensystem beschrieben, so daß ∫ d^{d}x ĥ ≡ 0. Bei J_{NG} wie (70) erhält man so aus Gleichung (68)

∂_{t}h = − Δ^{2}h − ∇ · [(1 − (∇h)^{2}) ∇h] , (72)

und mit dem Ansatz (69)

∂_{t}h = − Δ^{2}h − ∇ · \frac{∇h}{1 + (∇h)^{2}} . (73)

Vereinfachend ist hier h statt ĥ geschrieben. Aus der Form der Gleichung ist immer ersichtlich, ob h eine dimensionsbehaftete oder -lose Größe ist.

6.2. Vergleich zur Kinetik der Phasenordnung

Bei geeigneter Wahl der Stromfunktion J_{NG}(∇h) läßt sich Gleichung (68) als von einem Ginzburg–Landau–Energiefunktional getriebene Relaxation auffassen. Viele theoretische Vorhersagen und Näherungsmethoden lassen sich damit von der Theorie der Kinetik von Phasenordnungsprozessen [3] übernehmen. Diese Gemeinsamkeit rechtfertigt vor allem eine wesentliche Hypothese für Gleichung (68), die Irrelevanz von Rauschen für das Verhalten der Oberfläche auf Skalen oberhalb der Volumenkorrelationslänge in Dimensionen d > 1. Phasenordnung findet im späten Stadium effektiv bei T = 0 statt, dem Fixpunkt für alle Temperaturen unterhalb

6.2. VERGLEICH ZUR KINETIK DER PHASENORDNUNG 69

T_{c} [3]. Man erwartet also, daß auch für die Vergröberung der Oberflächenstrukturen (mounds) auf großen Skalen Fluktuationen im Diffusionsstrom und im auftreffenden Teilchenstrahl irrelevant sind.

Als wohlberechtigte Annahme und hilfreich in der theoretischen Auswertung hat sich bei Phasenordnungsprozessen die Skalenhypothese erwiesen [3]. Neben der mikroskopischen Skala der Defekte, z.B. der Dicke der Phasengrenzschicht bei skalarem Ordnungsparameter oder der Ausdehnung der Wirbel bei einem O(n)–Modell in d = n Dimensionen, gibt es im System nur eine ausgezeichnete Längenskala L(t), die dem mittleren Abstand zwischen den Defekten entspricht. Korrelationsfunktionen, wie C(x, t) := ⟨ϕ(x+y, t) ϕ(x, t)⟩ haben Skalenform C(x, t) := f(|x|/L(t)), wobei die Länge L(t) ∼ t^{1/z} mit einer bestimmten Potenz in der Zeit zunimmt ([3], Abschnitt 2.1).

Bei Vergröberung von Oberflächenstrukturen beobachtet man sowohl in Experimenten als auch in Simulationen [64, 65] Skalenverhalten. Auch wenn die jeweiligen Exponenten von Material und Wachstumsbedingungen abhängen, treten Potenzgesetze und Skalenverhalten sowohl bei verschiedenen Halbleitern [66] als auch bei mehreren Metallen [67] auf. In die theoretische Modellierung geht die Skalenhypothese folgendermaßen ein: Für die gleichzeitige Korrelationsfunktion der Fluktuationen um die mittlere gewachsene Höhe der Oberfläche, die wie in (71) vom Höhenprofil abgezogen wird

G(x, t) := ⟨h(x+y, t) h(x, t)⟩ (74)

setzt man eine Skalenform an:

G(x, t) = w(t)^{2}g(|x|/ξ(t)). (75)

Auch hier wächst die laterale Längenskala mit einer Potenz in der Zeit,

ξ(t) ∼ t^{1/z}, (76)

zusätzlich nimmt aber auch die mittlere Breite zeitlich zu

w(t) ∼ t^{β}. (77)

Diese beiden Größen, sowie die sich daraus ergebende typische Steigung der Oberflächenstrukturen, der mounds,

s(t) := w(t)/ξ(t) ∼ t^{α} (78)

werden im weiteren Verlauf dieses Kapitels eine wichtige Rolle spielen.

70 6. STUFENBARRIEREN IM KONTINUUMSBILD

Neben diesen Gemeinsamkeiten werfen auch die Unterschiede zwischen beiden Problemklassen wichtige Fragen auf, deren Beantwortung zum weiteren Verständnis führt. Die Verwandtschaft beider Systeme ist leicht zu sehen. Zunächst betrachtet man statt (68) die Dynamik des lokalen Gradientenfelds u(x, t)=∇h(x, t)

∂_{t}u = −∇ ∇· [σ∇ ∇·γ̃ u + J_{NG}(u)] (79)

Hat die Neigungsabhängigkeit des Oberflächenstroms die Form eines Gradienten einer skalaren Funktion

J_{NG}(u) = −∇_{u}Ψ(u), (80)

kann man Gleichung (79) ähnlich der Modell-B–Relaxationsdynamik schreiben,

(81)  ∂_{t}u = −∇ ∇· \frac{δ}{δu}F_{Ψ}\{u\}  mit
F_{Ψ}(u) = ∫ d^{d}x [\frac{1}{2} ∑_{ijlm} γ̃_{lm} (∂_{i}u_{l})σ_{ij}(∂_{j}u_{m}) + Ψ(u)] .

Eine im linearen Anteil isotrope Gleichung wie (59), mit einem Nichtgleichgewichtsstrom wie in (69) ergibt sich also aus dem Funktional

F_{Ψ}\{u\} = ∫ d^{d}x [\frac{1}{2} ∑_{ij}(∂_{i}u_{j})^{2} − log (u^{2} + 1)] . (82)

Laterale und vertikale Längeneinheiten wurden gemäß (71) reskaliert. Ein Strom wie in (70) ergibt sich aus der Wahl

Ψ(u) = −\frac{1}{2}u^{2} (1 − \frac{1}{2}u^{2}) . (83)

Dies ist gerade das Flaschenboden- oder Mexikanerhut-Potential der Ginzburg-Landau–Theorie eines O(n)–Modells in Dimension d = n. Mit einer isotropen Abhängigkeit des Funktionals F_{Ψ} vom Gradienten ∇u wie in (82) entspricht Gleichung (81) fast der Relaxationsgleichung des O(n)– Modells bei T = 0 [3], mit dem wesentlichen Unterschied, daß hier der Differentialoperator ∇∇· statt ∇ · ∇ am Anfang der rechten Seite steht. Zwar ist ∇h ein wirbelfreies Gradientenfeld, nicht aber das Feld J_{NG}, das nichtlinear von ∇h abhängt, so daß die Vertauschung der ∇–Operatoren in Dimensionen d > 1 nicht trivial ist. In d = 1 ist Gleichung (81) mit (83) äquivalent zur erhaltenen Relaxation des ϕ^{4}–Modells, so daß sich alle Ergebnisse übertragen lassen [63].

Besonderes Augenmerk richtet sich auf den Exponenten 1/z, der das Anwachsen der lateralen Längenskala bestimmt. In Kontinuumsmodellen

6.3. LINEARE STABILITÄTSANALYSE 71

für Oberflächen vom Typ (68) ergibt sich immer z = 4, was in den folgenden Abschnitten näher begründet wird. Dieser Wert erscheint auch als Richtwert in vielen Experimenten [66, 67] und Monte–Carlo–Simulationen von Gittermodellen [64], wobei auf Abweichungen von z = 4 noch einzugehen ist. Derselbe Wert für den lateralen Exponenten tritt auch auf bei Modell-B–Relaxation mit einem Ordnungsparameter, der ein Kontinuum von Grundzuständen besitzt, also insbesondere für das erhaltene O(2)– Modell in d=2 [3].

Unterschiede bestehen für ein Potential Ψ(u) mit diskreten Minima. Beim Modell B bilden sich von Wänden getrennte Domänen der einzelnen Phasen, deren Krümmung wegen einer Gibbs-Thomson–Beziehung die Randwerte des chemischen Potentials für die Minoritätsphase(n) in den Domänen festlegt [3]. Bei skalarem Ordnungsparameter erhält man so das t^{1/3}–Wachstumsgesetz von Lifshitz und Slyozov bzw. Wagner [68, 69]. Simulationen der Relaxation eines Ginzburg-Landau–Modells mit zweidimensionalem Ordnungsparameter und diskreten Minima, des Uhrzeiger– oder clock–Modells, ergeben ebenfalls Domänenwachstum mit z =3 [70].

Nicht so beim analogen Oberflächenmodell (79). Auch mit diskreten Minima wachsen die lateralen Strukturen mit z = 4 [70]. Die Ursache dieses Unterschieds ist relativ leicht zu sehen, detaillierte Schlußfolgerungen daraus jedoch noch nicht gezogen. Den Domänen, in denen der Ordnungsparameter jeweils einen Minimalwert von Ψ annimmt entsprechen nun Flächen mit konstantem ∇h, also Seitenflächen der Wachstumspyramiden. Domänenwände sind die Kanten oder Knicke an den Pyramiden. Wie leicht einzusehen ist, können die Knicke nur bestimmte feste Orientierungen haben, nämlich solche, daß die Parallelkomponente von ∇h an ihnen stetig ist. Andernfalls existierte kein eindeutiges Höhenprofil h. Insbesondere können die Knicke also nicht gekrümmt sein und die Bewegung der Wände kann also nicht von der Krümmung getrieben sein.

6.3. Lineare Stabilitätsanalyse

Die in Kapitel 4 beschriebenen Instabilitäten einer ebenen Oberfläche treten auch im Kontinuumsmodell (68) mit (80) auf. Dieses erste Ergebnis des Kontinuumsansatzes soll hier im Detail untersucht werden, nicht zuletzt weil das Konzept einer (in)stabilen Flanke in den folgenden Abschnitten noch benötigt wird.

Bei einer nahezu gleichmäßig geneigten Oberfläche

h(x, t) = m · x + ϵ(x, t)

entwickeln sich die kleinen Störungen gemäß

∂_{t}ϵ(x, t) = −K Δ^{2}ϵ + ∑_{ij} \frac{∂^{2}Ψ(m)}{∂m_{i}∂m_{j}} \frac{∂^{2}ϵ(x, t)}{∂x_{i}∂x_{j}} . (84)

Oberflächen mit Orientierungen, bei denen Ψ eine negative Krümmung aufweist, sind also instabil. Als Beispiele dafür werden im folgenden die Wachstumsinstabilitäten aus Kapitel 4 in der Formulierung der Kontinuumsmodelle dargestellt.

In der Nähe von m = 0 ist in allen isotropen Ansätzen J_{NG} = αm (mit α > 0) und Ψ(m) = −(α/2)m^{2}. Gleichung (84) ist in diesem Fall ∂_{t}ϵ = [−αΔ − KΔ^{2}] ϵ und alle langwelligen Störungen (mit Wellenvektoren |k| < \sqrt{K/α}) wachsen exponentiell in der Zeit an. Es bilden sich Buckel, deren Ausdehnung der am meisten instabilen Wellenlänge ξ ≈ 2π\sqrt{2α/K} entspricht, und die das Äquivalent im Kontinuumsbild zu den zuvor beschriebenen Pyramiden (wedding cakes) sind.

Bleiben wir weiter bei der Näherung eines isotropen Stroms, wenn Ψ also nur von m = |m| abhängt. Dann hat Ψ(m) zwei Hauptkrümmungen, nämlich ∂^{2}Ψ/∂m^{2} in Richtung von m und m^{−1}∂Ψ/∂m in den übrigen d−1 Richtungen orthogonal zu m. Bezeichnen ∂_{∥} bzw. ∇_{⊥} die Gradienten in diesen Unterräumen und J(m) = m · J_{NG}/|m| die Komponente des Stroms in Richtung von m, so kann man (84) umformen zu

∂_{t}ϵ = −KΔ^{2}ϵ + (ν_{∥}∂_{∥}^{2} + ν_{⊥}∇_{⊥}^{2})ϵ , (85)

wobei ν_{∥} = −J′(m) und ν_{⊥} = −J(m)/m. Nimmt der Nichtgleichgewichtsstrom mit zunehmender Neigung ab, so ist der Vorfaktor ν_{∥} > 0 und die Oberfläche ist linear stabil gegen Fluktuationen, die entlang der Richtung von m variieren. Dies ist der Fall beim Stufenflußwachstum, bei dem eine Folge äquidistanter Stufen bei Neigungen |m| > a/l_{D} stabil ist, während sich für kleinere Neigungen auf den Terrassen Inseln bilden können und die Oberfläche instabil ist (siehe Kapitel 4 und [48]). Der phänomenologische Ansatz für J_{NG} wie (69) erfaßt diesen Effekt: J′(m) wechselt gerade bei m=a/l_{D} das Vorzeichen.

Der Vorfaktor ν_{⊥} vor ∇_{⊥}^{2}ϵ bestimmt die (In)Stabilität einer glatten Oberfläche gegenüber Fluktuationen senkrecht zu m. Ist der Strom bergauf gerichtet, so ist ν_{⊥} = −J(m)/m < 0, die Oberfläche h(x, t) = m · x also instabil. Umgekehrt wirkt ein bergab gerichteter Strom stabilisierend. Dies ist konsistent mit der elementaren Anschauung aus Kapitel 4: Aufwärts gerichteter Strom bewirkt die Bales-Zangwill–Instabilität [49], und Abwärtsstrom glättet die lateralen Fluktuationen der Kristallstufen.

6.3. LINEARE STABILITÄTSANALYSE 73

Der Ansatz (70) für J_{NG} beinhaltet für alle Steigungen ab einem bestimmten m_{c} die Positivität von ν_{∥} und ν_{⊥}, also Stabilität in beiden Richtungen. Wenn bei der Anfangskonfiguration überall |∇h(x, t=0)| < m_{c}, so wird sich aufgrund der Instabilitäten die Oberfläche aufsteilen, und die lokalen Steigungen nähern sich, wo immer dies möglich ist, an m_{c} an. Im Sinne des vorangegangenen Abschnitts entspricht dem die spontane Entmischung eines O(n)–Modells mit Potential (83) und einer Anfangsmagnetisierung, die überall |m|<1 aufweist. Da die so in vielen Experimenten und Simulationen beobachtete untere stabile Neigung oftmals nicht einer Flanke mit niedrigen Miller–Indizes entspricht, mithin “unnatürlich” im Vergleich zur Kristallsymmetrie ist, hat sich dafür die Bezeichnung “magic slope” eingebürgert.

Im Rahmen der linearen Stabilitätsanalyse in diesem Abschnitt haben wir Gleichungen erhalten, die gerade den linearen Anteil der in Kapitel 3.3 behandelten anisotropen KS–Gleichung enthalten. Es liegt also nahe, die Vorhersagen der Analyse dieser Gleichung auf Oberflächenwachstum unter MBE zu übertragen. Unter einer zusätzlichen Bedingung erhält man in der Näherung kleiner Neigungen gerade Gleichung (36): Diese unterscheidet sich von den in diesem Kapitel behandelten Gleichungen dadurch, daß sie die Masse auf dem wachsenden Kristall nicht erhält, ihr nichtlinearer Term ist nicht die Divergenz eines Oberflächenstroms. Allerdings unter Wachstumsbedingungen, bei denen Evaporation nicht mehr vernachlässigbar ist, muß die Gleichung (68) um einen Evaporationsterm ergänzt werden:

∂_{t}h = −K Δ^{2}h − ∇ · J_{NG} + F − R(∇h). (86)

Vereinfachend ist hier der Relaxationsterm in der Mullins–Näherung (59) gewählt und die Evaporationsrate R nur in ihrer Abhängigkeit von der Neigung der Oberfläche betrachtet. R hängt hauptsächlich von der Adatomdichte ab, also von der Terrassengröße und somit von der lokalen Neigung [43]. Eine Linearisierung wie in (84) erfordert eine Taylorentwicklung von R nach ∇h. Allerdings verschwindet der lineare Term bei einer Galileitransformation, wie sie in den Kapiteln 3.2 und 3.3 eingeführt wurde, so daß die niedrigste nichttriviale Ordnung zu Termen λ_{∥}(∂_{∥}h)^{2} + λ_{⊥}(∂_{⊥}h)^{2} führt. Eine Abschätzung der Verhältnisse ν_{∥}/ν_{⊥} (nach [43]) und λ_{∥}/λ_{⊥} (nach [29]) für Wachstum bei großen Neigungen ergibt, daß die so erhaltene anisotrope Kuramoto-Sivashinsky–Gleichung im in Kapitel 3.3 beschriebenen instabilen Bereich (siehe Bild 3.8) liegt [71].

Nach der hier untersuchten anfänglichen Entwicklung der Wachstumsinstabilitäten auf nahezu glatten Oberflächen wenden wir uns in den nächsten

74 6. STUFENBARRIEREN IM KONTINUUMSBILD

beiden Abschnitten ihrem Wachstum im Endstadium zu. Der nächste Abschnitt behandelt instabiles Wachstum auf anfänglichen Hochsymmetrieflächen für die Modelle (70) und (69), im letzten Abschnitt dieses Kapitels werden Besonderheiten der Instabilität auf vizinalen Oberflächen mit dem Nichtgleichgewichtsstrom (69) betrachtet.

6.4. Strukturvergröberung auf singulären Oberflächen

Zwei Wege zur Beschreibung des Spätzeitverhaltens sollen hier begangen werden. Ein Skalenansatz (75) mit Potenzgesetzen (76,77) für die eingehenden Längen führt zu Skalenrelationen in Form von Ungleichungen zwischen den auftretenden Exponenten. Ferner liefert eine vom sphärischen Modell oder dem Limes O(n → ∞) inspirierte Approximation [72, 73] ein Multiskalenform von G(x, t) (75), die “nahe” an einer Skalenform liegt mit Exponenten, die die soeben erwähnten Skalenungleichungen erfüllen.

Zur Abkürzung schreiben wir die reskalierten Gleichungen (72) und (73) in der gemeinsamen Form

∂_{t}h = −Δ^{2}h − ∇ · [f(∇h^{2}) ∇h] , (87)

wobei entweder f(m^{2}) = 1/(1+m^{1+γ}) (im folgenden als Modell (i) bezeichnet) oder f(m^{2}) = 1 − m^{2} (Modell (ii)). Modell (i) ist eine Verallgemeinerung von (73), bei der die Form der Abnahme des Nichtgleichgewichtsstroms mit der Neigung mit dem Parameter γ variiert. Er hat Einfluß auf den Exponenten β in w(t), und der Limes γ → ∞ soll mit Modell (ii) verglichen werden.

6.4.1. Skalenrelationen. Betrachten wir zunächst das zeitliche Anwachsen der mittleren Breite der Oberfläche, w(t) := ⟨h(x, t)^{2}⟩,

\frac{1}{2} ∂_{t}w^{2} = −⟨[Δh(x, t)]^{2}⟩ + ⟨∇h(x, t)^{2} f (∇h(x, t)^{2})⟩
(88)  =: −I_{1} + I_{2}.

Randterme treten bei den hier durchgeführten partiellen Integrationen nicht auf, wenn die mittlere Neigung m = 0 ist, und wenn bei einem endlichen System periodische Randbedingungen angenommen werden. Im nächsten Abschnitt wird ein ähnlicher Ansatz für vizinale Oberflächen gemacht. In jedem Fall ist in Gleichung (88) I_{1} ≥ 0. Ist in der Wahl des Modells der Strom immer aufwärts gerichtet, so ist auch an jedem Punkt f(∇h(x, t)^{2})≥0. Aber auch für Modell (ii) erhält man eine untere Schranke für I_{2}: Da die mittlere Breite zeitlich zunimmt, ist wegen I_{1} ≥0

0 ≤ \frac{1}{2} ∂_{t}w(t)^{2} ≤ I_{2} (89)

und daher

I_{1} ≤ I_{2}. (90)

Daraus erhält man sofort eine erste Ungleichung für den Exponenten β, der das Wachstum der Breite der Oberfläche beschreibt. Für Modell (ii) und Modell (i) mit γ ≥ 1 hat ∇h^{2} f(∇h^{2}) und damit auch I_{2} eine obere Schranke. Die Zunahme ∂_{t}w(t)^{2} ≤ const. und w(t) ≤ const. t^{1/2}. Ein von Stufenbarrieren erzeugter Nichtgleichgewichtsstrom macht ebene singuläre Oberflächen instabil, läßt ihre Breite aber nicht schneller als mit der Quadratwurzel der Zeit zunehmen.

Mit der Skalenhypothese in der Folge von (75) kann man abschätzen I_{1} ∼ (s/ξ)^{2} ∼ (w/ξ^{2})^{2} und ∂_{t}w^{2} ∼ w^{2}/t ∼ (sξ)^{2}/t. Bei Modell (i) erhält man außerdem I_{2} ∼ s^{2}f(s^{2}) ∼ s^{1−γ}. Für die oben definierten Skalenexponenten liefern die Ungleichungen (89) 2(α + 1/z) − 1 ≤ α(1 − γ) und (90) 2α − 2/z ≤ α(1 − γ). Beide Ungleichungen zusammen ergeben

\frac{1 + γ}{2} α ≤ \frac{1}{z} ≤ \frac{1}{2} − \frac{1 + γ}{2} α. (91)

Eine weitere Beschränkung besteht darin, daß die Flanken der Mounds stabil sein müssen gegen die laterale Bales-Zangwill–Instabilität [49]. Nach den Überlegungen zur linearen Stabilität (85) führen Fluktuationen mit einer Wellenlänge über 2π/\sqrt{|ν_{⊥}|} zur Instabilität. Zu gegebener Steigung der Flanken s ist dies bei Modell (i) −ν_{⊥} = f(s^{2}) ∼ s^{−(1+γ)}. Die Einschränkung für die laterale Größe der Hügel ist damit ξ ≤ 2π/\sqrt{|ν_{⊥}|} ∼ s^{(1+γ)/2}, bzw. für die Skalenexponenten

\frac{1}{z} ≤ \frac{1 + γ}{2} α. (92)

Damit wird die erste Ungleichung in (91) zur Gleichheit und die zweite ergibt

z ≥ 4, α = \frac{2}{z(1 + γ)} ≤ \frac{1}{2(1 + γ)}, β ≤ \frac{3 + γ}{4(1 + γ)}. (93)

Bei der Dynamik von Modell (ii) nähert sich die typische Steigung der Hügel an den stabilen Wert s=1 an. Nimmt man wieder ein Potenzgesetz an, diesmal s−1 ∼ t^{−α′} mit α′ > 0, so wird der Schätzwert für den letzten Term in (88) I_{2} ∼ s^{2}(1 − s^{2}) ≃ 1 − s^{2} ∼ t^{−α′}. Aus Ungleichung (89) erhält man damit 2/z − 1 ≤ −α′ und aus (90) −2/z ≤ −α′. Eine umgekehrte Abschätzung ergibt sich wie für Modell (i) aus ξ ≤ 2π/\sqrt{|ν_{⊥}|} mit |ν_{⊥}| = 1−s^{2} ∼ t^{−α′}, oder zusammengesetzt

\frac{1}{z} = \frac{α′}{2} = β ≤ \frac{1}{4}. (94)

Fassen wir nochmals die beiden Ergebnisse dieses Abschnitts zusammen. Neben der strengen Schranke für die Zunahme der Breite β ≤ 1/2 erhält man aus einer einfachen Skalenhypothese eine obere Schranke für die Zunahme der lateralen Längenskala ξ der Oberflächenstrukturen, 1/z ≤ 1/4 in beiden Modellen.

Das letzte Ergebnis läßt sich anschaulich so interpretieren. Man stellt zunächst einen Zusammenhang zwischen der lateralen Ausdehnung ξ und der Steigung der Seiten der Hügel s her: Die Bales-Zangwill–Instabilität beschränkt ξ gegenüber einer Funktion der Neigung. Umgekehrt muß die Krümmung der Oberfläche wegen (90) kleiner als eine Funktion der Neigung sein, die laterale Ausdehnung ξ also größer als eine untere Schranke. Da das Anwachsen der mittleren Breite beschränkt ist (89), und w ∼ sξ ergibt sich aus der Beziehung zwischen s und ξ eine Schranke für ξ und damit auch für s (93,94).

6.4.2. Sphärische Näherung. Betrachten wir im folgenden die Zeitabhängigkeit der Zweipunkt–Höhenkorrelationsfunktion (74)

∂_{t}G(x, t) = −2 Δ^{2}G(x, t) − 2 ∇ · ⟨h(0, t) f(∇h(x, t)^{2}) ∇h(x, t)⟩. (95)

Um eine geschlossene Gleichung zu erhalten ersetzen wir f(∇h^{2}) durch f(⟨∇h^{2}⟩) auf der rechten Seite der Gleichung. Diese Näherung entspricht der sphärischen Approximation oder dem Limes n → ∞ der Phasenordnungsdynamik [3], daher die Überschrift dieses Abschnitts. Das Argument von f kann man dann einfach ausdrücken

⟨∇h(x, t)^{2}⟩ = −Δ⟨h(x, t)^{2}⟩ = −ΔG(0, t), (96)

und Gleichung (95) schließt sich zu

∂_{t}G(x, t) = −2 Δ^{2}G(x, t) − 2 f(|ΔG(0, t)|) ΔG(x, t). (97)

Die betrachtete Dynamik (87) ist isotrop in der Substratebene. Gilt dies zudem für die Anfangsbedingungen h(x, t=0), hängt G(x, t) räumlich nur von |x| ab. Entsprechend ist der Strukturfaktor S(k, t), die Fouriertransformierte von G, eine Funktion von t und k =|k| und genügt wegen (95)

∂_{t}S(k, t) = −2 [k^{4} − f(a(t))k^{2}] S(k, t). (98)

Die Funktion a(t) ist definiert durch

a(t) = (2π^{d/2}/Γ(d/2))∫_{0}^{∞}dk k^{d+1}S(k, t), (99)

mit der Dimension d des Substratraums. Formal ist (98) lösbar

S(k, t) = S_{0}(k) exp [−2tk^{4} + 2k^{2}∫_{0}^{t}ds f(a(s))] . (100)

Die Anfangsbedingung S_{0}(k) spiegelt die Unordnung in der Anfangskonfiguration in Gleichung (87) wider. Sie besteht aus den Fluktuationen zu frühen Zeiten des epitaktischen Wachstumsprozesses, also der Konfiguration mit den ersten nukleierten Inseln bei einer Teilbedeckung von θ ≈ 0.1 … 0.3 über der ursprünglichen Substratfläche. Simulationen lassen folgenden Verlauf sinnvoll erscheinen [74]: S_{0}(k) besitzt ein ausgeprägtes Maximum bei einer (endlichen) Wellenzahl, die dem mittleren Inselabstand bzw. der Diffusionslänge l_{D} entspricht. Es verschwindet k → ∞ und nähert sich für k → 0 einem endlichen Wert an.

Zu späten Zeiten behält S(k, t) das Maximum, was ein bestimmtes k_{max} in der Nähe des Maximums im Exponenten von (100) auszeichnet. Die zugehörige Längenskala ξ bezeichnet den typischen Abstand benachbarter Buckel in der Oberflächenstruktur. Da für späte Zeiten k_{max} gegen Null geht, genügt es, S_{0}(k) in der Nähe von k = 0 zu betrachten. Tatsächlich liefert sogar S_{0}(k) ≡ S_{0}(0) = c den korrekten führenden Beitrag zu ξ, womit hier auf noch folgende ausführlichere Bemerkungen vorgegriffen wird. Der Wert von c hat keinen Einfluß auf den dynamischen Exponenten z. Um die Formeln bequem handzuhaben, wählt man S_{0}(k) = (2π)^{−d/2}, bzw. G(x, t= 0) = δ(x).

Aus der im Anhang A bestimmten asymptotischen Form von b(t) := ∫_{0}^{t} ds f(a(s)) erhält man nun die Zeitabhängigkeit der wachsenden Oberflächenstruktur. Fügt man das Ergebnis in den Ausdruck für den Strukturfaktor S(k, t) (100) ein, zeichnet sich für jedes t eine Wellenzahl aus, nämlich die Lage des Maximums von S(k, t).

k_{m}(t) = (\frac{1}{2} (\frac{d + 2}{4} + \frac{1}{γ} ) \frac{log t}{t})^{1/4} . (101)

Man interpretiert es als die Inverse einer ausgezeichneten Längenskala ξ ∼ (t/ log t)^{1/4}. Bis auf eine logarithmische Korrektur ergibt sich laterale Vergröberung der Strukturen mit einem dynamischen Exponenten z = 4. Die im letzten Abschnitt hergeleitete Ungleichung wird dadurch nicht verletzt.

78 6. STUFENBARRIEREN IM KONTINUUMSBILD

Im Gegensatz zum Ausgangspunkt des letzten Abschnitts (75) hat der hier erhaltene Strukturfaktor S(k, t) keine Skalenform S(k, t) = w^{2}k_{m}^{−d}S(k/k_{m}). k_{m} kann also nicht die einzige in dieser Lösung enthaltene Längenskala sein. Vielmehr hat die Lösung eine Multiskalenform

S(k, t) = L(t)^{φ(k/k_{m}(t))/d}, (102)

mit φ(x) = 2x^{2}−x^{4} und der zweiten Längenskala L(t) ∼ t^{( \frac{d+2}{4} + \frac{1}{γ})/d}. Anders als ξ ∼ k_{m}^{−1} hängt die zeitliche Zunahme von L(t) sehr wohl ab von γ, der Form der Stromfunktion f(s^{2}).

Die typische Steigung in den Flanken der Oberflächenbuckel ergibt sich direkt aus der Konstruktion unserer Lösung, a(t) = ⟨(∇h)^{2}⟩. In Modell (ii) nähert sich a(t) dem stabilen Wert s^{2} =1 an. Die führende Abweichung ist

a(t) = 1 − ḃ(t) ≃ 1 − \frac{1}{2}\sqrt{\frac{d + 2}{2}} (\frac{log t}{t})^{1/2} . (103)

Die Skalenungleichung für α′ wird nicht verletzt, sie wäre sogar bis auf eine logarithmische Korrektur als Gleichung erfüllt.

Bei Modell (i) steilen sich die Flanken der Hügel auf,

a(t) ≃ ḃ^{−2/(1+γ)} ≃ ( \frac{8}{\frac{d+2}{4} + \frac{1}{γ}} )^{1/(1+γ)} (\frac{t}{log t})^{1/(1+γ)} , (104)

im Einklang mit der oben hergeleiteten (93) Schranke α = 1/(2 + 2γ).

Die Form der Korrelationsfunktion, wie sie sich aus der sphärischen Näherung ergibt, sei noch näher erläutert. Mit den oben benutzten Anfangsbedingungen S_{0}(k)≡c ist der Strukturfaktor zu jeder Zeit und an jedem k analytisch (100). Die Höhenkorrelationsfunktion G fällt daher exponentiell ab, moduliert mit Oszillationen der Frequenz k_{m}(t).

Interessant ist die Vorhersage der Entwicklung langreichweitiger Korrelationen, wenn sie in der Anfangskonfiguration vorhanden waren. Ein Potenzabfall in G(x, t = 0) bedeutet eine Singularität in S(k, t). Liegt sie an einem Wert k_{0} > 0, ist also der Potenzabfall von G durch Ozsillationen moduliert, so bleibt die Singularität bestehen, wird aber zu späten Zeiten mit einem Faktor exp −tk_{0}^{4} unterdrückt, wenn k_{m}(t) ≪ k_{0}. Für extrem große |x| hat G(x, t) also einen sehr schwachen mit einem Potenzgesetz abfallenden “Schwanz”, nachdem es bis zu einem zeitlich zunehmenden x_{0} fällt exponentiell abgefallen ist.

Eine Singularität bei k_{0} = 0 wird dagegen nicht unterdrückt. In Gleichung (100) wird sie immer mit einem Faktor exp 0 = 1 multipliziert. Ein Potenzabfall von G ohne Oszillationen bleibt bestehen.

6.4. STRUKTURVERGRÖBERUNG AUF SINGULÄREN OBERFLÄCHEN 79

Die oben hergeleiteten Skalengesetze verändern sich nicht, wenn S_{0}(k) am Ursprung singulär ist. Sei z.B. S_{0}(k) ∼ k^{σ} für k → 0. Die Potenz σ kann in das Integral über den Phasenraum (99,126) einbezogen werden, und führt zu einer Verschiebung der effektiven Dimension d → d+σ. Dies verändert die Vorfaktoren vor den Potenzgesetzen (101),(103) und (104), nicht aber die Potenzen 1/z, α, α′ und β selbst.

Abbildung 6.1. Höhenkorrelationsfunktion G(x, t) von Modell (i) mit γ = 1 als Ergebnis einer Simulation und der sphärischen Näherung. Die Kurven sind für Zeiten t = 1000, 1100, …, 10000 so reskaliert, daß G(0, t) und die erste Nullstelle jeweils zu 1 gesetzt wurden. Die Oszillationen der sphärischen Näherung sind deutlicher ausgeprägt als bei der vollen Dynamik. Im Eisatz ξ(t), die erste Nullstelle von G als Funktion der Zeit. Die durchgezogenen Linien sind beste Anpassungen der Form (t−t_{0})^{1/z} bzw. ((t−t_{0})/ log(t−t_{0}))^{1/z}, mit variabel angepaßtem Nullpunkt der Zeitskala t_{0}.

In Bild 6.1 vergleichen wir die aus numerischen Integrationen zu zufälligen Anfangsbedingungen erhaltene Korrelationsfunktion G(x, t) mit der Lösung in sphärischer Näherung. Berechnet wurden beide für Modell (i) mit γ =1, zu Zeiten t = 1000, 1100, …, 10000, und in reskalierter Form gemeinsam dargestellt. Die erste Nullstelle und die Höhe am Ursprung wurden jeweils fixiert, G(x, t) = w(t)^{2}g(|x|/ξ(t)) zu g(0)=1 und g(1)=0. Die “Skalenform” der Näherungslösung oszilliert für große Entfernungen stärker und

80 6. STUFENBARRIEREN IM KONTINUUMSBILD

fällt etwas weniger schnell ab. Dennoch stimmen beide Kurven bemerkenswert gut überein, was die zunächst unkontrollierte Näherung von (95) zu (97) nochmals rechtfertigt.

Eine Beobachtung bietet einen Anlaß, das Skalenverhalten des nicht genäherten G genauer zu untersuchen. Die reskalierten Näherungslösungen liegen nahezu aufeinander, die logarithmischen Korrekturen zum Skalenverhalten scheinen sehr schwach. Es stellt sich die Frage, ob die “exakte” Lösung des Modells skaliert, ob die Streuung der Datenpunkte statistische oder systematische Abweichung ist. Ein sensibler Test wäre, aus der numerischen Lösung den Strukturfaktor und den skalenabhängigen Exponenten φ(k/k_{m}) zu bestimmen, was für die Phasenordnung des O(n)–Modells vorgenommen wird [75].

6.4.3. Weitere Näherungsverfahren. Die sphärische Approximation ist nur eines von mehreren möglichen Näherungsverfahren für die Dynamik der Zweipunktkorrelationen [76, 77, 78, 3]. Erarbeitet wurden sie für die Phasenentmischung mit erhaltener Magnetisierung (Modell-B– Dynamik), weshalb sie in den Oberflächenmodellen am leichtesten auf die Dynamik des Gradientenfeldes, in Gleichung (79), angewandt werden. Allerdings eignen sie sich nur für den Fall endlicher nicht divergierender Steigungen an den Flanken der wachsenden Hügel, Model (ii) aus dem letzten Abschnitt nach Gleichung (72).

Die gleichzeitige Zweipunkt-Korrelationsfunktion des Gradientenfelds

C(x, t) := ⟨u(x, t) · u(0, t)⟩

genügt analog zu 95) wegen (79) der Gleichung

\frac{1}{2} ∂_{t}⟨u(x, t) · u(0, t)⟩ = −∂_{i}∂_{j}[⟨u_{i}(0, t)Δu_{j}(x, t)⟩

+⟨u_{i}(0, t) f(u(x, t)^{2})u_{j}(x, t)⟩], (105)

wobei sich die Ableitungen ∂_{i} auf x beziehen, und über gleichlautende Indizes summiert wird. Der Unterschied in der Abfolge der Differentialoperatoren zwischen (81) und der entsprechenden Gleichung für Phasenordnung im O(n)–Modell, ∇∇· statt Δ spielt in der Tat für Gleichung (105) keine Rolle. Für den ersten Term auf der rechten Seite ist das leicht einsichtig, denn ∂_{i}u_{j} = ∂_{ij}^{2}h = ∂_{j}u_{i}. Zum zweiten Term überlegt man sich, daß für einen beliebigen Nichtgleichgewichtsstrom J(u) wegen der Translationsinvarianz der Korrelationsfunktion

\frac{∂}{∂x_{i}} \frac{∂}{∂x_{j}} ⟨u_{i}(y, t)J_{j}(u(x, t))⟩ = \frac{∂}{∂y_{i}} \frac{∂}{∂y_{j}} ⟨u_{i}(y, t)J_{j}(u(x, t))⟩,

6.4. STRUKTURVERGRÖBERUNG AUF SINGULÄREN OBERFLÄCHEN 81

gilt und damit dasselbe Argument wie für den ersten Term.

Die zeitliche Entwicklung der Korrelationsfunktion (105) gleicht also der entsprechenden Gleichung für magnetische Systeme

(106) \frac{1}{2} ∂_{t}C(x, t) = −Δ^{2}C(x, t) − Δ⟨u(0, t) · f(u(x, t)^{2})u(x, t)⟩.

Natürlich besteht zwischen beiden Modellklassen – denen mit einem Magnetisierungsfeld und denen mit Gradientenfeld eines Höhenprofils – ein Unterschied im Term höherer Ordnung, der jedoch verloren geht, wenn der zweite Term durch einen Ausdruck mit C(x, t) angenähert und so Gleichung (106) geschlossen wird. Die Ergebnisse in einer entsprechenden Näherung sind daher in beiden Modellklassen identisch.

Eine Verbesserung gegenüber der sphärischen Näherung ließe sich insofern erzielen, als die Lösung z.B. in Mazenkos Näherung [78] konventionelles Skalenverhalten aufweist, was numerischen Auswertungen wohl besser entspricht als das Multiscaling-Szenario (s.o.). Möglicherweise liegt auch die Näherungsform von C(x, t) näher an den Daten der numerischen Integration als in Bild 6.1. Allerdings bedeutet eine solche Verfeinerung der Approximation keine fundamentale Verbesserung in der Berechnung des Skalenverhaltens. Insbesondere tragen diese Näherungsformen nicht dazu bei, den im Abschnitt 6.2 angesprochenen Unterschied bei diskreten Minima in der Potentialfunktion des Nichtgleichgewichtsstroms Ψ(u) – ein Wachstumsgesetz ξ(t) ∼ t^{1/4} bei wirbelfreien Feldern des Ordnungsparameters gegenüber ξ(t) ∼ t^{1/3} im anderen Fall – theoretisch zu klären.

Die Renormierungsgruppenanalyse von Bray [3], die Aufschluß über Skalenexponenten geben könnte, ist auf das Wachstum der Oberflächenstrukturen nicht übertragbar. Ihr Ausgangspunkt ist die dynamische Gleichung (81) im Fourierraum

∂_{t}u_{i}(k, t) = k_{i}k_{j} \frac{δ}{δu_{j}} F_{Ψ}\{u\}. (107)

Anstelle von k^{2}δ_{ij} beim O(n)–Modell steht in Gleichung (107) k_{i}k_{j}. Als Matrix ist aber k_{i}k_{j} selbst für k ≠ 0 nicht invertierbar. δF_{Ψ}/δu ist daher nicht nur im infraroten Limes bei k = 0 singulär, weshalb das Renormierungsschema von [3] nicht angewandt werden kann.

Gleich dem Renormierungsansatz liefert eine elegante Abschätzung der Energiedissipation während des Phasenordnungsprozesses die Skalenexponenten [79]. In Verbindung mit der Skalenhypothese ergeben sich auch für die hier behandelten Probleme sinnvolle Werte der Exponenten [80]. Leider läßt sich aus einem Ansatz der Dissipation der freien Energiedichte keine

82 6. STUFENBARRIEREN IM KONTINUUMSBILD

scharfe Ungleichung wie aus (88) für die Breite w(t) gewinnen. Die Terme im entsprechenden Ausdruck tragen kein eindeutiges Vorzeichen.

Wie auch auf dem Gebiet der Kinetik der Phasenordungsprozesse stehen auch bei der Dynamik der Vergröberung von Oberflächenstrukturen noch viele Fragen offen. Einer ersten Näherung, wie sie in diesem Abschnitt präsentiert wurde, sollten noch viele Verbesserungen folgen. Wenden wir uns nun der Situation zu, die der Phasenordnung mit mittlerer Magnetisierung m ≠ 0, aber unter dem Gleichgewichtswert im Volumen, |m| < m_{0}, nach einem sogenannten “off critical quench” entspricht, nämlich einer vizinalen Oberfläche mit einer Neigung unterhalb der stabilen magic slope.

6.5. Instabilität auf vizinalen Oberflächen

Das zentrale Ergebnis der linearen Stabilitätsanalyse im vorletzten Abschnitt war die Instabilität geneigter Oberflächen mit aufwärtsgerichtetem Nichtgleichgewichtsstrom gegen laterale Fluktuationen und bei hinreichend großer Neigung (|m| > a/l_{D}) ihre lineare Stabililtät gegen Fluktuationen in den Stufenabständen. Die durch die Neigung aufgeprägte Anisotropie spielt also eine wesentliche Rolle in der Oberflächendynamik. Andererseits ergeben sich bei Wachstum auf einer singulären Oberfläche generisch Konfigurationen mit Buckeln, deren Flanken sich aufsteilen und die ursprüngliche Neigung m deutlich übertreffen können. Die Anfangsneigung scheint somit irrelevant, und die gleiche Endkonfiguration müßte auch von einer vizinalen Oberfläche ausgehend erreicht werden. In 1+1 Dimensionen ist ein solches Szenario untersucht worden, und man findet die Metastabilität des Stufenflußwachstums[44]. In 2+1 Dimensionen erweist sich die laterale Instabilität als Vorbereiterin zur Bildung von (im wesentlichen symmetrischen) Hügel. Die zunächst entstehenden entlang der Neigung verlaufenden Rippen sind ihrerseits instabil, brechen auf und ermöglichen so die Bildung der Buckel, die das Endstadium bestimmen. Dies gilt sowohl für Kontinuums- als auch für Gittermodelle, und die Beobachtung des Übergangs zwischen beiden Stadien erlaubt die Anpassung der Parameter in der Kontinuumsgleichung [81].

Früh- und Spätstadium des Wachstums sind in den Bildern 6.2 und 6.3 dargestellt, die Aufsicht (Bild 6.2) zeigt eine Rippenkonfiguration im Gittermodell. Im nächsten Bild 6.3 finden sich in Grauskalendarstellung zwischen Gitter- und Kontinuumsmodell vergleichend die Rippen und die daraus im Endstadium entstandenen Buckel. Das Rastertunnelmikroskop nachahmend ist wie in den ersten Kapiteln die Darstellung so angepaßt, daß Abbildung 6.2. Aufsicht auf eine Konfiguration des in diesem Kapitel vorgestellten SOS–Modells (Teilansicht eines Systems der Größe 300×300). Auf der ursprünglich gleichmäßig geneigten Fläche haben sich Rippen gebildet.

schwarze Farbe die tiefsten Stellen, weiße die höchsten Stellen der Oberfläche anzeigt. Um durch diese Kontrastierung wirklich die Struktur der Oberfläche aufzulösen, wurde die Anfangskonfiguration abgezogen: Bei Kontinuumsmodell wurde m·x, beim Gittermodell die anfängliche “Treppenkonfiguration” abgezogen, was noch an den senkrecht durchs Bild verlaufenden Kontrasten zu erkennen ist. An dieser Stelle verweisen wir nochmals auf Bild 1.1: Es zeigt einen Ausschnitt aus einer Platin–(111)–oberfläche mit einem vizinalen Flächenstück, einer sogenannten Makrostufe (macrostep), nach Deposition von weiteren fünf Monolagen Platin. Auf dem singulären Oberflächenstück sind Pyramiden gewachsen, auf dem geneigten Teil haben sich in Hangrichtung orientierte Rippen gebildet, die allerdings aufbrechen und Pyramiden weichen.

6.5.1. Skalenabschätzung. Wie für anfangs singuläre, so lassen sich auch für das Wachstum auf geneigten Oberflächen Skalenrelationen herleiten. Wegen der neigungsbedingten Anisotropie in der Substratebene müssen jedoch bei einem Skalenansatz wie (77) ff. die lateralen Längenskalen Abbildung 6.3. Vergleich des Früh- und Spätstadiums des Wachstums in Grauskalendarstellung (siehe Text) im Kontinuums- und Gittermodell. Isotrope Buckel (Hügel) treten an die Stelle der anfangs gebildeten Rippen.

ξ_{∥} und ξ_{⊥} und die zugehörigen Neigungen s_{∥} und s_{⊥} unterschieden werden. Die Irrelevanz der anfänglichen Neigung m für die Moundkonfigurationen zu späten Zeiten wird sich unter anderem darin zeigen, daß sich ξ_{∥} und ξ_{⊥} einander annähern. Ein erster Indikator dafür sind die Skalenrelationen für die Skalen entlang (∥) und senkrecht zu (⊥) m, die jeweils beide den Beziehungen für Modell (i) (93) und Modell (ii) (94) entsprechen.

6.5. INSTABILITÄT AUF VIZINALEN OBERFLÄCHEN 85

Um wie in (88) partiell integrieren zu können, benötigt man ein Höhenfeld mit periodischen Randbedingungen. Deshalb betrachtet man in einem gekippten Koordinatensystem h̃(x, t) := h(x, t)−m·x. Gleichung (87) schreibt sich in diesem Feld

∂_{t}h̃ = −Δ^{2}h̃ − ∇ · [f ((∇h̃+m)^{2}) (∇h̃+m)] , (108)

Die typische Aufsteilung in Längsrichtung, s_{∥}^{2} := ⟨(∂_{∥}h̃)^{2}⟩, wächst an wie

\frac{1}{2}∂_{t}s_{∥}^{2} = − ⟨(∂_{∥}Δh)^{2}⟩ + ⟨f(∇h^{2})(∂_{∥}∇h)^{2}⟩

+ 2 ⟨f′(∇h^{2})(∇h · ∂_{∥}∇h)^{2}⟩ (109)

=: − I_{1}^{∥} + I_{2}^{∥} − I_{3}^{∥}.

Offensichtlich ist I_{1}^{∥} > 0. Ebenso I_{3}^{∥} > 0, weil für die betrachteten Modelle (i) und (ii) f′(s^{2}) = df/d(s^{2}) < 0. I_{2}^{∥} ist im Modell ohne Abwärtsstrom (i) in jedem Fall positiv, muß aber auch in Modell (ii) postiv sein, damit s_{∥}^{2} zunimmt. , (Natürlich gilt dies nicht für eine stabile Orientierung |m| ≥ 1. Analog zu (89) und (90) erhält man hier die Ungleichungen

0 ≤ ∂_{t}s_{∥}^{2} ≤ I_{2}^{∥} (110)

und

I_{1}^{∥} ≤ I_{2}^{∥} sowie I_{3}^{∥} ≤ I_{2}^{∥}. (111)

Drückt man nun I_{i}^{∥} durch die Längenskalen ξ_{∥}, s_{∥} bzw. w aus, ergeben sich dieselben Skalenrelationen wie (93) und (94). Für s_{⊥}^{2} := ⟨(∂_{⊥}h̃)^{2}⟩ gilt dieselbe Rechnung und damit in beiden Richtungen (∥ wie ⊥) dieselben Skalenrelationen.

6.5.2. Numerische Ergebnisse. Aus der Zweipunkt-Korrelations-funktion von h̃

(112) G(x−y, t) := ⟨ (h(x, t)−m·x) (h(y, t)−m·y) ⟩

lassen sich in Simulationen alle oben bezeichneten Längenskalen gewinnen. Wie auf der isotropen Oberfläche ist auch hier w(t)^{2} = G(0, t). Weiter untersucht man G entlang den bezüglich m longitudinalen und orthogonalen Richtungen. Diese Einschränkungen werden bezeichnet mit G_{∥}(x, t) := G((x/|m|)m, t) und ihr Gegenstück mit G_{⊥}. Die lateralen Längenskalen ξ_{∥} (ξ_{⊥}) erhält man nun aus der Breite von G_{∥} (bzw. G_{⊥}). Besser als die erste Nullstelle von G_{∥(⊥)} eignet sich für unsere Zwecke die Breite, bei der die Korrelationsfunktion auf die Hälfte ihres Maximalwerts

86 6. STUFENBARRIEREN IM KONTINUUMSBILD

G_{∥(⊥)}(0, t) abgefallen ist. Analog zu (78) erhält man aus der Korrelationsfunktion die typischen Steigungen parallel und senkrecht zu m, s_{∥}:= w/ξ_{∥} und ebenso für s_{⊥}. Wir dürfen annehmen, daß sich diese Definitionen zu späten Zeiten nur durch einen unwesentlichen Faktor von den im letzten Abschnitt gebrauchten ⟨(∂_{∥(⊥)}h̃)^{2}⟩ unterscheiden.

6.5.2.1. Das Solid-on-Solid–Modell. Es ist sinnvoll, an dieser Stelle detaillierter auf die mikroskopische Modellierung im SOS–Modell einzugehen, wie sie zum Vergleich mit der Kontinuumsbeschreibung vorgenommen wurde. Die Simulationen wurden ausgeführt auf Substraten der Größe 150×150 und 300×300 mit periodischen Randbedingungen. In Richtung der Oberflächenneigung wird dabei ein “Sprung” von 25 bzw. 50 Gitterebenen berücksichtigt, weshalb man auch von “helikalen” Randbedingungen spricht. Die mittlere Neigung beträgt damit |m|=1/6. Sie liegt bei den hier präsentierten Simulationen deutlich über a/l_{D}, was aus der Inselbildung auf einer singulären Oberfläche zu 1/14.3 bestimmt wurde (siehe unten). Anfangs wächst also die Oberfläche im Stufenflußmodus.

Die Dynamik besteht aus zwei elementaren Prozessen, dem Auftreffen eines Teilchens auf der Oberfläche, h_{x} → h_{x}+a, und dem Hüpfprozeß entlang der Oberfläche, h_{x} → h_{x} −a und h_{x+e} → h_{x+e} +a, wobei x+e einen Nachbarplatz von x in der Substratebene des Gitters bezeichnet. Teilchen treffen auf mit Rate F, die Hüpfrate hängt ab von der Temperatur und der Energiebarriere E zwischen den benachbarten Plätzen x und x+e. Ein oft gewählter Ansatz [64] ist k(E, T) = k_{0} exp(−E/k_{B}T). Der Vorfaktor k_{0} ist die Rate, mit der die Adatome von den Phononen im Kristall zu Sprüngen angeregt werden. Sie kann sinnvoll abgeschätzt werden durch k_{0} =2k_{B}T/h [43, 64]. Die Sprungbarriere zwischen benachbarten Plätzen, E = E_{S} + nE_{N} + (n_{i} − n_{f})θ(n_{i} − n_{f})E_{B}, hängt über die Parameter E_{S}, E_{N} und E_{B} von n, der Anzahl der bindenden nächsten Nachbaratome in derselben Gitterebene, sowie der Anzahl der Nachbaratome in der nächsttieferen und nächsthöheren Ebene vor (i) und nach (f) einem Sprung ab. Der Term nE_{N} beschreibt also Arrhenius–artige Sprungraten, der Term mit E_{B} modelliert die erschwerte Annäherung an eine Stufenkante von oben. Eine detailliertere Beschreibung der Dynamik des SOS–Modells für MBE findet sich in [64], in der hier verwendeten Form in [81].

Die in die Energiebarriere eingehenden Parameter wurden experimentellen Beobachtungen angepaßt. Quantitativ ausgewertet wurde ein Parametersatz, der das Wachstum von Platin bei T = 400 K mit F = 1/40 Monolagen je Sekunde beschreibt: E_{S} = 0.75 eV, E_{N} = 0.18 eV und E_{B} = 0.15 eV. Dies entspricht relativ starken Stufenbarrieren, außer von der obersten

6.5. INSTABILITÄT AUF VIZINALEN OBERFLÄCHEN 87

Insel eines Buckels sind Abwärtssprünge von Adatomen nahezu vollständig unterdrückt.

6.5.2.2. Längen- und Zeitskalen in den Simulationen. Die sich aus den Simulationen ergebenden Längenskalen sind hier in einer Übersicht dargestellt. Die beiden Abbildungen in 6.4 zeigen die zeitliche Entwicklung von ξ_{∥} und ξ_{⊥}, sowie in der eingefügten Graphik die typischen Steigungen s_{∥} und s_{⊥}. Ergänzend findet man in 6.5 die zunehmenden Breiten der Oberfläche und deren logarithmische Zunahme mit der Zeit d log w(t)/d log t.

Abbildung 6.4. Die lateralen Längenskalen ξ_{∥} und ξ_{⊥} als Funktion der Zeit, links für die Kontinuumsgleichung (87) mit Nichtgleichgewichtsstrom nach Modell (i) und rechts für das im Text beschriebene SOS–Modell (Zeitskala in Monolagen). Zu späten Zeiten verschwindet in beiden Modellen der Unterschied zwischen ∥– und ⊥–Richtung, ξ_{∥} ∼ ξ_{⊥} ∼ t^{1/4}. Eingefügt die Entwicklung typischer Steigungen s_{∥(⊥)} an den Oberflächenstrukturen.

In beiden Modelltypen ist zunächst ξ_{∥} deutlich größer als ξ_{⊥}. Anfangs dominieren die entlang m orientierten Rippen. Dieser Unterschied ist beim Kontinuumsmodell wesentlich ausgeprägter als beim SOS–Modell. In der linearen Näherung (85) ist ν_{∥} > 0, im Anfangsstadium werden daher Fluktuationen entlang m gedämpft. Im Gegensatz dazu fluktuieren die Rippen beim SOS–Modell deutlich, siehe Bild 6.3, und begrenzen ξ_{∥} durch ihre Windungslänge.

Im Endstadium der Dynamik nähern sich in beiden Modellklassen orthogonale und longitudinale Längenskalen einander an. Die gemeinsame Längeskala ξ_{∥} =ξ_{⊥} =:ξ ∼t^{1/z} wächst mit einer Potenz, die nahe am theoretischen Wert z = 4 liegt. Unterschiedlich verhält sich bei beiden Modellen die Breite w(t). Das Kontinuumsmodell (Modell (i) mit γ = 1) nähert sich dem theoretischen Wert β = (3+γ)/(4γ+4) = 1/2, während die Breite der Abbildung 6.5. Die Breite der Oberfläche (Höhe der Hügel bzw. Rippen) als Funktion der Zeit (bzw. der aufgedampften Schichtdicke) und ihre logarithmische Ableitung, für beide Modellierungen wie in der vorhergehenden Abbildung.

SOS–Oberfläche mit β ≃ 0.35 zunimmt. Nach den Überlegungen des vorhergehenden und vor allem des ersten Abschnitts dieses Kapitels sehen wir, daß der Aufwärtsstrom beim SOS–Modell für starke Neigungen stärker als 1/|m| abnimmt. Eine stabile Neigung oder magic slope wird bei den Simulationen nicht erreicht. Die Neigungen s_{∥} ≃ s_{⊥} nehmen mit einer Potenz α ≃ 0.14 zu. Zwei einfache mikroskopische Erklärungen für eine Abwärtskomponente des Nichtgleichgewichtsstroms sind denkbar: Werden Adatome an den Stufenkanten deponiert, springen sie bevorzugt nach unten. Zudem können Stufen in (1,1)–Richtung mit der oben beschriebenen Sprungbarriere E entlang ihrer “Kerben” leichter überwunden werden als die “glatten” (1,0)–Stufen.

6.5.2.3. Die Übergangsphase. Auch in der Übergangsphase zwischen anisotropem und isotropem Wachstum lassen sich beide Modelle vergleichen. Insbesondere kann man hier die Parameter im Kontinuumsmodell zur Anpassung an das SOS–Modell eichen. Die Frage ist, wievielen dimensionslosen Zeit- und Längeneinheiten im Kontinuumsmodell die “natürlichen” Einheiten des Gittermodells, die Gitterkonstante und die Wachstumszeit einer Monolage, entsprechen. Dies kann z.B. über die Zeitskala τ erfolgen: Genau dann, wenn die Rippen aufbrechen, hat ξ_{∥} ein Minimum. Diese Zeitpunkte stimmen überein, wenn das Wachstum einer Monolage in der dimensionslosen Kontinuumsgleichung 120 Zeiteinheiten entspricht.

Aus der natürlichen Wahl der Zeitskala des Gittermodells erfolgt F =1, die Diffusionslänge l_{D} dagegen kann man aus der Inselbildung auf einer großen Hochsymmetriefläche ablesen. Bei obiger Parameterwahl ist l_{D} ≃

6.6. INSTABILES WACHSTUM, EINE SCHLUSSBETRACHTUNG 89

14.3 Gitterplätze. Eine Gitterkonstante im SOS–Modell wird daher durch l =l_{D}\sqrt{τ} ≃1.31 dimensionslose Längeneinheiten repräsentiert.

Vergleicht man nun mit dieser Reskalierung die Verläufe von ξ_{∥}(t) und ξ_{⊥}(t) in beiden Teilen von 6.4, so fallen merhrere Unterschiede auf. Wie bereits erwähnt ist zu frühen Zeiten ξ_{∥} im Kontinuumsmodell deutlich größer als im SOS–Modell. Die Breite der Rippen, ξ_{⊥} zu frühen Zeiten, stimmt in beiden Modellen dagegen gut überein.

Allerdings unterscheiden sich beide Modelle durch die Breite der Mounds im Endstadium: ξ_{∥(⊥)} ist im Kontinuumsmodell etwa doppelt so groß wie im Gittermodell. Diese Unterschied entsteht bereits in der Übergangsphase. Im linken Teilbild von 6.4 sieht man, daß ξ_{⊥} gerade dann besonders stark zunimmt, wenn ξ_{∥} abnimmt; wenn die Rippen aufbrechen, entstehen Buckel, die ihre Breite deutlich überragen. Nicht so beim SOS– Modell (rechtes Teilbild): ξ_{⊥} nimmt hier gleichmäßig zu, ohne den ausgeprägten Sprung beim Übergang.

Ein Vergleich zweier Konfigurationen in der Übergangsphase macht dies deutlich (Bild 6.6). Das Rippenmuster ist vereinzelt aufgebrochen und Buckel haben sich gebildet. Diese Buckel haben in beiden Modelltypen die Vorzugsrichtung m verloren. Allerdings sind sie im Kontinuumsmodell deutlich breiter als die zugrundeliegenden Rippen, während sie im SOS– Modell etwa deren Breite behalten. Eine Ursache dafür: Das hier verwendete Kontinuumsmodell hat einen isotropen Relaxationsterm K =σγ̃ in der Näherung nach Mullins (59). Auf einer vizinalen Oberfläche ist aber der Tensor der Biegesteifigkeit γ̃ anisotrop, insbesondere an den am stärksten geneigten Stellen, den Flanken der Bales-Zangwill–Rippen. Eine Krümmung der Oberfläche in Richtung der Neigung variiert die Abstände der Stufen, eine Krümmung senkrecht dazu aber deren Länge [82], weshalb die Oberfläche in dieser Richtung deutlich steifer ist. Ein Kontinuumsmodell mit neigungsabhängiger und anisotroper Biegesteifigkeit und entsprechend tensoriell neigungsabhängigem K müßte also diese Ausbuchtungen der entstehenden Buckel unterdrücken. Eine Abschätzung der Anisotropie und ihrer Auswirkung auf die Verbreiterung der Buckel findet sich im Anhang von [81].

6.6. Instabiles Wachstum, eine Schlußbetrachtung

Die Instabilität von Epitaxie auf Hochsymmetrieflächen ist experimentell und in Simulationen in den letzten Jahren eingehend untersucht worden Abbildung 6.6. Konfigurationen im Übergang zwischen Rippen- und Hügelmuster. Hügel haben sich vor allem an den Defekten im Ripenmuster gebildet. Ein Unterschied zwischen Gittermodell (oben) und Kontinuumsgleichung (unten) zeigt sich am Größenverhältnis der Buckel zur Rippenbreite (siehe Text).

[66, 67]. Der Versuch diese Phänomene theoretisch in Kontinuumsbeschreibung zu fassen führt zu Problemen, wie sie auch in der Theorie der Phasenordnungskinetik bestehen. In 1+1 Dimensionen sind die Lösungsansätze für Phasenentmischung direkt auf die Oberflächendynamik übertragbar [63], während in höheren Dimensionen fundamentale Unterschiede bestehen, die die Schwierigkeiten zusätzlich vergrößern [70]. In jenem Bereich erarbeitete Näherungsverfahren müssen neu auf diese Problemstellung angewandt

6.6. INSTABILES WACHSTUM, EINE SCHLUSSBETRACHTUNG 91

werden. Ein erster Versuch ist im Abschnitt 6.2 dieses Kapitels dokumentiert. Aus dem Skalenansatz (75) wurde jüngst für (in unserer Bezeichnung) Modell (i) mit γ = 1 ein Satz gewöhnlicher Differentialgleichungen für ξ(t) und w(t) abgeleitet, dessen Lösung mit unseren Ergebnissen übereinstimmt [80].

Hinreichend stark geneigte Oberflächen können wachsen, indem ihre Stufen gleichmäßig über die vor ihnen liegenden Terrassen streichen. Den gerichteten Materiefluß bewirken Stufenbarrieren, die Abwärtsspünge erschweren (siehe Kapitel 4, auch [48] und [29]). In Abschnitt 6.5 in diesem Kapitel wurde eine erste Beschreibung aufeinanderfolgender Instabilitäten des Stufenflußwachstums präsentiert [81]. Die Mäander-Instabilität der Stufen [49] führt zu einer Rippenstruktur, die ihrerseits den Weg zur Bildung dreidimensionaler Buckel bereitet. Im Endstadium nähert sich das Wachstum dem auf ursprünglich singulären Flächen an.

Dieser Übergang wurde in zwei fundamental verschiedenen Modellen untersucht, die für entgegengesetzte Grenzfälle geeignet sind. Oberflächen mit sehr starken Stufenbarrieren haben kurze Zeit- und Längenskalen der instabilen Wachstumsphänomene und werden durch ein Gittermodell gut beschrieben. Der Limes schwacher Stufenbarrieren führt zu einer großen Kapillaritätslänge. Die Oberfläche variiert kaum auf der Skala der Gitterkonstanten, weshalb die Kontinuumsbeschreibung angemessen ist (siehe Kapitel 5).

Dennoch liefern beide Modellierungen vergleichbare Ergebnisse, bedingt sogar quantitative Übereinstimmung. Wir dürfen daher annehmen, daß das beschriebene Phänomen robust und auch experimentell beobachtbar ist. Einen ersten Hinweis dafür bietet Bild 1.1 [2]

Die letzte Bemerkung zu den in diesem Kapitel behandelten Kontinuumsgleichungen leitet zugleich über in das folgende Kapitel. Die Gleichungen sind invariant unter der Vertauschung h → −h. Die Spitzen der entstehenden Mounds haben die gleiche Form wie die Gruben zwischen ihnen. In Simulationen von SOS–artigen Modellen oder in Experimenten beobachtet man jedoch eher abgerundete Gipfel und scharf eingeschnittene Täler. Solche Formen erhält man mit Gleichungen, die einen die ±h–Symmetrie brechenden Term enthalten, etwa im phänomenologischen Ansatz von J_{NG} einen Beitrag ∼ ∇(∇h)^{2}. Im nächsten Kapitel wird er mikroskopisch interpretiert. In einer Renormierungsrechnung wird er durch das Zusammenspiel von treibender Kraft F und vertikaler Gitterstruktur des Kristalls erzeugt.

92 6. STUFENBARRIEREN IM KONTINUUMSBILD

KAPITEL 7 Das Kristallgitter im Kontinuumsmodell

Es erscheint zunächst widersinnig, von einer Gitterstruktur in Kontinuumsmodellen zu sprechen. Dies ist richtig für die laterale Gitterstruktur, über die beim Übergang von diskretem zu Kontinuumsmodell gemittelt wird (siehe Kapitel 5). Die vertikale Gitterstruktur kann beim Mittelungsvorgang erhalten bleiben und so bestimmte Phänomene beschreiben. Ein Beispiel ist die Kontinuumstheorie zum Rauhigkeitsübergang, das Sinus-Gordon–Modell, das im nächsten Abschnitt in den Zusammenhang mit MBE–Wachstumsmodellen gestellt wird.

7.1. Die erhaltene Sinus-Gordon–Gleichung

Bei tiefen Temperaturen haben Kristallstufen eine große freie Energie β. Die Freie Energie eines mesoskopischen Flächenstücks der Größe l_{m} geht ein in die Oberflächenspannung (44) und das Funktional der Freien Oberflächenenergie F{h} (50). Ein in Richtung einer Hochsymmetriefläche des Kristalls orientiertes Flächenstück S_{m} durchziehen weniger Stufen, wenn es auf der Höhe einer Gitterebene h = na liegt. Andernfalls wechselt die Oberfläche zwischen benachbarten Gitterebenen hin und her, was die Stufendichte in S_{m} und auch seinen Beitrag zu F{h} erhöht. Dieser Effekt wird berücksichtigt durch einen Term

∫_{V} d^{d}x \sqrt{1 + ∇h(x, t)^{2}} V_{G}(h(x, t)),

wobei V_{G} eine Potentialfunktion mit der Periodizität des Gitters ist, die bei ganzzahligen Werten h/a ∈ ℤ Minima aufweist.

Nähert man wie in (54) die Oberflächenspannung durch einen isotropen Term an, und berücksichtigt nur die niedrigste Harmonische des Gitterpotentials, erhält man das Sinus-Gordon–Energiefunktional

F_{SG}\{h\} = ∫_{V} d^{d}x [ \frac{γ̃}{2} ∇h^{2} − V cos \frac{2π}{a} h ] . (113)

In 2+1 Dimensionen besitzt dieses Modell einen Phasenübergang in der Kosterlitz-Thouless–Universalitätsklasse, den Rauhigkeitsübergang [55, 83]. Die Höhendifferenzkoorelationsfunktion G_{2}(x, t) = ⟨(h(x, t)−h(0, t))^{2}⟩ nimmt im Limes t → ∞ oberhalb einer kritischen Temperatur T_{R} mit dem

94 7. DAS KRISTALLGITTER IM KONTINUUMSMODELL

Abstand |x| logarithmisch zu, während sie für T < T_{R} ab einer temperaturabhängigen Korrelationslänge beschränkt bleibt. Dasselbe universelle Verhalten findet sich bei einem exakt lösbaren zweidimensionalen Gittermodell für Oberflächen, dem BCSOS–Modell [84].

Vorgehend wie in Kapitel 5 beschrieben, erhält man eine Gleichung, die die Dynamik unter den Bedingungen von MBE beschreibt, mit äußerem Fluß und Erhalt der aufgetroffenen Masse, aber ohne Berücksichtigung der Stufenbarrieren,

(114) ∂_{t}h(x, t) = −KΔ^{2}h(x, t) + σV Δ sin \frac{2π}{a} h(x, t) + F (x, t).

Die Adatommobilität σ wurde isotrop angenommen, und K = σγ̃. Anders als im letzten Kapitel soll hier nicht das Endstadium von Strukturbildungsprozessen untersucht werden, weshalb in (114) die räumliche und zeitliche Fluktuation des auftreffenden Massenflusses berücksichtigt wird:

F (x, t) = F + η(x, t) (115)

Die Fluktuationen um den mittleren Fluß F sind gaußisch mit Mittelwert ⟨η⟩ = 0 und Korrelator ⟨η(k, t)η(k′, t)⟩ = 2Dθ(Λ−|k|) δ(t − t′)δ(k + k′) angesetzt. Die auftreffenden Teilchen haben endliche Ausdehnung a. Die Fluktuationen sind daher räumlich über die Ausdehnung eines Gitterplatzes korreliert, was durch den Cutoff im k–Raum Λ = 1/a berücksichtigt wird. Das Auftreffen der Teilchen kann als Poissonprozeß aufgefaßt werden, seine Fluktuationen werden auch als “Schrotrauschen” oder shot noise bezeichnet. Die Rauschstärke D ist daher proportional zu F, aus Dimensionsgründen D = a^{3}F. Selbstverständlich fluktuiert auch die Diffusionsbewegung auf der Oberfläche. Dies kann in (114) durch einen zusätzlichen Term berücksichtigt werden. Im Korrelator dieses Diffusionsrauschens steht statt dem D des Schrotrauschens ein k^{2}D_{diff}, weil es die Divergenz eines fluktuierenden Oberflächenstroms ist und die Masse auf dem Kristall erhält. Aus Dimensionsgründen ist jedoch Diffusionsrauschen im Vergleich zu Schrotrauschen irrelevant.

Es wird nun untersucht, inwiefern das Gitterpotential in vertikaler Richtung relevant ist und das Wachstumsverhalten beeinflußt. Die Methode dafür, eine störungstheoretische Renormierungsgruppenrechnung, wird im nächsten Abschnitt vorgestellt. Zwei wichtige Ergebnisse gewinnt man aus dieser Analyse: In einer mesoskopischen Gleichung, die sich ergibt, wenn Gleichung (114) über genügend viele laterale Gitterabstände gemittelt wird, also in der Gleichgung für ein “geglättetes” Höhenfeld, muß ein Term der

7.2. RENORMIERUNG NACH NOZIÈRES UND GALLET 95

Form Δ(∇h)^{2} stehen, der unter der Renormierung erzeugt wird. Im Abschnitt 7.4 wird er zudem heuristisch motiviert. Zweitens liefert diese Methode eine Vorhersage über die Längen- und Zeitskalen, die in das Abklingen des Lagenwachstums eingehen, wie im Abschnitt 4.3, insbesondere in 4.3.1.2 angedeutet ist. Heuristische Argumente und Simulationen [47] zeigen dieselben Ergebnisse wie die Interpretation der Renormierung in Abschnitt 7.3.

7.2. Renormierung nach Nozières und Gallet

7.2.1. Die Methode. Wegen der besonderen Struktur von Gleichung (114) verwenden wir eine Methode, die sich für allgemeine Formen der Nichtlinearität, nicht nur einfache Potenzen im zu untersuchenden Feld, eignet. Sie wurde von Nozières und Gallet für die dynamische Renormierung der nichterhaltenen Sinus-Gordon–Gleichung in 2+1 Dimensionen und die Theorie des Rauhigkeitsübergangs eingeführt [55, 82]. Dort findet sich auch eine detaillierte Beschreibung, die hier kurz widergegeben werden soll

Man betrachte eine Langenvingleichung

∂_{t}h = Lh + N (h) + η (116)

mit einem linearen Anteil Lh, einem nichtlinearen Term N (h) und gaußischem weissen Rauschen wie in (114) mit Cutoff Λ im Impulsraum. Die Renormierung von (114) läuft in folgenden Schritten ab:

• Über die Komponenten kurzer Wellenlänge des Rauschens δη wird gemittelt. Im k–Raum ist δη nur für Moden k mit (1−dl)Λ < |k| ≤ Λ von Null verschieden. Man definiert das gemittelte oder geglättete (coarse grained) Feld h̄ ≡ ⟨h⟩_{δη}.

• Gleichung (114) wird in zwei Anteile aufgespalten, einen für das geglättete Feld

∂_{t}h̄ = Lh̄ + ⟨N (h̄ + δh)⟩_{δη} + η̄

und einen zweiten für die Differenz δh ≡ h − h̄. Man nähert ⟨N (h̄ + δh)⟩_{δη} mit den (hoffentlich allen) relevanten Operatoren für eine Oberflächendynamik wie in Gleichung (114). Dazu berechnet man δh (bzw. seine Korrelationsfunktionen) in einer Störungsentwicklung nach N. Da N nicht von einfacher polynomialer Form ist, verwendet man eine Rayleigh-Schrödinger–Entwicklung. Sie ist die einzige unmittelbar zugängliche, wenn auch schlecht kontrollierbar.

• Zeit, laterale und vertikale Längen werden mit unterschiedlichen Exponenten reskaliert: x → (1 − dl)x, h → (1 − ζ dl)h und

96 7. DAS KRISTALLGITTER IM KONTINUUMSMODELL

t → (1 − z dl)t. Die Terme in (114) werden entsprechend reskaliert und erzeugen so die Renormierungsgruppen(RG)flußgleichungen von L und N.

Diese Schritte werden nun auf Gleichung (114) angewandt.

7.2.2. Anwendung auf die erhaltene Sinus-Gordon–Gleichung. Dazu betrachten wir (114) in einem mit der Wachstumsgeschwindigkeit F mitbewegten Bezugssystem.

∂_{t}h = −KΔ^{2}h − V̂ Δ sin \frac{2π}{2} (h − F t) + η. (117)

(V̂ ist hier eine Abkürzung für σV). Epitaktisches Wachstum beginnt auf einem atomar glatten Substrat. Die Anfangskonfiguration bei t=0 ist also h≡0. Fluktuationen werden vom Rauschen zu späteren Zeiten verursacht. Dies ist entscheidend für die Interpretation unserer Rechnung.

Mittelung über δη ergibt zwei gekoppelte Gleichungen, die bis zur jeweils niedrigsten Ordnung in der infinitesimalen Größe δh entwickelt werden.

(118)∂_{t}h̄ = −K Δ^{2}h̄ − V̂ Δ [sin \frac{2π}{a} (h̄ − F t) (1 − \frac{2π^{2}}{a^{2}} ⟨δh^{2}⟩)] + η̄
∂_{t}δh = −K Δ^{2}δh − V̂ Δ [cos \frac{2π}{a} (h̄ − F t) \frac{2π}{a} δh] + δη.

Die zweite Gleichung wird durch einen Störungsansatz gelöst

(119) δh^{(0)}(x, t) = ∫ d^{d}x′ ∫_{0}^{t} dt′ Γ(x − x′, t − t′) δη(x′, t′)
δh^{(1)}(x, t) = −\frac{2πV̂}{a} ∫ d^{d}x′ ∫_{0}^{t} dt′ Γ(x − x′, t − t′) ×

Δ′ [cos \frac{2π}{a} (h̄(x′, t′) − F t′) δh^{(0)}(x′, t′)] .

Im k–Raum ist der lineare Propagator Γ = exp −K(t − t′)|k|^{4}. Im folgenden bezeichnen h̄′, η′, … Größen, deren Argument die Integrationsvariable (x′, t′) ist, und Symbole ohne ′ die Größen bei (x, t). Ableitungen (Laplace-Operatoren) werden mit Δ und Δ′ bezeichnet, um diesem Unterschied Rechnung zu tragen. Das Faltungsintegral ∫ d^{d}x′ ∫_{0}^{t} dt′ wird mit ∫ abgekürzt.

Aus ⟨δh^{2}⟩=⟨δh^{(0)}^{2}⟩ + 2⟨δh^{(0)}δh^{(1)}⟩ + O(V̂^{2}) in die erste Gleichung von (118) eingesetzt ergeben sich die Korrekturterme zu h̄ zur Ordnung V̂^{2}.

⟨δh^{2}⟩ = ⟨δh^{(0)}^{2}⟩ − \frac{4πV̂}{a} ∫ (Δ′Γ) cos \frac{2π}{a} (h̄′ − F t′) C, (120)

7.2. RENORMIERUNG NACH NOZIÈRES UND GALLET 97

wobei C ≡ ⟨δh^{(0)}δh^{(0)}′⟩ die ungestörte Korrelationsfunktion der auszuintegrierenden Moden bezeichnet. Ihre Fouriertransformierte ist e^{−K(t−t′)|k|^{4}}D/(K|k|^{4}) [1 − e^{−2Kt|k|^{4}}].

Der Term niedrigster Ordnung ⟨δh^{(0)}^{2}⟩ ist eine Konstante, mit der die sinusförmige Kraft multipliziert wird, und liefert die Korrektur

δV̂ = −\frac{4π^{3}D}{Ka^{2}Λ^{2}} (1 − e^{−2KΛ^{4}t}) dl V̂ . (121)

Diese Gleichung gilt für den physikalisch relevanten Fall d=2 einer Kristalloberfläche. Allgemein stünde in d+1 Dimensionen statt 2π/Λ^{2} der Term S_{d}Λ^{d−4}, wobei S_{d} das Oberflächenmaß der d-dimensionalen Einheitssphäre bezeichnet.

Zur Auswertung des Beitrags in nächsthöherer Ordnung, nämlich von ⟨δh^{(0)}δh^{(1)}⟩,

−\frac{8π^{3}V̂^{2}}{a^{3}} Δ ∫ sin \frac{2π}{a} (h̄−F t) cos \frac{2π}{a} (h̄′−F t′) (Δ′Γ) C,

spalten wir das Produkt sin a cos b = (sin(a + b) + sin(a − b))/2 auf und vernachlässigen den Term mit sin(a+b). Er erzeugt höhere Harmonische der periodischen Kraft. Gleichung (121) zeigt, daß sie weniger relevant sind als die Grundfrequenz. Im verbleibenden Ausdruck unterteilen wir in Terme mit h̄ und mit F t

−\frac{4π^{3}V̂^{2}}{a^{3}} Δ ∫ [sin \frac{2π}{a} F (t′ − t) cos \frac{2π}{a} (h̄ − h̄′) (Δ′Γ) C
+ cos \frac{2π}{a} F (t′ − t) sin \frac{2π}{a} (h̄ − h̄′) (Δ′Γ) C].

Eine Entwicklung von Termen mit h̄−h̄′ nach Potenzen von ξ := x′−x unter Berücksichtigung der Symmetrie bei Integration über d^{2}ξ ergibt

=  \frac{4π^{5}V̂^{2}}{a^{5}} [Δ(∇h̄)^{2}] ∫ sin [\frac{2π}{a} F (t′ − t)] (Δ′Γ) C ξ^{2}
(122)  − \frac{2π^{4}V̂^{2}}{a^{4}} [Δ^{2}h̄] ∫ cos [\frac{2π}{a} F (t′ − t)] (Δ′Γ) C ξ^{2}.

Der zweite Term liefert eine Korrektur δK zum Koeffizienten des Relaxationsterms, der erste erzeugt einen nichtlinearen Term

−λΔ(∇h(x, t))^{2}, (123)

der zu Gleichung (114) hinzugefügt werden muß. Damit haben wir den gemittelten Anteil der Nichtlinearität ⟨N (h)⟩_{δη} bis zur Ordnung V̂^{2} entwickelt in Projektionen auf die relevanten Operatoren des Problems.

98 7. DAS KRISTALLGITTER IM KONTINUUMSMODELL

Der letzte Schritt besteht in der Reskalierung von Raum und Zeit. Um das Verhalten nah an den linearen Fixpunkten (K, V̂ = 0, λ = 0) zu untersuchen, wählen wir deren Skalenexponenten, in 2+1 Dimensionen z =4 und ζ = 2−d/2 = 1. Die treibende Kraft (bzw. Wachstumsgeschwindigkeit) F und die Gitterkonstante a werden ebenfalls reskaliert: dF/dl =(z−ζ)F und da/dl =−ζa.

Die folgenden “natürlichen” Einheiten eignen sich besser zur Auswertung des RG-Flusses. Die Zeit t wird ausgedrückt durch die Bedeckung oder Anzahl aufgedampfter Monolagen θ:=F t/a. Dies führt zu einem dimensionslosen Ausdruck für K, nämlich K:=K/(F a^{3})=Kt(θ=1)/a^{4}. Die Größe Kt(θ = 1) ist die vierte Potenz einer Länge, über die sich Korrelationen lateral ausbreiten während des Wachstums einer Lage. Beim schichtweisen Wachstum ist diese Länge der Abstand zweier benachbarter Nukleationszentren, d.h. K = (l_{D}/a)^{4}. Die Rauschstärke D wird nicht reskaliert und auch in einem System mit erhaltener Dynamik nicht renormiert. Man benutzt daher das dimensionslose Verhältnis der Stärke der Gitterkraft zum Rauschen U := a_{0}V̂ /D. Da für d = 2 der Rauhigkeitsexponent ζ = 1 ist, wächst der nichtlineare Operator Δ(∇h)^{2} beim Reskalieren an, was im dimensionslosen Produkt L = λ/(a^{2}F) ausgeglichen wird. Als Skalenparameter verwendet man κ := exp 4l = (A/a)^{4}, wobei A die Cutoff länge ist, unterhalb der alle Fluktuationen ausgeglättet wurden. Damit erhält man den Renormierungsgruppenfluß

(i) 4κ \frac{dU}{dκ} = [1 − 4π^{3} \frac{\sqrt{κ}}{K} (1 − e^{−2Kθ/κ})] U,
(124)  (ii) 4κ \frac{dK}{dκ} = \frac{U^{2}}{K} g(K/κ, θ),
(iii) 4κ \frac{dL}{dκ} = \frac{U^{2}}{K} f(K/κ, θ).

Die Funktionen f und g ergeben sich aus den Integralen in (122) und lauten

g(K/κ, θ) = 8π^{5} \frac{1 − e^{−2Kθ/κ}}{[(K/κ)^{2} + π^{2}]^{3}}
[K/κ(−(K/κ)^{4} − 4π^{2}(K/κ)^{2} + 5π^{4})
+(B cos(2πθ) + A sin(2πθ)) e^{−2Kθ/κ}]

f(K/κ, θ) = 16π^{6} \frac{1 − e^{−2Kθ/κ}}{[(K/κ)^{2} + π^{2}]^{3}}

[π(−(K/κ)^{4} − 8π^{2}(K/κ)^{2} + π^{4})

+(−A cos(2πθ) + B sin(2πθ)) e^{−2Kθ/κ}]

mit den Polynomen

A = 4πθ^{2}(K/κ)^{6} + (π − 8π^{3}θ^{2})(K/κ)^{4}
+8π^{3}θ(K/κ)^{3} + (8π^{3} − 4π^{5}θ^{2})(K/κ)^{2}
+8π^{5}θK/κ − π^{5}
B = −4θ^{2}(K/κ)^{7} + 4θ(K/κ)^{6} + (1 − 8π^{2}θ^{2})(K/κ)^{5}
+16π^{2}θ(K/κ)^{4} + 4(π^{2} − π^{4}θ^{2})(K/κ)^{3} + 12π^{4}θ(K/κ)^{2}
−5π^{4}K/κ.

Den Beitrag des zur Ordnung V̂^{2} erzeugten λ, d.h. die Rückkopplung von L auf die übrigen Parameter wird vernachlässigt. Wir werden im nächsten Abschnitt sehen, daß bei einer getriebenen Oberfläche eine Längenskala l̃ existiert, ab der das Kristallgitter keinen Einfluß auf die übrigen Parameter hat. Eine Nichtlinearität vom Typ (123) jedoch wird erst auf wesentlich größeren Skalen relevant, führt dann allerdings zu einem neuen Fixpunkt mit neuen Skalenexponenten [85]. Zunächst, und so weit, wie hier die Renormierung betrachtet werden soll, ist er von der Ordnung L ∼ U^{2}, seine Rückkopplung von der Ordnung L^{2} kann also gegen den Einfluß des Gitters vernachlässigt werden.

7.3. Dämpfung der Wachstumsoszillationen

Die Entwicklung von U unter dem Renormierungsgruppenfluß (124) entscheidet, ob die Oberfläche deutlich auf eine Gitterebene fixiert (“gepinnt”) ist, oder aber auf großen Skalen der Einfluß des Kristallgitters auf die Morphologie irrelevant ist. Die RG-Analyse des Rauhigkeitsübergangs benutzt dasselbe Prinzip [83, 55]. Oberhalb T_{R} verschwindet die Stärke des Pinningpotentials asymptotisch, unterhalb divergiert sie. Zunächst werten wir das Ergebnis (124) im Limes unendlicher Wachstumsdauer θ → ∞ aus, danach berücksichtigen wir, daß zu Zeiten t < ∞ der Einfluß der Fluktuationen nur über eine endliche Zeitspanne bestanden hat.

100 7. DAS KRISTALLGITTER IM KONTINUUMSMODELL

Die Flußgleichung für U in (124) unterscheidet sich in einem Punkt wesentlich von der nichterhaltenen Sinus-Gordon–Gleichung in 2 + 1 Dimensionen, die in [83, 55] untersucht wird. Sie ist nicht autonom, die Skalenvariable κ steht explizit in (124,i) und treibt asymptotisch U gegen Null. Im allgemeinen Fall stünde statt \sqrt{κ} die Potenz κ^{ζ/2}, wobei der Rauhigkeitsexponent ζ von der Dimension, der Art des Rauschens und des linearen Relaxationsmechanismus abhängt (siehe oben). Die einzige Größe, die im Korrekturterm δV̂ /V̂ für t → ∞ in Gleichung (121) skalenabhängig ist, ist die Gitterkonstante a = a_{0}κ^{−ζ/4}. Wenn ζ < 0, wächst U asymptotisch und die Oberfläche ist makroskopisch flach, wenn ζ >0, dominiert der Term in (124,i), der U vermindert, wie im hier untersuchten Beispiel. Eine allgemeine Bedingung für einen Rauhigkeitsübergang ist daher ζ = 0, wie es z.B. beim Sinus-Gordon–Modell in [83, 55] der Fall ist.

Unter den für MBE typischen Bedingungen ist die Depositionsrate F deutlich kleiner als die Diffusionsrate, in den Parametern von Gleichung (124) bedeutet dies K ≫ 1. Die relative Änderung von K ist also klein, d.h. in sehr guter Näherung nimmt U unter Renormierung zu bis zur Skala κ_{max} = K^{2}/(16π^{6}) (bzw. κ_{max} ∼ K^{2/ζ} im oben angedeuteten allgemeineren Fall), um dann sehr schnell auf Null abzufallen. Dies hat eine wichtige Konsequenz für Substrate endlicher Ausdehnung und führt zum Begriff der Schichtenkohärenzlänge (layer coherence length) l̃, die in [47] eingeführt wird. Auf Substraten, die kleiner sind als

l̃= a_{0} κ_{max}^{1/4} = a_{0} \frac{\sqrt{K}}{2π^{3/2}}

(oder l̃∼ a_{0}K^{1/2ζ} ∼ a_{0}K^{1/(4−d)}) bleibt der Einfluß des Gitters relevant. Die Oberfläche bleibt glatt und wächst lagenweise. Das RHEED–Signal, das die mittlere Breite wiedergibt, oszilliert beliebig lang.

Die Kohärenzlänge l̃ist aus Dimensionsüberlegungen mit einer Zeit oder Bedeckung verbunden θ̃ ≡ κ_{max}/K, die wir als Abklingzeit der Oszillationen erkennen werden. Dazu betrachtet man nun endliche Bedeckungen θ. Höhen an Punkten, die weiter als l_{θ} = a_{0}(Kθ)^{1/4} voneinander entfernt sind, sind nicht korreliert. Integriert man Fluktuationen jenseits der Skala l_{θ} aus, ergibt sich kein Renormierungsbeitrag außer trivialer Reskalierung der Parameter. Siehe dazu Gleichung (124,i): Für κ ≫ Kθ lautet sie 4κ dU/dκ = [1 − O(κ^{−1/2})] U ≃ U und U wird nur noch reskaliert wie κ^{1/4}. Nichttriviale Beiträge verschwinden.

Ist nun l_{θ} < l̃, so endet die Renormierung von Ubei großen Werten von U, also starker Pinningkraft, andernfalls ist der Gittereinfluß schwach und die Oberfläche rauh. Anders ausgedrückt: Bis zur Zeit (Bedeckung)

θ̃ = \frac{κ_{max}}{K} = \frac{1}{16π^{6}} K

(oder θ̃ ∼ K^{(2−ζ)/ζ} ∼ K^{d/(4−d)}) wächst die Oberfläche schichtweise und geht dann über zu rauhem Wachstum, in Übereinstimmung mit den Ergebnissen von [47].

Abbildung 7.1. Rastertunnelmikroskopische Aufnahmen einer wachsenden Fe(100)–Oberfläche bei verschiedenen Temperaturen [86]. Die Intensität des RHEED–Signals zeigt Oszillationen mit Maxima bei möglichst glatten Flächen. Niedrige Oberflächentemperatur bewirkt geringe Diffusionsrate, kleine Diffusionslänge l_{D} und ein schnelles Abklingen der Oszillationen (Bild a). Bei der höchsten Temperatur (c) ist das Lagenwachstum nahezu perfekt. Alle Rasterbilder wurden bei einer Bedeckung von fünf Monolagen aufgenommen (Pfeile).

102 7. DAS KRISTALLGITTER IM KONTINUUMSMODELL

7.4. Die erhaltene KPZ–Nichtlinearität

In den Korrekturtermen (122) steht auch ein Beitrag der Form ∼ −Δ(∇h)^{2}. Die heuristische Begründung eines solchen Terms ist leicht einzusehen. Die Adatomdichte ρ auf den Terrassen ist eine Funktion ihrer Ausdehnung und damit der Größe der Neigung (∇h)^{2}. Auf jeder Terrasse erfüllt ρ Gleichung (63). Es ergibt sich ein Dichteprofil mit konstanter mittlerer Krümmung F/σ. Die Adatomdichte auf einer Terrasse der linearen Ausdehnung l ist daher von der Größenordnung F l^{2}/σ.

Für große Neigungen ist l = a/|∇h|, während es für kleine Neigungen durch die Diffusionslänge l_{D} begrenzt ist. Eine interpolierende Formel ist

l^{2} = \frac{l_{D}}{1 + (l_{D}/a)^{2}(∇h)^{2}} .

Für kleine Neigungen entwickelt man ρ nach Potenzen von (∇h)^{2}:

ρ ≈ \frac{l^{2}F}{σ} ≈ \frac{F l_{D}^{2}}{σ} (1 − \frac{l_{D}^{2}}{a^{2}} (∇h)^{2} + …) .

Variationen in der Adatomdichte führen zu einem mesoskopischen Oberflächenstrom j = −σ∇ρ. (Wie bei der Kontinuumsgleichung (114) wird auch in dieser heuristischen Betrachtung der Effekt von Stufenbarrieren vernachlässigt.) Die Divergenz von Oberflächenströmen bedeutet eine lokale Änderung des Höhenfelds h, in diesem Fall erhält ∂_{t}h den zusätzlichen Beitrag −λ Δ(∇h)^{2} mit

λ = σ \frac{F l_{D}^{4}}{σa^{2}} ∼ F (\frac{σ}{F} )^{2/(d+1)} , (125)

wobei die Beziehung (40) zwischen l_{D} und dem Verhältnis von Diffusionszu Aufdampfrate σ/F eingesetzt wurde.

Der Koeffizient λ gehört zur Nichtgleichgewichtseigenschaft des Systems, er ist proportional zu F und verschwindet bei F = 0. Ohne den Fluß auftreffender Teilchen sind die diffundierenden Adatome im Gleichgewicht mit den Stufen und der unterliegenden Gitterebene. Es besteht keine “Quelle” für Adatome auf den Terrassenflächen, weshalb ihr Dichteprofil keine mittlere Krümmung aufweist. Die Abhängigkeit von der Terrassenausdehnung l und von (∇h)^{2} geht verloren.

An dieser Stelle soll auch auf die nichterhaltene Sinus-Gordon–Gleich-ung verwiesen werden. Ihre mittlere Wachstumsgeschwindigkeit ist neigungsabhängig [55]. In einer mesoskopischen Gleichung muß also eine lokale Wachstumsgeschwindigkeit v(|∇h|) berücksichtigt werden. In der Entwicklung für kleine Neigungen ergibt sich daraus ein Term der Form λ(∇h)^{2}

7.4. DIE ERHALTENE KPZ–NICHTLINEARITÄT 103

[6]. Unter einer der Rechnung in diesem Kapitel vergleichbaren Renormierung wird in der nicht erhaltenen Sinus-Gordon–Gleichung ein solcher Term erzeugt [56].

Im hier untersuchten Parameterbereich (K ≫ 1, dessen relative Änderungen vernachlässigt werden) findet man bei Auswertung des RG–Flusses (124) die Proportionalität zwischen λ und der Aufdampfrate, die soeben heuristisch begründet wurde. Mit L nimmt unter dem Fluß (124) das Verhältnis λ/F einen konstanten Grenzwert an. F ist proportional zu F_{0}, also ist auch λ∝F_{0}.

Von der Stärke des Gitterpotentials hängt λ quadratisch ab. Gleichung (124i) ist linear in U. Ist die Rückkopplung von U auf K vernachlässigbar, so ist daher das renormierte U proportional zu U_{0}. Ferner wird L von 0 ausgehend erzeugt, es ist also L ∝ U^{2} ∝ U_{0}^{2} und damit

λ ∝ F_{0}U_{0}^{2}.

Schwieriger ist die Abhängigkeit des renormierten λ bzw. des Grenzwerts von L für κ → ∞ von K. Die heuristischen Überlegungen zu Beginn dieses Abschnitts gelten für einen stationären Zustand, wenn die Oberfläche bereits weit gewachsen ist. Gleichung (125) ließe daher im Limes θ → ∞ erwarten, daß L ∼ K. Als Ergebnis einer (numerischen) Integration der Flußgleichungen ist L in Abhängigkeit von K in Abbildung 7.2 dargestellt. Die Integration wurde abgebrochen beim im letzten Abschnitt (7.3) beschriebenen Abfall des Gitterpotentials U. Die Stärke der Nichtlinearität L nimmt für κ → ∞ einen Grenzwert an, weil die rechte Seite von (124(iii)) hinreichend schnell verschwindet (sowohl für endliche θ als auch im hier vorher genommenen Limes θ → ∞).

Quantitativ entspricht das Ergebnis nicht der Vorhersage von Gleichung (125). Zwar nimmt L mit K ∼ l_{d}^{4} zu, aber nicht proportional, sondern nähert sich einer Konstanten an. Die RG–Analyse beruht auf einer Störungsentwicklung für schwaches Pinning an den Gitterebenen. Das Konzept der Adatomdichte auf Terrassenstücken, wie es in obige Herleitung einging, benötigt ein starkes Pinning, so daß Terrassen deutlich hervortreten. Dies ist jedoch mit dieser und vergleichbaren Methoden nicht zugänglich. Bemerkenswert ist jedoch die Gegenwart eines erhaltenen KPZ–Terms ∼ −Δ(∇h)^{2} auch in diesem stark vereinfachenden Ansatz. Sie deutet darauf hin, daß er bei erhaltener Oberflächendynamik eine ebenso universelle Rolle spielt wie der Term (∇h)^{2} in der Dynamik ohne Massenerhaltung [6].

Abbildung 7.2. Abhängigkeit des renormierten Parameters L = λ/(a^{2}F) von der Stärke des linearen Terms K = (l_{D}/a)^{4} in Gleichung (124). Die Stärke des nichtlinearen Terms −Δ(∇h)^{2} nimmt mit l_{D} zu, wenn auch nicht so wie nach Gl. (125) quantitativ erwartet.

KAPITEL 8 Zusammenfassung

Dynamik instabiler Oberflächen, dieser Themenbereich wurde in dieser Arbeit unter mehreren Gesichtspunkten und in verschiedenen Zusammenhängen behandelt.

Den physikalischen Hintergrund bildet die Morphologie wachsender Kristalloberflächen. Den Theorien zugrunde liegt eine nahezu unermeßliche Vielfalt an Experimenten, sowohl in den Techniken des Wachstums selbst und ihren Verfeinerungen, als auch in den Methoden zur Überwachung und Steuerung. Instabilitäten theoretisch zu verstehen dient zwei gegensätzlichen Zwecken: Beabsichtigt man stabiles Wachstum glatter Flächen, geht es darum, sie zu vermeiden. Umgekehrt besteht die Hoffnung, die Formen instabilen Wachstums gezielt zu nutzen und so mikroskopische Oberflächenstrukturen zu erhalten, sogenannte Quantendrähte und -punkte.

Für sich allein betrachtet stellt die theoretische Modellierung von Oberflächendynamiken ein weit entwickeltes Forschungsgebiet dar. Zwei ihrer Hauptbereiche werden berührt: Statistisch physikalische Methoden werden eingeführt, wie sie im Zusammenhang mit Problemen der kinetischen Aufrauhung entwickelt wurden. Angewandt werden sie zunächst in einem rein theoretischen Rahmen auf deterministische (chaotische) Gleichungen in ihrem Verhalten auf großen Skalen.

Ziel dieser Arbeit ist es, den experimentellen und den theoretischen Ansatzpunkt miteinander zu verbinden und in diesem Zusammenhang neue Forschungsergebnisse darzustellen. Sehr viel Wert gelegt wird auf die Herleitung von Kontinuumsgleichungen für Oberflächen wachsender Kristalle. Die allgemeine Struktur von Kontinuumsgleichungen erhält so einen anschaulichen Hintergrund. Der Herleitung voraus gehen eine Übersicht über die experimentelle Methode der Molekularstrahlepitaxie und die Phänomene instabilen Wachstums mit Stufenkantenbarrieren, sowie über die oben angesprochenen theoretischen Gesichtspunkte. Beide Darstellungen verbinden sich in den nachfolgenden Kapiteln in drei detaillierten Rechnungen: Behandelt werden in ihnen die Vergröberung von Oberflächenstrukturen

106 8. ZUSAMMENFASSUNG

(coarsening) auf singulären – entlang einer Hochsymmetrieebene des Kristalls – und vizinalen, relativ zu einer Symmetrierichtung verkippten Oberflächen. Die dritte Rechnung betrifft ein früheres Stadium epitaktischen Wachstums, den Übergang von schichtweisem zu rauhem Wachstum.

Dabei besteht der Anspruch, die Möglichkeiten und Grenzen deterministischer Kontinuumsgleichungen im Hinblick auf Strukturbildung an wachsenden Kristallen umfassend darzustellen. Insbesondere die “Schwächen” des Kontinuumsansatzes weisen auf weitere offene Fragen hin: Deuten die unterschiedlichen in Experimenten und Simulationen beobachteten Wachstumsgesetze nicht auf einen inneren Widerspruch der Skalenhypothese hin? Generell bildet Rauschen einen Mechanismus, der Skaleninvarianz erzeugt [12] (siehe auch Kapitel 2). Bei Molekularstrahlepitaxie liegen die Fluktuationen im auf die Oberfläche treffenden Teilchenstrom, dem sogenannten Schrotrauschen. Ungeklärt ist, ob sie das vom deterministischen Anteil der Dynamik erzeugte Skalenverhalten (Kapitel 6) überlagert, also fundamental ändert, oder nur transiente Korrekturen erzeugt.

Diese Fragestellungen gehen über das Gebiet der Oberflächendynamik hinaus. Wegen der Verwandtschaft dieser Gleichungen zu Problemen der Phasenordnung in magnetischen Systemen dürften dort Abweichungen von reinem Skalenverhalten und Einfluß von Rauschen eine ähnliche Rolle spielen. Wie bei der Behandlung rein deterministischer Gleichungen in Kapitel 6 darf man erwarten, daß es diese Verwandtschaft auch für verrauschte Gleichungen erlaubt, Ergebnisse zwischen beiden Problemklassen zu übertragen, ihre Unterschiede zu vergleichen und so zu verstehen.

ANHANG A Strukturfaktor in sphärischer Näherung

In den Gleichungen (98) und (99) besteht eine wechselseitige Abhängigkeit von a(t) und S(k, t): a ist ein Funktional von S, geht aber in dessen Berechnung ein. Eine Lösung des Strukturfaktors muß also einer Selbstkonsistenzbedingung genügen, die sich in folgende Form bringen läßt

\frac{db}{dt} = f (2/(2^{d/2}Γ(d/2)) ∫_{0}^{∞}dk k^{d+1} exp (−2tk^{4} + 2k^{2}b(t))) . (126)

Die Anfangsbedingung S_{0} ist gewählt wie im Text motiviert, und die neu eingeführte Funktion ist definiert als b(t) = ∫_{0}^{t} ds f(a(s)). Führt man das Integral aus, ergibt sich für b(t) die Differentialgleichung

\frac{db}{dt} = f (\frac{d}{2}(4t)^{−(d+2)/4} 2^{−d/2} D_{− \frac{d+2}{2}} (−b/\sqrt{t}) exp \frac{b^{2}}{4t}) , (127)

in der D eine parabolische Zylinderfunktion bezeichnet [87]. Gleichung (127) ist nicht explizit lösbar, aber für späte Zeiten divergiert b/\sqrt{t}, wie wir gleich zeigen werden. D läßt sich daher ausgehend von einer Sattelpunktsnäherung des Integrals in (126) annähern.

Gleichung (127) hat die Form

\frac{db}{dt} = f (t^{−(d+2)/4} F (\frac{b}{\sqrt{t}})) . (128)

Bliebe b/\sqrt{t} beschränkt, würde das Argument von f in (127) für große t gegen 0 gehen. Die Differentialgleichung für b(t) wäre damit db/dt ≃ f(0) = 1 im Widerspruch zur Beschränktheit von b/\sqrt{t}.

Die Näherung von (127) ist damit für große t (und b/\sqrt{t})

\frac{db}{dt} = f (\sqrt{\frac{π}{2^{d−1}}} (4t)^{−(d+2)/4}/Γ(d/2) (\frac{b}{\sqrt{t}})^{d/2} exp \frac{b^{2}}{2t}) , (129)

woraus sich ergibt, daß b/\sqrt{t} langsamer als jede Potenz in der Zeit zunimmt: Gälte b ∼ t^{1/2+ϵ}, folgte daraus db/dt ∼ t^{ϵ−1/2}. Andererseits nähme das Argument in f exponentiell mit einem dominierenden Term exp t^{2ϵ} zu.

Je nach Wahl der Form des Nichtgleichgewichtsstroms f ergibt sich ein unterschiedliche Stärke des führenden Logarithmus beim Anwachsen von B = b^{2}/(2t). Für Modell (ii) ist Gleichung (127)

t \frac{dB}{dt} = −B + \sqrt{2Bt} (1 − C_{(ii)} B^{d/4} t^{−(d+2)/4} exp B) , (130)

wobei C_{(ii)} eine beliebige Konstante bezeichnet. Keiner der Terme darf mit einer Potenz der Zeit anwachsen. Dies erfordert, daß der Term in Klammern asymptotisch verschwindet, und damit exp B ∼ t^{(d+2)/4}. Der führende Anteil von B ist daher

B ≃ \frac{d + 2}{4} log t. (131)

Für Modell (i) nützt man das asymptotische Verhalten f(s^{2}) ≃ (s^{2})^{−γ/2} bei großen s^{2} aus. Gleichung (127) wird zu

t \frac{dB}{dt} = −B + C_{(i)} t^{(γ/2)(d+2)/4+1/2} B^{−(γ/2)d/4+1/2} exp −(γ/2)B. (132)

Wie oben müssen sich auch hier die Potenzen von t im letzten Term aufheben, was nur unter der Bedingung

B ≃ (\frac{d + 2}{4} + \frac{1}{γ} ) log t. (133)

erfüllt ist.

Die Ergebnisse von Modell (i) und (ii) weisen eine Gemeinsamkeit auf, die sich mit den einführenden Überlegungen zum Nichtgleichgewichtsstrom ergänzt. Wird γ vergrößert, so fällt der Strom mit zunehmender Neigung stärker ab. Man möchte daher erwarten, daß der Limes γ → ∞ einer Nullstelle der Stromfunktion bei endlicher Neigung entspricht, also einer “magic slope”. Tatsächlich trifft dies auf die hier erhaltenen Ergebnisse der sphärischen Näherung zu: Der Term 1/γ in (133) trägt in diesem Grenzfall nichts bei, und man erhält das Ergebnis von Modell (ii) (131). Dasselbe trifft auch auf die Skalenrelationen (93,94) zu, nämlich z ≥4, β ≤1/4 und α = 0.

ANHANG B Allgemeine Erläuterungen

Einige Begriffe und Methoden sollen gesondert erläutert werden, weil auf sie in allen Kapiteln bezug genommen wird.

B.1. Verwendete Begriffe

Manche der in dieser Arbeit verwendeten Bezeichnungen haben sich im engeren Gebiet der statistischen Physik wachsender Oberflächen herausgebildet. Zum besseren Verständnis sollen sie hier definiert werden.

Unter Oberflächen werden Hyperflächen verstanden, also in einen (d+ 1)–dimensionalen Raum eingebettete d–dimensionale Objekte. Um auf die Eigendimension und die zusätzliche Kodimension aufmerksam zu machen, spricht man oft von d+1 Dimensionen.

Die Lage einer Oberfläche im Raum zur Zeit t, oft durch eine Funktion h(·, t) : ℝ^{d} ∋ (x, t) → h(x, t) ∈ ℝ ausgedrückt, bezeichnen wir als Höhenkonfiguration oder -profil.

Ein Zustand ist eine Wahrscheinlichkeitsverteilung auf dem Raum der Konfigurationen.

Uneinheitlich ist in der Literatur die Benennung der Wachstumsformen auf Oberflächen, oft gibt es nur die englischen Bezeichnungen, die hier “eigenmächtig” übersetzt sind:

Im Frühstadium des Wachstums auf einer atomar glatten Fläche entstehen zunächst Inseln aus zusammengelagerten Teilchen in der ersten Schicht über der Ausgangsfläche. Kleine Gebilde aus mehreren übereinanderliegenden Inseln werden wegen ihrer Form pyramids oder wedding cakes (Hochzeitstorten) genannt, in dieser Arbeit Pyramiden. Im späteren Stadium entstehen größere Buckel, deutlich größer als die in den Pyramiden noch sichtbare atomare Skala. Man nennt sie mounds, es gibt auch den Namen SLUG, super large unstable growth feature^{1}, vor allem für langgezogene Gebilde – wir wählen die Bezeichnung Hügel. Eine langgezogene Wellenoder Rinnenstruktur auf einer Oberfläche besteht aus sogenannten Rippen (englisch ripples).

1“Nacktschnecke” auf deutsch.

110 B. ALLGEMEINE ERLÄUTERUNGEN

Die Bezeichnung der Parameter in Wachstumsmodellen ist jeweils im Text definiert. Auf doppelte Verwendung einiger Benennungen, insbesondere zwischen den Kapiteln 6 und 7 sei hier hingewiesen, auch wenn sie sich im jeweiligen Zusammenhang selbst erklären. In mikroskopischen Modellen und bei Monte-Carlo–Simulationen bezeichnet man die Dynamik mit den Raten bestimmter Ereignisse: D ist die Sprungrate der Adatomdiffusion, F ist die Rate, mit der Teilchen auf jedem Gitterplatz landen. In den Kontinuumsmodellen bezeichnet F die mittlere Wachstumsgeschwindigkeit der Oberfläche, hat also die Dimension Länge/Zeit. D bezeichnet üblicherweise und so auch hier die Stärke einer zufälligen Rauschkraft.

Die in den Formeln verwendeten Bezeichnungen sind in nochmals in einer Tabelle aufgeführt:

In Kontinuumsmodellen von Oberflächen treten auf:

h(x, t)

Höhenprofil als Funktion von Substratkoordinate und Zeit.

η(x, t)

Rauschkraft, mit Stärke D.

F{h}

Freie Energie, Funktional der Konfiguration.

γ(∇h)

neigungsabhängige Oberflächenspannung.

γ̃

Steifigkeit der Oberfläche.

σ, σ

Adatommobilität, als tensorielle bzw. skalare Größe.

σ_{⊥}

Mobilität der Oberfläche gegen senkrechte Verschiebung.

μ(x, t)

chemisches Potential für Adatome.

J Oberflächenstrom. Definitionen für Gittermodelle von Oberflächen:

a

Gitterkonstante.

F

Aufdampfrate.

D

Sprungrate der Adatome.

l_{D}

Diffusionslänge.

l_{+} und l_{−} von Stufenbarrieren induzierte Längenskala.

B.2. NUMERIK 111

Zur Auswertung der statistischen Eigenschaften von Oberflächenkonfigurationen:

G(x, t)

Zweipunkt-Höhenkorrelationsfunktion ⟨h(x, t)h(0, t)⟩.

S(k, t)

Strukturfaktor, Fouriertransformierte von G.

G_{2}(x, t)

Korrelationsfunktion der Höhendifferenzen

⟨(h(x, t)−h(0, t))^{2}⟩.

ξ, ξ_{∥}, ξ_{⊥} laterale Längenskala, parallel und senkrecht zu einer ausgezeichneten Richtung.

w,s

Breite der Oberfläche und Steilheit .

z, ζ

dynamischer und statischer Rauhigkeitsexponent.

β, α Skalenexponenten w ∼t^{β}, s∼t^{α}.

B.2. Numerik

Kurz skizziert sei ebenfalls die numerische Methode zur Integration partieller Differentialgleichungen in 2+1 Dimensionen, die für verschiedenen Modelle in dieser Arbeit verwendet wurde. Es handelt sich um Anfangswertprobleme mit einer Dynamik erster Ordnung in der Zeit, die in allgemeinster Form folgende Gestalt haben

∂_{t}h = Lh + N [h]. (134)

Die rechte Seite besteht aus einem linearen und einem nichtlinearen Differentialoperator, die getrennt behandelt wurden. Für die Anfangskonfiguration wurde an jedem Gitterpunkt x= (i, j) unabhängig eine Höhe h_{0}(x) aus dem Intervall [−ϵ/2, ϵ/2) mit konstanter Wahrscheinlichkeitsdichte ausgewürfelt. Meist wurde ϵ=1/1000 gewählt.

Der lineare Operator besitzt hohe (vierte) Ableitungen, so daß zur besseren Stabilität ein semiimplizites Crank-Nicholson–Schema verwendet wurde. Ausgehend von einer Konfiguration h_{t} ergibt sich die Konfiguration zum Zeitpunkt t+τ zu

h_{t+τ}^{L} = (1 − \frac{τ}{2}L)^{−1} (1 − \frac{τ}{2}L) h_{t}

In Fourierdarstellung ist der Differentialoperator L diagonal, so daß h_{t+τ}^{L} mit Hilfe einer schnellen Fouriertransformation (FFT) einer Bibliotheksroutine (NAG) berechnet wurde. Es wurden daher Gitter mit 2^{n} × 2^{n}, n = 5, …, 10 Gitterpunkten verwendet.

Für den nichtlinearen Anteil wurde aus Zeitgründen der Integrationsschritt explizit berechnet,

h + τN [h],

112 B. ALLGEMEINE ERLÄUTERUNGEN

wobei die Schrittweite τ angepaßt wurde, um die Größe des zweiten Terms zu kontrollieren.

Die beiden Operatoren L und N wurden nach zwei Schemata miteinander verknüpft, einmal durch direkte Addition beider Beiträge ausgehend von h_{t}:

h_{t+τ} = h_{t+τ}^{L} + τN [h_{t}],

zum anderen, indem beide Schritte hintereinander ausgeführt wurden (method of fractional steps, [88]):

h_{t+τ} = h_{t+τ}^{L} + τN [h_{t+τ}^{L}].

Für hinreichend kleine Zeitschritte ergaben sich keine merklichen Abweichungen zwischen beiden Methoden.

Literaturverzeichnis

[1] T. Mann, Doktor Faustus, S. 26, Fischer Taschenbuch (1993).

[2] S. Esch, T. Michely und G. Comsa, unveröffentlicht.

[3] A.J. Bray, Advances in Physics, 43, 357 (1994).

[4] P.C. Hohenberg und B.I. Halperin, Rev. Mod. Phys. 49, 435 (1977).

[5] S.F. Edwards und D.R. Wilkinson, Proc. Roy. Soc. London A 381, 17 (1982).

[6] M. Kardar, G. Parisi und Y.C. Zhang, Phys. Rev. Lett. 56, 889 (1986); E. Medina, T. Hwa, M. Kardar und Y.C. Zhang, Phys. Rev. A 39, 3053 (1989).

[7] Y. Kuramoto, Chemical Oscillations Waves and Turbulence, Springer, Berlin (1984).

[8] G.I. Sivashinsky, Ann. Rev. Fluid Mech. 15, 179 (1983).

[9] K. Sneppen, J. Krug, H.M. Jensen, C. Jayaprakash und T. Bohr, Phys. Rev. A 46, R7351 (1992).

[10] C. Chow und T. Hwa, Physica D 84, 494 (1995).

[11] M. Rost und J. Krug, Physica D 88, 1 (1995).

[12] J. Krug, Origins of scale invariance in growth processes, Habilitationsschrift, Heinrich-Heine-Universität Düsseldorf (1995). Erschienen als Bericht Jül-3031 des Forschungszentrums Jülich, sowie in Advances in Physics 46, 139 (1997).

[13] D. Forster, D.R. Nelson und M.J. Stephen, Phys. Rev. A 16, 732 (1977).

[14] P.C. Hohenberg und B.I. Shraiman, Physica D 37, 109 (1989).

[15] G.I. Sivashinsky und D.M. Michelson, Prog. Theor. Phys. (Progress Lett.) 63, 2112 (1980).

[16] B. Nicolaenko, B. Scheurer und R. Temam, Physica D 16, 155 (1985); P. Collet, J.-P. Eckmann, H. Epstein und J. Stubbe, Commun. Math. Phys. 152, 203 (1993).

[17] V. Yakhot, Phys. Rev. A 24, 642 (1981).

[18] S. Zaleski, Physica D 34, 427 (1989).

[19] L. Golubović und R. Bruinsma, Phys. Rev. Lett. 66, 321 (1991); L. Golubović und R. Bruinsma, ebd. 67 2747 (1991).

[20] R. Cuerno und K.B. Lauritsen, Phys. Rev. E 52, 4853 (1995).

[21] I. Procaccia, M.H. Jensen, V.S. L’vov, K. Sneppen und R. Zeitak, Phys. Rev. A 46, 3220 (1992); V.S. L’vov und I. Procaccia, Phys. Rev. Lett. 69, 3543 (1992); V.S. L’vov, V.V. Lebedev, M. Paton und I. Procaccia, Nonlinearity 6, 25 (1993).

[22] F. Hayot, C. Jayaprakash und C. Josserand, Phys. Rev. E 47, 911 (1993); C. Jayaprakash, F. Hayot und R. Pandit, Phys. Rev. Lett. 71, 12 (1993).

[23] U. Frisch, Z.S. She und O. Thual, J. Fluid Mech. 168, 221 (1986).

[24] J. Krug, P. Meakin und T. Halpin-Healy, Phys. Rev. A 45, 639 (1992).

[25] M. Vergassola, B. Dubrulle, U. Frisch und A. Noullez, Astron. Astrophys. 289, 325 (1994).

[26] S. Toh und T. Kawahara, J. Phys. Soc. Japan 54, 1257 (1985).

[27] H. Chaté und P. Manneville, Phys. Rev. Lett. 58, 112 (1987).

[28] R. Cuerno und A.-L. Barabási, Phys. Rev. Lett. 74, 4746 (1995).

[29] J. Villain, Journal de Physique I 1, 1 (1991).

[30] D.E. Wolf, Phys. Rev. Lett. 67, 1783 (1991).

[31] O. Pierre-Louis und C. Misbah, Phys. Rev. Lett. 76, 4761 (1996).

[32] R.M. Bradley und J.M.E. Harper, J. Vac. Sci. Technol. A 6, 2390 (1988).

[33] A.Y. Cho und J.R. Arthur, Prog. Solid State Chem. 10 157 (1976).

[34] J.Y. Tsao, Materials Fundamentals of Molecular Beam Epitaxy, Academic Press, San Diego (1993).

[35] A. Zangwill, H.T. Dobbs und D.D. Vvedensky, Vorabdruck.

[36] D.E. Wolf in Scale Invariance, Interfaces, and Non-Equilibrium Dynamics, Hrsg.: A. McKane et al. Plenum Press, New York (1995).

[37] P. Politi, private Mitteilung.

[38] J. Villain und A. Pimpinelli, Physique de la Croissance Cristalline, Aléa–Saclay, Paris (1994).

[39] C. Ratsch, A. Zangwill, P. Šmilauer und D.D. Vvedensky, Phys. Rev. Lett. 72, 3194 (1994).

[40] M. Bott, T. Michely und G. Comsa, Surface Science 272, 161 (1992).

[41] G. Ehrlich und F.G. Hudda, Journal of Chemical Physics 44, 1039 (1966).

[42] E.W. Müller, Science 149, 591 (1965).

[43] W.K. Burton, N. Cabrera und F.C. Frank, Phil. Trans. Roy. Soc. (London) A 243, 299 (1951).

[44] J. Krug und M. Schimschak, J. Phys. I France 5, 1065 (1995).

[45] J. Krug, J. Stat. Phys. 87, 505 (1997).

[46] G. Comsa, private Mitteilung.

[47] H. Kallabis, L. Brendel, D.E. Wolf und J. Krug, Int. J. Mod. Phys. B, (im Druck).

[48] R.L. Schwoebel und E.J. Shipsey, Journal of Applied Physics 37, 3682 (1966).

[49] G.S. Bales and A. Zangwill, Phys. Rev. B 41, 5500 (1990).

[50] J. Krug und H. Spohn in Solids far from equilibrium, Hrsg.: C. Godrèche, Cambridge University Press (1991).

[51] E. Arzt, O. Kraft, W.D. Nix und J.E. Sanchez Jr., J. Appl. Phys. 76, 1563 (1994).

[52] J. Krug und H.T. Dobbs, Phys. Rev. Lett. 73, 1947 (1994).

[53] G.L. Eyink, J.L. Lebowitz und H. Spohn, J. Stat. Phys. 83, 385 (1996).

[54] G. Wulff, Z. Krist. Mineral. 34, 449 (1901).

[55] P. Nozières und F. Gallet, Journal de Physique 48, 353 (1987).

[56] M. Rost und H. Spohn, Phys. Rev. E 47, 3097 (1994); M. Rost und H. Spohn, Phys. Rev. Lett. 72, 784 (1994).

[57] W.W. Mullins, J. Appl. Phys. 30, 77 (1959).

[58] H. Spohn, J. de Physique I 3, 69 (1993).

[59] W. Selke und P.M. Duxbury, Phys. Rev. B 52, 17468 (1995).

[60] L. Onsager, Phys. Rev. 37, 405 (1931) und 38, 2265 (1931).

[61] H.B. Callen, Thermodynamics and an Introduction to Thermostatistics, John Wiley & Sons, New York (1985).

[62] D.E. Wolf und J. Villain, Europhys. Lett. 13, 389 (1990).

[63] J. Krug, M. Plischke und M. Siegert, Phys. Rev. Lett. 70, 3271 (1993).

[64] P. Šmilauer und D.D. Vvedensky, Phys. Rev. B 52, 14263 (1995).

[65] M. Siegert und M. Plischke, Phys. Rev. E 53, 307 (1996).

[66] G.W. Smith, A.J. Pidduck, C.R. Whitehouse, J.L. Glasperund J. Spowart, J. Cryst. Growth 127, 996 (1993); C. Orme, M.D. Johnson, K.-T. Leung, B.G. Orr, P. Šmilauer und D.D. Vvedensky, ebd. 150, 128 (1995); J.E. Van Nostrand, S.J. Chey, M.-A. Hasan, D.G. Cahill und J.E. Greene, Phys. Rev. Lett. 74, 1127 (1995); M.D. Johnson, C. Orme, A.W. Hunt, D. Graff, J. Sudijono, L.M. Sander und B.G. Orr, Phys. Rev. Lett. 72, 116 (1994).

[67] H.-J. Ernst, F. Fabre, R. Folkerts and J. Lapujoulade, Phys. Rev. Lett. 72, 112 (1994); K. Thürmer, R. Koch, M. Weber und K.H. Rieder, Phys. Rev. Lett. 75, 1767 (1995); J.A. Stroscio, D.T. Pierce, M. Stiles, A. Zangwill und L.M. Sander, Phys. Rev. Lett. 75, 4246 (1995).

[68] I.M. Lifshitz und V.V. Slyozov, J. Phys. Chem. Solids 19, 35 (1961).

[69] C. Wagner, Z. Elektrochem. 65, 581 (1961).

[70] M. Siegert in Scale Invariance, Interfaces, and Non-Equilibrium Dynamics, Hrsg.: A. McKane et al. Plenum Press, New York (1995).

[71] M. Rost und J. Krug, Phys. Rev. Letters 75, 3894 (1995).

[72] M. Rost und J. Krug, Phys. Rev. E 55, 3952 (1997).

[73] A. Coniglio und M. Zannetti, Europhys. Lett. 10, 575 (1989).

[74] H. Kallabis und P. Šmilauer, private Mitteilung.

[75] C. Castellano, private Mitteilung.

[76] T. Ohta, D. Jasnow und K. Kawasaki, Phys. Rev. Lett. 49, 1223 1982).

[77] K. Kawasaki, M.C. Yalabik und J.D. Gunton, Phys. Rev. A 17, 445 (1978).

[78] G.F. Mazenko, Phys. Rev. Lett. 63, 1605 (1989); G.F. Mazenko, Phys. Rev. B 42, 4487 (1990); G.F. Mazenko, Phys. Rev. B 43, 5747 (1991).

[79] A.D. Rutenberg und A.J. Bray, Phys. Rev. E 50, 1900 (1994).

[80] L. Golubović, Phys. Rev. Lett. 78, 90 (1997).

[81] M. Rost, P. Šmilauer und J. Krug, Surface Science 369, 393 (1996).

[82] P. Nozières in Solids far from equilibrium, Hrsg.: C. Godrèche, Cambridge University Press (1991).

[83] S.T. Chui und J.D. Weeks, Phys. Rev. B 14, 4978 (1978).

[84] H. van Beijeren und I. Nolden in Structure and dynamics of surfaces II , Hrsg.: W. Schommers und P. von Blanckenhagen, Springer (1990).

[85] H.K. Janssen, Vorabdruck.

[86] J.A. Stroscio, D.T. Pierce und R.A. Dragoset, Phys. Rev. Lett. 70, 3615 (1993).

[87] I.S. Gradshteyn and I.M. Ryzhik, Tables of Integrals, Series and Products, Section 3.462, Academic Press (1980).

116 LITERATURVERZEICHNIS

[88] W.H. Press, B.P. Flannery, S.A. Teukolsky und W.T. Vetterling, Numerical Recipies, Cambridge University Press (1986).

Danksagung

Vielen Personen ist für ihre Zusammenarbeit und Unterstützung beim Entstehen dieser Arbeit zu danken.

Allen voran danke ich J. Krug für die Themenvergabe, behutsame Führung und ständige Mitarbeit, und nicht zuletzt, weil er mir viele wissenschaftliche Kontakte ermöglicht hat.

Für anregende Diskussionen sei gedankt L.A.N. Amaral, T. Ala-Nissilä, A.-L. Barabási, A. Bray, H. Chaté, R. Cuerno, H. Dobbs, R. Grauer, J. Hager, A. Hanisch, T. Hwa, I. Janosi, H. Kallabis, H. Kinzelbach, M. Kraus, R. Lipowsky, S. Majaniemi, T. Newman, O. Pierre-Louis, P. Politi, M. Prähofer, N. Provatas, F. Rojas Iníguez, M. Schimschak, G. Schuckert, U. Schwarz, W. Selke, P. Šmilauer, H. Spohn, L.H. Tang, W. Wintz, D.E. Wolf und I. Yekutieli.

Im Arbeitsalltag und bei Schwierigkeiten am Rechnernetzwerk haben mich unterstützt U. Funk, C. Jansen, F. Jülicher, M. Kraus, R. Oberhage, H. Paffen, B. Sacha und M. Schimschak. Ihnen ein herzliches Dankeschön.

Für ihre Gastfreundschaft möchte ich danken dem Isaac Newton Institute in Cambridge, der Universität von Helsinki und dem ICTP Triest.

Finanziell gefördert wurde diese Arbeit vom Forschungszentrum Jülich, von der Universität GHS Essen und von der DFG im Sonderforschungsbereich 237.

Lebenslauf

Personalien Martin Rost Herwarthstraße 89 45138 Essen

26.6.1967 Geboren in München.

1973-1977 Grundschule in Leutershausen und Nußloch bei Heidelberg.

1977-1986 Kurfürst–Friedrich–Gymnasium in Heidelberg.

Abitur.

1986-1987 Militärdienst.

1987-1989 Physikstudium an der Ruprecht-Karls-Universität in Heidelberg. Vordiplom.

1989-1990 Physikstudium an der Université Paris–Sud in Orsay/Frankreich. Abschluß: Licence de Physique.

1990-1993 Physikstudium an der Ludwig-Maximilians-Universität in München. Diplomarbeit Beweglichkeit von Phasengrenzflächen unter Anleitung von Prof. Herbert Spohn. Diplom.

1993-1996 Promotionsstudium an der Universität zu Köln unter Anleitung von Dr. Joachim Krug am Institut für Festkörperforschung, Forschungszentrum Jülich.

seit 1.4.96 Promotionsstudium an der Universität Gesamthochschule Essen. unter Anleitung von Prof. Joachim Krug.

LITERATURVERZEICHNIS 119

Ich versichere, daß ich die von mir vorgelegte Disertation selbständig angefertigt, die benutzten Quellen und Hilfsmittel vollständig angegeben und die Stellen der Arbeit – einschließlich Tabellen, Karten und Abbildungen –, die anderen Werken im Wortlaut oder dem Sinn nach entnommen sind, in jedem Einzelfall als Entlehnung kenntlich gemacht habe. Die Bestimmungen der Promotionsordnung sind mir bekannt. Die von mir vorgelegte Dissertation ist von Prof. Dr. J. Krug betreut worden.