Anwendungen der Zeitabhängigen Dichtefunktionaltheorie und der Wigner-Maxwell Gleichungen in der Plasmonik für Simulationen im Ultrakurzzeitbereich vorgelegt von M. Sc. Mathias Wand Fakultät III- Naturwissenschaften der Universität Paderborn Department Physik Computational Nanophotonics Dissertation zur Erlangung des akademischen Doktorgrades Dr. rer. nat. genehmigte Dissertation Datum der mündlichen Prüfung: 13.01.2014 Promotionsausschuss: Vorsitzender: Prof. Dr. rer. nat. Thomas Zentgraf 1. Gutachter: Prof. Dr. rer. nat. Jens Förstner 2. Gutachter: Prof. Dr. rer. nat. Torsten Meier Zusammenfassung In dieser Arbeit werden Theorien und Simulationen entwickelt, mit denen die nichtlokalen und nichtlinearen optischen Eigenschaften metallischer Nanostrukturen untersucht werden können. Diese Nanostrukturen nden zur Zeit Anwendung im Bereich neuartiger Metamaterialien und Nanoantennen für den infraroten und den sichtbaren Spektralbereich. Eine detaillierte Analyse von Längenskalen zeigt, dass die klassischen Theorien zur Modellierung der Materie, speziell in den Oberächen der Strukturen, nicht angewendet werden können. Insbesondere wird versucht, die Licht-Materie Interaktion im Ultrakurzzeitbereich zu berechnen. Das Ergebnis solcher Berechnungen ergibt die zeitabhängige Polarisation der Materie und der damit verbundenen elektromagnetischen Felder über wenige hundert Femtosekunden. Damit lassen sich Anregungen durch ultrakurze Lichtpulse analysieren, wie sie bei Anrege-Abfrage Experimenten verwendet werden. Diese Experimente sollen letztlich durch Simulationen im Rechner durchgeführt werden können. Um dieses Ziel zu erreichen, wird eine Vielzahl moderner numerischer Methoden benötigt, um die relevanten Gleichungen für Licht und Materie mit angemessenem Rechenaufwand lösen zu können. Der Themenschwerpunkt dieser Arbeit liegt zunächst bei der Lösung von quantenmechanischen Vielteilchenproblemen in der Plasmonik. Dazu wird zum einen die Zeitabhängige Dichtefunktionaltheorie und zum anderen ein Ansatz, der die Dynamik der Dichtematrix über Wignerfunktionen beschreibt, auf Anwendbarkeit untersucht. Um optische Eigenschaften von Nanostrukturen berechnen zu können, muss dazu die selbstkonsistente Kopplung der Materiegleichungen an die Maxwellgleichungen berücksichtigt werden. Ein weiterer Themenschwerpunkt stellt die Erweiterung der Zeitabhängigen Dichtefunktionaltheorie um eine phänomenologische Dissipation zur Simulation von Streuprozessen dar. Diese wird auch benötigt, um Rückstreuungen von Ladungsdichtewellen in den begrenzt groÿen Nanostrukturen zu verhindern. Auÿerdem wird noch auf die Realisierung eines Finite-Dierenzen Verfahrens zur Berechnung der elektromagnetischen Potentiale in der Coulomb-Eichung eingegangen. Diese Eichung ndet in der Beschreibung von Optik in Festkörpern viele Anwendungen und es existieren bisher keine numerischen Lösungsverfahren zu diesem Problem. Zusammenfassend werden in dieser Arbeit verschiedene voll quantenmechanische Modellrechnungen für den Zeitbereich entwickelt und durchgeführt, welche zeigen, wo es in Nanostrukturen zur Entstehung von höheren Harmonischen und nicht-lokalen Eekten kommt. Alle Simulationen berechnen die optischen Eigenschaften der Strukturen auf mikroskopischer Ebene im Rahmen des Jellium-Modells. Dabei werden nicht-lokale Abschirmungseekte unmittelbar berücksichtigt. Mit Ausnahme der Simulationen zur dissipativen Dichtefunktionaltheorie werden keine phänomenologischen Parameter benötigt. 1 Abstract In this work theories and simulations are developed which can be used to study the nonlocal and nonlinear optical properties of metallic nanostructures. These nanostructures are currently used in the area of novel metamaterials and nanoantennas for the infra-red and visible spectral range. A detailed analysis of length scales shows that the classical theories for modelling the optical response of matter specically at the surfaces of the structures can not be applied. In particular, an attempt is made to calculate the light-matter interaction in the ultrafast time regime. The result of these calculations is given by the time-dependent polarization of matter and the associated electromagnetic elds over a few hundred femtoseconds. This allows the analysis of excitations caused by ultrashort light pulses which are used in pump-probe experiments. These experiments shall be carried out by computer simulations. To achieve this objective, a variety of modern numerical methods are required to solve the relevant equations of light and matter with reasonable computational complexity. A central topic of this thesis is the solution of quantum-mechanical many body problems in the eld of plasmonics. The time-dependent density functional theory and an approach, which describes the dynamics of the density matrix by Wignerfunctions, are examined for applicability. In order to calculate the optical properties of nanostructures, the self-consistent coupling of the equations for matter with Maxwell's equations must be considered. Another focal point is the extension of the time-dependent density functional theory by a phenomenological dissipation for the simulation of scattering processes. This is also required in order to prevent backscattering of charge density waves in small nanostructures. Moreover, the implementation of a nite dierence method for the calculation of electromagnetic potentials in the Coulomb gauge is discussed. This gauge is very common in the description of optics in solids and there are so far no numerical solution methods for this problem. In summary, several fully quantum mechanical model calculations in the time domain are developed and carried out, showing where it comes to the emergence of higher harmonics and nonlocal eects in nanostructures. All simulations calculate the optical properties of the structures at the microscopic level within the Jellium model. Nonlocal screening eects are ultimately taken into account. With the exception of simulations which involve dissipative density functional theory, no phenomenological parameters are required. 3 Inhaltsverzeichnis 1. Einleitung 1.1. Motivation.................................... 1.2. Experimente................................... 1.2.1. Ultraschnelle Nanooptik........................ 1.2.2. Messung reektierter Zweiter Harmonischer............. 1.2.3. Erzeugung Zweiter Harmonischer in magnetischen Metamaterialien 1.2.4. Messung von Oberächen- und Bulkbeiträgen zur SHG....... 1.3. Jellium-Modell................................. 1.4. Modelle für das Elektronengas......................... 1.5. Skalenanalyse und Gröÿenverhältnisse.................... 1.6. Vorgehensweise................................. 9 9 11 11 12 13 13 14 14 16 19 2. Physikalische Grundlagen 21 2.1. Plasmonik.................................... 21 2.1.1. Maxwell Gleichungen und Propagation elektromagnetischer Wellen 21 2.1.2. Dielektrische Funktion des freien Elektronengases.......... 22 2.1.3. Dispersionsrelation des freien Elektronengases............ 24 2.1.4. Reale Metalle.............................. 25 2.2. Dichtefunktionaltheorie............................ 26 2.2.1. Grundzustand.............................. 26 2.2.2. Zeitentwicklung............................. 31 3. Anwendungen der DFT auf Nanostrukturen 3.1. Metalllme................................... 3.1.1. Elektronischer Grundzustand..................... 3.1.2. Friedel Oszillationen.......................... 3.1.3. Bewegung der Metallelektronen im Lichtfeld............. 3.1.4. Erzeugung Zweiter Harmonischer................... 3.1.5. Mikroskopische Struktur der Fresnelfelder.............. 3.2. Nichtlokale Suszeptibilität der Elektronendichte............... 3.2.1. Denition der Antwortfunktion.................... 3.2.2. Numerische Berechnung der Antwortfunktion............ 3.2.3. Zusammenfassung der Berechnungsmethoden............ 3.2.4. Eigenschaften der Suszeptibilität im Metalllm........... 3.3. Metallische Nanodrähte............................ 3.3.1. Elektronischer Grundzustand..................... 3.3.2. Mikroskopische Ladungsdichte an Metallkanten und-ecken..... 3.3.3. Erzeugung Zweiter Harmonischer................... 35 35 35 36 37 40 42 44 44 45 47 48 51 52 53 53 5 Inhaltsverzeichnis 3.4. Metallische Nanopartikel............................ 3.4.1. Elektronischer Grundzustand..................... 3.4.2. Ladungsdichte in einer Metallkugel.................. 3.5. Multiskalensimulationen............................ 3.5.1. Multiskalenansatz für Oberächen.................. 54 55 56 56 58 4. Dissipative Zeitabhängige Dichtefunktionaltheorie 4.1. Motivation.................................... 4.1.1. Energieverlust durch Streuprozesse.................. 4.1.2. Bedeutung für Modellrechnungen................... 4.1.3. Vergleich von Verlusten durch Streuung und e/m-Abstrahlung... 4.2. Übersicht zu existierenden Ansätzen..................... 4.3. Methode von Neuhauser............................ 4.3.1. Herleitung für Ein-Teilchen Systeme................. 4.3.2. Dimensionsanalyse........................... 4.3.3. Anwendung auf Kohn-Sham Systeme................. 4.3.4. Energieerhaltung in zeitabhängigen Kohn-Sham Systemen..... 4.3.5. Energieverlustgleichung........................ 4.3.6. Alternative Formulierungen...................... 4.4. Analyse der Methoden............................. 4.4.1. Berechnung der Drude Streuzeit................... 4.4.2. Dämpfungseigenschaften........................ 4.4.3. Vorgegebene Energieabnahme..................... 4.4.4. Rückstreuung.............................. 4.4.5. Anmerkungen.............................. 61 61 61 62 62 66 67 67 69 70 70 72 74 77 78 78 80 81 83 5. Wigner-Maxwell Gleichungen 85 5.1. Allgemeine Formulierung in drei Dimensionen................ 85 5.1.1. Hamiltonoperator............................ 86 5.1.2. Impulsbasis............................... 87 5.1.3. Bewegungsgleichung der Kohärenzenmatrix............. 88 5.1.4. Formulierung im Wignerbild...................... 89 5.1.5. Gleichungssystem für ein Zweikomponentenplasma......... 91 5.1.6. Anwendung auf metallische Nanostrukturen............. 92 5.1.7. Phänomenologische Relaxationsterme................ 93 5.2. Nanodrähte................................... 94 5.2.1. Hamiltonoperator............................ 95 5.2.2. Eindimensionale Impulsbasis..................... 95 5.2.3. Coulombmatrixelemente........................ 97 5.2.4. Matrixelemente des Störpotentials.................. 98 5.2.5. Bewegungsgleichung der Kohärenzenmatrix............. 99 5.2.6. Observablen............................... 99 5.2.7. Eigenschaften der Kohärenzenmatrix................. 100 5.2.8. Transformation in das Wignerbild.................. 102 5.3. Wigner-Poisson System............................ 103 5.3.1. Herleitung................................ 103 6 Inhaltsverzeichnis 5.3.2. Fluidmodell für ein Quantenplasma.................. 105 6. Numerik 109 6.1. Lösungsverfahren für stationäre Kohn-Sham Gleichungen.......... 109 6.1.1. Problemübersicht............................ 109 6.1.2. Propagation in Imaginärzeit...................... 110 6.1.3. Orthogonalisierungsverfahren..................... 111 6.1.4. ITP-basierter Diagonalisierungs-Algorithmus............ 112 6.1.5. SCF Iterationsverfahren........................ 113 6.1.6. Adaptives Mixing Verfahren...................... 115 6.1.7. Bedeutung der Abschirmkonstante für die SCF-Schleife....... 116 6.2. Propagatoren für die Kohn-Sham Gleichungen................ 117 6.2.1. Problembeschreibung.......................... 117 6.2.2. Numerische Lösung der zeitabhängigen Ein-Teilchen Schrödingergleichung ................................ 118 6.2.3. Alternative Propagationsmethoden.................. 120 6.2.4. Extrapolation des Hamiltonoperators................. 123 6.3. Propagatoren für implizite Kohn-Sham Gleichungen............. 123 6.3.1. Problembeschreibung.......................... 123 6.3.2. Lösungsverfahren von Neuhauser................... 124 6.3.3. Diskretisierung der impliziten Kohn-Sham Gleichungen....... 125 6.3.4. Implizite Runge-Kutta Verfahren................... 126 6.3.5. Magnusentwicklung........................... 128 6.3.6. Analyse der Konditionszahl...................... 129 6.4. Methode der Finiten Volumen......................... 135 6.4.1. Problemübersicht............................ 135 6.4.2. Grundlagen............................... 136 6.4.3. Kurganov-Tadmor Methode...................... 138 6.4.4. Euler-Gleichungen........................... 139 6.4.5. Wigner-Gleichungen.......................... 140 6.5. Wigner-Maxwell Gleichungen......................... 141 6.5.1. Problemübersicht............................ 141 6.5.2. Präparation des Grundzustands.................... 143 6.5.3. Diskretisierung des Phasenraums................... 144 6.5.4. Zeitentwicklung............................. 145 k 6.5.5. Drude Bewegungsgleichung im-Raum............... 145 6.6. Lösungsverfahren für die Poisson-Gleichung................. 146 6.6.1. Problemübersicht............................ 146 6.6.2. FFT-Methode.............................. 147 6.7. Lösungsverfahren für elektromagnetische Potentiale in Coulomb-Eichung. 149 6.7.1. Problemübersicht............................ 149 6.7.2. Lösung der Wellengleichung...................... 151 7. Schlussfolgerung und Ausblick 155 7 Inhaltsverzeichnis A. Anhang 159 A.1. Zustandsdichten................................. 159 A.2. Gröÿen in der linearen Optik von Metallen.................. 159 A.3. Drude-Parameter für Edelmetalle....................... 160 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum...... 160 A.4.1. Stromdichteoperator im kontinuierlichen Impulsraum........ 160 A.4.2. Reibungsterm im kontinuierlichen Impulsraum........... 164 A.4.3. Stromdichteoperator im diskreten Impulsraum............ 168 A.4.4. Reibungsterm im diskreten Impulsraum............... 173 A.5. Nanodraht.................................... 174 A.5.1. Bewegungsgleichung der Kohärenzenmatrix(Herleitung)...... 174 A.5.2. Transformation in das Wignerbild(Herleitung)........... 176 A.6. Hartree-Energiefunktional........................... 180 Literaturverzeichnis 181 8 1. Einleitung 1.1. Motivation Die technologische Entwicklung von Nanofabrikationsprozessen ermöglicht seit einiger Zeit die Herstellung neuartiger metallischer Nanostrukturen welche in Metamaterialien[1] und Nanoantennen[2, 3] zum Einsatz kommen. Die Strukturen dienen zur Manipulation von Licht im sichtbaren und infraroten Spektralbereich und haben daher Abmessungen in der Gröÿenordnung der Lichtwellenlänge. Bei Metamaterialien(s. Abb. 1.1) geht es speziell darum, ein Medium mit vorgegebener eektiver Permittivität r ( ) und Permeabilität µ r ( ) zu schaen, was diverse technologische Anwendungen im Bereich der Optik ermöglicht[4]. Die Nanoantennen(s. Abb. 1.2) sollen die für Licht geeignete und daher entsprechend miniaturisierte Versionen der Radio- und Mikrowellenantennen sein: Deren Zweck besteht also ebenfalls darin, das Strahlungsfeld möglichst ezient mit einem Empfänger oder einem Sender zu koppeln, d.h. die elektromagnetische Strahlung auf einen Punkt zu konzentrieren oder von diesem Punkt in den freien Raum zu führen. Für spektroskopische Anwendungen dienen die Nanoantennen als Empfänger und Sender zugleich. Die Skalierung der Radioantennen in den Nanometerbereich zeigt vor allem, welche Besonderheiten auf diesen Längenskalen auftreten, die bei den(Zenti-)Meter groÿen Pendants nicht vorhanden sind: Im sichtbaren Spektralbereich können Plasmonresonanzen und Interbandübergänge angeregt werden. Vor allem liegt die Skintiefe des Lichtes bei Metallen in der Gröÿenordnung von den Strukturen selbst. Auch im Bulkbereich entstehen Polarisationsströme wodurch ohmsche Verluste stark an Bedeutung gewinnen. Bei makroskopischen Antennen wird dagegen nur die Oberäche angeregt, da die Strahlung an dieser fast vollständig reektiert wird. Wegen der geringen Skintiefe kann diese für Radiowellen als unendlich dünn betrachtet werden(s. Abb. A.1a und A.1b). Besonders die Beschreibung der Nanostrukturen mittels klassischer Elektrodynamik in der Form wie sie bei Radioantennen verwendet wird- wirft die Fragestellung auf, inwieweit diese auch hier angewendet werden darf: Aufgrund der hohen Lichtintensitäten, welche bei Metamaterialien und Nanoantennen auftreten, ist zwar die Annahme, dass die Quantisierung des Lichtfelds bedeutungslos ist, gesichert. Dagegen darf die Verwendung einer räumlich homogenen Permittivität ( ) in den Maxwellgleichungen mittlerweile als physikalisch unzureichend angesehen werden[5, 6], weil sich aufgrund der Längenskalenverhältnisse(s. Kap. 1.5) ortsabhängige, nichtlokale optische Eekte bemerkbar machen und die Homogenitätsannahme der Polarisierbarkeit in den Oberächenbereichen der Strukturen versagt. Dazu lässt sich auch eine sehr elementare Überlegung anstellen, die zeigt, dass gerade der Bulkbereich der Strukturen(in dem die Homogenitätsannahme noch am ehesten zutrit) bei Verkleinerung der Abmessungen immer unbedeutender wer9 1. Einleitung (a)(b)(c) (a) Abbildung 1.1.: Split-Ring Resonatoren(SRR) für den Mikrowellenbereich[Quelle: (b) NASA Glenn Research]. Elektronenmikroskopische Aufnahme eines (c) SRRs[1] für den optischen Bereich. Die nur wenige hundert Nanometer groÿen SRR aus Abb.(b) sind herstellungsbedingt auf einer ITO-Schicht mit Glas-Substrat angeordnet. Die Abbildung zeigt eine Elementarzelle des Metamaterials, welche in der Ebene periodisch fortgesetzt wird. Abbildung 1.2.: Elektronenmikroskopische Aufnahme einer Yagi-Uda Antenne für den optischen Spektralbereich[2]. den muss : Das Volumen skaliert mit x 3 und die Oberächen mit x 2 (s. Abb. 1.3). Dabei sei x= l/l 0 ein dimensionsloser Parameter, der durch(etwas willkürliche) Festlegung einer Referenzlänge l 0 darüber entscheided, wie extrem eine Struktur von Bulk- oder Oberächeneigenschaften bestimmt wird. Für die Nanostrukturen, deren Beschreibung diese Arbeit als Ziel hat, wird angenommen, dass diese in den Bereich x 1 gehören. Nichtlineare Eekte in Nanostrukturen, welche ebenfalls Gegenstand dieser Arbeit sind, lassen am Beispiel der Photoemission[7] sofort erkennen, dass diese eine quantenmechanische Berechnung der Elektronendichte erfordern. Die nichtlineare Antwortfunktion der Elektronendichte lässt sich in der klassischen Elektrodynamik über Suszeptibilitätstensoren höherer Ordnung[8] in die(makroskopischen) Maxwellgleichungen integrieren oder als selbstkonsistent zu berechnende Ströme und Ladungen an die VakuumMaxwellgleichungen koppeln. 1 Damit ist speziell die Unterscheidung von Bulk- und Oberächenbereich gemeint. Die Oberächen haben selbst eine Dicke in der Gröÿenordnung der Skintiefe. Die extrem kleinen Elemente der Nanoantenne in Abb. 1.2 bestehen demnach nur aus Oberächen. 10 1.2. Experimente Abbildung 1.3.: Die Oberächenbeiträge, welche mit x 2 skalieren(rote Kurve), gewinnen mit zunehmender Verkleinerung gegenüber den Bulkbeiträgen, welche mit x 3 skalieren(blaue Kurve), immer mehr an Bedeutung. Diese Überlegung dient nur als Denkanstoÿ im Bereich der Nanostrukturen und kann nicht als allgemeingültig angesehen werden: Die Grenze der Gültigkeit für x 1 zeigen z.B. Radioantennen, in denen die Strahlung aufgrund der Skintiefe nur Ladungen an der Oberäche anregt. Auf solchen makroskopischen Längenskalen gilt, dass die Suszeptibilität des Materials als homogen und lokal angesehen werden darf. Die andere Grenze der Gültigkeit für x 1 liegt bei Molekülen und Clustern, in denen die Begrie Oberäche und Bulk gar nicht mehr deniert sind. Das primäre Ziel dieser Arbeit besteht deshalb darin, ein Materialmodell für die Nanostrukturen zu nden, welches die Homogenitätsannahme vollständig umgeht und die optische Polarisierbarkeit als Lösung eines Vielteilchenproblems zu einer inhomogenen Ladungsverteilung liefert. Da die Skalenanalyse in Kap. 1.5 zeigen wird, dass auch im linearen optischen Bereich die Quanteneekte von Bedeutung sind, werden primär quantenmechanische Modelle betrachtet. Die Gültigkeit von klassischen Modellen muss sich durch Betrachtung bestimmter Grenzfälle erst noch zeigen. Die genaue Erkundung dieser Grenzfälle trägt nicht nur zum Verständnis der Optik in Nanostrukturen bei, sondern hilft auch dabei, künftige Modellrechnungen zu vereinfachen ohne deren Gültigkeit zu verletzen. 1.2. Experimente 1.2.1. Ultraschnelle Nanooptik Das neue Forschungsgebiet Ultraschnelle Nanooptik[9] umfasst eine Reihe von aktuellen wissenschaftlichen Fragestellungen, welche diese theoretische Arbeit motiviert haben: Es verbindet die Bereiche von Nanooptik(Photonik auf der Nanoskala, Plasmonik, Nano11 1. Einleitung Abbildung 1.4.: Geometrie zur Erzeugung eines 2 -Lichtstrahls durch Reektion eines 1 -Lichtstrahls an einer Metalloberäche(nach[11]). antennen, Nichtlineare Optik, etc.) und die ultraschnelle Laserspektroskopie(zeitaufgelöste Spektroskopie im Sub-Femtosekundenbereich, Dynamik ultraschneller Phänomene, kohärente Kontrolle, etc.) miteinander. In erster Linie geht es darum, optische Anregungen auf Zeitskalen im Femtosekundenbereich sowie auf Längenskalen im Nanometerbereich experimentell untersuchen und kontrollieren zu können. Diese Forschung ermöglicht die Entwicklung von neuartigen nanooptischen Geräten, wie z.B. Nanoantennen[2]. Auch für die Entwicklung von Metamaterialien, die im optischen Spektralbereich arbeiten[1, 4], können die Ergebnisse dieser Forschung genutzt werden. In[10] wird beschrieben, wie mittels Anrege-Abfrage(pump-probe) Methoden unter Verwendung neuester Femtosekundenlaser die Elektronenvorgänge(wie z.B. Augerprozesse) in Atomen und Molekülen experimentell im Zeitbereich beobachtet werden können. Die Entwicklung solcher experimenteller Methoden motiviert dazu, auch eine theoretische Beschreibung der Vorgänge im Zeitbereich zu entwickeln und durch Computersimulationen zugänglich zu machen. 1.2.2. Messung reektierter Zweiter Harmonischer Die Zweite Harmonische(SH) kann an Metalloberächen durch Messung der reektierten Strahlung(s. Abb. 1.4) erfolgen. In[11] werden zahlreiche experimentelle Ergebnisse zu diesem Versuch zusammengefasst: Es hat sich gezeigt, dass die SH-Strahlung extrem sensitiv auf die Beschaenheit der Oberäche reagiert und sogar eine Messung der Temperatur ermöglicht ! . Noch erstaunlicher ist das Ergebnis, dass sogar einzelne Atomlagen an der Oberäche die SH-Strahlung verändern können. Diese experimentelle Methode eignet sich daher für chemische Analysen und bietet für solche den Vorteil, dass sie nichtinvasiv ist und Messungen an Proben erlaubt, die sich in einem transparenten(z.B. üssigen) Medium benden. Diese experimentellen Resultate machen ganz klar deutlich, dass die theoretische Modellierung der Experimente eine genaue Erfassung der Oberächenmorphologie erfordert. Da in dieser Arbeit nur idealisierte Modelle von Metallen benutzt werden, ist ein Vergleich 2 engl.: second harmonic 3 Die Temperaturabhängigkeit ist deshalb erstaunlich, weil die Fermitemperatur des Elektronengases im Metall üblicherweise weit oberhalb des Schmelzpunktes liegt. 12 1.2. Experimente der hier angestellten Modellrechnungen mit experimentellen Ergebnissen nicht sinnvoll. Indem die Modellrechnungen schrittweise in Richtung realer Oberächen erweitert werden, können viele Erkenntnisse gewonnen werden, die auch ein Verständnis derzeitiger experimenteller Beobachtungen ermöglichen. In diesem Sinn stellt diese Arbeit einen wichtigen Grundstein für zukünftige Modellrechnungen dar. 1.2.3. Erzeugung Zweiter Harmonischer in magnetischen Metamaterialien In[1] wird berichtet, wie im Experiment die Erzeugung der Zweiten Harmonischen an magnetischen Metamaterialien beobachtet werden konnte. Das Metamaterial ist aus den Split-Ring Resonatoren, die in Abb. 1.1b und 1.1c dargestellt sind, aufgebaut. Diese Resonatoren kann man sich wie kleine LC-Schwingkreise vorstellen, bei denen der oene Ring die Spule und der Spalt den Kondensator bildet. Durch Einstrahlung von Licht mit horizontaler Polarisation kann ein zirkularer Stromuss angeregt werden, der das lokale Magnetfeld verstärkt. Bei diesem Experiment ist unklar, welche Mechanismen an der Erzeugung der beobachteten Zweiten Harmonischen beteiligt sind[12]. Der magnetische Anteil der Lorentzkraft q e ( v × B) liefert die einfachste Erklärung für dieses Phänomen, da die 1 -Anteile der Felder v und B zu 2 -Anteilen in der Kraft auf die Elektronen führen. In dieser Arbeit wird versucht, die Entstehung der Zweiten Harmonischen durch rein elektrische, nichtlineare Eekte an den Metall-Vakuum Grenzächen zu erklären. Wichtige theoretische Vorüberlegungen, welche maÿgeblich die Untersuchung dieser Eekte motiviert haben, wurden durch die Publikation von Rudnick und Stern[13] geliefert: Darin wurde insbesondere die Rolle der elektrischen Ströme entlang der Normalenrichtung in Oberächen von Metallen untersucht und gezeigt, dass diese quantenmechanisch berechnet werden müssen. 1.2.4. Messung von Oberächen- und Bulkbeiträgen zur SHG Die SH Oberächenpolarisation kann über den Oberächensuszeptibilitätstensor 2. Ordnung (s2f) phänomenologisch beschrieben werden: P sf (2 )= (s2f) : E( ) E( ) (1.1) Für die SH Bulkpolarisation, die auf magnetischen Dipol- und elektrischen Quadrupoltermen basiert[14], wird folgende phänomenologische Parameterisierung verwendet: P bulk (2 )= E() [ · E() ]+ [ E( ) · E() ]+ [ E( ) · ] E( ) (1.2) In[15] wird berichtet, dass es inzwischen möglich ist, die Komponenten des (s2f) -Tensors und die Bulkparameter experimentell mit einer Methode, die als two-phase beam secondharmonic generation bezeichnet wird, zu bestimmen. Die Autoren kommen zu dem Schluss, dass die Oberächenpolarisation die primäre Quelle für die gemessene SHStrahlung darstellt. Dieses Ergebnis stützt die angestrebte Vorgehensweise dieser Arbeit, sich primär auf die Entstehung der SH-Strahlung in Oberächenbereichen der metallischen Strukturen zu konzentrieren. 13 1. Einleitung Abbildung 1.5.: Jellium-Modell: Die metallische Nanostruktur setzt sich aus einer statischen, homogenen Ionendichte und einem Elektronengas zusammen. 1.3. Jellium-Modell Die metallische Bindung ist dadurch charakterisiert, dass es pro Atom ca. 1-2 ungebundene Elektronen gibt[16], die zu einem Gleichstrom beitragen können. Insgesamt steht eine groÿe Anzahl an frei beweglichen Leitungselektronen zur Verfügung. Für diese Elektronen näherungsweise erscheinen die Ionen als eine homogene Verteilung von positiven Ladungen, da das singuläre Potential der Atomkerne durch die(gebundenen) Valenzelektronen abgeschirmt wird. Jellium Diese Gegebenheiten führten zur Etablierung des-Modells[17], welches das Einfachste aller Modelle für Metalle darstellt: In Abb. 1.5 ist angedeutet wie sich damit metallische Nanostrukturen, wie z.B. die Split-Ring Resonatoren aus Abb. 1.1b modellieren lassen. Dazu wird die Geometrie der Struktur über eine statische Ionendichte n + ( r) vorgegeben und die Leitungselektronen werden als Elektronengas mit der Dichte n ( r, t) beschrieben. Dabei wird oengelassen, auf welchem theoretischen Niveau das Elektronengas beschrieben wird. Die wichtigsten Modelle werden in Kap. 1.4 besprochen. Die Ionendichte legt nicht nur die Geometrie, sondern auch die Art des Metalls über ihren Maximalwert n 0 fest. Dieser Maximalwert liegt üblicherweise überall im Inneren Wigner-Seitz der Struktur vor. Statt der Dichte n 0 wird beim Jellium-Modell der Radius r S benutzt: ( n 0 ) 1 = 4 3 r S3 (1.3) Es ist üblich diesen Parameter in atomaren Einheiten zu spezizieren, wodurch übliche numerische Werte im Bereich von 2 bis 6 liegen. Bei r S = 3 a 0 entspricht die Dichte n 0 ungefähr der Leitungselektronendichte von Gold. 1.4. Modelle für das Elektronengas Das Jellium-Modell aus Kap. 1.3 bildet in dieser Arbeit die Grundlage für die Modellierung von Festköpern. Für das Elektronengas werden folgende Modelle benutzt: Drude Modell 1.: Es handelt sich um ein(über hundert Jahre altes) klassisches Modell, welches die elektrische Leitfähigkeit von Metallen mikroskopisch zu erklären versucht[18]. Dieses Modell wird heutzutage immer noch benutzt und auch für die Optik von Metallen verwendet. In seiner Anwendung besteht das Modell letztlich 14 1.4. Modelle für das Elektronengas phänomenologischen darin, dass folgende klassische Bewegungsgleichung mit einem Streuzeitparameter für die Bewegung der Elektronen in einem elektrischen Feld E benutzt wird: m e ¨ r= m e r + q e E (1.4) Im Übergang von Punktladungen zum Kontinuumsmodell ¨ r j +( j · ) j t mit Stromdichtefeld j= q e v und Elektronendichte n 0 wird die nicht-lineare konvektive Ableitung weggelassen: j = j + q e2 n 0 E t m e (1.5) Letztlich wird die eigentliche Arbeit von Drude auf einen experimentell ermittelten Streuparameter und eine klassische Bewegungsgleichung reduziert. Diese newtonsche Bewegungsgleichung berücksichtigt vor allem die Trägheit der Elektronen, welche in schnell oszillierenden Feldern(wie in der Optik) berücksichtigt werden muss " . Im sichtbaren Spektralbereich kann das Drude Modell die optischen Eigenschaften von Metallen nur unzureichend beschreiben, da hier Interbandübergänge angeregt werden[20]. 2. Hydrodynamisches Modell : In der Terminologie dieser Arbeit wird damit die nicht-lineare und semiklassische Erweiterung des Drude Modells bezeichnet: j +( j · ) j= j + q e2 n 0 E+ 1 j × B p t m e m e m e (1.6a) + · j= 0 t (1.6b) Für den Druck p wird noch zusätzlich eine thermische Zustandsgleichung benötigt. Zur Beschreibung des Elektronengases wird hier die Zustandsgleichung des idealen Fermigases( p n 5/ 3 ) benutzt.(In Kapitel 5.9 von[17] wird beschrieben, wie sich die vollständige Navier-Stokes Gleichung auf ein Elektronengas anwenden lässt.) 3. Wechselwirkendes Elektronengas : Es ist ein rein quantenmechanisches Modell und wird durch folgenden Hamiltonoperator beschrieben: H ^ = N j=1 1 2 m e [ p ^ j q e A(^ r j , t)] 2 + 1 2 N N j=1 k=1 | ^ r j q e2 ^ r k | + N j=1 [ q e (^ r j , t)+ v ext (^ r j )] (1.7) Der Zustand dieses Systems ist durch eine antisymmetrische N -Teilchen Wellenfunktion ( r 1 , 1 ,..., r N , N , t) gegeben. 4 Bei Anwendung des Ohmschen Gesetzes J( t)= E( t) in der Elektrizitätslehre sind die Felder deutlich niederfrequenter und die Trägheit wird deshalb vernachlässigt. Auf die Bedeutung der Trägheit wird z.B. in[19] eingegangen. 15 1. Einleitung 1.5. Skalenanalyse und Gröÿenverhältnisse Ein wesentlicher Gegenstand dieser Arbeit ist es, die richtigen Längen- und Zeitskalen sowie die Verhältnisse zwischen diesen zu erkennen. Erst dadurch wird es möglich, eine theoretische Beschreibung zu nden, welche die relevanten Eekte korrekt wiedergeben kann und trotzdem nicht zu aufwendig ist, um das System beschreiben zu können. Häug besteht das Problem darin, dass die existierenden Modelle zu starke Vereinfachungen vornehmen und es zu erkennen gilt, welche dieser Vereinfachungen falsch sind. Ein solches Problem konnten Rudnick und Stern[13] lösen, indem sie das Verhältnis v Fermi / l em (1.8) betrachtet haben: Dabei ist v Fermi die Fermigeschwindigkeit der Metallelektronen, l em bezeichnet die Strecke, auf der sich das elektromagnetische Feld signikant ändert und die Frequenz des Lichtes. Der Quotient ist im Wesentlichen das Verhältnis von Periodendauer des Lichts T zu der Zeit T e , die ein Elektron für die Strecke l em benötigt. Die Strecke l em ist nahe der Oberäche sehr viel kleiner als im Inneren des Festkörpers, weil hier einerseits der Groÿteil der Strahlung auf der Länge der Skintiefe reektiert wird und andererseits die Oberäche eine steile Potentialbarriere darstellt # . In der Metalloberäche nicht in Richtung der Normalen gilt deshalb T T e wie im Inneren des Festkörpers: Stattdessen gilt, dass die Elektronen eine signikante Strecke(d.h. l em ) während einer Periode des Lichts zurücklegen können. Dieses Kriterium zeigt, dass die Berechnung der Oberächenpolarisierbarkeit ein quantenmechanisches Problem darstellt $ . Quantenmechanische Vielteilcheneekte führen insbesondere dazu, dass die Polarisierbarkeit(neben der zeitlichen) auch eine räumliche Nichtlokalität besitzt. In der Honung, auf ähnlich bedeutende Gröÿenverhältnisse zu stoÿen, wurden für diese Arbeit noch folgende weitere Gröÿen betrachtet % : · Verhältnis von Lichtwellenlänge zur Gröÿe der Nanostruktur: Wenn die Wellenlänge wesentlich gröÿer als die Struktur ist, können Retardierungseekte vernachlässigt werden. Die Maxwellgleichungen können dann in quasistatischer oder sogar elektrostatischer Näherung berücksichtigt werden. · Skintiefe: Diese kann benutzt werden, um die Dicke der Oberäche abzuschätzen und um damit zu entscheiden, wie groÿ das Verhältnis zwischen Bulk- und Oberächenanteilen einer Struktur ist(s. Abb. 1.3). · Verhältnis von Skintiefe zu Wellenlänge im Vakuum: s. Abb. A.1b. · Kräfte-Verhältnis zwischen treibender elektrischer Kraft durch ein externes Lichtfeld und der Kraft, welche von der Barriere des eektiven Potentials V eff am Rand 5 In diesem Zusammenhang erscheint es sinnvoll, die Strecke l em auf das eektive, zeitabhängige 6 Potential der Elektronen zu beziehen anstatt nur auf das elektromagnetische Feld. Die Quantenmechanik wird relevant, wenn sich das Potential auf der DeBroglie-Wellenlänge des Elektrons ändert[21]. Das Kriterium von Rudnick und Stern, welches das Verhältnis 1.8 benutzt, ist im 7 Prinzip identisch damit und hat die gleiche Konsequenz. Die Zusammenfassung der Gröÿen für die Charakterisierung des Elektronengases orientiert sich an [22]. 16 1.5. Skalenanalyse und Gröÿenverhältnisse des Festkörpers verursacht wird: | F em | = | q e | E 0 | F eff || V eff | (1.9) Dieses Verhältnis kann benutzt werden, um abschätzen zu können, ab welcher Feldstärke E 0 der Lichtquelle nichtlineare Eekte an den Oberächen des Festkörpers relevant werden. · Abschirmlänge in einem entarteten Plasma: F = 2 0 E F 3 n 0 q e2 (1.10) Diese Länge charakterisiert die Wellenlänge von Friedel-Oszillationen, die bei der Abschirmung von Überschussladungen in Plasmen auftreten. Diese Friedel-Oszillationen sind an den Metalloberächen zu erwarten, wo das Elektronengas den ionischen Hintergrund abschirmt. · Entartungsparameter: = T F /T (1.11) Für Metalle kann aufgrund der hohen Fermitemperatur T F immer von T= 0 K ausgegangen werden. Das Elektronengas ist also(maximal) entartet.(Deshalb ist hier nur die Fermi Abschirmlänge 1.10 und nicht die Debye Abschirmlänge relevant.) · Mittlere Coulomb-Wechselwirkungsenergie zwischen zwei Elektronen, deren Abstand von der Teilchendichte n 0 abhängt: U int = q e2 n 10/ 3 0 (1.12) · Kopplungsparameter(Verhältnis von Coulomb Energie U int und kinetischer Energie K= k B T F ): = U int K = 2 m e q e2 (3 2 ) 2/ 3 0 2 n 0 1/ 3 (1.13) An diesem Verhältnis kann man sehen, dass sich ein Elektronengas mit hoher Teilchendichte wie ein ideales Gas verhält. In[17](Kap. 1.3.3) wird dieser Zusammenhang am Hamiltonoperator(eines unendlich ausgedehnten Jellium-Festkörpers mit wechselwirkenden Elektronen) und der potentielle Anteil mit demonstriert: r S 1 . Der kinetische Anteil skaliert mit r S 2 Für den Kopplungsparameter in Gold & gilt Au 20 , so dass Kollisionen zwischen nicht den Elektronen oensichtlich vernachlässigt werden können. Allerdings ist nach[23] auf Grund des Pauli-Prinzips die Kollisionsrate zwischen den Elektronen sehr gering, da für T T F nur wenige Elektronen ihren Quantenzustand ändern können. Damit wäre ein kollisionsfreies Modell für das Elektronengas zu rechtfertigen. 8 Teilchendichte n 0 10 28 m 3 17 1. Einleitung · Nach[23] zeigen ausserdem folgende Zeitkonstanten für ein typisches Metall, dass Kollisionen im Ultrakurzzeit-Bereich vermutlich unerheblich sind: p = 10 16 s= 2 p 1 r = 10 14 s ee = 10 10 s (1.14a) (1.14b) (1.14c) Die Zeitkonstante p entspricht der Periodendauer der Plasmaoszillationen(Plasmafrequenz p ). Die Relaxationszeit r ist die Zeitkonstante mit der das System wieder in das thermische Gleichgewicht zurückkehrt(z.B. nach Abschalten eines elektrischen Feldes, welches einen Strom erzeugt hat). Wie r abgeschätzt werden kann, wird in[23] nicht erklärt. In[24](Kap. 1.2) wird r als Relaxationszeit des freien Elektronengases bezeichnet, die durch nicht näher spezizierte Kollisionen bedingt ist und als Streuzeit im Drude-Modell benutzt wird. Die Elektron-Elektron Streuzeit ee kann nach[23] über die Dicke der Fermischale und der Energie-Zeitunschärfe Relation abgeschätzt werden. · In unmittelbarem Zusammenhang mit den Zeitkonstanten 1.14a-1.14c steht die Frage, wie lange es dauert, damit sich eine Nichtgleichgewichtsverteilung von elektrischen Ladungen im Leiter wieder in den Gleichgewichtszustand begibt: Die Abschätzung der Dauer über die Kontinuitätsgleichung t + · J= 0 und dem Ohmschen Gesetz J= E liefert gänzlich falsche Ergebnisse ' , da sich dieser Vorgang nur durch Berücksichtigung der gesamten Maxwell-Gleichungen beschreiben lässt. Nach[25] besteht dieser Relaxationsvorgang nämlich aus drei Teilen: Erst werden die Ladungen aus dem Inneren des Leiters verdrängt, dann die(damit verbundenen) Ströme und elektromagnetischen Felder und letztlich verlieren die Oberächenladungen an kinetischer Energie durch elektromagnetische Abstrahlung und Streuprozesse. Die gesamte Relaxationszeit ist nach dieser Betrachtung auch von der Geometrie des Leiters abhängig. Für den Spezialfall eines Hohlzylinders mit Deckeln wird in [25] auch eine Formel hergeleitet. Aus dieser Analyse der Längen- und Zeitskalen ergeben sich im Wesentlichen folgende Konsequenzen: Es kann die elektrostatische Näherung für die Maxwellgleichungen benutzt werden, da in dieser Arbeit nur Strukturen, mit einer maximalen Gröÿe von ca. 5 nm im optischen Lichtfeld betrachtet werden. Die Strukturen haben einen vernachlässigbar kleinen Bulkanteil und eine Trennung von Bulk- und Oberäche ist nicht sinnvoll. Stattdessen wird für die gesamte Struktur die gleiche theoretische Beschreibung benutzt. Die Zeitkonstante ee (s. Gl. 1.14c) suggeriert, dass Streuprozesse zwischen den Elektronen vernachlässigt werden können. Nur die Multiskalensimulation in Kap. 3.5 macht hier eine Ausnahme: Diese unterscheidet zwischen Bulk- und Oberächenanteilen, und da hier ein ganzes Array von Split-Ring Resonatoren(wie in Abb. 1.1b) simuliert wird, müssen die Maxwellgleichungen ohne Näherung berücksichtigt werden. 9 Zur Gültigkeit des Ohmschen Gesetzes: s. Fuÿnote 4. 18 1.6. Vorgehensweise 1.6. Vorgehensweise Die Vorgehensweise in dieser Arbeit wird durch aktuelle Experimente im Bereich der Plasmonik(Kap. 1.2), die verschiedenen Modelle für Metalle(Kap. 1.4) und die Längenund Zeitskalenanalyse(Kap. 1.5) begründet: Die Experimente motivieren dazu, Simulationen der mikroskopischen Vorgänge im Zeitbereich zu beschreiben. Von den Metallmodellen kann nur das rein quantenmechanische Modell(Gl. 1.7) aufgrund der Ergebnisse der Skalenanalyse benutzt werden. Klassische und semi-klassische Modelle sind nur in Grenzfällen, bei denen die Oberächen eine untergeordnete Rolle spielen, anwendbar. Es müssen Theorien zur Lösung des quantenmechanischen Modells(Gl. 1.7) verwendet werden, die das Vielteilchenproblem so weit vereinfachen, dass die resultierenden Gleichungen einerseits technisch lösbar sind und andererseits noch die Quantennatur der Teilchen möglichst genau wiedergeben. Für dieses spezielle Problem wird in dieser Arbeit die(Zeitabhängige) Dichtefunktionaltheorie(Kap. 3 und 4) verwendet, da diese Theorie genau den genannten Anforderungen gerecht wird. Insbesondere können Verbesserungen an den Lösungen dadurch erreicht werden, dass man das sogenannte xc -Funktional durch eines ersetzt, welches besser für das spezielle System geeignet ist. Mögliche Verbesserungen bei der Genauigkeit des xc -Funktionals können in zukünftigen Arbeiten untersucht werden. Neben der Dichtefunktionaltheorie wird auch ein Dichtematrix-basierter Formalismus, der die Zeitentwicklung der Wigner-Verteilungsfunktion beschreibt, zur Lösung des quantenmechanischen Modells betrachtet(Kap. 5). Mit diesem Formalismus kann die Zeitentwicklung von statistischen Ensembles beschrieben werden. Davon wird allerdings kein Gebrauch gemacht, weil bei Metallen eine Beschreibung mit T= 0 K völlig ausreichend ist(s. Erläuterung zu Gl. 1.11) und somit nur ein reiner Quantenzustand(d.h. der elektronische Grundzustand als Startzustand) betrachtet werden muss. Hier soll speziell untersucht werden, ob die Gleichungen numerisch gelöst werden können und eine interessante Alternative zur Dichtefunktionaltheorie darstellen. Zeitbereich Die Gleichungen werden explizit im(und nicht im Frequenzbereich ) formuliert. Die Verwendung von optischen Suszeptibilitätstensoren(und auch von Antwortfunktionen des Elektronengases) wird vollständig umgangen. Alle Simulationen berücksichtigen auch einen Teil des Vakuums ausserhalb der Struktur. Damit wird weiter umgangen, dass an Materialgrenzen Randbedingungen gemacht werden müssen, welche im Bereich der nicht-lokalen Optik immer wieder für Diskussionen gesorgt haben[26], da deren Korrektheit unklar ist. Stattdessen wird ein Anfangswertproblem gelöst, bei dem die Randbedingungen im Vakuum liegen und sich die Elektronen an den Oberächen ohne den Einuss künstlicher Randbedingungen bewegen können. Die so erzielten Ergebnisse sollten gerade in dem Oberächenbereich wesentlich plausibler sein als solche, die mit Randbedingungen an den Oberächen hergeleitet wurden. Eine weitere Stärke der Formulierung im Zeitbereich ist, dass die Gleichungen ohne Störungstheorie hergeleitet werden und somit speziell für die Simulation der nicht-linearen Eigenschaften von Metallen geeignet sind. Es können zeitlich beliebig geformte Pulse mit beliebig hoher Feldstärke zur Anregung der Strukturen(genau wie im Experiment) 10 In der nicht-linearen Optik[8] ist die Formulierung im Frequenzraum üblich. 19 1. Einleitung benutzt werden, ohne die Gültigkeit der Gleichungen in Frage zu stellen . Der Nachteil der Methodik besteht darin, dass nur sehr kleine Strukturen(wie z.B. Teile von Nanoantennen < 10 nm ) beschrieben werden können und die Verwendung der technisch nicht-trivialen Multiskalensimulationen(Kap. 3.5) für die Simulation von gröÿeren Strukturen(wie z.B. Metamaterialien) sich als unumgänglich herausstellen wird. 11 Eine Ausnahme betrit die gemachte Langwellennäherung in Kap. 3, nach der die Wellenlänge des Lichtes deutlich gröÿer als die Gröÿe der Struktur sein muss. 20 2. Physikalische Grundlagen 2.1. Plasmonik Hier werden elementare Grundlagen aus dem Bereich der Plasmonik[16, 24] zusammengefasst und dabei auf Aspekte, die in dieser Arbeit besonders wichtig sind, genauer eingegangen. 2.1.1. Maxwell Gleichungen und Propagation elektromagnetischer Wellen Die Modelle für Metalle(s. Kapitel 1.3 und 1.4) sind rein mikroskopische Modelle, bei denen die Ladungs- und Stromdichten durch räumlich stetige Funktionen gegeben sind. Für diese Modelle besteht keine Notwendigkeit, die Ladungen und Ströme in gebunden und frei aufzuteilen, da keine gebundenen Ladungen vorhanden sind. Die Lichtausbreitung wird vollständig durch die mikroskopischen Maxwellgleichungen[27] beschrieben: 1 · E= 0 · B= 0 × E= B t × B= E µ 0 J+ µ 0 0 t (2.1a) (2.1b) (2.1c) (2.1d) Die Kopplung zwischen den elektromagnetischen Feldern und der Materie wird hier dadurch zum Ausdruck gebracht, dass man die Ladungs- und Stromdichte als Funktional des elektrischen und magnetischen Feldes schreibt: J= J[ E, B] = [ E, B] (2.2a) (2.2b) Da die Maxwellgleichungen die Kontinuitätsgleichung implizieren und die Ladungsdichte durch Gl. 2.1a gegeben ist, wird nur ein Gesetz für die Stromdichte benötigt. Im linearen Fall lässt sich die Gl. 2.2a in Form eines räumlich und zeitlich nichtlokalen Ohmschen Gesetzes schreiben: J( r, t)= d 3 r dt ( r, r r , t t ) E( r , t ) (2.3) Die Leitfähigkeit ist ein Tensor 2. Stufe, dessen Komponenten von einer absoluten und einer relativen räumlichen Koordinate sowie einer relativen zeitlichen Koordinate abhängen. In isotropen Medien reduziert sich die Leitfähigkeit auf ein Skalar. In homogenen 1 Die ist Ionendichte darf am Rand des Festkörpers unstetig, hier, dass es keine atomistischen Dichtefunktionen ds.ihn.ds:t u ( f r e ) na = rt ig j a q u e f N 3 ( u r ll ab r f j a ) llen. Entscheidend 21 2. Physikalische Grundlagen Medien entfällt die Abhängigkeit von der absoluten Raumkoordinate. Wie in Kapitel 1 erläutert wurde, trit die Homogenitätsannahme in Nanostrukturen nicht zu und daher kann auf die absolute Raumkoordinate für deren Beschreibung nicht verzichtet werden. Aus den Rotationsmaxwellgleichungen 2.1c und 2.1d folgen die inhomogenen Wellengleichungen für das elektrische und magnetische Feld, in denen die Stromdichte als Quellterm auftritt: 1 2 E J × × E+ c 02 t 2 = µ 0 t (2.4a) × × B+ 1 c 02 2 B t 2 = µ 0 × J (2.4b) Dabei ist c 0 = 1/ µ 0 0 die Lichtgeschwindigkeit. Mit Ausnahme der Multiskalensimulation in Kap. 3.5 kann in allen anderen Fällen angenommen werden, dass die Distanz c 0 / (wobei die Frequenz des externen Lichtfeldes ist) viel gröÿer ist als die Abmessung der Struktur ist. Deshalb ist für diese Systeme die elektrostatische Näherung der Maxwellgleichungen gerechtfertigt: · E 1 0 × E 0 (2.5a) (2.5b) Alle Felder sind dabei weiterhin zeitabhängig. Das elektrische Feld kann hier über ein elektrostatisches Potential ( r, t) , welches die Poissongleichung 2 ( r, t)= 1 ( r, t) 0 (2.6) erfüllt, in guter Näherung beschrieben werden. Die Multiskalensimulation verwendet die vollständige als auch die elektrostatische Form der Maxwellgleichungen. 2.1.2. Dielektrische Funktion des freien Elektronengases In Gl. 2.3 wurde die Leitfähigkeit als lineare Antwortfunktion der Stromdichte J auf das Feld E deniert. Das freie Elektronengas ist isotrop und homogen. In diesem Fall lautet die Relation 2.3: J( r, t)= dt ( t t ) E( r, t ) (2.7) Nach dem Faltungstheorem ergibt die Fouriertransformation dieser Relation in der Zeit: J( r, )= ( ) E( r, ) (2.8) In der Optik ist es üblich, die dielektrische Funktion r ( ) (bzw. die Suszeptibilität = r 1 ) als Antwortfunktion zu verwenden . Diese hat folgenden Zusammenhang mit 2 nur Bei den makroskopischen Maxwellgleichungen ist es üblich durch die dielektrische Funktion die gebundenen Ladungen und deren Polarisationsstrom zu berücksichtigen. Die freien Ladungen werden separat über die Leitfähigkeit einbezogen[21]. Hier wird von dieser Konvention abgewichen und keine Unterscheidung bei den Ladungen und Strömen gemacht, wie unmittelbar in Gl. 2.9 zum Ausdruck kommt. 22 2.1. Plasmonik der Leitfähigkeit: i( ) r ( )= 1+ 0 ( )= 0 i [ r ( ) 1] (2.9) Für das freie Elektronengas kann diese Funktion aus der klassischen Drude Bewegungsgleichung 1.4 hergeleitet werden, indem die Auslenkung r eines Elektrons mit der Polarisationsdichte P= q e n 0 r in Verbindung gesetzt wird. Als treibenden Term setzt man in die Bewegungsgleichung ein harmonisch oszillierendes Feld E( t)= E 0 e it ein. Insgesamt erhält man folgenden Zusammenhang zwischen Polarisation und elektrischem Feld: P= n 0 q e2 m e ( 2 + i) E (2.10) Über den Zusammenhang P= 0 E= 0 ( r 1) E kann die dielektrische Funktion abgelesen werden: r ( )= 1 2 P2 + i (2.11) Dabei ist P die Plasmafrequenz: P = n 0 q e2 0 m e (2.12) Auch wenn dieses Modell eine zu starke Vereinfachung darstellt, um adäquat die Optik von metallischen Nanostrukturen zu beschreiben, liefert diese Antwortfunktion qualitatives Verständnis und eine Menge an wichtigen Informationen um Abschätzungen durchzuführen. Die Einschränkungen des Modells sind nach der gegebenen Herleitung oensichtlich: 1. Dynamische Abschirmeekte in den Oberächenbereichen werden falsch beschrieben, weil in Gl. 2.10 angenommen wird, dass die Polarisation nur vom lokalen elektrischen Feld abhängig ist.(Grund für die Nichtlokalität: s. Kap. 1.5). Das führt eld enhancement zu Abweichungen beim im Vergleich zu quantenmechanischen Modellen[6, 7]. 2. Es kann weder die Photoemission noch das Auftreten von Tunnelströmen beschreiben. 3. Als lineares Modell kann es nicht die Entstehung von Höheren Harmonischen erklären. In Anhang A.2 sind die wichtigsten Formeln für dieses Modell zusammengefasst. 23 2. Physikalische Grundlagen 2.1.3. Dispersionsrelation des freien Elektronengases Unter Berücksichtigung der zeitlichen Nichtlokalität der Stromdichte, die durch die Leitfähigkeit in Gl. 2.7(oder durch die Dielektrische Funktion r ) beschrieben wird, lautet die inhomogene Wellengleichung 2.4a für das elektrische Feld: × × E( r, t)+ 1 c 02 2 E( r, t 2 t) = = µ 0 µ 0 t dt dt ( t t ) E( r, t ) ( t t ) E( r, t ) t Mittels Fouriertransformation der Wellengleichung in der Zeit( t i ) folgt: × × E( r, ) 2 c 02 E( r, )= µ 0 ( i) ( ) E( r, ) Als nächstes wird die ebene, monochromatische Welle E( r, )= E 0 ( ) e i k · r eingesetzt: { i k × [ i k × E( r,) ] } 2 c 02 E( r, )= iµ 0 ( ) E( r, ) { k[ k · E( r,) ] k 2 E( r, } ) 2 c 02 E( r, )= iµ 0 ( ) E( r, ) Für transversale Wellen ist k · E= 0 . Mit Gl. 2.9 folgt dann: [ 2 c 02 ] k 2 = iµ 0 ( )= 2 c 02 [1 r () ] Somit lautet die gesuchte Dispersionsrelation für transversale Wellen: k 2 2 c 02 r ( )= 0 (2.13) Zur Vereinfachung wird nun von einem idealen Metall( = 0 ) ausgegangen. Nach Gl. 2.11 ergibt sich dann folgende Dispersionsrelation: r ( )= 1 P2 2 k 2 2 c 02 + 2 P2 c 02 2 =0 1 k( )= c 0 2 P2 (2.14) Daran erkennt man, dass im Frequenzbereich < P die Wellenzahl k imaginär ist und keine Wellenausbreitung möglich ist. Diese Dispersionsrelation ist auch in Abb. 2.1 grasch dargestellt: Hier sieht man an der Gruppengeschwindigkeit d/dk , dass diese grundsätzlich kleiner als die Lichtgeschwindigkeit ist und für k 0 verschwindet. Letzteres beschreibt eine longitudinale, kollektive Schwingung des Elektronengases, die mit der Plasmafrequenz P oszilliert ! . 3 Die Plasmafrequenz kann anschaulich über die elektrostatische Kraft von Oberächenladungen in einer Metallschicht hergeleitet werden[16, 24]. 24 : P 2.1. Plasmonik 2 1.5 1 0.5 0 0 0.5 1 1.5 2 kc/ P Abbildung 2.1.: Plot der Dispersionsrelation für das freie Elektronengas 2.14(rote Kurve). Die blaue Linie zeigt die Relation für Licht: = ck . 2.1.4. Reale Metalle Die wichtigsten Metalle für plasmonische Anwendungen im nahen Infrarot- und sichtbaren Spektralbereich sind Gold und Silber[24]. Alle Modelle in dieser theoretischen Arbeit befassen sich ausschlieÿlich mit dem Elektronengas, welches die Leitungselektronen in der metallischen Bindung bilden(s. Kap. 1.4). Es soll daher kurz erläutert werden, welche Unterschiede zwischen den idealisierten Jellium-Modellen der Metalle und den realen Metallen bestehen: 1. Das Auftreten von Interbandübergängen ab einer bestimmten Mindestenergie der Photonen stellt den primären Unterschied dar. Bei Gold treten diese bereits im nahen Infrarotbereich auf[20]. 2. In realen Metallen streuen die Elektronen untereinander sowie an Phononen, Gitterdefekten und Verunreinigungen. Die mittlere freie Weglänge ist daher temperaturabhängig und kann auÿerdem in Nanostrukturen, deren Abmessungen klein genug sind, auch von deren Geometrie abhängen. Die Interbandübergänge können über das lorentzsche Oszillatormodell[27] in die Modellrechnungen integriert werden. Das trägt erheblich zu einer qualitativen und quantitativen Verbesserung der Modellrechnungen bei und stellt eine Voraussetzung dar, um die Ergebnisse von Theorie und Experiment überhaupt sinnvoll vergleichen zu können. Eine weitere Verbesserung des Modells betrit die mittlere freie Weglänge der Elektronen: In [28] wird eine Möglichkeit gezeigt, wie diese abhängig von der Form der Nanostruktur berechnet werden kann. 25 2. Physikalische Grundlagen 2.2. Dichtefunktionaltheorie Die Dichtefunktionaltheorie(DFT) wurde in den 1960er Jahren entwickelt, um die elektronischen Vielteilchenprobleme bei Molekülen, Polymeren und Festkörpern auf Grundlage der Quantenmechanik lösen zu können um letztlich deren damit verbundene Eigenschaften zu verstehen und erforschen zu können. An der Entwicklung der Grundlagen dieser Theorie waren Hohenberg, Kohn und Sham beteiligt[29, 30]. Die wichtigste Aussage dieser Theorie besteht darin, dass anstelle der hoch-dimensionalen VielteilchenwellenN funktion eines wechselwirkenden-Elektronensystems die Kenntnis der Grundzustandsdichte ausreicht um die physikalischen Eigenschaften des Systems im Grundzustand zu beschreiben. Die Eigenschaften lassen sich als Funktionale der Teilchendichte formulieren, was auch den Namen dieser Theorie erklärt. Die hier präsentierte Abhandlung über die Grundlagen der Dichtefunktionaltheorie orientiert sich an der Darstellung in[17], Kapitel 7. Aus Gründen der Übersichtlichkeit wird die Spinkoordinate in den Wellenfunktionen ausgelassen. 2.2.1. Grundzustand Die Hohenberg-Kohn Theoreme Das System von N wechselwirkenden Elektronen wird durch den Hamiltonoperator H ^ = T ^ + V ^ ee + V ^ ext = N j=1 2 2 m e 2 j + N j=1 j ) an der Mitte des Metalllms( x= 0 ) auntegriert: T N= dt j > ( x= 0, t) 0 (4.38) Die genaue Wahl einer solchen Gröÿe unterliegt einer gewissen Willkür, jedoch erfüllen diese alle die Eigenschaft, dass N gegen Null gehen muss, wenn man die Dämpfungsstärke gegen unendlich gehen lässt. Das Ergebnis der Abhängigkeit von N und der Dämpfungsstärke a 0 ist in Abb. 4.11 gezeigt. Wie man gut erkennen kann, acht die Kurve mit zunehmender Dämpfungsstärke immer mehr ab. Um die Gröÿe N um mehrere Gröÿenordnungen zu reduzieren(wie es für eine Entkopplung der beiden Oberächen nötig wäre), müsste ein so hoher Wert für a 0 gewählt werden, dass es in der numerischen Simulation zu erheblichen Genauigkeitsproblemen kommt. 82 4.4. Analyse der Methoden Abbildung 4.9.: Ein Metalllm wird durch einen Gausspuls angeregt. In dem räumlich und zeitlich begrenzten Intervall wird die reinieÿende Teilchenstromdichte integriert(gestrichelte Linie). 4.4.5. Anmerkungen Die Neuhauser-Methode führt eine gravierende Veränderung in die zeitabhängigen KohnSham Gleichungen ein und ebenso für Ein-Teilchensysteme in die zeitabhängige Schrödingergleichung: Weil im Reibungsterm H ^ f ( t) die Gröÿe d z /dt benötigt wird, und diese implizit vom Zustand des Systems abhängt, werden die Gleichungen. Diese strukturelle Änderung an den Gleichungen ist aber nur dadurch bedingt, dass diese kein abgeschlossenes System mehr beschreiben. Das eigentliche System mit Hamiltonoperator H ^ 0 wird einem externen Feld ausgesetzt, welches basierend auf dem Zustand des Systems so gewählt wird, dass es dem System Energie entzieht. Eine ungewöhnlich erscheinende Eigenschaft dieser Dämpfungsmethode ist die, dass sich die Dynamik auch nach langer Dämpfungszeit vollständig umkehren lässt: Normalerweise würde man von einem Vielteilchensystem erwarten, dass es durch die Streuprozesse, welche die Dämpfung bewirken, an Gedächtnis verliert und der ursprüngliche Zustand nicht mehr durch Zeitumkehr erreicht werden kann. Ein Vergleich mit einem klassischen Gas macht das ebenfalls deutlich: Wenn sich dieses in einem Kolben ausdehnt, wird es sich nicht wieder dadurch zusammenziehen, indem man alle mikroskopischen Impulse der Gasatome umkehrt. 83 4. Dissipative Zeitabhängige Dichtefunktionaltheorie Abbildung 4.10.: Der Plot zeigt die Stromdichte in der Mitte des Metalllms in Abhängigkeit der Zeit für eingeschaltete Dämpfung(blaue Kurve) und ausgeschaltete Dämpfung(rote Kurve). Das gezeigte Intervall umfasst 20 fs . Nach ca. 10 fs ist keine signikante Abnahme in der Stromdichte bei der blauen Kurve mehr zu erkennen. Dieses Phänomen wurde bereits von Neuhauser beobachtet aufgrund dessen er zwei verschiedene Drude Streuzeiten zur Charakterisierung der Dämpfung verwendet hat(s. Kap. 4.4.1). N [a.u.] 0.016 0.014 0.012 0.01 0.008 0.006 0.004 0.002 00 50 100 150 200 250 300 a0[a.u.] Abbildung 4.11.: Die Teilchenanzahl(pro Fläche) aus Gl. 4.38( T= 20 fs ) für verschiedene Dämpfungsparameter a 0 : Die roten Kreuzchen stellen Simulationsergebnisse dar, wohingegen die blaue Linie nur ein angetteter Spline ist. 84 5. Wigner-Maxwell Gleichungen In diesem Teil der Arbeit soll ein quantenmechanisches Modell des Elektronengases entwickelt werden, welches Austausch- und Korrelationswechselwirkungen berücksichtigt. Dazu wird das Elektronengas als statistisches Ensemble über die Dichtematrix beschrieben. Die resultierenden Gleichungen können in ihrer allgemeinsten Form die Dynamik eines Zweikomponentenplasmas(bestehend aus Elektronen und Ionen) in einem elekWigner-Maxwell tromagnetischen Feld beschreiben. Die Herleitung dieser Gleichungen in der hier verwendeten Form wurde bereits in einer unveröentlichten Arbeit von W. Hoyer durchgeführt[56]. In den folgenden Kapiteln wird versucht die Gleichungen auf ein System mit reduzierter Dimensionalität(d.h. speziell einen Nanodraht) anzuwenden, da schnell klar wird, numerisch dass die Gleichungen für drei Raumdimensionen nur schwer lösbar sind[57] und daher für die Anwendungen in dieser Arbeit keine praktische Alternative zur Dichtefunktionaltheorie darstellen. Des weiteren werden die Gleichungen für elektrostatische Wigner-Poisson Felder formuliert, was zu der Gleichung führt, welche in der Literatur bereits unter physikalischen[22, 23, 58] als auch numerischen[59] Aspekten untersucht wurde. Abschlieÿend wird speziell auf die Arbeit[22] eingegangen, in der gezeigt wird, wie sich Quantenkorrekturen in einem klassischen Fluidmodell des Elektronengases herleiten lassen. Die resultierende Gleichung wird als Quanten Euler-Gleichung bezeichnet und stellt eine Alternative zur(Kohn-Sham) Dichtefunktionaltheorie dar. 5.1. Allgemeine Formulierung in drei Dimensionen Die Herleitung der Bewegungsgleichung für das Zweikomponentenplasma erfolgt in mehreren Schritten: Zunächst wird der Hamiltonoperator des Systems in zweiter Quantisierung hergeleitet. Die Quantenfeldoperatoren werden in die Impulsbasis transformiert. Dann wird die Heisenbergsche Bewegungsgleichung für die Kohärenzenmatrix aufgestellt und das dabei auftretende Hierarchieproblem in Hartree-Fock Näherung gelöst. Zuletzt wird die Dichtematrix in das Wignerbild transformiert und das zentrale Ergebnis dieser Herleitung wird die Bewegungsgleichung für die Ein-Teilchen Wignerverteilungsfunktion [60] darstellen. Mit diesen Gleichungen wird die Austauschwechselwirkung exakt beschrieben. Die Berücksichtigung der Korrelation ist in den Gleichungen allerdings derzeit noch oen gelassen. Es sei darauf hingewiesen, dass die Herleitungen in den folgenden Kapiteln keine Zwischenschritte enthalten und aus der Arbeit[56] übernommen wurden. Allerdings ist das selbst entwickelte Modell im Kapitel 5.2, welches auf der Arbeit von Hoyer basiert, mit detaillierter Herleitung gegeben. 85 5. Wigner-Maxwell Gleichungen 5.1.1. Hamiltonoperator Der Hamilton-Operator beschreibt ein System aus zwei Sorten von fermionischen Teilchen, welche über den Index an den Operatoren berücksichtigt werden( = e für Elektronen und = i für Ionen). Die Spinkoordinate s hat zwei Einstellmöglichkeiten s= ± 1/ 2 . Die Teilchen benden sich in einem elektromagnetischen Feld, welches in Coulomb-Eichung über die Potentialfelder A( r, t) und ( r, t) beschrieben wird. Für die Potentiale wird die Coulomb-Eichung verwendet: · A( r, t)= 0 (5.1) Der Hamiltonoperator für die kinetische Energie in minimaler Kopplung lautet zunächst allgemein: H^ min = ,s d 3 r^ ( r, s) 1 2 m ( p ^ q A( r, t)) 2 ^ ( r, s) (5.2) Durch die Eichbedingung kann der Operator in Klammern vereinfacht werden: H^ min = ,s d 3 r^ ( r, s) 1 2 m ( p ^ 2 q A( r, t) · p ^ + q 2 A 2 ( r, t))^ ( r, s) (5.3) Die Coulombwechselwirkung innerhalb einer Teilchensorte wird durch folgenden Hamiltonoperator beschrieben: H^ C = 1 2 d 3 r d 3 r ^ ( r, s)^ ( r , s ) V( | r r | )^ ( r , s )^ ( r, s) s,s (5.4) Für die Coulombwechselwirkung zwischen Elektronen und Ionen lautet der Hamiltonoperator: H^ eCi = d 3 r d 3 r ^ e ( r, s)^ i ( r , s ) V( | r r | )^ i ( r , s )^ e ( r, s) s,s (5.5) Der Hamiltonoperator des Zweikomponentenplasmas lautet insgesamt: H^= H^ kin + H^ A · p + H^ A 2 + H^ eC + H^ iC + H^ eCi (5.6) Dieser Operator hängt vom elektromagnetischen Feld ab, welches von den Maxwellgleichungen beschrieben wird, deren Quellterme sich aus den Ladungen und Strömen des Plasmas ergeben. In der Originalarbeit von W. Hoyer wurde auch das Lichtfeld quantisiert, was aber für die Optik metallischer Nanostrukturen in dieser Arbeit keinen Sinn macht: Im Hinblick auf die nichtlinearen optischen Eigenschaften dieser Strukturen ist klar, dass Felder mit hoher Feldstärke auftreten und die Quantisierung des Lichtfeldes daher nicht beachtet werden muss. 86 5.1. Allgemeine Formulierung in drei Dimensionen 5.1.2. Impulsbasis Die Feldoperatoren haben die folgende Entwicklung in der Impulsbasis: ^ ( r, s)= 1 V e i k · r ^a ,s, k k ^ ( r, s)= 1 V e i k · r ^a ,s, k k (5.7) (5.8) Hier wird implizit angenommen, dass das Plasma ein endliches Volumen V füllt, welches wiederum periodisch den ganzen Raum füllt. Die Operatoren ^a,^a erfüllen die fermionischen Antikommutator-Relationen: [^a ,s, k ,^a ,s , k ] + = k,k s,s , [^a ,s, k ,^a ,s , k ] + = 0 [^a ,s, k ,^a ,s , k ] + = 0 (5.9a) (5.9b) (5.9c) Der Wellenvektor k wird ab jetzt als Verbundindex deniert, in dem auch die Spinquantenzahl enthalten ist: k, s k=( k x , k y , k z , s) Der Operator aus Gl. 5.3 besteht aus drei Teilen, welche in dieser Basis folgendermaÿen lauten: H^ kin = k ^a , k ^a , k H^ A · p = , k J k · A q ^a , k+ q/ 2 ^a , k q/ 2 , k, q H^ A 2 = q 2 2 m k, q, q A q · A q ^a , k+ q ^a , k+ q (5.10) (5.11) (5.12) Dabei wird mit k = 2 k 2 2 m das Matrixelement des kinetische Energieoperators T ^ = p^ 2 / 2 m , mit J das Stromdichtematrixelement J k = q m k (5.13) und mit A q ( t) die Fouriertransformierte des Vektorpotentials A( r, t) bezeichnet. Die Operatoren für die Coulombwechselwirkung(Gln. 5.4 und 5.5) lauten in dieser Basis: H^ C = 1 2 V q ^a , k ^a , k ^a , k + q ^a , k q k, k , q H^ eCi = V q ^a e, k ^a i, k ^a i, k + q ^a e, k q k, k , q (5.14) (5.15) 87 5. Wigner-Maxwell Gleichungen Das Coulombmatrixelement V q ist gegeben durch: V q = 1 0 V | q e | q 2 (5.16) Es gilt für die Ladung der Teilchen | q e | = q i . Damit steht nun der gesamte Hamiltonoperator(Gl. 5.6) auch in der Impulsbasis zur Verfügung, auf der die restliche Herleitung in den folgenden Kapiteln basieren wird. 5.1.3. Bewegungsgleichung der Kohärenzenmatrix Die Kenntnis der Kohärenzen ^a , k ^a , k ( t) erlaubt es, zentrale Gröÿen wie Teilchen- und Stromdichte des Zweikomponentenplasmas zu berechnen. Mit Hilfe der Heisenberggleichung kann die Bewegungsgleichung für diese zentralen Gröÿen hergeleitet werden: i O^=[O^, H^] t Mit O^=^a , k ^a , k und dem Hamiltonoperator (5.17) H^= H^ kin + H^ A · p + H^ A 2 + H^ eC + H^ iC + H^ eCi folgt für die Zeitableitung der Erwartungswerte der Kohärenzen i t ^a , k ^a , k =( k k ) ^a , k ^a , k + A p · ( J k ^a , k+ p ^a , k ) J k ^a , k ^a , k p p 2 qm 2 V p q, ( p A ^a , kp+ q · () A p ^a , k+ p p ^a , k - ^a , k ^a , k p p ) ^a , l ^a , l+ q ^a , k - ^a , k ^a , l ^a , l+ q ^a , k q + q, l V q ( ^a , k+ q ^a ¯ , l ^a ¯ , l+ q ^a , k ) ^a , k ^a ¯ , l ^a ¯ , l+ q ^a , k q q, l (5.18) Als nächstes wird eine Hartree-Fock Faktorisierung auf die 4-Punkt Erwartungswerte angewendet, wobei allerdings die dadurch fehlenden 2-Teilchen Korrelationen als nicht näher spezizierter Term ... mit deniert werden: ^a , k 1 ^a ¯ , k 2 ^a ¯ , k 2 ^a , k 1 = ^a , k 1 ^a ¯ , k 2 ^a ¯ , k 2 ^a , k 1 S + ^a , k 1 ^a ¯ , k 2 ^a ¯ , k 2 ^a , k 1 (5.19a) ^a , k 1 ^a , k 2 ^a , k 2 ^a , k 1 S = ^a , k 1 ^a , k 1 ^a , k 2 ^a , k 2 - ^a , k 1 ^a , k 2 ^a , k 2 ^a , k 1 ( 5.19b) ^a , k 1 ^a ¯ , k 2 ^a ¯ , k 2 ^a , k 1 S = ^a , k 1 ^a , k 1 ^a ¯ , k 2 ^a ¯ , k 2 (5.19c) Mit Hilfe der fouriertransformierten Teilchendichte n q = 1 V ^a , l ^a , l+ q l (5.20) 88 5.1. Allgemeine Formulierung in drei Dimensionen lässt sich die Bewegungsgleichung der Kohärenzenmatrix formulieren: i t ^a , k ^a , k = + ( k k ) ^a , k ^a , k A p · ( J k ^a , k+ p ^a , k ) J k ^a , k ^a , k p p + q 2 2 m A p p, p V V q ( n q ( · A p ^a , k+ p p ^a , k - ^a , k ^a , k p+ p () n q ¯ ) ^a , k ^a , k q - ^a , k+ q ^a , k ) q k, k + k ,, k k ,, k ¯ (5.21) Die letzten drei Terme beschreiben Austausch( ) und Korrelation( ): k, k = V q ( ^a , k ^a , l ^a , l q ^a , k q - ^a , k+ q ^a , l+ q ^a , l ^a , k ) k ,, k ¯ = q, l V q ( ^a , k ^a ¯ , l ^a ¯ , l+ q ^a , k q ^a , k+ q ^a ¯ , l ^a ¯ , l+ q ^a , k ) q, l (5.22) (5.23) 5.1.4. Formulierung im Wignerbild Die Dichtematrix ^ , welche das statistische Ensemble von Zweikomponentenplasmen beschreibt, kann über den gesamten Satz von N reduzierten Wignerfunktionen dargestellt werden. Diese sind als Analogon zu den reduzierten Verteilungsfunktionen der klassischen statistischen Mechanik zu betrachten. Allerdings können die Wignerfunktionen auch neQuasiverteilungen gative Werte annehmen[58] und werden daher auch als bezeichnet. Die Ein-Teilchen Wignerfunktion hängt folgendermaÿen mit der Dichtematrix zusammen: f k ( r)= e i q · r ^a , k q/ 2 ^a , k+ q/ 2 q (5.24) ^a , k ^a , k = 1 V e i( k k) · r f k+ k ( r) d 3 r 2 (5.25) Die Teilchen- und Geschwindigkeitsdichtefelder können sehr einfach aus der Wignerverteilung bestimmt werden: n ( r)= 1 V f k ( r) k v k an ( r)= 1 V k k m f k ( r) v k in ( r)= 1 V k k q A( r) m f k ( r) (5.26) (5.27) (5.28) 89 5. Wigner-Maxwell Gleichungen Das kinetische Geschwindigkeitsdichtefeld tritt in der Kontinuitätsgleichung der Teilchendichte auf: t n ( r)+ · v k in ( r)= 0 (5.29) Mit einer umfangreichen Rechnung(s. Anhang C von[56]), bei der eine Gradientenentwicklung der Wignerfunktion vorgenommen wird, lässt sich zeigen, dass die Bewegungsgleichung der Kohärenzenmatrix 5.21 im Wignerbild folgende Form hat: () t f k ( r) kin t f k ( r) A · p = = r · k m f k ( r) (5.30a) q m [ r · ( A( r) f k ( r)) ( r ( A( r) · k)) · ] ( k f k ( r)) (5.30b) t f k ( r) A 2 = 2 q 2 m ( r | A( r) | 2 ) · ( k f k ( r)) (5.30c) t f k ( r) C, L = q ( r ( r)) · ( k f k ( r)) (5.30d) t f k ( r) C, HF = | q | ([ k k ( r)] · [ r f k ( r)] [ r k ( r)] · [ k f k ( r)]) (5.30e) t f k ( r) C, Korr = 1 q | q | [ r U( r r )] · [ k f ^ k ( r)^ ( r ) ] d 3 r (5.30f) Die Gradientenentwicklung hat keinen Einuss auf den Term in Gleichung 5.30a. Alle anderen Terme sind Näherungen. In der Coulomb-Eichung gilt für das skalare Potential ( r) in Gl. 5.30d: ( r)= U( r r )( n i ( r ) n e ( r )) d 3 r (5.31) Die Energierenormierung in Gl. 5.30e ist durch eine Faltung im k -Raum gegeben: k ( r)= U k k f k ( r) k (5.32) Für das Coulombmatrixelement im Orts- und Fourierraum gilt: U( r)= 1 4 0 q i , r U q = 1 0 V q i q 2 (5.33) Im Korrelationsterm(Gl. 5.30f) ist der Operator f ^ durch den Operator in den Erwartungswertklammern aus Gl. 5.24(inklusive Summation und dem Vorfaktor e i q · r ) deniert. Der Operator ^ ist der Ladungsdichteoperator und setzt sich aus den Teilchendichteoperatoren zusammen: ^ ( r)= q i ( n^ i ( r) n^ e ( r)) (5.34) 90 5.1. Allgemeine Formulierung in drei Dimensionen 5.1.5. Gleichungssystem für ein Zweikomponentenplasma Im letzten Abschnitt wurden die Herleitung der Wigner-Gleichungen für ein Zweikomponentenplasma nach[56] skizziert. Hinzugenommen werden müssen noch die MaxwellGleichungen, welche in Gl. 5.31 nur unvollständig berücksichtigt wurden, da durch diese Gleichung nur das longitudinale elektrische Feld in die Dynamik der Wignerverteilung eingeht. Das vollständige Wigner-Maxwell Gleichungssystem umfasst die Bewegungsgleichung für f k ( r) , die Eichbedingung und die Wellengleichung für das Vektorpotential sowie die Poissongleichung für das skalare Potential: t f k ( r)= t f k ( r) kin + t f k ( r) C, L + t f k ( r) A · p + t f k ( r) A 2 + t f k ( r) C, HF + t f k ( r) C, Korr · A( r)= 0 2 A( r)= 1 c 2 2 t 2 A( r) µ 0 j ( T) ( r) 2 ( r)= 1 ( r) 0 (5.35a) (5.35b) (5.35c) (5.35d) (5.35e) (5.35f) (5.35g) (5.35h) Dabei bezeichnet j (T) den transversalen Anteil der Stromdichte. Die Auistung der Terme 5.35a-5.35e hat folgenden Hintergrund: Bei der Anwendung dieser Gleichungen auf ein bestimmtes System kann es sinnvoll sein, nicht alle Terme zu berücksichtigen. In Systemen, wo beispielsweise die Retardierung keine Rolle spielt, kann die Termgruppe 5.35c aus der Bewegungsgleichung für f k ( r) und die Gleichungen 5.35f-5.35g aus dem Gleichungssystem gestrichen werden. Ebenso kann durch Ein/Ausschalten der Terme ermittelt werden, welche Bedeutung diese für die Dynamik des Plasmas in einem bestimmten System haben. 91 5. Wigner-Maxwell Gleichungen 5.1.6. Anwendung auf metallische Nanostrukturen Die Wigner-Maxwell Gleichungen sollen nun verwendet werden, um die gesamten Metallelektronen einer metallischen Nanostruktur im Lichtfeld zu beschreiben. Die Herleitung dieser Gleichungen im Fourierraum impliziert, dass die zu beschreibenden Strukturen im periodisch realen Raum angeordnet sind. Diese Periodizität kommt über den Basiswechsel bei den Feldoperatoren 5.7,5.8 und der(dazu passenden) Entwicklung des Vektorpotentials in ebenen Wellen A( r) A q ins Spiel. Um eine Isolation der Strukturen von ihren Nachbarn zu erzielen, muss das Zellvolumen V , in dem je eine Struktur enthalten ist, entsprechend groÿ gewählt werden, bis beispielsweise die longitudinalen elektrischen Felder über die Distanz im Vakuum genügend abgeklungen sind und keinen Einuss auf die Bewegung der Elektronen in der Nachbarzelle mehr haben. Den Eekt der transversalen Felder auf diese Weise zu verhindern ist zwar auf gleiche Weise möglich, aber da diese nur mit 1 /r im Raum abfallen, nur durch sehr groÿe Zellvolumen erreichbar. Dies kann in der numerischen Lösung dieser Gleichungen zu einem zentralen technischen Problem werden. Dabei wird dann auch klar deniert, was unter einem groÿen Zellvolumen zu verstehen ist, da verfügbarer Speicherplatz und Rechenzeit die Grenzen festlegen. Für die metallischen Nanostrukturen soll hier wieder das Jellium-Modell zum Einsatz kommen, d.h. die Ionendichte n i ( r) ist zeitlich konstant und legt die Form der Struktur fest. In den Wignergleichungen kann deshalb der Index weggelassen werden, da ab jetzt nur noch das Elektronengas von diesen Gleichungen beschrieben wird. Entsprechend gilt nun: e ¯ i Wie oben bereits angedeutet, können die Wigner-Maxwell Gleichungen in unterschiedlicher Komplexität zum Einsatz kommen. Ein Ziel dieser Arbeit besteht darin zu ermitteln, welche Terme unter welchen Umständen weggelassen werden können. Im Hinblick auf die numerische Umsetzung dieser Gleichungen werden verschiedene Ausbaustufen betrachtet, die in Tabelle 5.1.6 aufgeführt sind. Als Relaxation wird hier ein Prozess, der die Elektronen abbremst und asymptotisch im elektronischen Grundzustand endet, bezeichnet. In Festkörpern spielt sich dieser Prozess in Form von Streuung der Elektronen untereinander und mit Gitterphononen ab. Der Korrelationsterm C,Korr hat die Aufgabe, die Streuung der Elektronen untereinander zu beschreiben. Da im Jellium-Modell allerdings kein Kristallgitter und somit auch keine Gitterphononen existieren, kann dieser Relaxationsprozess nur unvollständig von diesen Gleichungen beschrieben werden. Für die Berechnung der Grundzustandsdichte ist ein solcher Relaxationsmechanismus aber unverzichtbar: Diese Dichte lässt sich nämlich extrem einfach durch Lösung der Gleichungen im Zeitbereich ermitteln, indem zunächst eine ungefähre Startdichte der Elektronenverteilung vorgegeben wird(z.B. n e ( r, t= 0)= n i ( r) ist eine mögliche Startdichte). Anschlieÿend wartet man, bis sich ein Endzustand der Dichte eingestellt hat. Für eine Simulation dieses Vorgangs wäre ein phänomenologisch beschriebener Relaxationsmechanismus(wie die Drude-Streuung in Metallen) ideal, da die Dämpfungskonstante künstlich hoch gewählt und somit die nötige Simulationsdauer stark reduziert werden 92 5.1. Allgemeine Formulierung in drei Dimensionen Stufe 1 2 3 4 5 6 Terme kin+ C,L kin+ C,L+ Relaxation kin+ C,L+ C,HF+ Relax. kin+ C,L+ C,HF+ A · p + A 2 kin+ C,L+ C,HF+ C,Korr (alle Terme) Anwendung Einuss statischer Felder untersuchen Berechnung der Grundzustandsdichte Einuss der Austauschwechselwirkung auf die Grundzustandsdichte Elektronendynamik im Lichtfeld ( A 2 -Term vernachlässigbar?) Berücksichtigung von Streuprozessen zwischen den Elektronen Validierung von und Vergleich mit anderen Modellen Tabelle 5.1.: Unterschiedlich komplexe Ausbaustufen der Wigner-Maxwell Gleichungen kann. Der Endzustand muss natürlich unabhängig von diesem Vorgang und dessen Parametrisierung sein. Ansätze dazu werden im nächsten Abschnitt(5.1.7) gegeben. Diese Methodik wurde ebenfalls bei der Berechnung der Grundzustandsdichte im Hydrodynamikmodell(s. Abb. 3.2) verwendet. 5.1.7. Phänomenologische Relaxationsterme Der erste sinnvolle Ansatz besteht darin kleine Abweichungen der Verteilungsfunktion vom Grundzustand zu betrachten: f k ( r, t)= f k (0) ( r)+ f k ( r, t) (5.36) Im ungestörten System muss diese Abweichung zeitlich asymptotisch gegen Null gehen. Daher kann in die Bewegungsgleichung für f einfach ein Term der Form f k ( r, t), (5.37) wobei eine Dämpfungskonstante ist, hinzugefügt werden. Leider kann diese Methode oensichtlich nicht dafür verwendet werden, die Verteilungsfunktion des Grundzustands f (0) überhaupt erst zu ermitteln. Wenn diese Verteilung aber erst mal vorliegt, bietet der Term 5.37 die einfachst denkbare Möglichkeit an eine gedämpfte Dynamik für die Verteilungsfunktion zu gelangen. Ein weiterer Ansatz zur Dämpfung basiert auf der Annahme, dass im ungestörten und gedämpften System das Geschwindigkeitsdichtefeld asymptotisch zeitlich gegen Null gehen muss. Einen nicht verschwindenden Beitrag zum Geschwindigkeitsfeld an einem Ort r erhält man immer dann, wenn die Verteilungsfunktion im k -Raum an diesem Ort unsymmetrisch ist. Die Idee besteht daher darin, eine symmetrische Verteilung asymptotisch zu erzwingen, indem die Bewegungsgleichung für f um folgenden Term erweitert wird: [ f k ( r, t) f k ( r, t)] (5.38) Wenn eine Asymmetrie in der Verteilung vorliegt, wird diese exponentiell zeitlich mit Zeitkonstante 1 abgebaut. 93 5. Wigner-Maxwell Gleichungen Abbildung 5.1.: Nanodraht mit Radius a . Der Querschnitt zeigt das Betragsquadrat der Wellenfunktion 0 ( r ) , die zum niedrigsten Subband gehört. Die Elektronen können sich durch die Beschränkung auf das unterste Subband in der yz -Ebene nicht bewegen. 5.2. Nanodrähte In diesem Abschnitt sollen die Wigner-Maxwell Gleichungen für metallische Nanodrähte entwickelt werden. Der Draht erstrecke sich entlang der x -Achse und besteht aus periodisch zusammengesetzten Elementen der Länge L , die durch eine Ionendichte n i ( r) charakterisiert sind. Die Bewegung in der Querschnittebene zum Draht ist quantisiert. Zur Vereinfachung des Modells wird nur das energetisch niedrigste Subband betrachtet, was letztlich dazu führt, dass keine Bewegung in der Ebene möglich ist. Die Gleichungen, die in den nächsten Unterkapiteln hergeleitet werden, sollen vor allem in numerischen Simulationen benutzt werden können. Um die Komplexität möglichst gering zu halten, werden hier nur statische Coulombfelder berücksichtigt und das Vektorpotential des Lichtfeldes wird ausser Acht gelassen. Um elektronische Anregungen im Draht zu erzeugen, dient ein eindimensionales Störpotential v P im Draht. Dessen Wirkung entspricht dem eines elektrostatischen Feldes entlang des Drahtes mit einer so schwachen Zeitabhängigkeit, dass Retardierungseekte nicht beachtet werden müssen. In diesem Zusammenhang sollten Wigner-Poisson Gleichungen, die in den nächsten Kapiteln hergeleitet werden, besser als Gleichungen bezeichnet werden. Es besteht ein wesentlicher Unterschied zwischen den hier betrachteten Nanodrähten und jenen, welche in Kapitel 3.3 über die DFT beschrieben werden: Hier wird nur die Bewegung der Elektronen entlang des Drahtes beschrieben, wohingegen die DFTRechnungen nur die Bewegung in der Querschnittsebene berücksichtigt haben. Lässt man den Radius a des Drahtes gegen Unendlich gehen, sollten sich die Berechnungen mit den DFT-Rechnungen von Metalllmen(s. Kap. 3.1) vergleichen lassen. 94 5.2. Nanodrähte 5.2.1. Hamiltonoperator Der Hamiltonoperator für den Nanodraht ist ähnlich wie der in Gl. 5.6 aufgebaut. Die Bestandteile sind selbsterklärend: 2 H^ kin H^ C = = d 3 r^ ( r, s) 2 ^( r, s) (5.39a) 1 s d 3 r 2 m e d 3 r ^ ( r, s)^ ( r , s ) U( | r r | )^( r , s )^( r, s) (5.39b) 2 s,s H^ eCi = d 3 r d 3 r ^ ( r, s) n i ( r ) U( | r r | )^( r, s) (5.39c) s, s H^ P = d 3 r^ ( r, s) v P ( r, t)^( r, s) (5.39d) s H^= H^ kin + H^ eC + H^ eCi + H^ P (5.40) Der Draht ist zunächst wie alle physikalischen Systeme dreidimensional. Im nächsten Kapitel wird gezeigt, wie die Symmetrie des Systems zur Vereinfachung der Gleichungen genutzt werden kann. 5.2.2. Eindimensionale Impulsbasis Es wird zunächst ein vollständiger Satz von orthonormalen Funktionen benötigt, mit denen der Nanodraht besonders einfach beschrieben werden kann. Die Zylindersymmetrie und die unendliche( L -periodische) Ausdehnung entlang der x -Achse legen folgende Form der Wellenfunktionen nahe: n,k, ( x, r , s)= 1 L n ( r ) e ikx s, , d 2 r n 1 n 2 = n 1 ,n 2 (5.41) Das Zahlentupel n enthält bei Zylindersymmetrie die Quantenzahlen ( l, m) und die Quank tenzahl charakterisiert die ebene Welle entlang des Drahtes. Für diese Funktionen können z.B. die Eigenfunktionen des eektiven Potentials des Drahtes gewählt werden. Die Quantenfeldoperatoren können in dieser Basis entwickelt werden: ^( r, s)= 1 L n,k n ( r ) e ikx ^a n,k,s ^ ( r, s)= 1 L n,k n ( r ) e ikx ^a n,k,s (5.42a) (5.42b) Die Entwicklung kann nun in die Operatoren 5.39a-5.39d eingesetzt werden. Dabei wird von der bereits erwähnten Beschränkung auf das unterste Subband Gebrauch gemacht. Die folgende Rechnung zeigt exemplarisch, wie diese im Fall des kinetischen Anteils H^ kin 95 5. Wigner-Maxwell Gleichungen aussieht: 2 H^ kin = s d 3 r^ ( r, s) 2 ^( r, s) 2 m e 2 1 = 2 m e L s n 1 ,k 1 n 2 ,k 2 () d 3 r n 1 e ik 1 x ( ik 2 ) 2 e ik 2 x n 2 + e ik 2 x r 2 n 2 ^a n 1 ,k 1 ,s ^a n 2 ,k 2 ,s Beschränkung auf unterstes Subband( n 1 , n 2 0 ): H^ kin = = 2 2 m e 2 2 m e 1 L 1 L s s ( ) d 3 r 0 e ik 1 x ( ik 2 ) 2 e ik 2 x 0 + e ik 2 x r 2 0 k 1 ,k 2 [ dxe ik 1 x ( ik 2 ) 2 e ik 2 x d 2 r 0 0 k 1 ,k 2 ] ^a k 1 ,s ^a k 2 ,s + dxe ik 1 x e ik 2 x d 2 r 0 r 2 0 ^a k 1 ,s ^a k 2 ,s = 2 2 m e 1 L s k 1 ,k 2 [ ( ik 2 ) 2 L k 1 ,k 2 + L k 1 ,k 2 T ] ^a k 1 ,s ^a k 2 ,s = 2 2 m e ( k 2 sk + T )^a k,s ^a k,s = k ^a k,s ^a k,s sk mit k := 2 [ k 2 T ] 2 m e Der orthogonale Anteil 0 der Basisfunktionen führt zu einer Konstante T , die in die Gröÿe k integriert wurde. Wenn man statt des Drahtes einen Metalllm beschreiben will, gilt 0 = 1 und die Konstante T verschwindet. Für die restlichen Operatoren 5.39b-5.39d zeigt sich ebenfalls, dass die Berücksichtigung des orthogonalen Anteils 0 über Korrekturen an den Matrixelementen, welche aus Herleitungen für ein eindimensionales System stammen, möglich ist. Letztlich bedeutet das, dass für den Nanodraht anstelle der Entwicklung 5.42a-5.42b für die Feldoperatoren direkt die eindimensionale Impulsbasis gewählt werden kann: ^( x, s)= 1 L k e ikx ^a k,s ^ ( x, s)= 1 L k e ikx ^a k,s (5.43a) (5.43b) Die Spinquantenzahl s wird ab hier gelegentlich im Verbundindex k untergebracht, wenn es der Übersichtlichkeit dient: ( k, s) k 96 5.2. Nanodrähte Da sich das System L -periodisch entlang der x-Achse fortsetzt, sind folgende diskrete Werte für die Quantenzahl k erlaubt: {} k 2 n, n Z L (5.44) Durch Einsetzen der Feldoperatoren 5.43a-5.43b lassen sich nun die eindimensionalen Versionen der Operatoren 5.39a-5.39d formulieren: H^ kin = k ^a k ^a k k mit k := 2 ( k 2 T ) 2 m e H^ C = 1 2 U q ^a k ^a k ^a k + q ^a k q k,k ,q H^ eCi = n i ( q) U q ^a k ^a k+ q H^ P ( t)= k,k ,q v P ( q, t)^a k ^a k+ q k,q (5.45a) (5.45b) (5.45c) (5.45d) Die Form des Drahtes(d.h. der 0 -Anteil) steckt in den Matrixelementen U q (s. Kap. 5.2.3). Bei dem Matrixelement v P fällt dagegen der 0 -Anteil raus(s. Kap. 5.2.4). 5.2.3. Coulombmatrixelemente Es wird von der dreidimensionalen Yukawa-Form der Coulombwechselwirkung ausgegangen : v( r, )= q e2 e r 4 0 r (5.46) Die folgende Herleitung zeigt nur die wesentlichen Schritte(siehe z.B. Anhang 1 in[17]). Mit | n, k, wird ein Zustand bezeichnet, dessen Wellenfunktion die Form aus Gleichung 5.41 hat. Wegen der Beschränkung auf das unterste Subband wird zur Abkürzung | k, := | 0, k, deniert. Unter Berücksichtigung der Impulserhaltung lässt sich das Matrixelement U q über folgenden Ausdruck berechnen: U q ( )= q e2 1; 4 0 k+ q, | 2; k q, | e r^ 12 r^ 12 | 2; k, | 1; k, , r^ 12 = | ^ r 1 ^ r 2 | (5.47) An dieser Stelle kann von der Fouriertransformation von e r 12 /r 12 Gebrauch gemacht werden: e r 12 = r 12 d 3 q 2 3 q 2 4 + 2 e i q · ( r 1 r 2 ) (5.48) 1 Die Verwendung einer abgeschirmten Coulombwechselwirkung hat sich bei den DFT-Rechnungen als nützliches Werkzeug erwiesen(s. Kap. 6.1.7) und soll daher hier ebenfalls berücksichtigt werden können. 97 5. Wigner-Maxwell Gleichungen Das Matrixelement ist dann gegeben durch: U q ( )= q e2 4 0 d q (2 ) 2 q 2 + 4 q 2 + 2 | F( q ) | 2 F( q )= d 2 r | 0 ( r ) | 2 e i q · r (5.49) Die Gröÿe | F( q ) | 2 , welche die Dichte | 0 ( r ) | 2 im Querschnitt beschreibt, wird als Formfaktor bezeichnet. Für einen Draht mit Radius a und gaussförmiger Wellenfunktion im Querschnitt(s. Abb. 5.1) erhält man: 0 ( r )= 2 a 2 e r 2 /a 2 U q ( )= q e2 e ( q 2 + 2 ) a 2 Ei( [ q 2 + 2 ] a 2 ) 4 0 (5.50a) (5.50b) Bei Ei( x) handelt es sich um die Exponential-Integral Funktion. 5.2.4. Matrixelemente des Störpotentials Das Störpotential sei im Ortsraum nur von x und t abhängig und durch eine L -periodische Funktion f p ( x, t) gegeben: v p ( r, t)= f p ( e x · r, t) (5.51) Um zu zeigen, dass dadurch der 0 -Anteil der Wellenfunktion keinen Einuss auf das Matrixelement v P hat, wird nochmals bei dem Operator in Gl. 5.39d angesetzt: H^ P = d 3 r^ ( r, s) v P ( r, t)^( r, s) = s 1 L d 3 r 0 ( r ) e ik 1 x f p ( e x · r, t) 0 ( r ) e ik 2 x ^a k 1 ,s ^a k 2 ,s = 1 s k 1 k 2 L dx e i[ k 1 k 2 ] x f p ( x, t)^a k 1 ,s ^a k 2 ,s = s k 1 k 2 v p ( k 1 k 2 , t)^a k 1 ,s ^a k 2 ,s = s k 1 v P ( q k 2 , t)^a k ^a k+ q ( Verbundindexschreibweise ) k,q Das Matrixelement v p ( q, t) ist dabei als folgender Fourierreihen-Koezient gegeben: v p ( q, t)= 1 L L/ 2 L/ 2 f p ( x, t) e iqx dx (5.52) In einer numerischen Simulation des Nanodrahtes würde man vermutlich f p im Ortsraum vorgeben wollen, so dass man die Matrixelemente v p ( q, t) mittels(diskreter) Fouriertransformation berechnen muss. 2 Def.: Ei( x)= x dt e t t , s.[61] 98 5.2. Nanodrähte 5.2.5. Bewegungsgleichung der Kohärenzenmatrix Die Herleitung der Bewegungsgleichung ist in Anhang A.5.1 gegeben. Die Gleichung lautet in Hartree-Fock Näherung: i t ^a k ^a k ( k k ) ^a k ^a k + k,k { U q L[ n e ( q) n i ( q)]+ v P ( q, t) } ( ^a k ^a k q - ^a k+ q ^a k ) (5.53) q Der Austauschterm k,k wird in Gl. A.32 deniert. 5.2.6. Observablen Für ein eindimensionales System wie dem Nanodraht besteht die Möglichkeit, die Kohärenzenmatrix anstelle der Wignerverteilung in einer numerischen Simulation zu betrachten. Damit können z.B. die Näherungen, welche bei der Transformation in das Wignerbild (Kap. 5.2.8) gemacht werden, untersucht werden. Daher werden in diesem Abschnitt einige der wichtigsten Observablen aufgeführt, welche ansonsten aus den Momenten der Wignerverteilung berechnet werden würden: · Teilchendichte: n^( r) = ^ ( r)^( r) = e iqx n^( q) q ( ) = 1 e iqx L ^a l ^a l+ q ql (5.54) · Stromdichte: ^j P ( q) = m ( k+ q 2 ) ^a k q ^a k k (5.55) Da kein Vektorpotential A vorhanden ist, ergibt sich nur ein paramagnetischer Strom und der diamagnetische Anteil entfällt. · Gesamtenergie: E ges = H ^ = T ^ + V ^ = k ( 2 k 2 2 m e ) ^a k ^a k + V ^ e, e + V ^ e, i (5.56) Die potentielle Energie innerhalb des Elektronengases ist gegeben durch eine Summe von 4-Punkt Erwartungswerten. Diese werden in Hartree-Fock Näherung ap99 5. Wigner-Maxwell Gleichungen Abbildung 5.2.: Ausschnitt aus der Kohärenzenmatrix ^a i ^a j . Die Indizes i, j stehen für Zeilen- und Spaltenindex. Mit der Hilfsgröÿe q ( k) (s. Gl. 5.60) bewegt man sich auf der q -ten Nebendiagonalen. Mit q> 0 bewegt man sich innerhalb der unteren- und mit q< 0 auf der oberen Dreiecksmatrix. proximiert: V ^ e, e = = = 1 2 U q ^a k,s ^a k ,s ^a k + q,s ^a k q,s s,s k,k ,q 1 2 U q ^a k,s ^a k ,s ^a k + q,s ^a k q,s S s,s k,k ,q 1 2 () U q ^a k,s ^a k q,s ^a k ,s ^a k + q,s - ^a k,s ^a k + q,s ^a k ,s ^a k q,s s,s k,k ,q L 2 U q n^( q) n^( q) 1 2 U q ^a k,s ^a k + q,s ^a k ,s ^a k q,s (5.57) q q s,s k,k Die potentielle Energie zwischen beiden Teilchensorten ist gegeben durch: V ^ e, i = U q L n^( q) n i ( q) q (5.58) 5.2.7. Eigenschaften der Kohärenzenmatrix Die Eigenschaften der Kohärenzenmatrix sind für numerische Simulationen, welche direkt bei der Bewegungsgleichung 5.53 dieser Matrix ansetzen, besonders relevant. Als erstes soll untersucht werden, wie man eine Teilchendichte n ( x)= n^( x) mit verschwindender Stromdichte vorgeben kann, da eine solche Konguration als Anfangsbedingung dienen kann. Das Verschwinden der Stromdichte wird durch folgende Gleichung ausgedrückt: 0= ^ j P ( q) = m ( k+ q ) 2 a k q a k q k (5.59) Dazu wird eine Hilfsgröÿe deniert: q ( k):= ^a k q ^a k (5.60) 100 5.2. Nanodrähte Der Index q gibt die Nebendiagonale an entlang der man sich mit dem Argument k bewegt. Die Bedingung 5.59 für Erwartungswert der Stromdichte lautet unter Verwendung dieser Hilfsgröÿe: 0= ( q ) m k+ 2 q ( k)( k k+ q/ 2) k = m k q ( k q/ 2) [ k ] = m k q ( k q/ 2) k q ( k q/ 2) k> 0 k> 0 (5.61) Die Hilfsgröÿe q ( k) muss für die Forderung ^ j P ( q) = 0 daher folgende Bedingung erfüllen: q ( q/ 2+ )= q ( q/ 2 ) (5.62) In Abbildung 5.2 benden sich die Matrixelemente mit = 0 auf der rot gestrichelten Diagonalen. Anschaulich liegt die Bedingung ^ j P ( q) = 0 dann vor, wenn die q -te Nebendiagonale symmetrisch bezüglich des Matrixelementes ist, dass sich auf dieser Diagonalen bendet. Für die Fouriertransformierte der Ladungsdichte(s. Gl. 5.54) gilt: n^( q) = 1 L ^a l q ^a l = 1 L q ( l) ll (5.63) Die Summe über die q -te Nebendiagonale ergibt die Ladungsdichte. Damit ist nun klar, wie im Prinzip ein Zustand mit verschwindender Stromdichte und vorgegebener Ladungsdichte in der Kohärenzenmatrix vorgegeben werden kann. In einer zeitabhängigen Simulation des Nanodrahtes würde man zur Zeit t= 0 mit einem elektronischen Grundzustand starten, der die Eigenschaften ^ j P ( q) = 0 als solcher besitzen muss. Das Problem, diesen Zustand zu bestimmen, lässt sich mathematisch über ein Variationsproblem formulieren: E ges [ n e ( x)]= 0, Nebenbed.: n e ( x) dx= N el = N ion (5.64) Die Energie E ges ist nach Gl. 5.56 zu berechnen. Die Variation enthält dabei zwei Variationsschritte: Zunächst wird eine Ladungsdichte, die mit der Nebenbedingung kompatibel ist, ausgewählt. Dann werden die q -ten Nebendiagonalen(inkl. Hauptdiagonale) unter Berücksichtigung der Bedingung 5.63 so variiert, dass die Energie minimal ist. Für die Numerik ist die Symmetrie der Kohärenzenmatrix des weiteren noch wichtig, damit keine redundanten Informationen gespeichert werden müssen. Für den Erwartungswert eines Operators b^ gilt b^ = b^ und daher gilt für die Matrixelemente: (^a k ^a l ) = ^a l ^a k = ^a k ^a l (5.65) Die Matrix ist also hermitesch und es muss nur eine Dreieckshälfte gespeichert werden. Die Besetzungen ^a k ^a k (d.h. die Diagonalelemente) müssen reellwertig sein und aufgrund des Pauli-Prinzips die Bedingung 0 ^a k ^a k ^a k,s ^a k,s 1 erfüllen. 101 5. Wigner-Maxwell Gleichungen 5.2.8. Transformation in das Wignerbild Die Herleitung ist in Kapitel A.5.2 gegeben. Dazu werden die Hin- und Rücktransformationen für die Wignerfunktion benötigt: f( x, k)= e iqx ^a k q/ 2 ^a k+ q/ 2 q ^a k ^a k = 1 L e i( k k) x f( x,[ k+ k ]/ 2) dx (5.66a) (5.66b) Die Bewegungsgleichung der Wignerverteilung lautet insgesamt: f( x, k)= f( x, k)+ f( x, k)+ f( x, k)+ f( x, k) t t kin t H t F t P k f( x, k) m e x + q e ( x) · f( x, k) x k + | q e | ([ k k ( x)] · [ x f( x, k)] [ x k ( x)] · [ k f( x, k)]) 1 + x v P ( x) · k f( x, k) (5.67) Aufgrund der Gradientenentwicklung(Gl. A.38) im Hartree-, Fock- und Störterm handelt es sich um eine Näherung ! der Bewegungsgleichung. 3 Das Fehlen des Korrelationsterms stellt natürlich eine weitere Näherung dar. 102 5.3. Wigner-Poisson System 5.3. Wigner-Poisson System In Kapitel 5.3.1 wird die Herleitung der eindimensionalen Wigner-Poisson Gleichungen skizziert. Dieses fundamentale System von Gleichungen stellt eine starke Vereinfachung der Wigner-Maxwell Gleichungen dar, bei der Austausch und Korrelation der Teilchen sowie die Retardierung vernachlässigt werden. In Kapitel 5.3.2 wird das Quantenmechanische Hydrodynamik-Modell(QHD) aus[22, 23] vorgestellt und dessen Stärken und Schwächen im Vergleich zu den Wigner-Poisson Gleichungen und DFT-Rechnungen erläutert. 5.3.1. Herleitung Ein Ensemble von M Systemen mit N fermionischen Teilchen wird durch eine Dichtematrix ^ beschrieben: M ^ ( t)= p | ( t) ( t) | =1 (5.68) ( r 1 ,..., r N , t)= r 1 ,..., r N | ( t) (5.69) Die N -Teilchen-Wignerverteilung hat 2 N+ 1 Argumente: f N ( x 1 , v 1 ,.., x N , v N , t)= N ( m e 2 ) N M p d N s ( exp im e N i=1 ) v i s i =1 × ( x 1 + s 1 / 2,..., x N + s N / 2, t) × ( x 1 s 1 / 2,..., x N s N / 2, t) (5.70) Diese hat folgende Normierungseigenschaft: d N x d N v f N ( x 1 , v 1 ,.., x N , v N , t)= N (5.71) Um Systeme von identischen Fermionen zu beschreiben, müssen die N -Teilchen Wellenfunktionen die entsprechende Antisymmetrie aufweisen. Über die Schrödingergleichung kann nun die Dynamik der Verteilung hergeleitet werden: i t = 2 2 m e N j=1 2 x 2 j + V( x 1 ,.., x N ) (5.72) Die potentielle Energie soll speziell folgende Form haben, wie sie für Systeme mit Coulombwechselwirkung vorliegt: V( x 1 ,.., x N )= W( | x i x j | ) i F (Thomas-Fermi Abschirmlänge, s. Gl. 1.10) gültig. Daher zeigen die Ergebnisse für die elektronische Grundzustandsdichte eines Metalllms(s. Abbildungen 1 und 2 in[65]) auch keine Friedel-Oszillationen (s. Kap. 3.1.2). · Das Phänomen der Landau-Dämpfung, welches in den Wigner-Poisson Gleichungen noch enthalten ist[59], fehlt in der makroskopischen Beschreibung. Als zentrale Verbesserung des Modells wird in[65] noch das Hinzufügen eines Dämpfungsterms, wie er bereits in Gl. 1.6b enthalten ist, vorgeschlagen. 107 6. Numerik 6.1. Lösungsverfahren für stationäre Kohn-Sham Gleichungen 6.1.1. Problemübersicht Die zeitunabhängigen Kohn-Sham Gleichungen haben aus numerischer Sicht folgende Struktur: { H ^ KS [ n ] j ( r)= j j } N ( r) j=1 N n ( r)= w j | j ( r) | 2 j=1 (6.1a) (6.1b) Bei H ^ KS [ n ] handelt es sich um einen linearen Operator, welcher allerdings in nichtlinearer Form von der Teilchendichte n ( r) abhängig ist. Die Lösung dieser Gleichungen hat die Eigenschaft, dass die Eingabedichte n in Gl. 6.1a mit der aus den Orbitalen resultierenden Dichte in Gl. 6.1b übereinstimmt. Dieses Problem wird numerisch folgendermaÿen gelöst: Für eine vorgegebene Eingabedichte n stellt die Gleichung 6.1a ein gewöhnliches Eigenwertproblem dar. Eine Diskretisierung mittels Finiter Dierenzen macht aus dem Operator H ^ KS [ n ] eine endlich groÿe Matrix und aus der Orbitalfunktion j einen Vektor: { H KS [ n ] y j = j y j } N j=1 (6.2) Es werden die N niedrigsten Eigenwerte und Eigenfunktionen der Matrix H KS [ n ] benötigt. Dieses numerische Eigenwertproblem kann im Prinzip mit einer breiten Palette von Algorithmen[66] gelöst werden. Die Matrix H KS [ n ] verfügt aber über bestimmte Eigenschaften, welche nur den Einsatz ganz spezieller Methoden erlaubt. Eine der derzeit geeignetsten Methoden wird in Kapitel 6.1.2 beschrieben. Wenn das Eigenwertproblem 6.1a für feste Eingabedichte gelöst werden kann, besteht die nächste Aufgabe darin, genau die Dichte zu nden, welche beide Gleichungen 6.1a und 6.1b gleichzeitig erfüllt. Ein geeignetes Verfahren wird in Kapitel 6.1.5 präsentiert. Beide Teilprobleme sind Gegenstand aktueller Forschung, wie an den Referenzen in den nachfolgenden Kapiteln zu sehen sein wird. Hier gilt es einen Kompromiss zwischen hoher Ezienz und angemessener Einarbeitungszeit zu nden. 1 Hohe Ezienz ist meistens gleichbedeutend mit der Verwendung neuester Algorithmen, die dann aber nur für ganz spezielle Probleme geeignet sind. 109 6. Numerik 6.1.2. Propagation in Imaginärzeit Die Matrix H KS [ n ] ist reellwertig, dünn-besetzt und symmetrisch. Die Anzahl von Zeilen und Spalten ist gleich der Anzahl von Gitterpunkte des diskretisierten Simulationsraums. In der Praxis liegt diese Anzahl meistens in der Gröÿenordnung 10 2 ... 10 6 . Für 10 2 ... 10 3 können zur Diagonalisierung noch Verfahren benutzt werden, welche auch für dicht-besetzte Matrizen entwickelt wurden. Für gröÿere Matrizen muss ausgenutzt werden, dass diese dünn-besetzt sind und nur eine kleine Anzahl der niedrigsten Eigenwerte und Eigenvektoren benötigt werden. In dieser Arbeit wird das Eigenwertproblem mittels Propagation in Imaginärzeit (ITP) gelöst[67]: Dazu wird die zeitabhängige Schrödingergleichung in Imaginärzeit = it formuliert: = H ^ (6.3) Der Operator H ^ ist der zeitunabhängige Hamiltonoperator, zu dem das Eigenwertproblem H ^ j = j j gelöst werden soll(d.h. hier: H ^ = H ^ KS [ n ] ). Der Propagationsoperator in Imaginärzeit lautet entsprechend: U ^ ( )= e H ^ / (6.4) Die Funktionsweise dieser Methode lässt sich über die Spektralzerlegung des Operators H ^ erläutern: H ^ = j | j j | ( j < j+1 ) j=0 (6.5) (Zur Vereinfachung wird hier ein rein diskretes Spektrum ohne Entartung betrachtet.) Wendet man nun den Propagator 6.4 auf einen Zustand | mit diversen spektralen Anteilen an, erhält man für den Grundzustand | 0 : | = c 0 | 0 + c 1 | 1 +... (6.6) U ^ ( ) | = c 0 e H ^ / | 0 + c 1 e H ^ / | 1 +... = c 0 e 0 / | 0 + c 1 e 1 / | 1 +... (6.7) Da 0 der kleinste Eigenwert ist, verschwindet dieser spektrale Anteil von | am langsamsten. Zerlegt man das Propagationsintervall in N kleine Intervalle = N und normiert nach jedem Zeitschritt den resultierenden Vektor, erhält man ein Extraktionsverfahren, welches den Spektralanteil mit kleinstem Eigenwert aus | in normierter Form für N hervorbringt. 2 engl.: imaginary time propagation 110 6.1. Lösungsverfahren für stationäre Kohn-Sham Gleichungen Damit bleiben noch zur vollständigen Lösung des Eigenwertproblems 6.2 die Fragen oen, wie der Propagationsoperator 6.4 numerisch realisiert werden kann und wie man nicht nur den Grundzustand von H ^ sondern die ersten N Eigenzustände berechnen kann: In[67] wird die Splitoperator-Methode SO4(s. Gl. 6.43) vorgeschlagen, um den Operator U ^ ( ) zu approximieren. Die Taylorreihenentwicklung 4. Ordnung(s. Gl. 6.38) hat sich allerdings in dieser Arbeit für alle Eigenwertprobleme als völlig ausreichend erwiesen. Wenn es auf aller höchste Performance ankommt, sollten die Analysen aus[67] berücksichtigt werden ! um die entscheidenden Optimierungen am verwendeten Algorithmus einzubauen. Die Berechnung der ersten N Eigenzustände erfordert ein spezielles Orthogonalisierungsverfahren, welches im nächsten Kapitel beschrieben wird. 6.1.3. Orthogonalisierungsverfahren Um die ersten N Eigenwerte und Eigenvektoren einer Matrix H mit der ITP-Methode zu berechnen, macht man folgenden Ansatz: Man wählt N Startvektoren | (00) ,... | ( N 0) , welche jeweils über alle spektralen Anteile der Matrix H verfügen. Das erzielt man dadurch, dass man die zugehörigen Vektordarstellungen dieser Vektoren mit Zufallszahlen vorbesetzt. Im Prinzip geht man dann einfach so vor, dass diese Vektoren jeweils alle um einen kleinen Zeitschritt propagiert und anschlieÿend untereinander orthogonalisiert und normiert werden: | ~ ( j 0) = U ^ ( ) | ( j 0) Orth o.+ Norm. | ( j 1) (6.8) Man erzeugt nach diesem Schema aus jedem Satz von Vektoren | ( jk ) einen neuen Satz normalen | ( jk +1) von ortho Vektoren. Naheliegender Weise würde man hier das Orthogonalisierungsverfahren von Gram-Schmidt " benutzen. Es ergeben sich nach genügend Iterationen die gesuchten Eigenvektoren(aus denen auch die Eigenwerte als Erwartungswert des Operator H ^ berechnet werden können): lim k | ( jk ) = | j (6.9) Die Erfahrung mit der Gram-Schmidt Methode für die Systeme in Kapitel 3 hat aber gezeigt, dass es zu erheblichen numerischen Genauigkeitsproblemen bei der Orthogonalisierung kommt, und diese umso mehr zunehmen, je mehr Vektoren man berechnen möchte. Dieses Problem wird auch in[67] beschrieben und ein spezielles Orthogonalisierungsverfahren von Aichinger und Krotscheck[70] benutzt, das auf der Diagonalisierung der Überlappmatrix M basiert. Die Matrixelemente M ij der Überlappmatrix M sind gegeben durch: M ij = ( ik ) U ^ ( ) U ^ ( ) ( jk ) = ~ ( ik +1) ~ ( jk +1) (6.10) 3 Die dort vorgestellte ITP-Implementation mit adaptiven Zeitschritten erweist sich sogar als noch 4 ezienter als der legendäre ARPACK -Algorithmus: Implicitly Restarted Lanczos Method[68]. s. z.B. Kapitel 1.7 in[69]. 111 6. Numerik (Der Unterschied zwischen und ~ ist an Gl. 6.8 ersichtlich.) Für diese Matrix wird ein Eigenwertproblem aufgestellt: M ij c ( jn ) = m n c ( in ) j (6.11) Die Koezienten c ( jn ) bilden den n -ten Eigenvektor c ( n) zum Eigenwert m n . Die Diwesentlich mension der Überlappmatrix ist kleiner als die Anzahl der Gitterpunkte des diskretisierten Simulationsraums. Diese Matrix kann daher mit generischen Algorithmen (beispielsweise aus der lapack -Bibliothek) für dicht-besetzte Matrizen diagonalisiert werden. Mit dem Ergebnis { m n , c ( n) } dieser Diagonalisierung wird nun ein neuer Satz von Vektoren | ( jk +1) (als Ergebnis des Schemas in Gl. 6.8) generiert: | ( jk +1) = 1 m j c i ( j) | ~ ( ik +1) i (6.12) Mit einer Rechnung lässt sich schnell zeigen, dass diese Vektoren orthonormal sind: ( ik +1) | ( jk +1) = 1 m i m j c ( ni ) c ( lj ) ~ ( nk +1) | ~ ( lk +1) n,l = 1 m i m j c ( ni ) c ( lj ) M nl n,l = 1 m i m j c ( ni ) c ( lj ) M nl nl = mm ij m j c ( ni ) c ( nj ) n = m j m i ij = ij = m j c ( nj ) (6.13) Im nächsten Kapitel wird der fertige Algorithmus zur Diagonalisierung der Matrix H KS unter Verwendung dieses speziellen Orthogonalisierungsverfahrens zusammengefasst. 6.1.4. ITP-basierter Diagonalisierungs-Algorithmus Zur Lösung des Eigenwertproblems in Gl. 6.2 wurde für alle Ergebnisse aus Kapitel 3 dieser Algorithmus verwendet: 1. Startwellenfunktionen | ( j 0) bestimmen(Vektoren mit Zufallswerten initialisieren) 2. In Imaginärzeit propagieren: | ~ ( jk ) = U ^ ( ) | ( jk ) 3. Orthogonalisieren: a) Überlappmatrix(Gl. 6.10) berechnen b) Überlappmatrix diagonalisieren(Gl. 6.11) c) Eigenvektoren aufsteigend nach den zugehörigen Eigenwerten sortieren 112 6.1. Lösungsverfahren für stationäre Kohn-Sham Gleichungen Abbildung 6.1.: Schleife zur Berechnung einer selbstkonsistenten Elektronendichte für die Kohn-Sham Gleichungen 6.1a-6.1b. Die Konvergenzgeschwindigkeit wird von der Wahl der Startdichte und dem Dichtekorrekturverfahren entscheidend bestimmt. d) Neuen Satz von Vektoren nach Gl. 6.12 berechnen 4. Weiter bei Schritt 2 bis Konvergenz der Vektoren eintritt: ( jk +1) j ( k) < (6.14) Eine sehr wichtige praktische Erfahrung mit diesem Algorithmus besteht darin, dass dessen Konvergenz erheblich verbessert wird, wenn man ein paar mehr Vektoren # berechnen lässt, als für das Problem erforderlich ist. Die Eigenvektoren mit niedrigsten Eigenwerten konvergieren am schnellsten(wie man Gl. 6.7 erkennen kann). Dieses Konvergenzverhalten bringt noch einen praktischen Vorteil mit sich: Die Eigenvektoren mit den niedrigsten Eigenwerten können, sobald Konvergenz bei diesen eingetreten ist, eingefroren werden: Diese brauchen vom Algorithmus während der Iteration nicht mehr weiter aktualisiert werden. 6.1.5. SCF Iterationsverfahren Die stationären Kohn-Sham Gleichungen 6.1a-6.1b werden iterativ gelöst(s. Abb. 6.1): Um die gesuchte selbstkonsistente Dichte n ( r) zu nden, wird eine Startdichte n (0) ( r) vorgegeben, welche eine physikalisch naheliegende Näherung der tatsächlichen Dichte darstellt $ . Aus dieser Dichte wird das eektive Potential V eff [ n ] berechnet. Dann wird die Kohn-Sham Eigenwertgleichung 6.1a mit diesem Potential gelöst. Die resultierenden Orbitale werden benutzt, um die zugehörige Elektronendichte dieser Orbitale nach Gl. 6.1b zu berechnen. Der zentrale Schritt, welcher die Konvergenzgeschwindigkeit dieses Verfahrens(neben der Wahl der Startdichte) entscheidend bestimmt ist folgender: Aus der Dichte n ( k) ( r) muss mit einem geeigneten Korrekturverfahren eine neue Dichte n ( k+1) ( r) 5 Für die 2D/3D-Systeme aus Kapitel 3 wurden 5-20 Vektoren zusätzlich verwendet. Bei den 1D6 Systemen können statt ITP auch Verfahren aus der lapack -Bibliothek verwendet werden. In Kapitel 3 wurde für die gezeigten Systeme jeweils die Ionendichte n + ( r) als Startdichte n (0) ( r) verwendet. 113 6. Numerik ermittelt werden, welche der tatsächlichen Dichte n ( ) ( r) immer näher kommt. Alle diese Verfahren müssen die Eigenschaft haben, nach unendlich vielen Iterationen die korrekte Dichte zu liefern: lim k n ( k) ( r)= n ( ) ( r) (6.15) Weil die Diagonalisierung der Kohn-Sham Matrix H KS [ n ] extrem viel Rechenzeit beansprucht und die Berechnung des eektiven Potentials eine ebenfalls Rechenzeit intensive Lösung der Poissongleichung(s. Kapitel 6.6) beinhaltet, wurden viele Verfahren entwickelt[7073] % , die eine möglichst geschickte Dichtekorrektur erlauben um die Iterationsanzahl so gering wie möglich zu halten. Diese Verfahren unterscheiden sich nicht nur in ihrer Konvergenzgeschwindigkeit, sondern insbesondere auch in ihrer Stabilität. Das stabilste und gleichzeitig auch langsamste simple mixing Verfahren ist das-Verfahren[73]. Die neue Dichte n ( k+1) ( r) wird über einen Mischparameter (0, 1) der vorherigen Dichte beigemischt: [] n ( k+1) ( r)=[1 ] n ( k) ( r)+ F n ( k) ( r) (6.16) Die Abbildung F ordnet der Eingabedichte n diejenige Dichte F [ n ] zu, welche sich aus den Orbitalen ergibt, wenn die Kohn-Sham Gleichungen mit dem eektiven Potential V eff [ n ] gelöst werden. Die gesuchte Dichte erfüllt also die Gleichung n = F [ n ] . Das Residuum R( r) ist somit gegeben durch: R( r)= n ( r) - F [ n ]( r) (6.17) Damit kann Gleichung 6.16 umgeschrieben werden: [] n ( k+1) ( r)= n ( k) ( r) n ( k) ( r)+ F n ( k) ( r) = n ( k) ( r) R( r) (6.18) Geometrisch bedeutet das, dass man sich bei diesem Iterationsverfahren immer um ein kleines Stück des Residuenvektors weiterbewegt, bis dieser klein genug geworden ist und die Iteration abgebrochen werden kann: R( r) < (6.19) Ein alternatives Abbruchkriterium, welches nicht das gesamte Dichtefeld benutzt, besteht darin, die Konvergenz der Kohn-Sham Eigenwerte j zu analysieren(dieses Kriterium kann z.B. in dem octopus -Softwarepaket[74] benutzt werden). Die Ergebnisse in Kapitel 3 wurden alle mit einem Abbruchkriterium vom Typ 6.19 erzielt, d.h. speziell: max { | n ( k+1) ( r) n ( k) ( r) | } < r n 0 (6.20) 7 Einige dieser Referenzen stellen Verfahren vor, in denen es allgemein um nichtlineare Gleichungen geht. 114 6.1. Lösungsverfahren für stationäre Kohn-Sham Gleichungen Dabei ist der Simulationsraum, n 0 die maximale Ionendichte und die relative Dichteänderung bei deren Unterschreitung abgebrochen wird. Für wurden Werte von 10 12 bis 10 6 verwendet. Hier muss häug ein Kompromiss zwischen Genauigkeit und Rechenzeit simple mixing gefunden werden, da das-Verfahren eine konstante Anzahl an Schritten benötigt um das Residuum um jeweils eine Gröÿenordnung zu verkleinern. 6.1.6. Adaptives Mixing Verfahren Es ist nicht Ziel dieser Arbeit gewesen, eine Verbesserung der Selbstkonsistenzschleife zu entwickeln, da auf diesem Gebiet schon genügend Methoden existieren. Jedoch wurde simple mixing folgende naheliegende Verbesserung des Verfahrens implementiert, die sich an dem Konzept der adaptiven Zeitschritte zur Integration von Bewegungsgleichungen anlehnt. Folgende Modikationen sind dafür notwendig: · Vor jedem Schritt wird die Elektronendichte n ( k) ( r) zwischengespeichert. · Jeder Schritt besteht aus 2 N substep Unterschritten: Es werden N substep simple mixing-Schritte mit dem momentanen Mischparameter durchgeführt. Zusätzlich werden N substep unabhängige simple mixing- Schritte mit 2 durchgeführt, die ebenfalls auf der Eingabedichte n ( k) ( r) basieren. · Die beiden resultierenden Dichten werden über das Kriterium 6.20 verglichen und die Dichte mit besserem Konvergenzwert wird als Ausgabedichte n ( k+1) ( r) benutzt. · Wenn die 2 -Dichte bessere Konvergenz zeigte, wird der Mischparameter verdoppelt und ansonsten halbiert. · Wird eine signikante Vergröÿerung des Residuums festgestellt, muss die zwischengespeicherte Elektronendichte wiederhergestellt werden. Dieser Fall impliziert, dass der Mischparameter bereits einmal halbiert wurde. Der Mischparameter muss in diesem Fall nochmals halbiert werden, weil sonst zweimal die selbe erfolglose Rechnung mit identischem Mischparameter durchgeführt wird. Eine Implementation dieses Verfahrens sollte die Skalierung des Mischparameters frei einstellbar lassen, da der Faktor 2 für manche Systeme zu groÿ ist. Ausserdem hilft es, während der Iteration das Konvergenzverhalten über die letzten 5-10 Schritte zu betrachten, um oszillierendes oder stagnierendes Verhalten zu erkennen. Wenn dieses Verhalten beobachtet wird, muss die Schleife abgebrochen werden, da das vorgegebene simple mixing Konvergenzziel nicht mehr erreicht wird & . Das-Verfahren ist dann nicht in der Lage, die Kohn-Sham Gleichungen zu lösen ' . Die Ergebnisse für eine Verbesserung des Konvergenzverhaltens sind in Abbildung 6.2 illustriert. Die Erfahrung mit dem Verfahren hat gezeigt, dass es sich gut eignet, um einen 8 Es handelt sich dabei um praktische Erfahrungswerte mit dem Konvergenzverlauf. 9 Einige 3D-Systeme, wie die Kugel aus[42], zeigten dieses Verhalten. Hier hat es geholfen, die Geometrie oder die Diskretisierung leicht zu verändern. 115 6. Numerik Abbildung 6.2.: Das Diagramm zeigt das Residuum(6.20) nach jedem Iterationsschritt bei konstantem und adaptivem Mischparameter . Das adaptive Verfahren hat ca. 150 Schritte weniger benötigt. Die CPU-Zeit hat bei diesen Rechnungen für = 0. 03 14. 2 s und bei adaptiven 11. 5 s betragen. Der kleine Peak in der roten Kurve zeigt, dass bei dem adaptiven Verfahren der Mischparameter zwischendurch zu groÿ gewählt wurde und reduziert werden musste, um divergierendes Verhalten zu verhindern. optimalen Wert für den Mischparameter zu erhalten, den man dann in das = const. Verfahren einsetzen kann. Um die benötigte CPU-Zeit gegenüber dem nicht-adaptiven Verfahren zu verringern, muss leider viel mit den Parametern experimentiert werden. 6.1.7. Bedeutung der Abschirmkonstante für die SCF-Schleife Die Berechnung des Hartree-Anteils V H [ n ] des eektiven Potentials beinhaltet die Lösung der Poissongleichung: [ 2 + 2 ] ( r)= 1 ( r) 0 (6.21) Die elektrostatische Abschirmkonstante hat erheblichen Einuss auf die Stabilität der SCF-Schleife : Speziell bei eindimensionalen Systemen, wie dem Metalllm aus Kapitel künstli3, konnte die SCF-Schleife nicht zur Konvergenz gebracht werden, wenn keine che elektrostatische Abschirmung vorhanden war. Es ist leicht nachvollziehbar, weshalb der eindimensionale Fall besonders anfällig für die Instabilität ist: Das Coulombpotential einer Ladung am Ort x ist proportional zu | x x | (entsprechend dem Feld einer unendlich ausgedehnten Flächenladung). Jede kleinste Änderung der Ladungsdichte n ( x) , welche die SCF-Schleife vornimmt, hat im Simulationsraum überall gleich groÿe Auswirkungen auf das elektrostatische Feld. In Anlehnung an das Konzept der Konditionszahl 10 s. Anhang von[75]. 116 6.2. Propagatoren für die Kohn-Sham Gleichungen für lineare Gleichungssysteme(s. Kapitel 6.3.6) kann man hier bei diesen nichtlinearen Gleichungen von einem schlecht konditionierten Problem sprechen, da kleine Änderungen an der Eingabedichte groÿe Wirkung auf die Ausgabedichte haben. In dieser Arbeit konnte dieses Problem durch folgende Vorgehensweise gelöst werden [76]: Für die Jelliumsysteme aus Kapitel 3 wurde zunächst in der Gröÿenordnung des Fermivektors k F gewählt. Dann wurden die Kohn-Sham Gleichungen gelöst. Die resultieerneut rende Ladungsdichte n wird dann als Startdichte n (0) ( r) benutzt, um die Kohn-Sham Gleichungen mit einer verkleinerten Abschirmkonstante zu lösen. Diese Schritte wurden so oft wiederholt, bis die Elektronendichte unabhängig von der Abschirmkonstante war. Eine alternative Methode wird im Anhang von[75] und in Kapitel 2 von[32] beschrieben. 6.2. Propagatoren für die Kohn-Sham Gleichungen Eine umfassende Beschreibung von Propagatoren für die zeitabhängigen Kohn-Sham Gleichungen ist in[77] zu nden. In diesem Kapitel wird ein Überblick über das Thema gegeben und es werden die Methoden, welche speziell in dieser Arbeit benutzt werden, beschrieben. 6.2.1. Problembeschreibung Die zeitabhängigen Kohn-Sham Gleichungen für ein N -Elektronen System lauten: { i t j ( r, t)= ( 2 2 m e 2 + V eff [ n ]( r, ) t) j ( r, } N t) j=1 (6.22) N n ( r, t)= w j | j ( r, t) | 2 j=1 (6.23) Diese haben folgende strukturelle Eigenschaften, die für die Problemanalyse von zentraler Bedeutung sind: 1. Es handelt sich mathematisch gesehen um N Ein-Teilchen Schrödingergleichungen. 2. Der Hamiltonoperator H ^ KS ( t) ist zeitabhängig, weil auch die Teilchendichte n ( r, t) und eventuell vorhandene externe Störpotentiale zeitabhängig sind. 3. Die Gleichungen sind über das eektive Potential miteinander gekoppelt. Die Kopplung der Gleichungen ist dabei als nebensächliches Problem anzusehen, da die gröÿte Schwierigkeit darin besteht, die zeitabhängige Ein-Teilchen Schrödingergleichung ezient zu propagieren(s. Kapitel 6.2.2). 117 6. Numerik 6.2.2. Numerische Lösung der zeitabhängigen Ein-Teilchen Schrödingergleichung Die Ein-Teilchen Schrödingergleichung lautet: i ( r, t)= H ^ ( t) ( r, t) t (6.24) Das Lösen dieser Gleichung erfordert die Vorgabe von Randbedingungen im Raum und einer Anfangsbedingung in der Zeit: ( , t)= 0 ( r, 0)= 0 ( r) (6.25) (6.26) Mit wird der Simulationsraum bezeichnet. Die Wellenfunktion wird auf einem regulären, kartesischen Gitter in diesem Gebiet dargestellt: j,k,l ( t)= ( j x e x + k y e y + l z e z , t) (6.27) Da bei der numerischen Verarbeitung nur endlich viele Gitterpunkte gespeichert werden können, ist auch das Gebiet nur endlich groÿ. Die Randbedingung 6.25, welche einem unendlich tiefen Potentialtopf entspricht, wurde gewählt, weil sich diese besonders einfach implementieren lässt und für alle Anwendungen in dieser Arbeit geeignet ist. Der Operator H ^ ( t) beinhaltet den Operator der kinetischen Energie, der im Ortsraum ein skalierter Laplaceoperator ist. Dieser Operator wird auf dem Gitter mittels Finiter Dierenzen dargestellt[66, 78]. Üblicher Weise verwendet man einen 3- oder 5-Punkt Stempel. Der Hamiltonoperator kann letztlich auf dem Gitter in Form einer endlichen Matrix H ( t) dargestellt werden. Der räumliche Anteil der Schrödingergleichung sei damit nun ein gelöstes Teilproblem. Das verbleibende Teilproblem besteht in der Propagation der Wellenfunktion in der Zeit: In einer Simulation muss aus der Wellenfunktion zur Zeit t diese zu einem Zeitpunkt t+ t berechnet werden. Dieses Problem kann formal über den Zeitentwicklungsoperator U ^ beschrieben werden: ( t+ t)= U ^ ( t+ t, t) ( t) (6.28) U ^ ( t 1 , t 0 )= ( T ^ D exp i t 1 ) H ^ ( t) dt t 0 (6.29) Bei dem Operator T ^ D handelt es sich um den dysonschen Zeitordnungsoperator[79] und der Operator T ^ D exp wird als zeitgeordneter Exponent bezeichnet[77]. Für die Numerik muss der Operator 6.29 auf eine endliche Anzahl an elementaren Rechenoperationen, die auf die diskretisierte Wellenfunktion j,k,l ( t) anzuwenden sind, runter gebrochen werden. Das setzt eine Approximation dieses Operators voraus, deren Fehler die Eigenschaft haben soll, für t 0 möglichst schnell gegen Null zu gehen. An dieser Stelle muss ein Kompromiss gemacht werden: Je besser die Approximation ist, 118 6.2. Propagatoren für die Kohn-Sham Gleichungen umso gröÿere Zeitschritte können zwar gemacht werden, jedoch steigt der Rechenaufwand pro Zeitschritt ebenfalls an . Die systematische Approximation des Operators 6.29 wird als Magnus-Entwicklung [80, 81] bezeichnet. Das Entscheidende an dieser Entwicklung ist, dass sich der Zeitenteinfacher wicklungsoperator als Exponent eines Operators ^ schreiben lässt, der in einer Reihe entwickelt werden kann: {} U ^ ( t 1 , t 0 )= exp ^( t 1 , t 0 ) (6.30) ^( t 1 , t 0 )=^ k ( t 1 , t 0 ) k=1 (6.31) Die Operatoren dieser Reihenentwicklung können nach folgendem(rekursivem) Schema unter Verwendung eines weiteren Operators S ^ und der Bernoulli-Zahlen B j berechnet werden: ^ k ( t 1 , t 0 )= k 1 j=0 B j j! t 1 t 0 S ^ kj ( ) d S ^ 10 ( )= 1 H ^ ( ) i S ^ k 0 ( )= 0( k> 1) S ^ kj ( )= k j [ ^ m ( t 1 , t 0 ), ] S ^ kj -1 m ( )(1 j k 1) m=1 (6.32) (6.33) (6.34) (6.35) Die Integration in Gl. 6.32 muss mit einer numerischen Quadraturformel durchgeführt werden. Um eine Approximation des Magnusoperators ^ zur Ordnung 2 n zu erhalten, muss die Reihe in Gl. 6.31 nach dem n -ten Term abgebrochen werden und die auftretenden Integrale mit einer Quadraturformel n -ter Ordnung ausgewertet werden[77]. In dieser Arbeit wurden die Operatoren in der zweiten und vierten Ordnung benutzt. Das Ergebnis für diese beiden(praktisch wichtigen) Fälle ist ebenfalls in[77] dokumentiert. Für einen Zeitschritt t mit t 0 = t und t 1 = t+ t lauten die Operatoren: ^ M(2) ( t+ t, t) ^ M(4) ( t+ t, t) = = 1 H ^ ( t+ t/ 2) t i i t 2 [ H ^ ( 1 )+ ] H ^ ( 2 ) 3( t) 2 12 2 [ H ^ ( 2 ), ] H ^ ( 1 ) (6.36) (6.37) 1, 2 = t+[(1/ 2) 3/ 6] t Von diesen Operatoren muss jetzt nur noch der einfache Exponent aus Gl. 6.30 gebildet werden und auf eine Wellenfunktion angewendet werden. Bei der Methode der Finiten Dierenzen liegt ein solcher Operator in Form einer Matrix vor. Um den Exponenten einer Matrix zu berechnen, existieren zahlreiche Algorithmen 11 In[77] sind Erfahrungswerte für dieses Problem ausführlich dokumentiert. 119 6. Numerik [82, 83]. Allerdings ist diese Vorgehensweise, nämlich diesen Exponenten der Matrix explizit zu berechnen, technisch kaum oder gar nicht möglich: Die resultierende N × N Matrix ist dicht besetzt und N ist dabei die Anzahl an verwendeten Gitterpunkten, welche durchaus in der Gröÿenordnung 10 3 ... 10 6 liegen kann. prozedural Stattdessen werden Operatoren der Form exp( A ^) grundsätzlich nur auf Wellenfunktionen angewendet, d.h. der Exponentialoperator wird durch ein Polynom approximiert, dessen Auswertung im Wesentlichen nur die Rechenoperation A ^ wiederholt benötigt. Dazu gibt es folgende Möglichkeiten[77]: · Der Exponentialoperator wird in einer Taylorreihe entwickelt: exp( A ^) = 1 A ^ ! =0 (6.38) Die prozedurale Auswertung erfolgt nach dem Horner-Schema[84]. · An Stelle der Monombasis kann die Reihe über Chebyshev-Polynome dargestellt werden: exp( A ^) = c T ( A ^) =0 (6.39) Die prozedurale Auswertung erfolgt nach dem Clenshaw-Algorithmus[85]. Die Formel zur Berechnung der Entwicklungskoezienten c ist in[77, 86] zu nden. Bei dem Operatorpolynom T ( A ^) muss beachtet werden, dass der Operator in seiner Spektralzerlegung keine Eigenwerte ausserhalb des Intervalls [ 1, 1] haben darf, weil die Chebyshev-Polynome nur auf diesem Intervall deniert sind. · In[87] wird eine Krylov-Unterraum Methode beschrieben, mit der sich exp( A ^) ebenfalls berechnen lässt. Nach der Analyse der Methoden in[77] ist diese Methode die ezienteste von den hier Genannten. Unter anderem lässt sich mit der Methode auch der Fehler der Approximation leicht abschätzen. Die Methoden unterscheiden sich letztlich darin, wie viele Matrix-Vektor Multiplikationen (speziell A ^ ) zur Propagation pro Zeitintervall durchzuführen sind, um eine vorgegebene Mindestgenauigkeit einzuhalten. Für die gezeigten Beispiele in[77] benötigte die KrylovMethode bis zu Faktor 2 weniger Operationen. Für alle Anwendungen in dieser Arbeit ist die Taylorreihenentwicklung 4. Ordnung vollkommen ausreichend gewesen, weil der theoretisch mögliche Performancegewinn durch die beiden anderen Methoden deren Implementationsaufwand nicht gerechtfertigt hat. Die Verwendung der 4. Ordnung wird ausserdem in[77] empfohlen. 6.2.3. Alternative Propagationsmethoden Im letzten Kapitel wurde der Zeitentwicklungsoperator mit einer Magnusentwicklung approximiert. Diese liefert nicht nur eine Möglichkeit, den zeitgeordneten Exponenten 120 6.2. Propagatoren für die Kohn-Sham Gleichungen zu vereinfachen, sondern bietet auch den groÿen Vorteil, dass sie in jeder Ordnung den unitären Charakter des Zeitentwicklungsoperators beibehält. Von den populären Standardverfahren wie der Crank-Nicholson(CN) Methode[66] und den Allzweck-Integratoren wie dem expliziten Runge-Kutta(RK) Verfahren wird in[77] (ohne nähere Begründung) abgeraten. Mögliche Gründe sind aber recht oensichtlich: Die CN-Methode verwendet zwar eine unitäre Approximation des Zeitentwicklungsoperators, allerdings muss bei der Methode in jedem Zeitschritt ein lineares Gleichungssystem gelöst werden. Die Form des Gleichungssystems wird im Wesentlichen vom Hamiltonoperator des Systems bestimmt. Beschränkt man sich auf eindimensionale Systeme, dann stellt das Gleichungssystem kein numerisches Problem dar, weil die Finite Dierenzen Form des Hamiltonoperators eine tridiagonale Matrix ergibt und daher in O ( N) (mit N als Anzahl der Gitterpunkte) gelöst werden kann[66]. Für mehrdimensionale Systeme ist diese Methode ungeeignet, da das Gleichungssystem diese Eigenschaft verliert. Die populären RK-Verfahren 2. und 4. Ordnung wurden in dieser Arbeit an verschiedenen Stellen getestet. Dabei zeigte sich immer wieder, dass in 2. Ordnung die Normierung der Wellenfunktion nicht erhalten blieb und die Simulation schnell instabil wurde. Das Verfahren 4. Ordnung zeigte keine signikanten Abweichungen in der Norm der Wellenfunktion(bedingt durch die hohe Ordnung des Verfahrens) und blieb stabil. zeitumkehrbar Die RK-Verfahren haben den Nachteil, dass die Lösung nicht ist. Propagiert man den Zustand ( t) um t und anschlieÿend um t ergibt sich ein anderer Zustand ~ ( t) als der Ausgangszustand: ( t) t ( t+ t) t ~ ( t) = ( t) Die korrekte Lösung der zugrunde liegenden Gleichungen hat aber diese Eigenschaft der Zeitumkehrbarkeit. Deshalb sollten nur Methoden benutzt werden, welche diese berücksichtigen. Das RK-Verfahren 4. Ordnung kann allerdings aufgrund seiner hohen Genauigkeit diese Unzulänglichkeit kompensieren. Eine wichtige Klasse von Propagatoren bilden die Split-Operator Methoden[88]: Man nehme zunächst an, dass der Zeitschritt t klein genug ist, damit die Zeitabhängigkeit des Hamiltonoperators vernachlässigbar ist. Der Zeitentwicklungsoperator vereinfacht sich dadurch erheblich[79]: U ^ ( t)= () exp i H ^ t = ( exp i [ T ^ ] + V ^ ) t (6.40) Die Idee besteht darin, den Exponenten zu faktorisieren. Nach der Baker-CampbellHausdor Relation müssen die Operatoren im Exponenten dafür miteinander kommutieren, was für T ^ und V ^ i. A. nicht der Fall ist. In der Praxis hat sich aber gezeigt, dass der Kommutator vernachlässigbar ist, wenn t klein genug gewählt wird. Es wird auÿerdem immer eine symmetrische Faktorisierung gewählt, bei der ein Operator zweimal auftritt: ( )()() U ^ SO2 ( t)= exp i t T ^ exp i t V ^ exp i t T ^ ( 2 )()( 2 ) (6.41) U ^ SO2 ( t)= exp i t V ^ 2 exp i t T ^ exp i t V ^ 2 (6.42) 121 6. Numerik Die Approximationen 6.41-6.42 haben die Ordnung O ( t 2 ) . In 4. Ordnung muss der Kommutator berücksichtigt werden : ( )() U ^ SO4 ( t)= exp i t V ^ exp i t T ^ ( 6 { exp 2 i V ^ 1 [ V ^ , 2 [ V ^ , ]] T ^ } t 2 t ) ( 3 ) 48 () exp i t T ^ exp i t V ^ 26 (6.43) Der Split-Operator Ansatz zielt vor allem darauf ab, den Operator der kinetischen Enerexakt gie T ^ =( i ) 2 / 2 m e (d.h. ohne Finite Dierenzen) anwenden zu können: Dafür muss die Wellenfunktion vor jeder Anwendung des Exponentialoperators mit T ^ -Term fouriertransformiert werden. In der Fourierbasis besteht die Wirkung des Operators T ^ bzw. exp(const. · T ^) nur noch in einer gitterpunktweisen Multiplikation. Anschlieÿend muss eine inverse Transformation durchgeführt werden um den Operator V ^ im Ortsraum auf exakt die gleiche eziente Weise(d.h. punktweise Multiplikation mit exp(const. · V ^ ) ) anwenden zu können. Diese Methode impliziert aufgrund der Fouriertransformation das Vorhandensein von periodischen Randbedingungen. Meistens liegen aber Wellenfunktionen vor, welche mit einer Dirichletrandbedingung ( )= 0 bestimmt wurden. Praktisch stören diese Unterschiede bei den Randbedingungen nicht, wenn die Wellenfunktionen weit genug im Inneren des Simulationsraums lokalisiert sind. In der Arbeit von Sugino und Miyamoto[89] wird gezeigt, wie sich die Split-Operator Methode ! auf noch höhere Ordnung erweitern lässt: Dabei wird vor allem die Zeitabhängigkeit des Hamiltonoperators berücksichtigt, die in den Gln. 6.41-6.43 vernachlässigt wird. Dazu ist es erforderlich, vom Hamiltonoperator H ^ ( t) das zeitabhängige Potential V( t) im Intervall [ t, t+ t] zu interpolieren " . Als Interpolationsmethode wird das Railway-Curve Schema # verwendet $ : V( s)= ( s t t ) 2 × [ 3 V( t)+ tV ( t)+ s t t [2 V( t)+ tV ] ( t)] t t (6.44) + ( s t ) 2 [ 3 V( t+ t) tV ( t+ t) s t [2 V( t+ t) tV ( t+ ] t)] t t Die Interpolationsformel erfordert leider auch die Kenntnis der Zeitableitung von der zu interpolierenden Gröÿe % . Der Rechenaufwand kann dadurch gegenüber anderen Interpolationsverfahren zwar deutlich höher sein, aber dafür ist durch Verwendung dieser 12 Der Term in dieser Form ist aus[67](Gl. 4) zum Thema Imaginärzeitpropagation entnommen. Die 13 Ersetzung i t ergibt Gl. 6.43. Speziell wird die Suzuki-Trotter Zerlegung[90] verwendet, welche ein Produkt aus fünf U ^ SO2 14 Operatoren mit unterschiedlicher Schrittweite t verwendet. Die Extrapolation des Hamiltonoperators von der Zeit t zu der Zeit t+ t wird im nächsten Kapitel beschrieben. 15 s.[91] 16 Die angegebene Interpolationsformel in[89](Gl. 14) enthält einen Vorzeichenfehler: Die korrekte Form 17 ist in Gl. 6.44 gegeben. Für Anwendungen in der zeitabhängigen DFT ist es an dieser Stelle erforderlich, vom Hartree- und xc-Potential die Zeitableitung zu bilden. Dazu wird in[89] eine Formel(Gl. 16) hergeleitet. 122 6.3. Propagatoren für implizite Kohn-Sham Gleichungen Zeitumkehrbarkeit Formel die gewährleistet, was entscheidend zur numerischen Stabilität der Propagationsmethode beiträgt. Dadurch können letztlich auch gröÿere Zeitschritte gewählt werden, so dass der zusätzliche Rechenaufwand kompensiert wird. Dieses Schema kommt in dieser Arbeit bei den Propagatoren für die impliziten KohnSham Gleichungen zum Einsatz. 6.2.4. Extrapolation des Hamiltonoperators In Kapitel 6.2.2 wurde bei der Magnusentwicklung des Zeitentwicklungsoperators bereits angedeutet, dass der Hamiltonoperator für den Zeitschritt t t+ t auch zu Zeitpunkten innerhalb dieses Intervalls bekannt sein muss(s. Gln. 6.36-6.37). Das ist unproblematisch, wenn der zeitabhängige Anteil des Hamiltonoperators nur aus einem Störpotential V P ( t) besteht, das zu jeder Zeit t per Formel berechnet werden kann. Bei den Kohn-Sham Gleichungen ist das denitiv nicht so einfach möglich, weil ja das eektive Potential von der Lösung der Kohn-Sham Gleichungen selbst abhängt. Für den allgemeinen Fall kann immer folgendes Selbstkonsistenz-Schema verwendet werden: 1. Näherungsweise Berechnung des Potentials V( t+ t) mittels einer Extrapolationsmethode, welche evtl. auch Werte von vorangegangenen Zeitschritten benutzen kann. 2. Interpolation von V( t t t+ t) mit Hilfe der Gl. 6.44. 3. Zeitschritt mit Hilfe der im letzten Schritt gewonnenen Zwischenwerte durchführen. 4. V( t+ t) erneut berechnen. Die Schritte 2-4 müssen solange wiederholt werden, bis das Potential V( t+ t) konvergiert ist. 6.3. Propagatoren für implizite Kohn-Sham Gleichungen Die impliziten Kohn-Sham Gleichungen treten in der dissipativen zeitabhängigen Dichtefunktionaltheorie(s. Kapitel 4) auf. Die Nomenklatur Implizite Kohn-Sham Gleichung ist derzeit(noch) nicht in der Literatur zu nden und wird nur in dieser Arbeit benutzt. 6.3.1. Problembeschreibung Die zeitabhängigen, impliziten Kohn-Sham Gleichungen für ein N -Elektronen System lauten: { i t j ( r, t)= ( 2 2 m e 2 + V eff [ n ]( r, t)+ H ^ f [ KS , ) KS ]( t) j ( r, } N t) j=1 (6.45) KS ( r 1 ,..., r N , t)= 1 | 1 · · · N | (Slaterdeterminante, s. Gl. 4.15) N! (6.46) 123 6. Numerik N n ( r, t)= w j | j ( r, t) | 2 j=1 (6.47) Diese haben folgende strukturelle Eigenschaften, die für die Problemanalyse von zentraler Bedeutung sind: 1. Es handelt sich mathematisch gesehen um N implizite Ein-Teilchen Schrödingergleichungen, bei denen jeweils die Zeitableitung der Wellenfunktion auf beiden Seiten der Gleichung auftritt. 2. Der Hamiltonoperator H ^ ( t)= H ^ KS ( t)+ H ^ f ( t) ist zeitabhängig, weil die Teilchendichte n ( r, t) , das Dämpfungsfeld im Reibungsterm H ^ f ( t) und eventuell vorhandene externe Störpotentiale zeitabhängig sind. 3. Die Gleichungen sind über das eektive Potential miteinander gekoppelt. Die Kopplung der Gleichungen ist dabei als nebensächliches Problem anzusehen, da die gröÿte Schwierigkeit darin besteht, die impliziten, zeitabhängigen Ein-Teilchen Schrödingergleichungen ezient zu propagieren. Für dieses Problem existiert in der Literatur derzeit nur der von Neuhauser verwendete Lösungsansatz[42] der im nächsten Kapitel beschrieben wird. 6.3.2. Lösungsverfahren von Neuhauser Zur Beschreibung des Lösungsverfahrens genügt es, sich auf eine Ein-Teilchen Schrödingergleichung mit Reibungsterm zu beschränken: i t ( r, t)= ( 2 2 m e 2 + V( r, t)+ H ^ f ) [] ( t) ( r, t) (6.48) Der zugehörige Zeitentwicklungsoperator lautet: U ^ ( t 1 , t 0 )= ( T ^ D exp i t 1 []) H ^ 0 ( t)+ H ^ f ( t) dt t 0 (6.49) Es wird nun angenommen, dass H ^ 0 ( t) für einen Zeitschritt t nur eine geringe Zeitabhängigkeit besitzt und die folgende Faktorisierung des Operators U ^ erlaubt: U ^ ( t)= e i H ^ 0 t/ 2 U ^ f ( t) e i H ^ 0 t/ 2 (6.50) Dieser Operator ist(trotz dieser Notation) von der absoluten Zeit t abhängig und wird folgendermaÿen auf eine Wellenfunktion ( t) angewendet: Die beiden äuÿeren Exponentialoperatoren werden nach dem Split-Operator SO2-Schema(s. Gl. 6.42) zerlegt. Nachdem der erste Exponentialoperator angewendet wurde, wird der Propagator U ^ f für den Reibungsterm angewendet: () U ^ f ( t)= exp i H ^ f [ ] t 1 i H ^ f [ ] t (6.51) 124 6.3. Propagatoren für implizite Kohn-Sham Gleichungen Im Operator H ^ f ist eine Gröÿe der Form d z /dt enthalten(s. Gln. 4.13a-4.13c), welche von abhängt und verantwortlich dafür ist, dass die Schrödingergleichung hier eine implizite Form hat. Dieses Problem wird im Neuhauser-Ansatz einfach dadurch gelöst, indem eine Finite Dierenz mit dem Wert von z des letzten Zeitschritts gebildet wird: d z dt z( t) - z( t t) t (6.52) Nachdem der Operator U ^ f auf angewendet worden ist, muss noch der letzte Exponentialoperator in Gl. 6.50 angewendet werden, wobei das Zeitargument des Hamiltonoperators nun t+ t/ 2 lautet. Diese Methode wird von den Autoren in[42] nur als funktionierend aber weder als besonders ezient noch besonders genau beschrieben. Gerade die Näherungen in den Gln. 6.51 und 6.52 werfen die Frage auf, wie schnell die Ergebnisse für t 0 konverexplizit gieren. Da die Methode ist, kann es gut sein, dass sie selbst bei Zeitschritten mit 0. 1 × t noch schneller ist, als ein implizites Verfahren mit Schrittweite t - und dabei die gleichen, auskonvergierten Ergebnisse liefert. 6.3.3. Diskretisierung der impliziten Kohn-Sham Gleichungen Die zeitabhängigen Kohn-Sham Orbitale können in einem Vektor y( t)=[ 1 ( t),..., N ( t)] zusammengefasst werden. Jedes einzelne Orbital j kann man sich ebenfalls als einen zeitabhängigen Vektor vorstellen, der durch die Darstellung der Wellenfunktion auf einem kartesischen Gitter(s. Gl. 6.27) gegeben ist. Der Vektor y( t) besteht also aus Anzahl N der Orbitale × Anzahl der Gitterpunkte. Der Hamiltonoperator auf der rechten Seite H der Kohn-Sham Gleichungen wird durch eine Matrix dargestellt. In dieser diskretisierten Form lautet das System von Gleichungen 6.45: i y ( t)= { H KS ( t)+ H f [ y, y ]( t) } y (6.53) nicht Zunächst wird der Fall betrachtet, dass der Reibungsterm H ^ f vorhanden sei: Der Kohn-Sham Hamiltonoperator ist eine Summe von Ein-Teilchen Operatoren, die j jeweils auf das-te Teilchen wirken. Die Matrixdarstellung entspricht daher einer direkten Summe identischer Matrizen: N N H ^ KS ( t)= H ^ ( j) ( t) H KS ( t)= H ( t) j=1 j=1 (6.54) Das ergibt eine blockdiagonale Struktur: H KS ( t)= . . . (6.55) Jetzt wird der Reibungsterm H ^ f ( t) in diese Betrachtung mit eingebracht: Dieser beinhaltet einen Operator Z ^ (s. Gln. 4.13a-4.13c), der sich als Summe von Ein-Teilchen 125 6. Numerik Operatoren z^ ( j) schreiben lässt: N N Z ^ = z^ ( j) Z = z j=1 j=1 (6.56) Als nächstes wird die Zeitableitung vom Erwartungswert Z benötigt: d Z = d y, Z y dt dt = d dt N j , z j j=1 N [] = j , z j + j , z j j=1 (6.57) Durch diese Summe geht die Blockdiagonalform 6.55 verloren: Jede Zeile der Matrix H f ist vollständig gefüllt, weil Z von allen Orbitalen gleichermaÿen abhängig ist. An dieser Stelle wird nochmals auf die Abhängigkeit des Operators H ^ f [ KS , KS ]( t) bzw. H f [ y, y ]( t) eingegangen: Der Reibungsterm darf in beliebiger(nicht-)linearer Form von den Orbitalen & abhängen. Die Zeitableitung der Orbitale geht dagegen nur linear in den Reibungsterm ein, wie man an Gl. 6.57 sehen kann. Bei der Gleichung 6.53 handelt es sich also um ein lineares Gleichungssystem für die unbekannte Zeitableitung y . 6.3.4. Implizite Runge-Kutta Verfahren In diesem Abschnitt wird gezeigt, wie sich das System von impliziten Kohn-Sham Gleichungen 6.45-6.47 mit Hilfe von impliziten Runge-Kutta Verfahren lösen lässt, ohne die Näherungen in den Gln. 6.51 und 6.52 dafür benutzen zu müssen. Als Vorarbeit dafür wurde im letzten Abschnitt die Struktur des Gleichungssystems 6.53 analysiert. Dieses System lösen zu können, reicht nicht aus, um an die Lösung y( t) für ein endliches Zeitintervall zu gelangen. Dazu muss ein numerisches Integrationsverfahren benutzt werden. Ziel ist es, den Zeitschritt y( t) y( t+ t) bzw. mit Zeitschrittindex y n y n+1 durchzuführen. Ein s -stuges Runge-Kutta Verfahren verwendet dafür folgenden Ansatz: s y n+1 = y n + b i k i i=1 s k i = t f( t n + c i t, y n + a ij k j ), j=1 i= 1,..., s (6.58) (6.59) 18 In[42] wird beispielsweise angedacht, dass man das Feld a( q) im Drude-Reibungsterm mit der lokalen Teilchendichte in unmittelbarer Umgebung von q in Verbindung bringt. 126 6.3. Propagatoren für implizite Kohn-Sham Gleichungen Die charakteristischen Koezienten a ij , b i und c i werden unter dem Begri ButcherTableau in der gängigen Literatur[78, 92] aufgeführt. Die Wahl der Koezienten entscheidet über die Ordnung des Verfahrens und darüber, ob es sich um ein explizites oder implizites Runge-Kutta Verfahren handelt. Bei den expliziten Verfahren ist nur die untere Dreieckshälfte der Koezientenmatrix a ij von Null verschieden. Bei den impliziten s Verfahren bilden die Gleichungen in 6.59 ein i. A. nicht-lineares Gleichungssystem für die unbekannten Inkremente k i . Um das Verfahren auf die impliziten Kohn-Sham Gleichungen anzuwenden, wird zunächst folgende Identikation der involvierten Vektoren vorgenommen: Der Vektor y n = [ 1 ( t n ),..., N ( t n )] enthält die diskretisierten Orbitale zur Zeit t n , die s Vektoren k i = [ (1 i) ,..., N ( i) ] enthalten die Inkremente, aus denen der Zeitschritt linearkombiniert wird(s. Gl. 6.58). Die vektorielle Funktion f ist dann folgendermaÿen deniert: f= f( t, y, y ) = 1 i { H KS [ y]( t)+ H f [ y, y ]( t) } y (6.60) Hier wurde zusätzlich das Argument für die Zeitableitung von y aufgeführt, von der die rechte Seite der impliziten Kohn-Sham Gleichungen abhängt. Als Argument wird hier jeweils das i -te Inkrement eingesetzt. Das Gleichungssystem lautet dann: s s k i = t f( t n + c i t, y n + a ij k j , k i / t) j=1 i=1 (6.61) Die Anzahl der Unbekannten im Gleichungssystem beträgt: Anzahl der Stufen s × Anzahl der Gitterpunkte × Anzahl der Orbitale Das führt bei praktischen Anwendungen auf eine so groÿe Anzahl von Unbekannten, dass man nur mit iterativen Verfahren arbeiten kann. Dieses Gleichungssystem ist nicht-linear bezüglich der Unbekannten k i : Im Zusammenhang mit Gl. 6.57 wurde zwar gesagt, dass die Zeitableitung der Orbitale nur linear auftritt, allerdings treten hier die Inkremente auch im zweiten Argument von f (in der Summation) auf. Dieses Argument geht nichtlinear in f ein, da aus diesem die Teilchendichte und aus dieser wiederum die Funktionale des Vielteilchensystems berechnet werden. Zum Lösen der nichtlinearen Gleichungssysteme bei impliziten Runge-Kutta Verfahren sind i. A. Newtonverfahren üblich. Das setzt allerdings die Berechnung des Gradienten der Funktion f voraus. Für die Funktion f kann man dafür zunächst versuchen, von der nicht-diskretisierten Form 6.45 auszugehen und einen analytischen Ausdruck der Funktionalableitung der rechten Seite zu bestimmen. An diesem Ansatz schreckt vor allem der hohe Aufwand bei der Berechnung des Gradienten(sowohl analytisch als auch numerisch) ab. Jede Änderung an der rechten Seite der impliziten Kohn-Sham Gleichungen erfordert eine Anpassung in der Herleitung des Gradienten. Das stört besonders, wenn man verschiedene Reibungsterme H ^ f ( t) testen möchte. Damit bleiben nur Lösungsverfahren übrig, welche ohne Gradienten auskommen: Glücklicherweise hat sich bei allen Anwendungen in dieser Arbeit gezeigt, dass das selbe 127 6. Numerik einfache Iterationsschema, wie es zur Lösung der stationären Kohn-Sham Gleichungen zum Einsatz kommt, auch hier verwendet werden kann. Dabei nutzt man aus, dass als Startvektor der Iteration immer das Ergebnis des letzten Zeitschritts gewählt werden kann, wodurch die Anzahl an Iterationen drastisch reduziert werden kann. 6.3.5. Magnusentwicklung Eine ebenfalls erfolgreich getestete Möglichkeit zur Lösung der impliziten Kohn-Sham Gleichungen 6.45-6.47 besteht darin, ein explizites Lösungsverfahren mit einer Selbstkonsistenzschleife zu kombinieren: Die expliziten Verfahren, welche die Magnusentwicklung verwenden, sind in Kapitel 6.2.2 bereits beschrieben worden. Wenn man diese Verfahren hier anwenden möchte, wird man gleich zu Beginn mit folgendem Problem konfrontiert: Der Zeitentwicklungsoperator ist aufgrund des Reibungsterms von der Kohn-Sham Wellenfunktion selbst abhängig, so dass man für einen Zeitschritt die Gleichung in der ungewöhnlichen Form KS ( t+ t)= U ^ [ KS , KS ]( t+ t, t) KS (6.62) schreiben müsste. Hier ergeben sich bereits Schwierigkeiten, wie man diesen Operator mathematisch sauber denieren kann. Um dieses Problem zu umgehen wird einfach angenommen, dass der Hamiltonoperator bekannte nicht eine Zeitabhängigkeit auf dem Intervall [ t, t+ t] besitzt, die von der Zeitableitung des Zustands abhängt. Dann lässt sich der Zeitentwicklungsoperator in bekannter Form auch für das dissipative System formulieren: U ^ ( t, t+ t)= ( T ^ D exp i t+ t dt []) H ^ KS ( t )+ H ^ f ( t ) t (6.63) Hier wurde der Operator H ^ f ( t) ohne das Argument KS angegeben. Bis auf das noch zu klärende Problem mit der Zeitabhängigkeit der Operatoren im Exponenten, kann hier die Magnusentwicklung aus Kapitel 6.2.2 auf den Operator 6.63 angewendet werden. Das Problem mit der Zeitabhängigkeit wird durch folgenden Algorithmus gelöst: 1. Auf das Kohn-Sham System den Operator exp[ iH ^ ( t) t] anwenden, um aus den Orbitalen n ( t) den Orbitalsatz jn=0 ( t+ t) zu generieren. Der Orbitalsatz 0 n ( t+ t) stellt eine erste Näherung des Zeitschrittes dar, der für Interpolation genutzt werden kann. 2. Das Kohn-Sham Potential V e j ff=0 ( t+ t) berechnen. 3. Wiederhole bis Konvergenz von V e j ff bzgl. des Iterationsschrittes j eintritt: a) Interpoliere die gröÿen V eff und z /t im Intervall [ t, t+ t] an den Stützstellen t j , welche für den Magnusoperator benötigt werden. Die Gröÿe z /t kann auch mit Hilfe der Tangenten an den Stützstellen t j der Interpolationsfunktion aus z berechnet werden. b) Hamiltonoperatoren H ^ KS und H ^ f initialisieren ' . 19 Elektrostatisches Potential, eektives Potential, etc. berechnen. 128 6.3. Propagatoren für implizite Kohn-Sham Gleichungen c) Wende exp(^ M ) auf den Orbitalsatz n ( t) an, um jn+1 ( t+ t) zu generieren. 4. Der konvergierte Orbitalsatz jn ( t+ t) stellt das Ergebnis des Zeitschritts dar. Das Feld V eff ( r, t) und die Gröÿe z( t) der vorherigen Zeitschritte können zur Verbesserung der Interpolation im Intervall [ t, t+ t] hinzugezogen werden. Wie sich gezeigt hat, ist von Interpolationsverfahren, welche nicht die Zeitumkehrbarkeit berücksichtigen(wie z.B. die Lagrangeinterpolation), abzuraten: Die numerische Stabilität steigt hier erheblich, wenn die Railway-Curve Interpolation(s. Gl. 6.44 in Kapitel 6.2.3) verwendet wird. Diese Interpolationsformel erfordert die Kenntnis der Zeitableitung von den zu interpolierenden Gröÿen bei t und t+ t : Zur Lösung dieses Problems sollten die Gln. 15 und 16 aus[89] herangezogen werden. 6.3.6. Analyse der Konditionszahl In diesem Kapitel werden Eigenschaften der impliziten Kohn-Sham Gleichungen 6.45-6.47 analysiert. Die Betrachtungen beschränken sich hier allerdings auf ein Ein-Teilchensystem mit der von Neuhauser vorgeschlagenen Drude-Dissipation. Wie in folgenden Unterkapiteln gezeigt wird, kann ein kompaktes lineares Gleichungssystem zur Bestimmung der Zeitableitung der Wellenfunktion aufgestellt werden. An diesem wird die Konditionszahl untersucht. Deren Bedeutung in der Numerik[93] wird hier nochmals kurz zusammengefasst: Für ein Gleichungssystem A x= b , das durch die Koezientenmatrix A deniert ist, kann die Lösung x als Ausgabe einer Funktion und die rechte Seite b als deren Eingabe betrachtet werden. Die Konditionszahl gibt nun an, wie stark sich eine kleine Änderung der Eingabe b auf die Ausgabe x( b) auswirkt. Bei linearen Gleichungssystemen ist eine groÿe Konditionszahl schlecht , weil kleine Fehler in b zu groÿen Fehlern in x führen. Zur Denition der Konditionszahl geht man von einem Fehler e in b aus, der zur Lösung als A 1 e beiträgt. Nun betrachtet man das Verhältnis zwischen relativen Fehler in der Lösung(Ausgabe) zu Fehler in der Eingabe: A 1 e / A 1 b e / b Darüber kommt man auf die Denition der Konditionszahl einer Matrix A : ( A )= A 1 · A (6.64) Für hermitesche Matrizen kann das Verhältnis vom gröÿten zum kleinsten Eigenwert genommen werden: ( A )= max ( A ) min ( A ) Motivation Eine wichtige Frage bei der impliziten Ein-Teilchen Schrödingergleichung ist, welcher Zusammenhang zwischen der Magnitude des Dämpfungsparameters a 0 und der Konditionszahl des Gleichungssystems besteht. Es liegt die Vermutung nahe, dass die Konditionszahl 129 6. Numerik sich umso mehr verschlechtert, je höher man die Dämpfung wählt. Dieser Aspekt ist im Zusammenhang mit der adaptiven Dämpfung(Kap. 4.4.3) wichtig, weil hier der Parameter a 0 sehr groÿ werden kann, um die vorgegebene Energieabnahme zu ermöglichen. In Simulationen zeigte sich, dass dieser drei bis vier Gröÿenordnungen durchläuft. Die Kenntnis der Konditionszahlabhängigkeit sollte insbesondere auch erlauben zu verstehen, weshalb die Simulationen mit vorgegebener Energieabnahme(s. Abb. 4.7) plötzlich extrem ungenau werden. Alle diese Aspekte können auch an einem Ein-Teilchensystem untersucht werden: Als Modellsystem kann das eindimensionale, eektive Potential des Metalllms(s. Kap. 3.1) verwendet werden. Als Startwellenfunktion kann eine der Eigenfunktionen dieses Potentials verwendet werden. Genau wie bei dem Vielteilchensystem wird auch dieses System durch ein zeitabhängiges Störpotential V P ( x, t) angeregt. Formulierung in diskreter Basis Die gedämpfte Ein-Teilchen Schrödingergleichung lautet: i ( r, t) t = [] H ^ 0 + H ^ f [, ]( t)+ V P ( r, t) ( r, t) H ^ f [, ]( t)= a 0 j( r, t) · J ^( r) d 3 r t (6.65) (6.66) Die Diskretisierung dieser Gleichung im Ortsraum ergibt ein sehr groÿes, lineares Gleichungssystem für die unbekannte Zeitableitung . In einer Raumdimension kann die Matrix des Gleichungssystems bereits 100... 1000 Zeilen und Spalten haben . Zur vereinfachten Analyse der Konditionszahl wird hier deshalb die Gleichung 6.65 von der kontinuierlichen Ortsraumbasis in eine diskrete Basis transformiert. Als Basis sollen die Eigenzustände {| j } des Hamiltonoperators H ^ 0 dienen. Dadurch ergibt sich auch der Vorteil, dass man besseren Einblick die Dynamik des Systems bzgl. der involvierten Quantenzustände erhält: Wie sich in den folgenden Gleichungen zeigt, wird diese durch die Symmetrieeigenschaften der Eigenzustände bestimmt. Die Wellenfunktion wird zunächst durch einen Koezientenvektor c( t) dargestellt: | ( t) = c l ( t) | l l (6.67) Die Zustände | n sind normiert und orthogonal. Durch Einsetzen in zeitabhängige Schrödingergleichung erhält man: [ ] i c l ( t) | l = H ^ 0 + H ^ f [ c ( t)]+ V ^ P ( t) c l ( t) | l ll (6.68) Um die Bewegungsgleichung für einen bestimmten Koezienten zu erhalten, wird mit 20 Hier kann man ggf. den Vorteil ausnutzen, dass die Matrix dünn besetzt ist. Diese Möglichkeit wurde in dieser Arbeit durch die Formulierung in einer endlichen, diskreten Basis nicht benötigt. 130 6.3. Propagatoren für implizite Kohn-Sham Gleichungen einem Bra-Vektor multipliziert: i c n ( t)= n | H ^ 0 c l ( t) | l + c l ( t) n | H ^ f [ c ( t)] | l l l + c l ( t) n | V ^ P ( t) | l l = E n c n ( t)+ c l ( t) H n f l [ c ( t)]+ c l ( t) V n P l ( t) ll Matrixelemente des Reibungsterms: H n f l [ c( t)]= d 3 q A[ c( t), c ( t)]( q) n | J ^( q) | l = d 3 q A[ c( t), c ( t)]( q) J nl ( q) (6.69) (6.70) Matrixelemente des Stromdichteoperators: J nl ( q)= d 3 q n ( q ) J ^( q) l ( q ) = d 3 q n ( q ) p^ ( q q)+ ( q 2 m e q) p^ l ( q ) = 1 2 m e d 3 q n ( q ) p^ ( q q) l ( q )+ 1 2 m e d 3 q n ( q ) ( q q) p^ l ( q ) 1 = 2 m e d 3 q [ p^ n ] ( q ) ( q q) l ( q )+ 1 2 m e d 3 q n ( q ) ( q q) p^ l ( q ) = 1 2 m e [ p^ n ] ( q) l ( q)+ 1 2 m e n ( q) p^ l ( q) = i 2 m e [ l ( q) n ( q) n ( q) l ( q)] Dämpfungsfeld(hier wird a( q)= a 0 gesetzt): A[ c( t), c ( t)]( q)= a 0 t j[ c( t)]( q) = a 0 t ( t) | J^( q) | ( t) = a 0 t c n ( t) c l ( t) J nl ( q) nl = i 2 m a 0 t c n ( t) c l ( t)[ l ( q) n ( q) n ( q) l ( q)] (6.71) nl In Gl. 6.69 einsetzen: i c ( t)= E c ( t)+ c l ( t) a 0 d 3 q t c n ( t) c m ( t) J nm ( q) · J l ( q) l nm + c l ( t) V n P l ( t) (6.72) l 131 6. Numerik Die nachfolgenden Betrachtungen konzentrieren sich auf den Reibungsterm. Zur Übervorerst sicht wird deshalb ab hier der Störpotentialterm weggelassen: i c ( t)= E c ( t)+ a 0 c l ( t) t [ c n ( t) c m ( t)] d 3 q J nm ( q) · J l ( q) l n,m = E c ( t)+ a 0 c l ( t)[ c n ( t) c m ( t)+ c n ( t) c m ( t)] d 3 q J nm ( q) · J l ( q) l n,m = E c ( t)+ a 0 c l ( t)[ c n ( t) c m ( t)+ c n ( t) c m ( t)] nlm (6.73) l n,m In der letzten Zeile wurde der Tensor durch das Integral in der Zeile darüber deniert. Die Elemente von werden nun für eine Basis aus den rein reellen Eigenfunktionen n ermittelt: nlm = = ( i ) 2 i 2 2 m 2 e 4 m 2e d 3 q d 3 q[ m ( q) n [ m ( q) n ( q) ( q) · l n ( q) m ( q) ( q) ( q)] · [ l ( q) n ( q) m ( q) · ( q) ( q) l ( q)] l ( q) ( q) = i 2 2 4 m 2e m ( q) n ( q) d 3 q[ m ( q) l · ( q) l ( q)+ ( q) n ( q) · n ( q) m ( q) · ( q) l ( q) n ( q) l ( q) m ( q) ( q)] · ( q) m ( q) ( q) n ( q) · l ( q)+ n ( q) ( q) m ( q) · l ( q)] I Zur Abkürzung für die Integrale wird die Gröÿe eingeführt, deren Indizes aus der letzten Gleichung sofort ersichtlich werden: 2 nlm = 4 m 2e [ I mln I nlm I mnl + I nml ] Die Parität der Eigenfunktionen erlaubt es, anhand der Indizes von I zu entscheiden, ob das Integral verschwindet oder nicht. Innerhalb der ersten beiden und letzten beiden Indexpaare von I liegt eine Symmetrie bei Vertauschung vor. Die Gleichung 6.73 repräsentiert ein gekoppeltes Gleichungssystem für die Komponenten des Koezientenvektors c . Die Dreifachsumme kann noch weiter vereinfacht werden, was für die eziente numerische Simulation der Gleichungen wichtig ist. Dazu werden folgende Eigenschaften vom -Tensor ausgenutzt: mln = nlm nlm = nlm nln = 0 nm = 0 mln = mln (6.74a) (6.74b) (6.74c) (6.74d) (6.74e) Die Summation in Gl. 6.73 kann über die Beziehungen 6.74c-6.74d sofort eingeschränkt werden: i c ( t)= E c ( t)+ a 0 c l ( t)[ c n ( t) c m ( t)+ c n ( t) c m ( t)] nlm l n,m l = n = m 132 6.3. Propagatoren für implizite Kohn-Sham Gleichungen Weil sich das Vorzeichen vom-Tensor bei Vertauschung zweier Indizes(jeweils oben oder unten) nach den Beziehungen 6.74a-6.74b ändert, kann die Doppelsumme über ( n, m) folgendermaÿen umgeschrieben werden: i c ( t) = E c ( t)+ a 0 c l ( t)[ c n ( t) c m ( t)+ c n ( t) c m ( t) c m ( t) c n ( t) c m ( t) c n ( t)] nlm l l = n,m n>m = E c ( t)+ a 0 c l ( t)[ c n ( t) c m ( t) c n ( t) c m ( t)+ c n ( t) c m ( t) c n ( t) c m ( t)] nlm l = l n n > ,m m = E c ( t)+ 2 i a 0 [Im { c n ( t) c m ( t) } + Im { c n ( t) c m ( t) } ] c l ( t) nlm (6.75) n,m n>m l l = Nimmt man nun noch das seit Gl. 6.73 vernachlässigte Störpotential wieder mit hinzu, erhält man als Ergebnis für die Ein-Teilchen Schrödingergleichung in Eigenbasis von H ^ 0 folgende Bewegungsgleichung der Koezienten c : i c ( t)= E c ( t)+ c l ( t) V P l ( t) l + 2 i a 0 [Im { c n ( t) c m ( t) } + Im { c n ( t) c m ( t) } ] c l ( t) nlm n,m n>m l l = (6.76) 133 6. Numerik Ergebnisse Um den Zusammenhang zwischen der Dämpfungskonstante a 0 und der Konditionszahl zu ermitteln wurde die zeitabhängige Ein-Teilchen Schrödingergleichung 6.65, 6.66 in diskreter Basis für das in Abbildung 4.4 gezeigte System gelöst. Das System wurde(wie in Abb. 4.4 dargestellt) durch einen Puls über das Potential V P angeregt. Während der Propagation wurde die Dämpfungskonstante schrittweise erhöht, wie im Text zu Abb. 6.3 erläutert wird. Diese schrittweise Vergröÿerung der Gröÿenordnung von a 0 ndet, basierend auf Beobachtung, auch bei der adaptiven Regulierung(s. Kap. 4.4.3) dieses Parameters statt. Die Ergebnisse in Abb. 6.3 bestätigen den vermuteten Zusammenhang zwischen a 0 und ( A ) . Abbildung 6.3.: Im oberen Diagramm ist die Anzahl an Iterationen dargestellt, die benötigt wurde um das Gleichungssystem des impliziten Runge-Kutta Verfahrens(s. Kap. 6.3.4) iterativ zu lösen. Im unteren Diagramm ist die reziproke Konditionszahl 1 des Gleichungssystems 6.76 jeweils gegen die Simulationszeit t aufgetragen. Nach jeweils 5 fs wurde die Dämpfungskonstante um Faktor 10 erhöht. Die schwarzen, senkrechten Linien markieren die Intervalle: a 0 = 0, 1, 10, 100, 1000, 10000 134 6.4. Methode der Finiten Volumen 6.4. Methode der Finiten Volumen Die Methode der Finiten Volumen(FVM ) wird zur Lösung hyperbolischer, partieller Dierentialgleichungen genutzt[94, 95]. Diese Gleichungen beschreiben häug bestimmte Erhaltungsgröÿen wie z.B. Masse, Energie oder Impuls(d.h. es handelt sich um Kontinuitätsgleichungen). Diese treten vor allem im Bereich der Fluiddynamik auf, in dem sich die FVM als ein Standardverfahren etabliert hat. Um solche Gleichungen numerisch lösen zu können, bietet diese Methode die Möglichkeit, die partielle Dierentialgleichung durch Umformung mittels des Divergenztheorems in ein äquivalentes und exaktes System von gewöhnlichen Dierentialgleichungen zu übersetzen. Diese beschreiben dann die zeitliche Entwicklung der Mittelwerte von einzelnen Zellen(nite Volumen) des Simulationsraums. nicht Die Erhaltungsgröÿen können dabei verletzt werden. Durch den FVM-Ansatz wird das ursprüngliche Problem durch zwei neue Probleme ersetzt: An den Zelloberächen muss ein Integral über den Fluss numerisch berechnet werden und die Gleichungen müssen mit einem numerischen Integrator in der Zeit propagiert werden. Die Berechnung des Flusses stellt die gröÿte technische Herausforderung dar und beeinusst maÿgeblich die Eigenschaften der Lösung: Ein typischer Fehler in numerischen FVM-Lösungen ist das Auftreten von unphysikalischer Diusion, welche durch das Schema zur Berechnung des Flusses verursacht wird. Es sei nochmals betont, dass die Erhaltungsgröÿen auch bei schlechten Approximationen der Flüsse nicht verletzt werden können, wie im Kapitel 6.4.2 gezeigt wird. Die FVM ndet in dieser Arbeit Anwendung bei der Beschreibung des Elektronengases über ein semiklassisches Hydrodynamikmodell(s. Kap. 6.4.4) und bei der Lösung der Wigner Gleichungen(s. Kap. 6.5.5 und 6.4.5). Zur Berechnung der Flüsse an den Zelloberächen kommt ein relativ neuartiges Verfahren zum Einsatz(s. Kap. 6.4.3), welches über einen geringen Diusionsfehler verfügt[96]. Wissenschaftliches Neuland wird bei der Anwendung dieses Verfahrens auf die Wigner-Gleichungen betreten. 6.4.1. Problemübersicht Euler-Gleichungen Bei der Verwendung eines hydrodynamischen Modells des Elektronengases in metallischen Strukturen müssen folgende Probleme gelöst werden: 1. Der elektronische Grundzustand n 0 ( r) muss ermittelt werden. 2. Die hydrodynamischen Gleichungen müssen die Metall-Vakuum Grenzächen korrekt wiedergeben und müssen daher auch im Vakuumbereich gelöst werden. Hier besteht das Problem, dass einige Gröÿen undeniert sein können, wie z.B. das Geschwindigkeitsfeld, welches sich als Quotient von Impuls- und Teilchendichte berechnen lässt. 21 engl.: Finite Volume Method 22 Die Elektronendichte tritt an Oberächen aus dem Festkörper aus. Es ist daher falsch, die Bewegungsgleichungen für das Elektronengas nur innerhalb des Festkörpers zu lösen. 135 6. Numerik 3. An der Metall-Vakuum Grenzäche hat die Elektronendichte einen groÿen Gradienten(s. Abb. 3.2), welcher spezielle Anforderungen an das Schema zur Berechnung der Zelloberächenüsse und die Diskretisierung stellt. Für die Berechnung der Teilchendichte im Grundzustand kann man eine Zeitpropagation durchführen: Man erweitert die Bewegungsgleichung um einen viskosen Term, der das Elektronengas abbremst und simuliert, wie sich eine geratene Startteilchendichte (z.B. die Ionendichte n ( r, t= 0)= n + ( r) ) in einen stationären Zustand mit Kräftegleichgewicht entwickelt. Dieses herrscht im Grundzustand zwischen dem Druck und der elektrostatischen Anziehung. Wigner-Gleichungen Die Verwendung der Wigner-Gleichungen zur Beschreibung des Elektronengases im Bereich der Plasmonik ist bisher wissenschaftlich gemieden worden: Die Notwendigkeit einen sechs-dimensionalen Phasenraum zu diskretisieren, ist ein wesentlicher Grund, weshalb alternative physikalische Beschreibungen gesucht werden[23, 65]. Nichtsdestotrotz soll hier versucht werden, diese Gleichungen direkt zu lösen, um später eine Möglichkeit zu haben, vereinfachte Modelle besser auf ihre Gültigkeit überprüfen zu können. Es müssen dabei folgende Probleme beachtet werden: 1. Die Wignerverteilung des elektronischen Grundzustandes ist unbekannt und muss auf gleiche Weise wie bei dem Hydrodynamik-Modell im Zeitbereich ermittelt werden. Allerdings ist unklar, wie man in diese Gleichungen einen viskosen Term einbaut. In Kapitel 5.1.7 wurde ein Vorschlag gemacht, wie sich eine im k -Raum symmetrische Verteilungsfunktion durch Dämpfung von Asymmetrien erzielen lässt. Diese Methode wurde allerdings in der Praxis noch nie getestet. 2. Die Anfangsverteilung, wie sie in Kapitel 6.5.2 vorgeschlagen wird, bedingt sowohl bei der Ladungsdichte im Ortsraum eine Stufe an der Metall-Vakuum Grenze als auch eine Stufe in der Verteilungsfunktion am Fermivektor im Impulsraum(an festem Ort). Das Vorhandensein von Unstetigkeiten in der Verteilungsfunktion haben die Verwendung des Kurganov-Tadmor Schemas(s. Kap. 6.4.3) zur Berechnung des Flusstensors motiviert. 6.4.2. Grundlagen Die Gleichungen, für welche die Methode der Finiten Volumen gedacht ist, haben folgende Form: u+ · f( u)= q+ Randbedingungen t (6.77) Die Komponenten des Vektors u sind orts- und zeitabhängige Felder, welche dem durch diese Gleichung ausgedrückten Erhaltungsgesetz unterliegen. Die Felder beschreiben den Zustand des Systems, weshalb man u auch als Zustandsvektor bezeichnen kann. Bei f 136 6.4. Methode der Finiten Volumen handelt es sich um den sogenannten Fluÿtensor, der eine Funktion des Zustands u darstellt. Bei einem Zustandsvektor mit K Komponenten und einem d -dimensionalen Raum hat der Fluÿtensor die Form einer K × d Matrix. Auf der rechten Seite der Gleichung hyperbolisch steht ein optionaler Quellterm q . Die Gleichung ist, wenn die Gradienten u f 1 ,..., u f d (bei denen es sich um d verschiedene k × k -Matrizen handelt) diagonalisierbar sind und über rein reelle Eigenwerte verfügen[94]. Um die Gleichung 6.77 in ein System von gewöhnlichen Dierentialgleichungen umzuwandeln, wird zunächst der Simulationsraum in N Zellen j zerlegt: N = j j=1 ( j untereinander disjunkt ) (6.78) Man bildet dann über jede dieser Zellen das Volumenintegral der Gleichung 6.77 und teilt durch deren Volumen: 1 | j | j u d 3 r t + 1 | j | j · f( u) d 3 r= 1 | j | q d 3 r j (6.79) Im ersten Integral kann die Zeitableitung vor das Integral gezogen werden, und bei dem zweiten Integral kann das Divergenztheorem angewendet werden: d 1 dt | j | j u d 3 r+ 1 | j | j f( u) · n ^ d= 1 | j | q d 3 r j (6.80) Das erste und letzte Integral kann durch die Zellmittelwerte ersetzt werden: d 1 dt u ¯ j + | j | f( u) · n ^ d= q ¯ j j (6.81) Das Oberächenintegral kann bei einfacher Geometrie der Zellen auf sinnvolle Weise als Summe der Flüsse F durch Teilächen ! geschrieben werden: 1 | j | j f( u) · n ^ d= 1 | j | F ( kj ) k (6.82) Der K -dimensionale Vektor F ( kj ) beschreibt den Fluss durch die k -te Oberäche des j -ten Volumenelementes. Die zentrale Gleichung der Finite Volumen Methode lautet: d dt u ¯ j + 1 | j | k F ( kj ) = q ¯ j (6.83) Diese Gleichung ist eine gewöhnliche Dierentialgleichung in der Zeit und enthält keine Näherungen. Erst bei der numerischen Berechnung der Flüsse F ( kj ) werden Näherungen erforderlich. Weil aber zwei benachbarte Zellen identische Flüsse an der gemeinsamen Oberäche haben, kann selbst eine schlechte Näherung bei der Berechnung der Gröÿen F ( kj ) zu keiner Verletzung von Erhaltungsgröÿen führen: Der Verlust einer Zelle ist immer gleich dem Zuwachs einer benachbarten Zelle. 23 Bsp.: Oberächen eines Quaders. 137 6. Numerik 6.4.3. Kurganov-Tadmor Methode Die Kurganov-Tadmor Methode[96] wird benötigt um die Flüsse in der Gleichung 6.83 zu berechnen. Diese Methode hat gegenüber den wesentlich einfacheren Methoden, wie z.B. UDS " und CDS # (s.[95]), den Vorteil nur einen geringen numerischen Diusionsfehler aufzuweisen und keine unphysikalischen Oszillationen in den Lösungen zu erzeugen. Vor allem kann die Methode benutzt werden, wenn die Lösung groÿe Gradienten aufweist. Um die Methode zu beschreiben wird nun von Gleichung 6.83 ausgegangen und angenommen, dass der Fluss in y - und z -Richtung verschwindet und die Oberächen im j -ten Volumen in x -Richtung den Abstand x j haben. Die Gleichung kann dann in folgender Form geschrieben werden: d u j dt + 1 x j [] f( u j+1/ 2 ) f( u j 1/ 2 )= 0 (6.84) Man benötigt zunächst eine Flussbegrenzungsfunktion $ ( r j ): R K R K , deren Argument r j eine Art Rauigkeit von u bei der Zelle j beschreibt: r j := u j u j 1 u j+1 u j Numerisch darf diese Gröÿe komponentenweise die Werte ± inf annehmen. Die Implementation der Funktion ( r j ) muss darauf vorbereitet sein. Die Berechnung der Flüsse erfordert weitere Hilfsgröÿen: u Lj +1/ 2 = u j + 0. 5 ( r i )( u j+1 u j ) u Rj +1/ 2 = u j+1 0. 5 ( r j+1 )( u j+2 u j+1 ) u Lj 1/ 2 = u j 1 + 0. 5 ( r j 1 )( u j u j 1 ) u Rj 1/ 2 = u j 0. 5 ( r j )( u j+1 u j ) Geometrisch betrachtet handelt es sich bei diesen Gröÿen um lineare Extrapolationen des Zustandsvektors u auf der linken ( L) bzw. rechten ( R) Seite der Zelloberäche bei j ± 1/ 2 . Die maximale Steigung in der Extrapolation wird durch die Funktion begrenzt. Die Flüsse in Gleichung 6.84 werden nun wie folgt berechnet: f( u j 1/ 2 ) f( u j+1/ 2 ) = = 1 2 1 2 {[ f {[ f ( u Rj 1/ 2 ) ( u Rj +1/ 2 ) + + f f ] ( u Lj 1/ 2 ) ] ( u Lj +1/ 2 ) a j 1/ 2 a j+1/ 2 [ u Rj 1/ 2 [ u Rj +1/ 2 ]} u Lj 1/ 2 ]} u Lj +1/ 2 Die Gröÿen a j ± 1/ 2 bezeichnen jeweils die betragsmäÿig gröÿten Eigenwerte der JacobiMatrix der Flusstensoren zwischen den Zellen j und j ± 1 : a j ± 1/ 2 = max[ ( f( u j )), ( f( u j ± 1 ))] (6.85) Dabei steht ( A) für den Spektralradius der Matrix A . Die Eigenwerte geben die lokale Ausbreitungsgeschwindigkeit im Raum an. Um die Gleichungen in der Zeit zu propagieren, kann aus dem maximalen a j -Wert auf eine Obergrenze für die Zeitschrittgröÿe entsprechend der Courant-Bedingung a t/ x< C max geschlossen werden. 24 upwind dierence scheme 25 central dierence scheme 26 ux limiter function In der englischen Fachliteratur entsprechend als bekannt. Diese Funktionen müssen problemabhängig gewählt werden und sind der Literatur zu entnehmen. 138 6.4. Methode der Finiten Volumen 6.4.4. Euler-Gleichungen Die Euler-Gleichungen der Fluiddynamik lauten ganz allgemein(s. z.B.[95]): + · ( v)= 0 t ( v) + · ( v ( v))+ p= 0 t E + · ( v( E+ p))= 0 t (6.86a) (6.86b) (6.86c) Diese Gleichungen stehen für Massen-, Impuls- und Energieerhaltung. Die Felder Massendichte , Impulsdichte v (Geschwindigkeitsfeld v ) und Energiedichte E bilden den Zustandsvektor u : u= v x v y v z (6.87) E Die Energiedichte E kann noch weiter aufgeteilt werden: E= u+ 1 2 m ( v x 2 + v y 2 + v z 2 )+ v ext ( r, t) (6.88) Diese Beiträge sind: Innere-+ kinetische-+ potentielle Energiedichte. Um die Gleichungen zu vervollständigen werden noch eine kalorische und eine thermische Zustandsgleichung benötigt. Hier dient als Modelluid für das Elektronengas ein Ideales Fermigas bei T= 0 K [21, 62]. Der Druck p ist dabei physikalisch durch das Pauli-Prinzip bedingt. Für den Druck p und die innere Energiedichte u= dU/dV gilt in Abhängigkeit der Teilchendichte n= /m e : p= 2 U 3 V = (3 2 ) 2/ 3 5 m e | q e | 2 n 5/ 3 dU u== 2 dV 3 2 (3 2 ) 5/ 3 10 2 m e n 5/ 3 (6.89) (6.90) Der Flusstensor hat im Fall der Eulergleichungen folgende Form: f( u)= v x p+ v x 2 v x v y v x v z v x ( E+ p) v y v x v y p+ v y 2 v y v z v y ( E+ p) v z v x v z v y v z p+ v z 2 =( f x , f y , f z ) v z ( E+ p) (6.91) Die Abbildung 3.2 zeigt ein Ergebnis für den Grundzustand, welches mittels eindimensionaler Finite Volumen Methode für das ideale Fermigas im Metalllm berechnet wurde. Das Ergebnis wurde auf folgende Weise generiert: 139 6. Numerik 1. Das Kurganov-Tadmor Schema war zum Zeitpunkt dieser Rechnung noch nicht implementiert worden. Daher wurden die Flüsse mittels eines einfachen Zentraldifferenzenschemas(CDS, s.[95]) berechnet. Weil dieses Schema(im Gegensatz zu dem KT-Schema) mit einer unstetigen Startteilchendichte der Form n ( r, t= 0)= n + ( r) (welche bei dem Metalllm stufenförmig an der Oberäche abfällt) nicht klar kommt, wurde eine geglättete Startdichte benutzt. 2. Es wurden periodische Randbedingungen gewählt, weil diese am einfachsten die Reexionen am Rand des Simulationsraumes verhindern und für das Ergebnis unerheblich sind. 3. Die Zeitintegration der Gl. 6.83 erfolgte mit einem Runge-Kutta Verfahren 4. Ordnung. Damit das Fluid von einem elektromagnetischen Feld in Bewegung versetzt werden kann, muss in die Kontinuitätsgleichung für die Impulserhaltung 6.86b noch ein Quellterm eingefügt werden. Auf gleiche Weise kann auch eine Senke für den Impuls ergänzt werden, welche die Reibung beschreibt. Die Gleichung lautet dann: ( v) t + · ( v ( v))+ p= q e [ E+ v × B] ( v) (6.92) Auf der rechten Seite steht die Lorentzkraftdichte und ein Dämpfungsterm mit Dämpfungskonstante . Letzterer wird insbesondere für die Berechnung der Grundzustandsdichte benötigt. 6.4.5. Wigner-Gleichungen Hier wird Vorarbeit für das Kapitel 6.5 geleistet: Die Wigner-Gleichungen 5.30a-5.30e werden durch elementare Umformungen, welche auf Anwendung der Produktregel des Nablaoperators basieren, in die Form einer Kontinuitätsgleichung 6.77 gebracht % . Das Ergebnis dieser Umformungen lautet: t f ( r, k) kin = m r · [ k f ( r, k)] (unverändert)(6.93a) t f ( r, k) A. p = q m [ k · ( f ( r, k) r [ A( r) · k])] + q m [ r ·{ A( r) f ( r, k) } ] t f ( r, k) A 2 = 2 q 2 m k · [ r { f ( r, k) | A( r) | 2 ] } (6.93b) (6.93c) t f ( r, k) H = q k · [ E( r) f ( r, k)] t f ( r, k) F = | q | ( r · [ f ( r, k) k ( k, r)] k · [ f ( r, k) r ( k, r)]) (6.93d) (6.93e) 27 Der Korrelationsterm 5.30f wird ausser Acht gelassen, weil dessen genaue Realisierung noch gar nicht klar ist. 140 6.5. Wigner-Maxwell Gleichungen In dieser Form lässt sich der Zustandsvektor u und der Flusstensor f( u) ablesen. Der Zustandsvektor ist einfach nur durch die Wignerfunktionen f gegeben: () u( r, k)= f e ( r, k) f i ( r, k) (6.94) Der Flusstensor hat entsprechend der zwei Komponenten von u und wegen des sechsdimensionalen Phasenraums die Form einer 2 × 6 -Matrix: ( f( u)= F r e F k e ) F r i F k i (6.95) Dabei sind F r und F k jeweils dreidimensionale Zeilenvektoren, welche sich auf den Ortsbzw. Impulsraumanteil des Phasenraums beziehen: F r ( r, k)= m k f ( r, k) q m A( r) f ( r, k) | q | f ( r, k) k ( k, r) F k ( r, k)= q m f ( r, k) r [ A( r) · k] 2 q 2 m r { f ( r, k) | A( r) | 2 } + q E( r) f ( r, k)+ | q | f ( r, k) r ( k, r) (6.96a) (6.96b) 6.5. Wigner-Maxwell Gleichungen Die numerische Lösung der Wigner-Maxwell Gleichungen 5.35a-5.35h stellt ein anspruchsvolles Problem dar, für das es kein universelles Lösungsverfahren gibt. Ein Versuch, diese mit Finiten Dierenzen Verfahren zu lösen, ist in[57] dokumentiert. Die Gleichungen für Licht(Maxwell) und Materie(Wigner) werden hier getrennt behandelt, da die selbstkonsistente Kopplung der Gleichungen vorerst als ein unbedeutendes Detail erscheint & . In diesem Kapitel wird es ausschlieÿlich um das Lösen der Wignergleichungen 5.35a5.35e für ein vorgegebenes elektromagnetisches Feld gehen. Für das Lösen der Maxwellgleichungen in Potentialform 5.35f-5.35h wird ein FDTD-Verfahren in Kapitel 6.7 vorgeschlagen, welches auch auf andere Licht-Materie Systeme angewendet werden kann, in denen die Materiegleichungen eine Coulomb-Eichung der elektromagnetischen Potentiale erfordern. 6.5.1. Problemübersicht Mit Hilfe der Abbildung 6.4 lässt sich leicht beschreiben, welches Resultat die numerische Simulation der Wignergleichungen bringen soll: Hier wird die Zeitentwicklung einer klassischen Verteilungsfunktion im Phasenraum illustriert, welche als Analogon zur quantenmechanischen Wignerverteilung betrachtet werden kann. Zur Zeit t 0 bende sich das 28 Bei der Kopplung der Gleichungen ist erfahrungsgemäÿ mit numerischen Stabilitätsproblemen zu rechnen. 141 6. Numerik Abbildung 6.4.: Zeitentwicklung der Verteilungsfunktion im Phasenraum. Bei klassischen Verteilungsfunktionen bleibt das blau markierte Volumen, welches sich aus den Mikrozuständen ( p, q) , die mit der Präparation kompatibel sind, zusammensetzt, während der Zeitentwicklung konstant[62]. Bei der Wignerverteilung gilt das nur für das Integral über den ganzen Raum. System im Grundzustand und wird anschlieÿen durch externe Felder angeregt. Dabei durchläuft es verschiedene Zustände bis es nach einer Zeit t 2 durch Relaxation wieder in den Grundzustand, den es bei der Zeit t 0 gehabt hat, zurückkehrt. Aus der Trajektorie im Phasenraum lassen sich Erwartungswert b( t) , Streuung b( t) und höhere Momente aus der Verteilungsfunktion zu Observablen b berechnen. Dabei handelt es sich z.B. um die Teilchen- und Stromdichte. Die Berechnung solcher Gröÿen ist das Ziel der Simulation. Für die Initialisierung der Simulation muss zunächst durch makroskopische Präparation die Verteilungsfunktion zur Zeit t 0 ermittelt werden. Dazu müssen alle mikroskopischen Zustände, die mit der Präparation kompatibel sind, mit einer entsprechenden Gewichtung versehen werden. In Abbildung 6.4 sind diese Zustände als klassische Mikrozustände ( p, q) in zusammenhängenden Gebieten dargestellt. Bei der Wignerverteilung sind die Mikrozustände durch Vielteilchenwellenfunktionen | gegeben. Als nächstes muss die Zeitentwicklung der Verteilungsfunktion berechnet werden. Speziell für dieses Problem wurde zuvor in Kapitel 6.4.5 gezeigt, dass sich die WignerGleichungen 5.35a-5.35e über einen Flusstensor f( u) in die Form der Kontinuitätsgleichung 6.77 bringen lassen: t ( f e ( r, k) f i ( r, k) )( + · F r e F k e ) F r i F k i =0 (6.97) Wenn man von den Schwierigkeiten absieht, welche die Berechnung des Flusstensors beinhaltet ' , fällt die Gleichung 6.97 in die Klasse von hyperbolischen partiellen Differentialgleichungen für die es unzählige FVM-basierte Lösungsansätze in der Numerik gibt. Die Vielfalt der Lösungsansätze hat den Grund, weil so einfache Gleichungen wie die Advektionsgleichung als auch die hochkomplizierten Navier-Stokes Gleichungen in die 29 Man denke dabei insbesondere an die Energierenormierung ( r, k) . 142 6.5. Wigner-Maxwell Gleichungen Abbildung 6.5.: Der Impulsraumanteil der Wignerverteilung gleicht im Anfangszustand einer Fermikugel. Die(normierte) Funktion fällt an der Oberäche der Kugel von Eins auf Null ab. Die blau markierte Schale soll den Bereich k min < k< k max darstellen. Ausserhalb von diesem Bereich bleibt die Verteilungsfunktion konstant(Annahme, s. Text). Form 6.97(evtl. mit Quellterm) gebracht werden können und es unzählige Anwendungen für die Gleichungen gibt. Jede Anwendung bringt besondere Ansprüche an bestimmte Teilaspekte mit sich, wie z.B. die Realisierung von Randbedingungen oder die Diskretisierung des Simulationsraumes. Für die Wignergleichungen werden diese in den folgenden Kapiteln beschrieben. 6.5.2. Präparation des Grundzustands Im Jellium-Modell wird von einer statischen Ionendichte n + ( r) ausgegangen. Daher reduziert sich die Gleichung 6.97 auf den elektronischen Anteil: () t f e ( r, k)+ · F r e F k e =0 (6.98) Zukünftig kann der Index e an der Verteilungsfunktion daher weggelassen werden: f e f Als Anfangszustand wurde die Teilchendichte n ( r, 0)= n + ( r) gesetzt und angenommen, dass das Geschwindigkeitsfeld v( r, 0) überall verschwindet. Dazu muss der k -Raumanteil der Ein-Teilchen Wignerverteilung f( r, k) an jedem Ort r mit einer Fermikugel(s. Abb. 6.5) initialisiert werden, deren Radius die Bedingung 1 f( r, k)= n + ( r) V k (6.99) erfüllt. 143 6. Numerik Abbildung 6.6.: Der Impulsraum kann in kartesischen oder sphärischen Koordinaten diskretisiert werden. Das Innere der Fermikugel braucht nicht mit gespeichert werden, da die Funktion dort konstant Eins ist. Ebenso muss fernab der Fermikugel die Verteilung nicht gespeichert werden, da sie dort konstant Null ist. Die Farbwerte geben die Mittelwerte der Wignerverteilung in den einzelnen Zellen(nite Volumen) wieder. In der Mitte der Schale (d.h. bei k= k F ) fällt die Verteilungsfunktion auf Null ab. 6.5.3. Diskretisierung des Phasenraums Der Phasenraum, auf dem die Ein-Teilchen Wignerverteilung f( r, k) deniert ist, hat sechs Raumdimensionen. Dieser kann für realistische Berechnungen nicht einfach durch N Gitterpunkte je Dimension aufgelöst werden, da der Speicherbedarf mit N 6 steigt. Stattdessen wird hier von einer unbewiesenen Vermutung Gebrauch gemacht, durch die zwar die Diskretisierung des Phasenraums erheblich komplizierter wird, aber gleichzeitig auch das Speicherplatzproblem gelöst werden kann: In jeder Fermikugel des k -Raumes gibt es einen Radius k min unterhalb dem die Wignerverteilung f konstant den Wert f 0 ( r) behält. Genauso existiert ein Radius k max oberhalb dem die Verteilung konstant(Null) bleibt. Diese Annahme kann nicht auf alle Systeme zutreen und ausserdem kann sie während der Zeitentwicklung verletzt werden. Für metallische Systeme wird angenommen, dass diese Aussage stimmt und sich die Wignerverteilung nur innerhalb der Kugelschale k min < k< k max ändert. Im k -Raum muss dadurch nur noch die Kugelschale diskretisiert werden, wie in Abbildung 6.6 gezeigt wird. Naheliegender Weise kann man hier neben einem kartesischen Gitter auch ein sphärisches Gitter in Betracht ziehen. Welches besser geeignet ist, müssen numerische Simulationen erst noch zeigen. Diese Lösung zur Einsparung von Speicherplatz macht aber auch spezielle Randbedingungen an den Rändern der Kugelschale im k -Raum notwendig. Ein erster Versuch besteht darin, hier folgende Dirichletrandbedingungen zu verwenden: f( r, | k | = k min )= f 0 ( r) f( r, | k | = k max )= 0 (6.100) (6.101) 144 6.5. Wigner-Maxwell Gleichungen Diese Randbedingungen sind zwar sehr einfach zu implementieren, führen aber zu unphysikalischen Reexionen am Rand der Kugelschale. 6.5.4. Zeitentwicklung Die Berechnung der Flüsse durch die Oberächen von Volumenelementen(s. Abb. 6.6) sollte bei den Wignergleichungen über die in Kapitel 6.4.3 beschriebene Kurganov-Tadmor Methode durchführbar sein. Die Zeitentwicklung kann dann durch ein explizites Runge-Kutta Verfahren berechnet werden. In dieser Arbeit konnte die hier vorgeschlagene Methode allerdings nicht mehr getestet werden. In der Literatur wurde für das Wigner-Poisson System eine numerische Methode[59] gefunden, welche die Gleichung für eine Raumdimension(d.h. in einem zweidimensionalen Phasenraum) durch eine doppelte Fouriertransformation löst. Die Anwendung einer Finiten Volumen Methode auf die Wignergleichungen ist vermutlich neu. Unklar ist, ob es überhaupt korrekt ist, eine Finite Volumen Methode auf eine Quasiverteilungsfunktion wie der Wignerverteilung anzuwenden. 6.5.5. Drude Bewegungsgleichung im k-Raum Die klassische Drude Bewegungsgleichung für ein Elektron im elektrischen Feld m e v = q e E( t) m e v (6.102) kann durch einen Spezialfall der Wignergleichung 5.35b reproduziert werden. Um die klassische Drude Bewegungsgleichung zu reproduzieren, wird ein fester Ort r der Verteilungsfunktion betrachtet und alle anderen Terme aus der Bewegungsgleichung für f gestrichen: f( k) t = q e ( r ( r) · ( k f( k)) = q e E( r) · ( k f( k)) = q e k · [ E( r) f( k)] (6.103) Um eine phänomenologische Relaxation mit der gleichen Wirkung wie der des klassischen Terms m e v zu erhalten, werden nur Abweichungen f aus der Anfangsverteilung f 0 betrachtet: f( k, t)= f( k, t) f 0 ( k) (6.104) Indem man diesen Term in die Gleichung 6.103 einsetzt, kann dann für ein Finites Volumen die Bewegungsgleichung(in der Form der Gleichung 6.83) formuliert werden: d dt f ¯ j = q e 1 | j | s F j ( s) f ¯ j (6.105) Für den Fluss ( q e /) E( r) f( k) (s. Gl. 6.103) durch die s -te Oberäche des j -ten Volumenelementes wird eine Näherung benötigt: F j ( s) E( r) · n ^ ( js ) f ¯( k ( j 1 ,s ) )+ 2 f ¯( k ( j 2 ,s ) ) (6.106) 145 6. Numerik Abbildung 6.7.: Ergebnis der FVM Simulation zur Drude Bewegungsgleichung 6.103. Die blaue Kurve im Vordergrund zeigt den zeitlichen Verlauf des elektrischen Feldes, welches für eine Beschleunigung sorgt. Die rote Kurve zeigt den Geschwindigkeitsverlauf, der aus der Verschiebung der Fermikugel resultiert. Die schwarzen Kreuzchen zeigen die Übereinstimmung mit der numerischen Lösung der Drude Bewegungsgleichung 6.102, die mit einem Standardverfahren zur Zeitintegration erhalten wurde. Mit n ^ ( js ) wird die Oberächennormale bezeichnet und f ¯( k ( j 1 ,s ) ) , f ¯( k ( j 2 ,s ) ) bezeichnen die f Mittelwerte der Verteilungsfunktion in den beiden Volumina, welche sich die Oberäche teilen. Der Dämpfungsterm f ¯ j in Gl. 6.105 ist künstlich hinzugefügt worden. Das Ergebnis der Simulation nach dieser speziellen Finite Volumen Methode, bei der die Mittelwertbildung 6.106 zur Berechnung des Flusses benutzt wurde, ist in Abbildung 6.7 dargestellt. Als Finite Volumen kamen sowohl die kartesischen als auch die sphärischen Volumen aus Abbildung 6.6 zum Einsatz. Für kleine Anregung durch das Feld E konnte kein Unterschied in den Ergebnissen festgestellt werden ! . 6.6. Lösungsverfahren für die Poisson-Gleichung 6.6.1. Problemübersicht Die Poissongleichung[27] tritt in den Kohn-Sham Gleichungen bei der Berechnung des Hartree-Potentials V H [ n ] und bei den Wigner-Maxwellgleichungen im Zusammenhang mit dem longitudinalen Anteil des elektromagnetischen Feldes(Gl. 5.31) auf: [ 2 + 2 ] ( r)= 1 ( r)+ RB 0 (6.107) Die Randbedingungen können je nach System völlig verschieden sein: 30 Die maximale Auslenkung der Fermikugel aus der Ruhelage im Verhältnis zur Dicke der Schale(s. Abb. 6.6) legt fest, ob es sich um eine kleine oder groÿe Anregung des Systems handelt. 146 6.6. Lösungsverfahren für die Poisson-Gleichung · Bei elektronischen Strukturrechnungen an einzelnen Atomen und Molekülen würde man oene Randbedingungen wählen: lim( | r | )= 0 | r | (6.108) · Für die Simulation von periodisch auf einem Substrat angeordneten Split-Ring Resonatoren[1] wären periodische Randbedingungen in zwei Dimensionen und eine oene Randbedingung in der dritten Dimension sinnvoll: ( x ± w, y, z)=( x, y, z) ( x, y ± h, z)=( x, y, z) lim( x, y, z)= 0 | z | (6.109a) (6.109b) (6.109c) · Für Festkörper werden rein periodische Randbedingungen benötigt: ( r+ R )=( r)( R : Position der Elementarzelle mit Index ) (6.110) Um die Gleichung 6.107 numerisch zu lösen, werden üblicherweise Finite Dierenzen oder Finite Elemente Verfahren verwendet[78], bei denen die Poissongleichung in ein lineares Gleichungssystem umgewandelt wird. Für die Lösung der dabei entstehenden Gleichungssysteme existiert eine gewaltige Fülle von Methoden[69] deren Anwendbarkeit u. a. von der Gröÿe des Gleichungssystems abhängt. Die Anzahl der Unbekannten in gängigen Anwendungen deckt mindestens den Bereich von 10 2 bis 10 7 ab. Die Gleichung kann für den rein periodischen Fall 6.110 über eine Fourierreihenentwicklung besonders einfach und numerisch extrem ezient gelöst werden(s. Kapitel 6.6.2). Die Realisierung von oenen Randbedingungen ist numerisch gesehen am schwierigsten, da der Simulationsraum nur endlich groÿ ist, aber die gesuchten Lösungen sich auf einen unendlich groÿen nicht-periodischen Raum beziehen. Eine einfache, aber nicht sehr elegante Methode zur Realisierung von oenen Randbedingungen besteht darin, rein periodische oder Dirichletrandbedingungen( ( )= 0 ) zu wählen und die Ladungen so weit vom Rand zu entfernen, dass die Lösung nicht mehr vom Abstand zum Rand abhängt. oene In[97] wird gezeigt, wie sich auch Randbedingungen mit einer schnellen Fourier Methode realisieren lassen. 6.6.2. FFT-Methode Die Möglichkeit, die Poissongleichung mittels schneller Fouriertransformation(FFT ! ) technische lösen zu können, bietet folgende entscheidende Vorteile gegenüber Verfahren, welche ein lineares Gleichungssystem lösen müssen: 1. Die Berechnung der Lösung erfordert eine feste Anzahl an Rechenoperationen, d.h. es wird kein iteratives Verfahren wie beim Lösen von linearen Gleichungssystemen benötigt. 31 engl: Fast Fourier Transform 147 6. Numerik N 2. Für die schnelle diskrete Fouriertransformation eines-dimensionalen Vektors werden nur O ( N log N) Rechenoperationen benötigt. Die Methode eignet sich also auch für groÿe zwei- und dreidimensionale Probleme, bei denen die Anzahl der Gitterpunkte N in der Gröÿenordnung 10 4 ... 10 6 liegen kann. fftw 3. Zur Durchführung der schnellen Fouriertransformation kann auf die-Bibliothek zurückgegrien werden, in der die ezientesten Verfahren implementiert sind. Die Bibliothek bietet auch Möglichkeiten zur Parallelisierung mit OpenMP und MPI. Für die Herleitung der Gleichungen der FFT-Methode genügt es, die eindimensionale Poissongleichung zu betrachten: [ 2 z 2 ] 2 ( z)= 4 ( z) ( 0 =(4 ) 1 in a.u.)(6.111) Das Potential und die Ladungsdichte werden in einer Fourierreihe auf dem Intervall [0, L] entwickelt: [ 2 z 2 ] 2 c n e ik n z n= = c n [( ik n ) 2 2 ] e ik n z n= = 4 r n e ik n z n= 2 n k n := L Durch Koezientenvergleich ergibt sich folgende Bedingung für die Koezienten c n : c n = 4 r n k n 2 + 2 (6.112) Für ladungsneutrale Systeme ist der Koezient r 0 = 0 . In diesem Fall muss auch der Koezient c 0 verschwinden, da gilt: c 0 2 = 4 r 0 (6.113) Der Koezient c 0 , der nur eine physikalisch unbedeutende additive Konstante zum Potential darstellt, ist oensichtlich bei nicht-neutralen Systemen nur dann deniert, wenn es eine elektrostatische Abschirmung gibt. In dieser Arbeit werden nur elektrisch neutrale Systeme betrachtet und daher gilt c 0 = r 0 = 0 . Die numerische Umsetzung dieser Methode verwendet eine Diskretisierung des Simulationsraumes mit einem regulären, kartesischen Gitter. Die Fourierintegrale zur Berechnung der Koezienten r n werden in der Diskreten Fouriertransformation[66, 98] durch Riemann-Summen auf dem Gitter ersetzt. Insgesamt sind folgende drei Schritte nötig, um die Poissongleichung auf diese Weise zu lösen: 1. Diskrete Fouriertransformation der Ladungsdichte durchführen 2. Koezienten c n nach Gl. 6.112 berechnen 3. Inverse Diskrete Fouriertransformation auf die Koezienten c n anwenden, um die Potentialfunktion zu erhalten 148 6.7. Lösungsverfahren für elektromagnetische Potentiale in Coulomb-Eichung 6.7. Lösungsverfahren für elektromagnetische Potentiale in Coulomb-Eichung Die Formulierung von quantenmechanischen Bewegungsgleichungen berücksichtigen das elektromagnetische Feld über den sogenannten Minimale Kopplungsterm im Hamiltonoperator: H ^ = 1 2 m e ( p ^ q e A) 2 + q e (6.114) Die Eichung der Potentiale kann frei gewählt werden. Die Motivation für die Wahl der Coulomb-Eichung besteht darin, dass die Wignergleichungen 5.35a-5.35e speziell in dieser Eichung hergeleitet wurden. Diese Eichung wird sehr häug bei der Beschreibung von Optik in Festkörpern verwendet und es scheint derzeit noch kein numerisches Lösungsverfahren zu existieren, um elektromagnetische Potentiale in dieser Eichung zu simulieren. In diesem Kapitel wird daher versucht ein Verfahren zu entwickeln, welches die Potentialform der Maxwellgleichungen speziell in der Coulomb-Eichung lösen kann: 2 A 1 c 2 2 A t 2 = µ 0 J t 2 = 1 0 · A= 0 (6.115a) (6.115b) (6.115c) Grob gesagt muss dazu einmal eine vektorielle Wellengleichung 6.115a für das VektorpoA tential und eine Poissongleichung 6.115b für das longitudinale elektrische Feld gelöst werden. Das Vektorpotential darf die Eichbedingung 6.115c nicht verletzen. Die Verwendung der Lorentzeichung[27] erscheint einfacher, weil bei dieser auch für das skalare Potential eine Wellengleichung gelöst werden muss und somit nur eine Art von Problem zu lösen ist. Ein numerisches Verfahren für die Lorentzeichung existiert bereits und ist in[44] beschrieben. 6.7.1. Problemübersicht Die Entwicklung einer neuen numerischen Methode zur Lösung der Maxwellgleichungen ist eine sehr komplexe Aufgabe. Bei den Maxwellgleichungen tritt zusätzlich das Problem auf, dass ein kompliziertes Modell für die Materie(welches über nichtlineares und/oder nichtlokales Antwortverhalten verfügen kann) simultan mit den Gleichungen gelöst werden muss. Bei der Entwicklung einer Methode müssen folgende Eigenschaften der numerischen Methode berücksichtigt werden: · Konsistenz: Für t 0 und/oder x 0 gehen die diskretisierten Gleichungen über in die exakten Gleichungen. · Stabilität: Die numerische Methode ist stabil, wenn sich kleine Fehler in der numerischen Lösung(von Zeitschritt zu Zeitschritt) nicht verstärken. · Konvergenz: Die numerische Lösung der Dierenzengleichung konvergiert gegen die Lösung der exakten Dierentialgleichung für immer feinere Gitterauösungen. 149 6. Numerik · Erhaltungsgesetze: In Abwesenheit von Quellen und Senken dürfen Erhaltungsgröÿen von der numerischen Methode nicht verletzt werden. · Einhaltung von Schranken: Schranken für physikalische Gröÿen dürfen nicht verletzt werden. Z.B. darf eine Teilchendichte nicht negativ werden. · Realisierbarkeit: Man sollte garantieren können, dass das Modell und die zugehörigen Gleichungen überhaupt eine physikalisch sinnvolle Lösung besitzen. Dieses Problem ist nicht numerischer Natur. · Genauigkeit: Es gibt drei Arten von Fehlern, welche die Genauigkeit der numerischen Lösung charakterisieren: 1. Modellierungsfehler: Beschreibt das mathematische Modell das reale System hinreichend genau? Falls nicht, können auch auskonvergierte Lösungen qualitativ vom experimentell beobachteten Verhalten abweichen. 2. Diskretisierungsfehler: Das ist der Unterschied zwischen der exakten Lösung der Dierentialgleichung und der algebraischen Dierenzengleichung. 3. Konvergenzfehler: Das ist der Unterschied zwischen der exakten Lösung eines algebraischen Gleichungssystems und der numerisch-iterativ ermittelten Lösung. Diese Auistung von Eigenschaften orientiert sich an der Darstellung in[95], wo es speziell um Simulationen in der Fluiddynamik geht. Da die Maxwellgleichungen auch eine Kontinuitätsgleichung für elektrische Ladungen implizieren, sind die oben genannten Aspekte (Erhaltungsgesetze, Einhaltung von Schranken) auch für das Maxwellproblem relevant. Weil die Komplexität dieser Aufgabe das Hauptproblem darstellt, wird ein pragmatischer Lösungsansatz gewählt: Die Wellengleichung für das Vektorpotential wird über das bekannte FDTD-Verfahren[41, 45] gelöst, welches bereits ausgiebig numerisch untersucht wurde. Folgende Probleme konnten in dieser Arbeit nicht mehr betrachtet werden: Die Poissongleichung sollte mit unterschiedlichen Randbedingungen gelöst werden können ! und es sollte überprüft werden, dass die Retardierung in der Coulomb-Eichung[27, 99] in der Simulation korrekt wiedergegeben wird. Letzteres Problem ist vermutlich sehr kompliziert, weil die winkelabhängige Gitterdispersion des FDTD-Verfahrens berücksichtigt werden muss um das longitudinale Feld(welches über die Poissongleichung beschrieben wird) mit dem transversalen Feld(welches über die Wellengleichung beschrieben wird) korrekt zu vereinen. In diesem Zusammenhang sei noch erwähnt, dass der Quellterm J t in der Vektorpotentialwellengleichung 6.115a nur den transversalen Anteil der Stromdichte J darstellt. Die Stromdichte, wie sie vom Materiemodell ausgegeben wird, setzt sich aber aus longitudinalen und transversalen Anteilen zusammen: J= J l + J t (6.116a) · J t = 0 (6.116b) × J l = 0 (6.116c) 32 Es werden vor allem oene Randbedingungen in mindestens einer Raumrichtung benötigt, um beispielsweise Transmissionsspektren von Nanostrukturen(s. Kap. 1.2.3) berechnen zu können. 150 6.7. Lösungsverfahren für elektromagnetische Potentiale in Coulomb-Eichung Der longitudinale Anteil muss also raus gerechnet werden. Dieser steht in folgendem Zusammenhang mit dem elektrostatischen Potential : 1 ( ) µ 0 J l = c 2 t (6.117) Um die Felder in longitudinale und transversale Anteile zu zerlegen, könnten sich die distributionsartigen Projektionsoperatoren(Gln. 15 und 19) in[99] als hilfreich erweisen. 6.7.2. Lösung der Wellengleichung Die vektorielle Wellengleichung 6.115a wird hier zunächst als System von drei unabhängigen, skalaren Wellengleichungen aufgefasst. Die folgende Betrachtung reduziert sich deshalb auf ein skalares Feld A( r, t) , welches einer inhomogenen Wellengleichung genügt: 2 A 1 c 2 2 A t 2 = J (6.118) (Der Quellterm wurde zur Vereinfachung der Gleichung umdeniert: In dieses J muss eine Komponente von µ 0 J t eingesetzt werden.) Diese Gleichung soll nun mit einem Finite Dierenzen Verfahren numerisch gelöst werden. Dazu wird als erstes diese Gleichung, welche zweiter Ordnung sowohl in Raum als auch in Zeit ist, in ein System von Gleichungen erster Ordnung umgeschrieben: t A = c · v c 2 J( t ) dt t t 0 v = c A t (6.119a) (6.119b) Für diese Umformung muss ein vektorielles Hilfsfeld v eingeführt werden. Der Quellterm steht nun in einem Zeitintegral !! . Die Felder A und v werden in Analogie zum Yee-Schema [45] im Simulationsraum auf zwei zueinander versetzten Gittern angeordnet, wie in Abb. 6.8 gezeigt wird: Die beiden Gitter werden jeweils durch die rot und blau markierten Diskretisierungspunkte gebildet. Die Gitterabstände seien x , y und z . Der Yee-Cube mit den Indizes ( j, k, l) bendet sich am Ort r ( j,k,l) = j x e x + k y e y + l z e z mit den Komponenten v x ( j,k,l) v y ( j,k,l) , v z ( j,k,l) A ( j,k,l) 33 Der Quellterm muss die Bedingung J( t< t 0 )= 0 erfüllen. 151 6. Numerik Abbildung 6.8.: Der Simulationsraum wird mit sogenannten Yee-Cubes diskretisiert: Dazu werden die Felder A, v x , v y , v z auf regulären, kartesischen Gittern dargestellt, welche den gezeigten Versatz zueinander besitzen. Die ausgefüllten Punkte stellen eine Elementarzelle von Diskretisierungspunkten dar, welche im Raum periodisch fortgesetzt wird. Eine solche Zelle beinhaltet die drei Komponenten des Hilfsfeldes v (rote Punkte) und ein Gitterpunkt, der zum A -Feld gehört(blauer Punkt). Die Dierentialoperatoren in den Gln. 6.119a-6.119b können auf diesem speziellen Gitter durch Zentraldierenzen, welche 2. Ordnung genau sind, dargestellt werden: · v ( j,k,l) v x ( j+1 ,k,l) v x ( j,k,l) + v y ( j,k+1 ,l) v y ( j,k,l) + v z ( j,k,l+1) v z ( j,k,l) (6.120a) x y z 1 v x ( j,k,l) A ( j 1 ,k,l) A ( j,k,l) c t x 1 v y ( j,k,l) A ( j,k,l) A ( j,k 1 ,l) c t y (6.120b) (6.120c) 1 v z ( j,k,l) A ( j,k,l) A ( j,k,l 1) c t z (6.120d) Die Komponenten des Gradienten vom A -Feld existieren nur an den Orten der v -Feldkomponenten, weshalb die Gln. 6.120b-6.120d mit der zugehörigen v -Feldkomponente formuliert wurden. Die Komponenten des Yee-Cubes sind nicht nur räumlich sondern auch zeitlich versetzt angeordnet: Auf der diskretisierten Zeitachse t n = n t existiert das A -Feld nur zu den Zeitpunkten t n und das v -Feld nur zu den Zwischenpunkten t n+1/ 2 . Dadurch können für die Zeitableitung ebenfalls Zentraldierenzen 2. Ordnung verwendet werden: A t ( t n+1/ 2 )= A( t n+1 ) A( t n ) + O ( t 2 ) t v t ( t n )= v( t n+1/ 2 ) v( t n 1/ 2 ) + O ( t 2 ) t (6.121a) (6.121b) 152 6.7. Lösungsverfahren für elektromagnetische Potentiale in Coulomb-Eichung Auf der linken Seite ersetzt man nun die Zeitableitungen A t und v t durch die räumlichen Dierenzenterme 6.120a-6.120d gemäÿ den Gleichungen 6.119a-6.119b. Die resultierenden Gleichungen können nach A( t n+1 ) und v( t n+1/ 2 ) aufgelöst werden, wodurch eine explizite Leapfrog Integration der Gleichungen in der Zeit möglich ist(-Verfahren: s.[41]). Es werden als Anfangsbedingungen die Felder A( t 0 ) und v( t 1/ 2 ) benötigt. Als letztes muss noch der Quellterm in Gl. 6.119a diskretisiert werden. Das Integral t Q( t)= c 2 J( t ) dt (6.122) t 0 muss zu den Zeiten t n+1/ 2 deniert sein, da es an der gleichen Stelle wie das v -Feld in der Zeitableitung vom A -Feld auftritt(s. Gl. 6.119a). Indem das Feld J( t) zu den Zeiten t n deniert wird, kann das Integral als einfache Riemann-Summe geschrieben werden: n Q( t n+1/ 2 )= c 2 t J( t k )+ O ( t 2 ) k=0 (6.123) Die Ortsabhängigkeit vom Feld J wurde bisher nicht explizit erwähnt: Das Feld muss an den gleichen Gitterpunkten im Raum wie das A -Feld deniert werden(da es als Quellterm in der Wellengleichung steht). Die bisherige Darstellung hat sich auf eine Komponente der vektoriellen Wellengleichung 6.115a bezogen. Im Prinzip kann man nun für jede Komponente A x , A y , A z entsprechende Hilfsfelder v x , v y , v z einführen und mit dem beschriebenen Verfahren das A -Feld in der Zeit propagieren. Hier stellt sich die Frage, ob tatsächlich zwölf Felder benötigt werden, um drei skalare Wellengleichungen zu lösen: Wie man vom Yee-Algorithmus zur Lösung der Maxwellgleichungen, welche ja die Wellengleichungen implizieren, weiÿ, werden dafür nur die sechs Felder E und B benötigt. Um herauszunden, welche Modikationen an bereits existierenden FDTD-Implementationen für Maxwellgleichungen nötig sind, um damit die Vektorpotentialwellengleichung 6.115a lösen zu können, werden nochmals die Wellengleichungen für die E - und B -Felder (mit Quellterm J ) betrachtet: 2 E µ 0 0 2 E t 2 = J µ 0 t 2 B µ 0 0 2 B t 2 = µ 0 × J t (6.124a) (6.124b) Das E -Feld würde man als A -Feld benutzen und das B -Feld übernimmt die Rolle der Hilfsfelder v x , v y , v z . Das E -Feld hat allerdings die Zeitableitung von µ 0 J als Quelle. Die Zeitintegration der Quelle(Gl. 6.123) wird also weiterhin benötigt. Sonstige Änderungen sind am Yee-Algorithmus nicht erforderlich. Es muss allerdings noch beachtet werden, dass die Komponenten der E - und J -Felder nicht am selben Ort sondern auf dem Yee-Cube verteilt liegen. Das kann für Anwendungen, in denen die Rotation des Vektorpotentials benötigt wird von Vorteil sein, weil der Yee-Cube extra für diese Berechnung gedacht ist. Andernfalls müssen Mittelwerte gebildet werden, welche alle Komponenten der Felder auf einen gemeinsamen Punkt im Yee-Cube konzentrieren. 153 7. Schlussfolgerung und Ausblick Für die Berechnung der optischen Eigenschaften von Nanostrukturen wurde im Rahmen der Dichtefunktionaltheorie zunächst der elektronische Grundzustand der Nanostrukturen durch Lösen der stationären Kohn-Sham Gleichungen berechnet: Diese Aufgabe konnte erfolgreich für Strukturen durchgeführt werden, deren Symmetrien ein-, zweiund dreidimensionale numerische Berechnungen erfordern. Die eindimensionalen Berechnungen am Metalllm haben gezeigt, dass bereits bei einer Filmdicke von nur 5 nm die elektronische Teilchendichte an den Oberächen mit der Teilchendichte eines Metallhalbraumes sehr gut übereinstimmt(s. Abb. 3.1). Die Teilchendichte ist durch FriedelOszillationen charakterisiert, welche im Inneren des Metalllms(bzw. Halbraums) verschwinden. Diese Beobachtung kann genutzt werden, um die Dicke einer Metalloberäche auf mikroskopischer Ebene abzuschätzen, auf der sich deutliche Abweichungen von idealisierten Grenzächen der makroskopischen Theorien zeigen. Ebenfalls konnte mit Hilfe von zweidimensionalen Rechnungen an Nanodrähten die mikroskopische Teilchendichte an Metallkanten und-ecken mit 90 -Prol berechnet werden(s. Abb. 3.10 u. 3.11). Die dreidimensionalen Berechnungen erlauben es in analoger Weise auch die Teilchendichte an Metallspitzen zu ermitteln. Bei diesen Modellrechnungen macht sich allerdings der hohe Rechenaufwand stark bemerkbar, wenn man versucht die Abmessungen der Struktur solange zu erhöhen bis die Dichte an den Spitzen unabhängig von den Abmessungen wird. Solche Berechnungen können in einer zukünftigen Arbeit durchgeführt werden. Dazu muss allerdings eine Rechenumgebung genutzt werden, in der 10 3 ... 10 4 Kohn-Sham Orbitale ezient parallel verarbeitet werden können. Mit den Simulationen im Zeitbereich konnten Ergebnisse zu den nicht-linearen und nicht-lokalen Eekten in den Nanostrukturen erzielt werden: An den Oberächen des Metalllms konnten die Stromdichten der zweiten und dritten Harmonischen beobachtet werden. Deren räumliche Lokalisierung und deren Skalierungsverhalten(s. Abb. 3.4 u. 3.6) hat sich als physikalisch plausibel erwiesen, so dass davon auszugehen ist, dass die Berechnungsmethode und deren sehr umfangreiche Implementation korrekte Ergebnisse liefert. Die Berechnung von SH-Strömen wurde ebenfalls an diversen metallischen Nanodrähten durchgeführt(s. Abb. 3.12) um zu untersuchen, welche geometrischen Formen eine besonders hohe SH-Stromdichte erzeugen. Um aussagekräftige Ergebnisse zu erhalten, muss in zukünftiger Arbeit die Strukturgröÿe solange erhöht werden, bis die Stromdichte unabhängig von dieser wird. In diesem Zusammenhang müssen auch oene Randbedingungen(statt den derzeit verwendeten periodischen Randbedingungen) in Betracht gezogen werden, welche technisch im Allgemeinen schwierig zu realisieren sind. Die Berechnung der nicht-lokalen Suszeptibilität der Elektronendichte des Metalllms im linearen Regime ist ebenfalls gelungen: Diese ist durch eine starke Ortsabhängigkeit 1 Die Beobachtung von noch höheren Harmonischen ist durch Auftreten von numerischem Rauschen begrenzt. 155 7. Schlussfolgerung und Ausblick charakterisiert, welche im Inneren des Metalllms abnimmt. Dadurch konnte eine Art Bulkbereich um das Zentrum des Metalllms ausgemacht werden. Alle Berechnungen danachträglich zu unterliegen allerdings einer an die Teilchendichte angefügten, sehr starken Dämpfung, ohne die diese Berechnung nicht durchführbar ist: Wenn man die Dämpfung weglässt, ergeben sich unendlich lange andauernde Dichtewellen, die sich durch die Nanostrukturen bewegen und an den Grenzächen reektiert werden. Eine Antwortfunktion kann für ein solches System(mit unendlich zurückreichendem Gedächtnis) nicht berechnet werden. Dieser Aspekt hat auch die Entwicklung von Methoden zur dissipativen Dichtefunktionaltheorie innerhalb dieser Arbeit motiviert: Als Grundlage wurde die Methode von Neuhauser verwendet und unter verschiedenen Aspekten weiterentwickelt und untersucht. Zu den Weiterentwicklungen zählen zum einen die Lösungsverfahren für implizite KohnSham Gleichungen(s. Kap. 6.3) und zum anderen die alternativen Formulierungen des Reibungsterms(dessen Herleitung im Impulsraum sich als äuÿerst aufwendig erwiesen hat). Wie bereits in der Originalarbeit von Neuhauser konnte das Phänomen beobachtet werden, dass die Dämpfungsezienz nach wenigen Femtosekunden stark abnimmt. Um diesem Phänomen nachzugehen, wurden die impliziten Kohn-Sham Gleichungen für eine Basis mit geringer Dimensionalität(bestehend aus den Eigenfunktionen des eektiven Potentials im Grundzustand) formuliert. Die Gleichungen erlauben in dieser Form einen viel besseren Einblick in die Dynamik des Systems, da Eigenschaften, wie z.B. die Symmetrien von Matrixelementen bestimmter Operatoren, direkt sichtbar werden. Es konnte damit erfolgreich gezeigt werden, dass ein wesentliches Problem beim Lösen der impliziten Kohn-Sham Gleichungen darin besteht, dass sich die Konditionszahl mit Fortschreiten der Simulation immer weiter verschlechtert. Die Ursache für die Abnahme der Dämpfungsezienz lässt sich in einer zukünftigen Arbeit über diesen Basisansatz vermutlich ebenfalls ergründen. Ein sehr nützliches und praktisches Ergebnis für die Verwendung der Methode von Neuhauser besteht in der Erkenntnis, dass sich statt des Stromdichteoperators genauso gut der Teilchendichteoperator benutzen lässt: Der Rechenaufwand wird dadurch(besonders in zwei- und dreidimensionalen Simulationen) geringer ohne dass es einen qualitativen Unterschied in den Ergebnissen gibt. Als Alternative zur Dichtefunktionaltheorie wurden die Wignergleichungen aus einer unveröentlichten Arbeit von W. Hoyer untersucht: Es konnte in dieser Arbeit gezeigt werden, dass sich diese Gleichungen in die Form einer Kontinuitätsgleichung im sechsdimensionalen Phasenraum bringen lassen. Die Wignerfunktion kann man sich für den Grundzustand des Systems(d.h. die Elektronen einer Nanostruktur) so vorstellen, dass man an jedem Punkt im Ortsraum eine sphärische Fermiverteilungsfunktion im Impulsraum hat. Die Annahme, dass sich die Funktion in ihrem Impulsraumanteil nur an der Oberäche dieser Fermikugel zeitlich verändert, lässt hoen, dass die numerische Verarbeitung dieser Wignerfunktion durch eziente Speicherung möglich ist. Weil die Funktion über groÿe Gradienten im Ortsraum(an Grenzächen von Nanostrukturen) und im Impulsraum(an der Oberäche der Fermikugel) verfügt, wurde nach einem numerischen Lösungsverfahren für Kontinuitätsgleichungen gesucht, deren Lösung diese Eigenschaften haben darf. Die Finite-Volumen Methode von Kurganov und Tadmor(s. Kap. 6.4.3) hat sich als möglicher Kandidat erwiesen. Um die Komplexität von physikalischer und numerischer Seite möglichst gering zu halten, wurden in dieser Arbeit die Gleichungen für 156 einen Nanodraht hergeleitet und diverse Vereinfachungen vorgenommen(z.B. Beschränkung auf das niedrigste Subband und auf elektrostatische Potentiale). Dadurch ist man zu den eindimensionalen Wigner-Poisson Gleichungen(für einen zweidimensionalen Phasenraum) gelangt, welche in der Literatur bereits ausgiebig beschrieben wurden. Bereits diese Gleichungen zeigen interessante Phänomene, wie z.B. die Landau-Dämpfung, welche mit einem neu zu entwickelnden numerischen Lösungsverfahren erst einmal nachvollzogen werden müssen. In einer zukünftigen Arbeit könnte man den Einsatz von Finite-Volumen Methoden(z.B. der von Kurganov und Tadmor) für diese Aufgabe untersuchen, da es derzeit in der Literatur noch keine Lösungsverfahren zu geben scheint, die diese für Wignergleichungen verwenden. In der Literatur ist allerdings bereits beschrieben worden, wie sich aus den Wigner-Poisson Gleichungen die Quanten Euler Gleichung herleiten lässt(s. Kap. 5.3.2). Die Lösungen dieser Gleichung sind zwar nur auf Längenskalen oberhalb der Fermiwellenlänge gültig(und zeigen z.B. daher auch keine Friedel-Oszillationen), berücksichtigen aber trotzdem Quanteneekte bei einem Rechenaufwand, der mit dem vom Lösen der klassischen Euler Gleichung vergleichbar ist. In einer zukünftigen Arbeit sollte speziell auf die Bedeutung des Bohm-Potentials(s. Gl. 5.89) eingegangen werden, da im Vorhandensein dieses Terms der wesentliche Unterschied zu derzeitigen hydrodynamischen Modellrechnungen im Bereich der Plasmonik liegt. Ein weiterer Aspekt dieser Arbeit bestand darin das bekannte FDTD-Verfahren zur Simulation elektromagnetischer Felder auf die elektromagnetischen Potentiale, speziell in der Coulomb-Eichung, zu übertragen: Es konnte gezeigt werden, dass bereits existierender FDTD-Code nur in der Berechnung der Quellterme verändert werden muss, um die Wellengleichung für das Vektorpotential zu lösen. Ein noch ungelöstes Problem stellt die Behandlung des skalaren Potentials dar, welche als Lösung der Poissongleichung vorliegt: Das Vektorpotential unterliegt einer winkelabhängigen Gitterdispersion, die nicht zu der Lösung des skalaren Potentials passt. Es ist sehr schwierig abzuschätzen, wie sich ein solcher Fehler auf die numerische Stabilität auswirkt. Insgesamt stellt diese Arbeit einen viel versprechenden Ausgangspunkt für die theoretische Untersuchung der optischen Eigenschaften von metallischen Nanostrukturen dar: Die Ergebnisse der Berechnungen auf Basis der Dichtefunktionaltheorie können bereits als qualitativ korrekt angenommen werden. Um auch mit experimentellen Resultaten vergleichen zu können, müssen allerdings noch Details, z.B. bei den verwendeten Randbedingungen und der Dissipation, angepasst und verbessert werden. 157 A. Anhang A.1. Zustandsdichten System Partikel Connement Zustandsdichte 3D D( E)= 2 j ( E E j ) Draht 2D D( E) /L= 2 m e · E 1 2 Film 1D D( E) /A= m e 2 Tabelle A.1.: Zustandsdichte von Systemen mit unterschiedlichem Connement. D( E) /L ist die Zustandsdichte pro Länge des Drahtes und D( E) /A die Zustandsdichte des Films pro Fläche. Beim Draht und Film bezieht sich die Zustandsdichte auf ein Subband. A.2. Gröÿen in der linearen Optik von Metallen Die Tabelle A.2 ist so konzipiert, dass die fundamentalen Parameter eines Metalls, welche in Tabellenwerken nachgeschlagen werden können, ganz oben stehen und in die Gleichungen darunter eingesetzt werden können um abgeleitete Gröÿen zu berechnen. 159 A. Anhang DC-Leitwert Elektronendichte Plasmafrequenz Streuzeit Suszeptibilität Dispersionsrelation AC-Leitwert Phasendierenz(E- und B-Feld) Skintiefe DC n 0 P = n 0 q e2 0 m e 1 = m e DC n 0 q e2 ( )= ( ) i 0 k( )= µ 0 0 [1+ () ] 2 ( )= 0 P2 i ( )= arg[ k() ] ( )= Im { k( ) } 1 Tabelle A.2.: Formelsammlung zur linearen Optik in Metallen. A.3. Drude-Parameter für Edelmetalle Gold Silber Platin n 0 [a 0 3 ] n 0 [m 3 ] 0. 00876 0. 008736 0. 002845 5. 9 · 10 28 5. 9 · 10 28 1. 9 · 10 28 r S [a 0 ] 3. 0 3. 0 4. 4 P [rad s 1 ] 1. 4 · 10 16 1. 4 · 10 16 7. 8 · 10 15 DC [S m 1 ] 4. 1 · 10 7 6. 1 · 10 7 5. 1 · 10 6 [s 1 ] 4. 1 · 10 13 2. 7 · 10 13 1. 1 · 10 14 Tabelle A.3.: Plasmafrequenz und Streuzeit aus[100]. r S [a 0 ] 2. 0 3. 0 4. 0 5. 0 6. 0 n 0 [a 0 3 ] 0. 0298 0. 0088 0. 00373 0. 0019 0. 0011 P [10 15 rad s 1 ] 25. 3 13. 8 8. 95 6. 40 4. 87 F [a 0 ] 6. 55 9. 82 13. 10 16. 37 19. 64 Tabelle A.4.: Teilchendichte, Plasmafrequenz und Fermi-Wellenlänge in Abhängigkeit vom Wigner-Seitz Radius. A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum A.4.1. Stromdichteoperator im kontinuierlichen Impulsraum Der Stromdichteoperator, der die Wahrscheinlichkeitsstromdichte im Impulsraum beschreibt, soll hier hergeleitet werden, da dieser in Lehrbüchern der Quantenmechanik normalerweise nicht zu nden ist . 1 In[42] wird die allgemeine Form eines Stromdichteoperators als J ^( q)= i[ H ^ ,( x ^ q)] (in a.u.) angegeben. Dieses Konstrukt ist auch im Impulsraum gültig und kann als alternative Form der hier 160 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum Skintiefe [m] Relative Skintiefe ( : ) 10 -1 10 -2 10 -3 10 -4 10 -5 10 -6 10 -7 10 -8 10 0 10 5 10 10 10 15 Frequenz[Hz] 10 20 10 10 10 5 10 0 10 -5 10 -10 10 0 10 5 10 10 10 15 Frequenz[Hz] 10 20 (a)(b) Abbildung A.1.: Frequenzabhängigkeit der Skintiefe in Gold:(a) Absolute Werte,(b) Verhältnis der Skintiefe zur Wellenlänge im Vakuum. Diese wurde mit den Formeln aus Tabelle A.2 und den Parametern in Tabelle A.3 berechnet. Ausgangspunkt ist die Schrödingergleichung im Impulsraum: i ~( p, t)= p 2 ~( p, t)+ dp V ~ ( p p , t) ~( p , t) t 2 m 2 (A.1) Die fouriertransformierten Gröÿen werden zur Kenntlichmachung der gewählten Konvention der Normierung explizit angegeben: ~( p, t)= 1 dx ( x, t) e ipx/ 2 V ~ ( p, t)= dx e ipx/ V( x, t) (A.2a) (A.2b) Die Verwendung zwei verschiedener Konventionen hat keinen besonderen Grund. Die Wahl einer unitären Transformation für die Wellenfunktion ist wegen der statistischen Interpretation jedoch naheliegend. Für die Rücktransformationen gilt: ( x, t)= 1 dp ~( p, t) e ipx/ 2 1 V( x, t)= dp V ~ ( p, t) e ipx/ 2 (A.3a) (A.3b) gezeigten Herleitung genutzt werden. 161 A. Anhang Analog zur Herleitung im Ortsraum erfolgt diese hier durch Betrachtung der Zeitableitung der Aufenthaltswahrscheinlichkeit des Teilchens im Intervall [ a, b] des Impulsraums: b P( a, b; t)= | ~( p, t) | 2 dp a b b dd P( a, b; t)= | ~( p, t) | 2 dp= | ~( p, t) | 2 dp dt dt t aa b =[ ~ ~ + ~ ~ ] dp a = b 1 [ ~ p 2 ~ + ~ V ^ ~ ~ p 2 ] ~ ~( V ^ ~) dp i 2 m e 2 m e a (A.4) d P( a, b; t)= b 1 [] ~ V ^ ~ ~( V ^ ~) dp dt i a = b [ dp ~ dp V ~ ( p p ) ~( p , t) ~ dp V ~ ( p p ) ~ ( p , ] t) i 2 2 a In Gleichung A.4 wurde der Operator der potentiellen Energie V ^ eingeführt. Auf einen ket-Vektor | angewendet, ergibt sich folgende Funktion im Impulsraum: ( p)= p | V ^ | = dp ( p p ) dp V ~ ( p p ) ~( p , t) 2 (A.5) Fortsetzung der vorherigen Rechnung: d P( a, b; t) dt = b dp 2 i 2 [] dp V ~ ( p p ) ~( p , t) ~ ( p, t) V ~ ( p p ) ~ ( p , t) ~( p, t) a b = dp 2 dp Im { V ~ ( p p ) ~( p , t) ~ ( p, t) } a = J( a, t) J( b, t) Die Stromdichte hat folgende Eigenschaft: lim J( p, t)= 0 p ± 162 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum Unter Ausnutzung dieser Eigenschaft lässt sich die Stromdichte berechnen: d dt P > ( p, t)= dp 2 dp Im { V ~ ( p p ) ~( p , t) ~ ( p , t) } p = J( p, t) (A.6) Der Term kann noch weiter umgeformt werden: d dt P > ( p, t)= dp 2 i 2 dp [ V ~ ( p p ) ~( p , t) ~ ( p , t) ~ ( p , t) ~( p , t)] p = dp 2 i 2 ~ ( p , t)[ V ~ ~]( p ) p dp 2 i 2 ~( p , t)[ V ~ ~ ]( p ) p Über das Faltungstheorem folgt: [ V ~ ~]( p )= = = dp V ~ ( p p ) ~( p , t) 2 dx e ip x/ V( x, t) ( x, t) 2 F [ V( x) ( x)]( p ) (A.7) Die letzten beiden Zeilen dienen der Kenntlichmachung der verwendeten Konvention in der Fouriertransformation mit Symbol F . Für die Stromdichte gilt nun: d dt P > ( p, t)= J( p, t)= i dp { ~ ( p , t) F [ V( x) ( x)]( p ) 2 3 p ~( p , t) F [ V( x) ( x)] ( p ) } = 2 dp Im { ~ ( p , t) F [ V( x) ( x)]( p ) } (A.8) 2 3 p Dieser Term sollte für numerische Berechnungen besonders gut nutzbar sein: Die Strominkrementell dichte im ganzen Raum kann durch Verschieben der Integrationsgrenze berechnet werden. Als nächstes wird der Stromdichteoperator mit Parameter p ermittelt, von dem J( p, t) der Erwartungswert ist: J( p, t) J ^( p) = ( t) | J ^( p) | ( t) (A.9) Wie man leicht nachrechnen kann, ist dieser Operator durch folgenden Term gegeben: J ^[ V ~ ]( p)= 1 [ ( p^ p) V ^ ] V ^ ( p^ p) i (A.10) 163 A. Anhang Darin kommt einmal der Operator p^ und die c-Zahl p vor. Eine Besonderheit stellt die Abhängigkeit von der potentiellen Energie V( x) mit der Fouriertransformierten V ~ ( q) (s. Gl. A.2b) dar, welche die Notation J ^[ V ~ ] motiviert. Das Matrixelement des Stromdichteoperators kann folgendermaÿen in Impulsbasis berechnet werden: | J ^[ V ~ ]( p) | = 1 2 i 2 dp dp [ V ~ ( p p ) ( p ) ( p ) p V ~ ( p p ) ( p ) ( p ) ] (A.11) Zum Vergleich werden die Stromdichteoperatoren für den Orts- und Impulsraum nochmals aufgeführt: Ortsraum: J ^( q)= p^ ( q x^)+ ( x^ q) p^ 2 m e Impulsraum: J ^( p)= ( p p^) V ^ V ^ ( p^ p) i A.4.2. Reibungsterm im kontinuierlichen Impulsraum Es wird die Wirkungsweise des Reibungsterms 4.29 auf eine Wellenfunktion im Impulsraum benötigt: p | H ^ f ( t) | = ds a( s) J( s, t) p | J ^( s) | t Nebenrechnung: p | J ^( s) | = 1 2 i [] 2 dp dp V ~ ( p p ) ( p p) ( p ) V ~ ( p p ) ( p p) ( p ) s = 1 2 i [] 2 dp dp V ~ ( p p ) ( p p) ( p ) V ~ ( p p ) ( p p) ( p ) s = 1 2 i 2 ( p s) dp V ~ ( p p ) ( p ) dp V ~ ( p p) ( p ) s Daraus folgt: p | H ^ f ( t) | 1 = 2 i 2 ds a( s) J( s, t) ( p s) t dp V ~ ( p p ) ( p ) dp V ~ ( p p) ( p ) p s = 1 2 i 2 J ds a( s) t dp V ~ ( p p ) ( p ) J ds a( s) dp V ~ ( p p) ( p ) t -- s 164 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum p = 1 2 i 2 dp V ~ ( p p ) ( p ) J ds a( s) t J ds a( s) dp V ~ ( p p) ( p ) t s p = 1 2 i 2 dp V ~ ( p p ) ( p ) J ds a( s) t J ds a( s) dp V ~ ( p p ) ( p ) t s (A.12) Die erste Faltung(erstes Integral in A.12) kann man umschreiben: dp V ~ ( p p ) ( p ) = dp dx e i[ p p ] x/ V( x, t) dx e ip x / ( x ) 2 = dx e ipx/ V( x, t) dx ( x )2 ( x x ) 2 = 2 F{ V( x, t) ( x) } ( p/) Das letzte Integral in A.12 kann durch eine Thetafunktion als Faltung geschrieben werden: dp V ~ ( p p ) ( p ) s = dp V ~ ( p p )( p s) ( p ) (A.13) = dp ( p s) dx e i[ p p ] x/ V( x, t) dx e ip x / ( x ) 2 = 1 dx dx dp ( p s) e i[ p p ] x/ ip x / V( x, t) ( x ) 2 -- = 1 dx e ipx/ V( x, t) dx ( x ) dp ( p s) e ip [ x x]/ 2 -- 165 A. Anhang Mit q= p s, p = q+ s, dq= dp substituieren: dp ( p s) e ip [ x x]/ = dq( q) e i( q+ s)[ x x]/ - = e is[ x x]/ dq( q) e iq[ x x]/ = e is[ x x]/ ( ( [ x x]/ )) i( x x ) Es folgt: dp V ~ ( p p ) ( p ) s = 1 2 dx dx e ipx/ V( x, t) ( x ) e is[ x x]/ ( i( x x) + ( [ x x]/ )) - () = 1 2 dx dx e i[ p s] x/ V( x, t) ( x ) e isx / i( x x) + ( x x) - Zum x -Integral: dx ... () = dx ( x ) e isx / i( x x) + ( x x) = i dx ( x ) e isx / x x + ( x) e isx/ Letztlich muss noch über x integriert werden: dx e i[ p s] x/ V( x, t) i dx ( x ) e isx / x x - = i dx e i[ p s] x/ V( x, t) dx ( x ) e isx / x x - + dx e i[ p s] x/ V( x, t) ( x) e isx/ + ( x) e isx/ = i dx e i[ p s] x/ V( x, t) dx ( x ) e isx / x x + dx e ipx/ V( x, t) ( x) -- = i dx e i[ p s] x/ { ( x ) } V( x, t) F x x ( s/)+ F{ V( x) ( x) } ( p/) 166 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum Das Integral mit der Singularität ist die Fouriertransformation des Produktes zweier Funktionen, das in die Faltung derer Fouriertransformierten übergeht: F { ( x ) x x } ( k)= F { ( x ) } F { x 1 x } ( k) f( k)= F{ ( x) } ( k)= dx e ikx ( x)= 2 ~( k) {} {} g( k)= F 1 x x ( k)= e ikx F 1 x ( k)= ie ikx sgn( k) ( f g)( k)= i 2 dk ~( [ k k ]) e ik x sgn( k ) Einsetzen in das Integral über x : dx e i[ p s] x/ V( x, t) F { ( x ) x x } ( s ) = i 2 dx e i[ p s] x/ V( x, t) dk ~( s k ) e ik x sgn( k )( r= s k ) - 2 = i dr ~( r) sgn( r s) V ~ ( p r, t) p | H ^ f ( t ) | p = 1 2 i 2 dp V ~ ( p p ) ( p ) J ds a( s) J ds a( s) t t dp V ~ ( p p ) ( p ) = 1 2 i 2 2 s p J J F{ V( x, t) ( x) } ( p/) ds a( s) ds a( s) · t t 1 i 2 dr ~( r) sgn( r s) V ~ ( p r, t)+ F{ V( x) ( x) } ( p/) 2 i p = 1 F{ V( x, t) ( x) } ( p/) J ds a( s) i 2 3 t 1 2 i 2 dr ~( r) V ~ ( p r, t) J ds a( s) t sgn( r s) - 1 1 F{ V( x) ( x) } ( p/) J ds a( s) 2 i 2 3 t 167 A. Anhang [] = 1 F{ V( x, t) ( x) } ( p/) J ds a( s)( p s) 1 i 2 3 t 2 1 2 i 2 dr ~( r) V ~ ( p r, t) J ds a( s) sgn( r s) t - Das ist der Term in Gleichung 4.30. A.4.3. Stromdichteoperator im diskreten Impulsraum Als erstes wird die Schrödingergleichung im periodischen Ortsraum benötigt. Dazu werden die Fourierreihen der Wellenfunktion und des Potentials in die Schrödingergleichung im Ortsraum eingesetzt: ( x, t)= c n ( t) e ik n x n V( x, t)= v n ( t) e ik n x n 2 2 i ( x, t) t = 2 m x 2 ( x, t)+ V( x, t) ( x, t) (A.14) Die Orthogonaliätsrelation lautet: dx e i[ k j k n ] x = | | jn (A.15) Die Schrödingergleichung für die Koezienten der Fourierreihe von ( x, t) lautet: i t c j ( t)= 2 2 m k j 2 c j ( t)+ v n ( t) c j n ( t) n (A.16) Normierungsbedingung: dx | ( x, t) | 2 = 1 dx e i[ k n + k n ] x c n ( t) c n ( t)= 1 n,n | c n ( t) | 2 n = 1 | | (A.17) Die Wahrscheinlichkeit, dass das Teilchen einen Impuls im Interval [ p a , p b ] hat, ist durch folgende Summation gegeben: b P( a, b; t)= | c n ( t) | 2 n= a 168 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum d P( a, b; t) dt = b [ c n c n + c n c n ] n= a [ ] b 1 = i 2 2 2 m k n 2 c n c n + c n v n c n n 2 m k n 2 c n c n c n v n c n n n= a n n = b 1 i [ c n v n c n n c n v n c n n ] = b 2 Im { c n v n c n n } n= a n = n= a n = = 2 b Im { c n [ v c]( n) } n= a Für die Faltung im diskreten Raum gilt: [ v c]( n)= = = || n 11 || 2 v n n d n x ce n d i x = k n e x | V 1 i[ | ( k 2 x n ) n k n d ] x xVd ( xx( ) x ) d | x d 1 x | e e ni k n ik e n i n k x n x e [ x ( i x xk n ) ] x V( x) ( x ) An dieser Stelle tritt der Dirac-Kamm auf: ( x nL)= 1 e ik n x L nn (A.18) Damit folgt: [ v c]( n)= 1 | | dx e ik n x V( x) dx ( x ) ( x x n L) n 1 = dx e ik n x V( x) ( x) | | =: FT[ V( x) ( x)]( n) ( n -ter Fourierkoezient)(A.19) Zur Wahrscheinlichkeitsstromdichte: d P( a, b; t)= J( a, t) J( b, t) dt Dieser Zusammenhang gilt allgemein und folgt bereits aus der Anschauung. Genau wie im nicht-periodischen Fall lässt sich nun eine Grenze gegen unendlich schieben und annehmen, dass die Stromdichte für unendlich groÿe Impulse verschwindet: lim n ± J( p n , t)= 0 169 A. Anhang Für die Stromdichte im diskreten Impulsraum gilt nun: J( p n 1/ 2 , t)= d dt P( p n , ; t)=: d dt P ( p n ; t) = 2 Im { c j [ v c]( j) } (wegen ± 1/ 2 , s. unten) j= n = 2 Im { c j FT[ V( x) ( x)]( k j ) } (A.20a) j= n J( p n+1/ 2 , t)= 2 n 1 Im { c j FT[ V( x) ( x)]( k j ) } j= (A.20b) Im Argument von J steht p n ± 1/ 2 : Dadurch wird sicher gestellt, dass die Summationsgrenzen auf der rechten Seite eindeutig sind(andernfalls wäre nicht klar, ob die Summen bei n oder bei n ± 1 anfangen bzw. enden). Über folgenden Term kann die Kontinuitätsgleichung im diskreten Raum hergeleitet werden: d dt | c n ( t) | 2 + d dt n 1 | c j ( t) | 2 + d dt | c j ( t) | 2 j = j = n +1 = d dt [ | c n ( t) | 2 + { 1 P ( p n , t) } + P ( p n+1 , ] t) = d dt | c n ( t) | 2 J( p n 1/ 2 , t)+ J( p n+1/ 2 , t) Auf der linken Seite steht insgesamt d dt j= | c j ( t) | 2 = 0 . Damit folgt die Kontinuitätsgleichung: d dt | c n ( t) | 2 = J( p n 1/ 2 , t) J( p n+1/ 2 , t) (A.21) Für den Erwartungswert des Stromdichteoperators gilt im kontinuierlichen Raum: J( p, t) J ^( p) = ( t) | J ^( p) | ( t) = dp ~ ( p, t) J ^( p) ~( p, t) Für den diskreten Impulsraum wird nun ein analoger Term gesucht: J( p n 1/ 2 , t)= c n ( t) [ J ^[ V ] ]( p) c n ( t) n= Die Schreibweise J ^[ V] bedeutet, dass dieser Operator von der potentiellen Energie V( x, t) bzw. v n ( t) abhängt. In diesem Impulsraum wirkt der Operator J ^[ V] auf einen unendlich dimensionalen Vektor von Koezienten c= { c n } n= . In Analogie zum nichtperiodischen Fall hat der Operator voraussichtlich folgende Form: J ^[ V]( p n 1/ 2 )= 1 i [ ( p^ p n ) V ^ ] V ^ ( p^ p n ) (A.22) 170 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum (Der Index n 1/ 2 wird am Erwartungswert dieses Operators in Verbindung mit Gleichung A.20a deutlich.) Die Form wird nun überprüft. Dazu wird die Wirkung der Theta-Funktion und des potentielle Energieoperators V ^ auf den Koezientenvektor c benötigt: f=( p^ p n ) g=( p^ p n ) g j e j j= = ( p j p n ) g j e j j= = g j e j j= n (A.23) Dabei wurde angenommen, dass gilt: ( x 0)= 1 . Wirkungsweise vom Operator V ^ : f= V ^ g=[ v g]( n) e n = v n j g j e n n nj (A.24) Damit können nun die Matrixelemente des Stromdichteoperators in der Impulsbasis bestimmt werden und letztlich überprüft werden, ob für das Diagonalelement die Gleichung A.20a resultiert: | J ^[ V]( p n ) | = | p l p l | J ^[ V]( p n ) | p m p m | l,m In der diskreten Impulsbasis seien nun f und g die Darstellungen der abstrakten Vektoren | und | : | := f, f m = p m | | := g, g m = p m | wobei := hier die Bedeutung wird dargestellt durch hat. Damit kann nun die Gleichung für das Matrixelement geschrieben werden als: | J ^[ V]( p n ) | = f, J ^[ V]( p n ) g mit dem komplexen Skalarprodukt x, y = n x n y n . 171 A. Anhang Betrachte nun die Wirkung von J ^[ V] auf g : 1 i [] ( p^ p n ) V ^ V ^ ( p^ p n ) g= 1 i ( p^ p n ) V ^ g 1 i V ^ ( p^ p n ) g = 1 i ( p^ p n ) v m j g j e m 1 i V ^ g j e j m,j j= n = 1 i v m j g j e m 1 i v m j ( e m e j ) g l e l m= n j m,j l= n = 1 i v m j g j e m 1 i v m j g l jl e m m= n j m,j l= n = 1 i v m j g j e m 1 i v m j g j ( j n) e m m= n j m,j = 1 i v m j g j e m m= n j= 1 i v m j g j e m m= j= n Für die Fourierkoezienten der potentiellen Energiefunktion gilt v n = v n : ...= 1 i v m j g j e m m= n j= 1 i v j m g j e m m= j= n = 1 i ( m n) v m j g j e m m= j= 1 i ( j n) v j m g j e m m= j= = 1 i [ ( m n) v m j g j ( j n) v j m g j ] e m m= j= (A.25) Also gilt für das Matrixelement(im 2. Term m und j vertauschen): f, J ^[ V]( p n ) g = 1 i [ ( m n) f m v m j g j ( m n) f j v m j ] g m m= j= (A.26) Diagonalmatrixelement(Erwartungswert) betrachten: | J ^[ V]( p n ) | = f, J ^[ V]( p n ) f =1 i [ ( m n) f m v m j f j ( m n) f j v m j ] f m m= j= 2 = Im { ( m n) f m f j v m j } m= j= Wegen Im[ a+ b]= Im[ a]+ Im[ b] gilt: ...= 2 m= Im ( m n) f m j v m j f j = 2 Im { f m [ v m= n f]( m) } 172 A.4. Mathematische Grundlagen für Reibungsterme im Impulsraum Nach Gl. A.19 gilt für die Faltung und somit letztlich für den Erwartungswert der Stromdichte: | J ^[ V]( p n ) | = 2 Im { f m FT[ V( x) ( x)]( k m ) } m= n = J( p n 1/ 2 , t) Das ist exakt das Ergebnis aus Gl. A.20a. Damit wurde die Form des Stromdichteoperators in Gl. A.22 für periodische Systeme veriziert. A.4.4. Reibungsterm im diskreten Impulsraum Im diskreten Impulsraum hat der Reibungsterm folgende Form: H ^ f ( t)= a( p n ) J( p n , t t) J ^[ V]( p n ) n (A.27) Als erstes muss die Wirkungsweise dieses Operators auf eine Wellenfunktion im diskreten Impulsraum hergeleitet werden. Dazu kann das Matrixelement aus Gl. A.26 benutzt werden: f, J ^[ V]( p n ) c = 1 i [ ( m n) f m v m j c j m= j= ( m n) f j v m j ] c m Dabei wird | durch den Vektor c dargestellt. Um nun die l -te Komponente von J ^[ V] c zu erhalten wird folgender Ansatz gemacht: ql e q , J ^[ V]( p n ) c q = 1 i [ ( m n) ml v m j c j ( m n) jl v m j ] c m m= j= 1 = i ( m n) ml v m j c j 1 i ( m n) jl v m j c m m,j= m,j= = 1 i ( l n) v l j c j 1 i j= ( m n) v m l c m m= 1 = i [ ( l n) v l j c j ( j n) v j l c j ] j= (A.28) 173 A. Anhang Insgesamt erhält man: c~ l = ql e q , H ^ f ( t) c q = a( p n ) J( p n , t t) ql e q , J ^[ V]( p n ) c nq = a( p n ) J( p n , t t) 1 i [( l n) v l j c j ( j n) v l j c j ] n j= = 1 i a( p n ) J( p n t , t) { ( l n)FT[ V( x) ( x)]( l) [ v ( n c)]( l) } = 1 i n a( p n ) J( p n t , t) n { } ( l n)FT[ V( x) ( x)]( l) ( q n) v l q c q q= Der Term in den geschweiften Klammern kann noch mit der Signumfunktion umgeschrieben werden: { ... } = 1 [ sgn+] ( l n)FT[ V]( l)+ 1 FT[ V]( l) 22 1 2 [ q= sgn+] ( q n) v l q c q 1 FT[ V 2 ]( l) = 1 [ 2 sgn+] ( l n)FT[ V]( l) 1 2 [ q= sgn+] ( q n) v l q c q Damit folgt der Term aus Gleichung 4.31: c~ l = 1 i n= a( p n ) J( p n , t t) { } ( l n)FT[ V( x) ( x)]( l) v l q c q q= n (A.29) A.5. Nanodraht A.5.1. Bewegungsgleichung der Kohärenzenmatrix(Herleitung) Ziel dieser Rechnung ist es, die Bewegungsgleichung 5.53 herzuleiten. Ausgangspunkt ist die Heisenberg-Bewegungsgleichung: i O^=[O^, H^] t (A.30) Die zu untersuchende Observable ist gegeben durch O^=^a k ^a k und der Hamiltonoperator ist in Gl. 5.40 angegeben: H^= H^ kin + H^ eC + H^ eCi + H^ P In den Kommutator A.30 werden die Anteile des Hamiltonoperators nacheinander eingesetzt: 174 A.5. Nanodraht · Kinetischer Anteil: [^a k ^a k , H^ kin ]=[^a k ^a k , k ^a k ^a k ] k =( k k )^a k ^a k · Coulomb-Wechselwirkungsanteil( e e ): [^a k ^a k , H^ eC ]= 1 2 K,K ,q U q [^a k ^a k ,^a K ^a K ^a K + q ^a K q ] [^a 1 ^a 2 ,^a 3 ^a 4 ^a 5 ^a 6 ] [^a k ^a k , H^ eC ] = = 23 ^a 1 ^a U 4 ^a q 5 ( ^a 6 ^a k+ q 2 ^a 4 l ^a^a 1 l+ ^a 3 q ^a^a 5 k ^a 6 + 15 ^a 3 ^a 4 ^a 6 ^a 2 ) 16 - ^a k ^a l ^a l+ q ^a k q ^a 3 ^a 4 ^a 5 ^a 2 q,l · Coulomb-Wechselwirkungsanteil( e i ): Die(starren) Ionen, die über die Dichte n i ( q) gegeben sind, haben das Potential U q n i ( q) und können wie das Störpotential v P ( q, t) (s. unten) behandelt werden. · Störpotential-Anteil: [^a k ^a k , H^ P ]= v P ( q, t)[^a k ^a k ,^a k ^a k + q ] k ,q [^a k ^a k ,^a k ^a k + q ]=^a k ^a k + q k ,k ^a k ^a k k,k + q [^a k ^a k , H^ P ]= v P ( q, t) ( ^a k ^a k + q ) ^a k q ^a k q Die Bewegungsgleichung hat damit zunächst folgende Form: i t ^a k ^a k = ( k k ) ^a k ^a k U q ( ^a k+ q ^a l ^a l+ q ^a k - ^a k ) ^a l ^a l+ q ^a k q + q,l [ v P ( q, t) U q n i ( q)] ( ^a k ^a k + q - ^a k q ^a k ) q Die Hartree-Fock Faktorisierung( Singlet Beitrag) wird in der zweiten Zeile angewendet: ^a 1 ^a 2 ^a 3 ^a 4 S = ^a 1 ^a 4 ^a 2 ^a 3 - ^a 1 ^a 3 ^a 2 ^a 4 (A.31) Bewegungsgleichung in Hartree-Fock Näherung: i t ^a k ^a k ( k k ) ^a k ^a k U q ( ^a k+ q ^a k ^a l ^a l+ q - ^a k+ q ^a l+ q ^a l ^a k q,l ) ^a k ^a k q ^a l ^a l+ q + ^a k ^a l+ q ^a l ^a k q + [ v P ( q, t) U q n i ( q)] ( ^a k ^a k + q - ^a k q ^a k ) q 175 A. Anhang In der zweiten und dritten Zeile kann man die Ladungsdichte n e ( q) einsetzen: i t ^a k ^a k ( k k ) ^a k ^a k U q ( ^a k ^a l+ q ^a l ^a k q - ^a k+ q ^a l+ q ^a l ^a k ) + q,l L U q n e ( q) ( ^a k ^a k q - ^a k+ q ^a k ) + q [ v P ( q, t) U q n i ( q)] ( ^a k ^a k + q - ^a k q ^a k ) q Der Fock-/Austauschterm steht in der zweiten Zeile: k,k = U q ( ^a k ^a l+ q ^a l ^a k q - ^a k+ q ^a l+ q ^a l ^a k ) = q,l U q ( ^a k ^a l ^a l q ^a k q - ^a k+ q ^a l+ q ^a l ^a k ) q,l (A.32) Durch Änderung der Indizes von l+ q und l in l und l q wurde am gesamten Ausdruck nichts geändert, da der Indexabstand erhalten bleibt. Der Rechenschritt dient nur dazu, um auf den gleichen Term wie in Gl. 5.22 zu kommen, der in[56] hergeleitet wurde. Die Zeilen mit den Teilchendichten n e ( q) , n i ( q) und dem Störpotential können noch zusammengefasst werden. Dann erhält man die gesuchte Bewegungsgleichung 5.53: i t ^a k ^a k ( k k ) ^a k ^a k + k,k { U q L[ n e ( q) n i ( q)]+ v P ( q, t) } ( ^a k ^a k q - ^a k+ q ^a k ) (A.33) q A.5.2. Transformation in das Wignerbild(Herleitung) Ausgangspunkt ist die Bewegungsgleichung der Kohärenzenmatrix(5.53) in Hartree-Fock Näherung und die Denition der Wignerverteilung(5.66a). Die Beiträge zur Dynamik sind folgende: f( x, k)= f( x, k)+ f( x, k)+ f( x, k) t t kin t C t P (A.34) Der Coulomb-Beitrag kann in einen Hartree- und einen Fock-Beitrag zerlegt werden: f( x, k) = f( x, k) + f( x, k) t C t H t F Der Hartree-Beitrag enthält hier auch die Elektron-Ion-Wechselwirkung. (A.35) 176 A.5. Nanodraht Kinetik-Beitrag: f( x, k) t kin = t e iqx ^a k q/ 2 ^a k+ q/ 2 kin q = e iqx 1 i [ H^ kin , ] ^a k q/ 2 ^a k+ q/ 2 q = e iqx 1 i ( k+ q/ 2 k q/ 2 ) ^a k q/ 2 ^a k+ q/ 2 = q e iqx 1 ( 2 [ k+ q/ 2] 2 i q 2 m e 2 [ k q/ 2] 2 2 m e ) ^a k q/ 2 ^a k+ q/ 2 = 2 im e e iqx q 2 kq ^a k q/ 2 ^a k+ q/ 2 = k im e e iqx q q ^a k q/ 2 ^a k+ q/ 2 = k i 2 m e x q e iqx ^a k q/ 2 ^a k+ q/ 2 = k f( x, k) (exakt) m e x (A.36) Hartree-Beitrag( e e , e i und i i ): f( x, k) = t H t e iqx ^a k q/ 2 ^a k+ q/ 2 H q = e iqx 1 i q [ H^ eCe , ] ^a k q/ 2 ^a k+ q/ 2 Hartree = e iqx 1 i () U p L[ n( p) n i ( p)] ^a k q/ 2 ^a k+ q/ 2 p - ^a k q/ 2+ p ^a k+ q/ 2 qp An dieser Stelle wird die Rücktransformation nach f benutzt: ^a k q/ 2 ^a k+ q/ 2 p = 1 e i( q p) x f( x , k p/ 2) dx L ^a k q/ 2+ p ^a k+ q/ 2 = 1 L e i( q p) x f( x , k+ p/ 2) dx Damit folgt: f( x, k)= t H t e iqx ^a k q/ 2 ^a k+ q/ 2 H q = e iqx 1 i U p [ n e ( p) n i ( p)] q,p · e i( q p) x [ f( x , k p/ 2) f( x , k+ ] p/ 2) dx (A.37) An dieser Stelle wird eine Gradientenentwicklung für f bezüglich k eingesetzt: f( x, k+ p/ 2) f( x, k p/ 2) p f( x, k) k (A.38) 177 A. Anhang Damit geht nun eine Näherung in die Bewegungsgleichung ein: f( x, k) t H U p i [ n e ( p) n i ( p)] e iqx L/ 2 e i( q p) x p f( x , k) dx k q,p L/ 2 = i U p [ n e ( p) n i ( p)] L/ 2 e iq( x x ) e ipx p f( x , k) dx k p q L/ 2 An dieser Stelle kann nun folgende Rechenregel angewendet werden, die für eine beliebige Funktion f( x) gilt: L/ 2 e iq( x x ) f( x ) dx = L/ 2 e iqx f( x ) dx e iqx q L/ 2 q L/ 2 = L c q e iqx ( c q : Fourierkoezient ) q = Lf( x) D.h. über das Integral werden einfach nur die Fourierkoezienten der Fourier-Reihe von f( x, k) bzgl. x berechnet. Die Summe darüber ergibt dann f( x, k) . Daher folgt: f( x, k) t H i U p L[ n e ( p) n i ( p)] e ipx p k f( x, k) p (A.39) Ziel ist es nun das elektrische Feld x ( x) einzubauen. Dazu wird zunächst p e ipx = i x e ipx verwendet: () 1 x L U p [ n e ( p) n i ( p)] e ipx f( x, k) k p L/ 2 = 1 x V( | x x | )[ n e ( x ) n i ( x )] dx k f( x, k) L/ 2 = q e ( x) · f( x, k) x k (A.40) In der ersten Zeile erkennt man die Fouriertransformierte der Faltung aus der zweiten Zeile, welche die Lösung der Poissongleichung für das Potential ( x) darstellt. Der Störterm, der durch das Potential v p ( x) (bzw. v p ( p) im Fourierraum) bedingt ist, kann analog zu dieser Rechnung, beginnend bei Gl. A.37 mit der Ersetzung U p [ n e ( p) n i ( p)] v p ( p) erfolgen. Um das Ergebnis für den Störterm zu erhalten, muss in Gl. A.40 deshalb nur die Ersetzung q e ( x) v p ( x) gemacht werden: f( x, k) t P 1 x v p · k f( x, k) (A.41) 178 A.5. Nanodraht Fock-Beitrag: f( x, k)= t F t e iqx ^a k q/ 2 ^a k+ q/ 2 F q = e iqx 1 i q [ H^ eCe , ] ^a k q/ 2 ^a k+ q/ 2 Fock = e iqx 1 i k q/ 2 ,k+ q/ 2 q ( k,k : s. Gl. A.32 ) = 1 i e iqx U p ( ^a k q/ 2+ p ^a l+ p ^a l ^a k+ q/ 2 q p,l ) ^a k q/ 2 ^a l ^a l p ^a k+ q/ 2 p = 1 i e iqx U p ( ^a k q/ 2+ p ^a l+ p ^a l ^a k+ q/ 2 q p,l ) ^a k q/ 2 ^a l+ p ^a l ^a k+ q/ 2 p Hier muss die Rücktransformation eingesetzt werden: ^a k ^a k = 1 L e i( k k) x f( x,[ k+ k ]/ 2) dx (A.42) Die Energierenormierung wird analog zur dreidimensionalen Form(Gl. 5.32) deniert: k ( x)= 1 | q e | U k k f( x, k ) k (A.43) Nach einer Gradientenentwicklung von f und erhält man: f( x, k) t F | q e | ( [ k k ( x)] · [ x f( x, k)] [ x k ( x)] · [ k f( x, ) k)] (A.44) Dieses Ergebnis konnte bis jetzt noch nicht detailliert nachgerechnet werden und erfordert eine Verikation. 179 A. Anhang A.6. Hartree-Energiefunktional In(zeitabhängigen) DFT-Simulationen ist eine genaue Aufschlüsselung der Beiträge zur Gesamtenergie wünschenswert. Die Beiträge sollten durch möglichst wenig redundante Rechenschritte erhalten werden können. Die Berechnung der Hartree-Energie ist ein gutes Beispiel dafür. Die folgende Rechnung zeigt in diesem Zusammenhang, wie sich die elektrostatische Energie W el zusammensetzt: W el = 1 2 d 3 r d 3 r ( r ) ( r) | r r | = 1 2 d 3 r( r) el ( r) (A.45a) = = q e 2 q e 2 d 3 r[ n ( r) n + ( r)] el ( r) [ d 3 r[ n ( r) n + ( r)] d 3 r q e n ( r ) | r r | d 3 r q e n + ( r ) | r r | ] = q e2 [ 2 d 3 r d 3 r n ( r) n ( r ) | r r | 2 d 3 r d 3 r n ( r) n + ( r ) | r r | + d 3 r d 3 r n + ( r) n + ( r ) | r r | ] = E H [ n ( r)]+ q e2 2 [ 2 d 3 r d 3 r n ( r) n + ( r ) | r r | + d 3 r d 3 r n + ( r) n + ( r ) | r r | ] = E H [ n ( r)]+ E ext [ n ( r)]+ W el, i i (A.45b) In der Simulation stehen die Felder und el unmittelbar zur Verfügung, so dass W el über die Integration in Gl. A.45a direkt berechnet werden kann. Der Beitrag E ext [ n ( r)] ist auf gleiche elementare Weise berechenbar. Der Beitrag W el, i i kann einmalig vorberechnet werden. Es müssen also nur zwei einfache Raumintegrale(einmal für W el und einmal für E ext [ n ] ) durchgeführt werden um an alle relevanten Gröÿen zu gelangen. Die HartreeEnergie ergibt sich entsprechend Gl. A.45b als Dierenz: E H [ n ( r)]= W el W el, i i E ext [ n ( r)] (A.46) Für die Berechnung der Hartree-Energie muss daher kein Doppelintegral durchgeführt werden. 2 Die Rechnung ist in atomaren Einheiten angegeben: 1/(4 0 )= 1 , q e = 1 180 Literaturverzeichnis [1] Matthias W. Klein, Christian Enkrich, Martin Wegener und Stefan Linden. Second Science harmonic generation from magnetic metamaterials., 313:502504, 2006. nature photonics [2] Lukas Novotny und Niek van Hulst. Antennas for light., 5:8390, 2011. [3] Alberto G. Curto, Giorgio Volpe, Tim H. Taminiau, Mark P. Kreuzer, Romain Quidant und Niek F. van Hulst. Unidirectional emission of a quantum dot coupled Science to a nanoantenna., 329(5994):930933, 2010. [4] Gunnar Dolling, Martin Wegener und Stefan Linden. Der falsche Knick im Licht. Phys. Unserer Zeit , 38:2429, 2007. [5] Jerey M. McMahon, Stephen K. Gray und George C. Schatz. Nonlocal optical rePhys. Rev. Lett. sponse of metal nanostructures with arbitrary shape., 103:097403, Aug 2009. [6] Emil Prodan Jorge Zuloaga und Peter Nordlander. Quantum plasmonics: Optical ACS Nano properties and tunability of metallic nanorods., 4(9):52695276, 2010. [7] D.C. Marinica, A.K. Kazansky, P. Nordlander, J. Aizpurua und A. G. Borisov. Quantum plasmonics: Nonlinear eects in the eld enhancement of a plasmonic Nano Letters nanoparticle dimer., 12(3):13331339, 2012. [8] Robert Boyd. Nonlinear Optics . Academic Press, 2008. Laser& [9] P. Vasa, C. Ropers, R. Pomraenke und C. Lienau. Ultra-fast nano-optics. Photonics Reviews , 3(6):483507, 2009. ISSN 1863-8899. Physik Journal [10] Markus Drescher und Ferenc Krausz. Schnappschüsse im Atom., 2:4550, 2003. [11] S. Janz und H. M. van Driel. Second-harmonic generation from metal surfaces. International Journal of Nonlinear Optical Physics , 2(1):142, 1993. [12] S. Linden, F. B. P. Niesler, J. Förstner, Y. Grynko, T. Meier und M. Wegener. Collective eects in second-harmonic generation from split-ring-resonator arrays. Phys. Rev. Lett. , 109:015502, Jul 2012. [13] Joseph Rudnick und E. A. Stern. Second-harmonic radiation from metal surfaces. Phys. Rev. B , 4:42744290, Dec 1971. 181 LITERATURVERZEICHNIS [14] N. Bloembergen, R. K. Chang, S. S. Jha und C. H. Lee. Optical second-harmonic Phys. Rev. generation in reection from media with inversion symmetry., 174: 813822, Oct 1968. [15] Xiang Wang, Francisco J. Rodriguez, Willem M. Albers, Risto Ahorinta, J. E. Sipe und Martti Kauranen. Surface and bulk contributions to the second-order nonlinear Phys. Rev. B optical response of a gold lm., 80:233402, Dec 2009. [16] Charles Kittel. Introduction to Solid State Physics . John Wiley& Sons, Inc., 8th edition, 2005. [17] Gabriele F. Giuliani und Giovanni Vignale. Quantum Theoy of the Electron Liquid . Cambridge University Press, 2008. Annalen der Physik [18] Paul Drude. Zur Elektronentheorie der Metalle., 306(3):566, 1900. Am. J. Phys. [19] Neil Ashby. Relaxation of charge imbalances in conductors., 43(6): 553555, 1975. Physical [20] P. B. Johnson und R. W. Christy. Optical constants of the noble metals. Review B , 6:43704379, 1972. [21] David J. Griths. Introduction to Electrodynamics . Prentice Hall, 3rd edition, 1999. [22] Fernando Haas. An introduction to quantum plasmas. Brazilian Journal of Physics , 41(4-6):349363, 2011. ISSN 0103-9733. [23] G. Manfredi und F. Haas. Self-consistent uid model for a quantum electron gas. Phys. Rev. B , 64:075316, Jul 2001. [24] Stefan A. Maier. Plasmonics: Fundamentals and applications . Springer, 2007. [25] Hans C. Ohanian. On the approach to electro- and magneto-static equilibrium. Am. J. Phys. , 51(11):10201022, 1983. [26] G. S. Agarwal, D. N. Pattanayak und E. Wolf. Electromagnetic elds in spatially Phys. Rev. B dispersive media., 10:14471475, Aug 1974. [27] John David Jackson. Classical Electrodynamics . John Wiley& Sons, 3rd edition, 1999. [28] Eduardo A. Coronado und George C. Schatz. Surface plasmon broadening for The Journal of arbitrary shape nanoparticles: A geometrical probability approach. Chemical Physics , 119(7):39263934, 2003. Phys. Rev. B [29] P. Hohenberg und W. Kohn. Inhomogeneous electron gas., 136:864, 1964. [30] W. Kohn und L. J. Sham. Self-consistent equations including exchange and correPhys. Rev. lation eects., 140:A1133, 1965. 182 LITERATURVERZEICHNIS Phys. Rev. [31] E. Wigner. On the interaction of electrons in metals., 46:10021011, Dec 1934. [32] Ansgar Liebsch. Electronic Excitations at Metal Surfaces . Springer, 1st edition, 1997. [33] Erich Runge und E. K. U. Gross. Density-functional theory for time-dependent systems. Phys. Rev. Lett. , 52:9971000, Mar 1984. [34] S. Kurth, G. Stefanucci, C.-O. Almbladh, A. Rubio und E. K. U. Gross. Timedependent quantum transport: A practical scheme using density functional theory. Phys. Rev. B , 72:035308, Jul 2005. [35] Adolfo G. Eguiluz. Self-consistent static-density-response function of a metal surPhys. Rev. B face in density-functional theory., 31:33033314, Mar 1985. [36] N. D. Lang und W. Kohn. Theory of metal surfaces: Charge density and surface energy. Physical Review B , 1:45554568, 12 1970. [37] P. M. Oppeneer T. Kampfrath, M. Battiato und M. Munzenberg. Terahertz spin Nat Nano current pulses controlled by magnetic heterostructures., 8(4):256260, 2013. [38] M. Anderegg, B. Feuerbacher und B. Fitton. Optically excited longitudinal plasmons in potassium. Phys. Rev. Lett. , 27:15651568, Dec 1971. [39] C.L. Phillips, J. Parr und E. Riskin. Signals, Systems, and Transforms . Pearson Education, 2011. ISBN 9780133002300. [40] A.D. Boardman. Electromagnetic Surface Modes . John Wiley& Sons Australia, Limited, 1982. ISBN 9780471100775. [41] Allen Taove und Susan C. Hagness. Computational Electrodynamics . Artech House, 3rd edition, 2005. [42] Daniel Neuhauser und Kenneth Lopata. Quantum Drude friction for timedependent density functional theory. The Journal of Chemical Physics , 129(13): 134106, 2008. [43] T. Weiland. A discretization method for the solution of Maxwell's equations for sixcomponent elds. Electronics and Communications AEUE , 31(3):116120, 1977. [44] Lingyi Meng, ChiYung Yam, SiuKong Koo, Quan Chen, Ngai Wong und GuanHua Chen. Dynamic multiscale quantum mechanics/electromagnetics simulation method. Journal of Chemical Theory and Computation , 8(4):11901199, 2012. [45] K. S. Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. IEEE Trans. Antennas Propagat. , 14:302307, 1966. 183 LITERATURVERZEICHNIS [46] William R. Frensley. Boundary conditions for open quantum systems driven far Rev. Mod. Phys. from equilibrium., 62:745791, Jul 1990. [47] Mathias Wand. Theoretische Modellierung der Erzeugung von höheren Harmonischen in plasmonischen Strukturen. Master's thesis, Universität Paderborn, 2009. [48] Mary B. James und David J. Griths. Why the speed of light is reduced in a Am. J. Phys. transparent medium., 60(4):309313, April 1992. [49] M. D. Kostin. On the Schrödinger-Langevin equation. The Journal of Chemical Physics , 57(9):35893591, 1972. Phy[50] A. Davidson. Damping in Schrödinger's equation for macroscopic variables. sical Review A , 41(6):33953398, 1990. [51] Kieron Burke, Roberto Car und Ralph Gebauer. Density functional theory of the Phys. Rev. Lett. electrical conductivity of molecular devices., 94:146803, Apr 2005. [52] Joel Yuen-Zhou, David G. Tempel, Cesar A. Rodriguez-Rosario und Alan AspuruGuzik. Time-dependent density functional theory for open quantum systems with Phys. Rev. Lett. unitary propagation., 104:043001, Jan 2010. [53] Xiao Zheng, ChiYung Yam, Fan Wang und GuanHua Chen. Existence of timedependent density-functional theory for open electronic systems: Time-dependent holographic electron density theorem. Phys. Chem. Chem. Phys. , 13:1435814364, 2011. [54] David G. Tempel und Alan Aspuru-Guzik. Relaxation and dephasing in open quantum systems time-dependent density functional theory: Properties of exact Chemical Physics functionals from an exactly-solvable model system., 391(1):130 142, 2011. ISSN 0301-0104. [55] Arno Schindlmayr. Private Kommunikation, 2009-2013. [56] W. Hoyer und S. W. Koch. Theory of the nonlinear optical response of metallic and plasma systems. Unveröentlicht, 2007. [57] Peter Kölling. Modellrechnungen zur nichtlinearen optischen Antwort von Metallen auf kurzen Zeitskalen. Master's thesis, Universität Paderborn, 2012. [58] R. Balescu. Statistical Mechanics Of Charged Particles . John Wiley& Sons, 1963. [59] Nam-Duk Suh, Marl R Feix und Pierre Bertrand. Numerical simulation of the quantum Liouville-Poisson system. Journal of Computational Physics , 94(2):403 418, 1991. Phys. [60] E. Wigner. On the quantum correction for thermodynamic equilibrium. Rev. , 40:749759, 1932. [61] M. Abramowitz und I.A. Stegun. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables . Applied mathematics series. Dover Publications, 1964. ISBN 9780486158242. 184 LITERATURVERZEICHNIS [62] Radu Balescu. Equilibrium and Non-Equilibrium Statistical Mechanics . Wiley, 1975. [63] P.A. Markowich. On the equivalence of the Schrödinger and the quantum Liouville equations. Math. Meth. Appl. Sci. , 11:459469, 1989. [64] D. Bohm und B.J. Hiley. The Undivided Universe: An Ontological Interpretation of Quantum Theory . Physics, philosophy. Routledge, 1993. ISBN 9780415065887. [65] N. Crouseilles, P.-A. Hervieux und G. Manfredi. Quantum hydrodynamic model Phys. Rev. B for the nonlinear electron dynamics in thin metal lms., 78:155412, Oct 2008. Numerical [66] W. H. Press, S. A. Teukolsky, W. T. Vetterling und B. P. Flannery. Recipes In C: The Art Of Scientic Computing . Cambridge University Press, 1992. [67] L. Lehtovaara, J. Toivanen und J. Eloranta. Solution of time-independent Schrödinger equation by the imaginary time propagation method. Journal of Computational Physics , 221(1):148 157, 2007. ISSN 0021-9991. [68] C. Yang R. B. Lehoucq, D. C. Sorensen. ARPACK Users' Guide: Solution of LargeScale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods . SIAM, 1998. [69] Yousef Saad. Iterative Methods for Sparse Linear Systems . SIAM, 2 edition, 2003. [70] M. Aichinger und E. Krotscheck. A fast conguration space method for solving local Kohn-Sham equations. Computational Materials Science , 34(2):188 212, 2005. ISSN 0927-0256. [71] C. G. Broyden. A class of methods for solving nonlinear simultaneous equations. Math. Comp. , 19:577593, 1965. [72] Peter Pulay. Convergence acceleration of iterative sequences. The case of SCF iteration. Chemical Physics Letters , 73(2):393398, 1980. [73] Stefan Bluegel. Density functional theory in practice. Technical Report A7, Institut für Festkörperforschung, Forschungszentrum Jülich GmbH, 2007. [74] George F. Bertsch M.A.L. Marques, Alberto Castro und Angel Rubio. octopus: a rst-principles tool for excited electron-ion dynamics. Comput. Phys. Commun. , 151:6078, 2003. [75] R. Nieminen M. Manninen und P. Hautojaervi. Electrons and positrons in metal Phys. Rev. B vacancies., 12(10):40124022, 1975. [76] M. Wand, A. Schindlmayr, T. Meier und J. Foerstner. Simulation of the ultrafast Phys. Status Solidi B nonlinear optical response of metal slabs., 248:887891, 2011. [77] Alberto Castro, Miguel A. L. Marques und Angel Rubio. Propagators for the time-dependent KohnSham equations. The Journal of Chemical Physics , 121(8): 34253433, 2004. 185 LITERATURVERZEICHNIS [78] Michael A. Celia und William G. Gray. Numerical Methods for Dierential Equations . Prentice Hall, 1992. Fundamental Concepts for Scientic and Engineering Applications. [79] Wolfgang Nolting. Grundkurs Theoretische Physik 5/1: Quantenmechanik- Grundlagen . Springer Verlag, 2004. [80] Wilhelm Magnus. On the exponential solution of dierential equations for a linear operator. Communications on Pure and Applied Mathematics , 7(4):649673, 1954. ISSN 1097-0312. [81] S. Blanes, F. Casas, J.A. Oteo und J. Ros. The magnus expansion and some of its Physics Reports applications., 470(56):151 238, 2009. ISSN 0370-1573. [82] R. Ward. Numerical computation of the matrix exponential with accuracy estimate. SIAM Journal on Numerical Analysis , 14(4):600610, 1977. [83] C. Moler und C. Van Loan. Nineteen dubious ways to compute the exponential of SIAM Review a matrix, twenty-ve years later., 45(1):349, 2003. [84] William George Horner. A new method of solving numerical equations of all orders, by continuous approximation. Philosophical Transactions of the Royal Society of London , 1:308335, 1819. Math. Comp. [85] C. W. Clenshaw. A note on the summation of Chebyshev series., 9: 118120, 1955. [86] R. Koslo H. Tal-Ezer. An accurate and ecient scheme for propagating the time J. Chem. Phys. dependent Schrödinger equation., 81(9):39673971, 1984. [87] J. C. Light Tae Jun Park. Unitary quantum time evolution by iterative Lanczos J. Chem. Phys. reduction., 85(10):58705876, 1986. [88] M.D Feit, J.A Fleck Jr. und A Steiger. Solution of the Schrödinger equation by a spectral method. Journal of Computational Physics , 47(3):412 433, 1982. ISSN 0021-9991. [89] Osamu Sugino und Yoshiyuki Miyamoto. Density-functional approach to electron Phys. Rev. B dynamics: Stable simulation under a self-consistent eld., 59(4): 25792586, 1999. [90] Masuo Suzuki und Takashi Yamauchi. Convergence of unitary and complex decompositions of exponential operators. Journal of Mathematical Physics , 34(10): 48924897, 1993. [91] H. Akima. A new method of interpolation and smooth curve tting based on local procedures. J. Assoc. Comput. Math. , 17:589602, 1970. [92] Gerhard Wanner Ernst Hairer, Christian Lubich. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Dierential Equations . Springer, 2006. 186 LITERATURVERZEICHNIS [93] E.W. Cheney und D. Kincaid. Numerical Mathematics and Computing . International student edition. Brooks/Cole, 2008. ISBN 9780495114758. [94] Randall J. Leveque. Finite Volume Methods for Hyperbolic Problems . Cambridge University Press, 2002. [95] J.H. Ferziger und M. Peri¢. Computational Methods for Fluid Dynamics . Springer London, Limited, 2002. ISBN 9783540420743. [96] Alexander Kurganov und Eitan Tadmor. New high-resolution central schemes for Journal of Comnonlinear conservation laws and convectiondiusion equations. putational Physics , 160:241282, 2000. [97] Luigi Genovese, Thierry Deutsch, Alexey Neelov, Stefan Goedecker und Gregory Beylkin. Ecient solution of poisson's equation with free boundary conditions. The Journal of Chemical Physics , 125(7):074105, 2006. [98] Bronstein, Semendjajew, Musiol und Mühlig. Taschenbuch der Mathematik . Harri Deutsch, 5th edition, 2001. American Journal [99] O. L. Brill und B. Goodman. Causality in the Coulomb gauge. of Physics , 35(9):832837, 1967. [100] R. W. Alexander M. A. Ordal, R. J. Bell und M. R. Querry. Optical properties of fourteen metals in the infrared and far infrared: Al, Co, Cu, Au, Fe, Pb, Mo, Ni, Appl. Opt. Pd, Pt, Ag, Ti, V, and W., 24:44934499, 1985. 187