IN LASER-DRIVEN DIATOMIC MOLECULES
DISSERTATION
zur Erlangung des akademischen Grades eines Doktors der Naturwissenschaften (Dr. rer. nat.)
im Fachbereich Naturwissenschaften der Universität Kassel
vorgelegt von
CHIRAG ASHOKKUMAR JHALA
Tag der Disputation:
15. Dezember 2011
Erstgutachter : Prof. Dr. Manfred Lein Zweitgutachter: Prof. Dr. Christiane Koch
Electron-Nuclear Correlation in Laser-Driven Diatomic Molecules c
2011 Chirag Ashokkumar Jhala All rights reserved
Printed in Hannover, 2011
Institut fu¨r Physik Universität Kassel Heinrich-Plett-Srtaße 40 D-34109 Kassel
Germany
This thesis was typeset with LATEX-system.
Zusammenfassung v
1 Introduction 1
2 Multicomponent density functional theory 7
2.1 Multicomponent density functional theory: An overview . . . 8
2.1.1 Coordinate transformation . . . 8
2.1.2 Defining the densities . . . 10
2.1.3 The Runge-Gross theorem for multicomponent system . . . . 10
2.1.4 The time-dependent Kohn-Sham scheme for multicomponent system . . . 12
2.2 The model system . . . 15
2.2.1 One-dimensional H+ 2 molecular ion . . . 16
2.3 Construction of the exact Kohn-Sham potential . . . 16
2.4 Dynamics on Born-Oppenheimer potential energy surface . . . 19
2.5 Dynamics in presence of an external field . . . 22
2.6 Discussion . . . 24
3 Multiconfiguration time-dependent hartree approach 27 3.1 The theory of MCTDH . . . 28
3.2 The model system . . . 29
3.3 Ground state calculation . . . 30
3.4 The H+ 2 molecule in an external field . . . 33
3.4.1 Dynamics in the presence of an external field with 800 nm wavelength . . . 33
3.4.2 Dynamics in the presence of an external field with 150 nm wavelength . . . 42
3.4.3 Discussion . . . 46
3.5 The H2 molecule in an external field . . . 48
3.5.1 Two-center interference in high-order harmonic generation from H2 molecules . . . 49
3.5.2 Lochfraß in H2 molecules . . . 52
3.5.3 Discussion . . . 59 iii
4 Effects of absorbing boundaries in the mean-field approximation 61
4.1 Derivation of equations of motion . . . 62
4.1.1 Hartree variational equations of motion . . . 63
4.1.2 Hartree TDDFT equations of motion . . . 64
4.2 Discussion of the equations . . . 64
4.3 Results . . . 65
4.4 Discussion . . . 67
5 Summary 69
Appendices 73
A Born-Oppenheimer approximation 75
B Derrivation of MCTDH equations of motion 77
Bibliography 81
Publications 89
Acknowledgments 91
Erklärung 93
Seit den Anfängen der Lasertechnologie im Jahr 1958 hat es der rasche Fortschritt in diesem Feld inzwischen ermöglicht, ultrakurze Laserpulse mit hoher Intensität zu erzeugen. Diese ultrakurzen Laserfelder können so stark sein, dass sie mit den Coulomb-Kräften, die die Elektronen mit dem Atomkern zusammenhalten, um die vorherrschende Rolle in der Dynamik des Systems konkurrieren. Wechselwirkung von intensiven Feldern mit Molekülen verursacht höchst nichtlineare Prozesse, die mit einem störungstheoretischen Ansatz nicht erklärt werden können. Stattdessen ist ein über die Störungstheorie hinaus gehender Ansatz zur Beschreibung der Mole-küldynamik nötig, der die Dynamik der Elektronen und Kerne vollständig quanten-mechanisch beschreibt. Die nichtlinearen Prozesse, die in Laser-Materie-Wechsel-wirkungen beobachtet werden, können mit Hilfe der zeitabhängigen Schrödinger-gleichung (TDSE) beschrieben werden. Die TDSE kann allerdings nur für Systeme mit sehr wenigen Freiheitsgraden exakt gelöst werden, da mit einem Anstieg der Zahl an Freiheitsgraden auch der numerische Aufwand, was die Anforderungen an Speicher und Rechenzeit betrifft, exponentiell anwächst. Dadurch wird die Aufgabe in der Praxis schnell unlösbar.
Um das Problem des exponentiellen Anstiegs zu umgehen, wurden in der Vergangen-heit bereits verschiedene Näherungsmethoden vorgeschlagen. Die Herangehensweise der Born-Oppenheimer(BO)-Näherung besteht darin, die Elektronenbewegung von der Kernbewegung zu entkoppeln, indem die stationäre Schrödingergleichung für die Elektronen parametrisch für jede Kernkonfiguration gelöst wird, während die Bewe-gung der Kerne zeitabhängig auf gekoppelten Potentialfächen stattfindet. Allerdings beschreibt die BO-Näherung Systeme, welche mit einem starken Feld wechselwirken, oftmals nicht akkurat.
In der Dichtefunktionaltheorie (DFT) ist die Dichte des Systems die wesentliche Größe. Alle Observablen sind hier Funktionale der Dichte. In den meisten der frühe-ren Rechnungen, die statische oder zeitabhängige DFT verwendet haben, wurden die Kerne entweder als statisch angenommen oder klassich behandelt. Wenn ein System starken Feldern ausgesetzt wird, sind die elektronischen Freiheitsgrade oft stark mit der Kernbewegung gekoppelt, und daher ist eine möglichst einheitliche Beschreibung beider Bewegungen wünschenswert. Die Multikomponenten-Dichtefunktionaltheorie
vi
(MCDFT) kann Elektronen und Kerne gleichzeitig und vollständig quantenmecha-nisch beschreiben.
In Kapitel 2 wurde ein Überblick über die zeitabhängige MCDFT (TD-MCDFT) gegeben. In der TD-MCDFT sind die wesentlichen Variablen die Elektronendich-te sowie die KerndichElektronendich-te. Zur Berechnung dieser DichElektronendich-ten wird ein SysElektronendich-tem nicht-wechselwirkender Elektronen herangezogen sowie ein System gekoppelter Kerne, die nicht mit Elektronen wechselwirken. Stattdessen treten in diesen sogenannten Kohn-Sham(KS)-Systemen effektive Potentiale auf: die KS-Potentiale. Die elektronischen KS-Gleichungen werden nicht parametrisch für feste Kernkonfigurationen gelöst. Die Information über die Kernkonfiguration geht vielmehr über die KS-Potentiale ein, da diese Funktionale der Kerndichte sind. Ebenso fließt jede Änderung der elek-tronischen Dichte in die KS-Potentiale ein, da diese auch Funktionale der Elektro-nendichte sind. Effekte, die eine Beschreibung jenseits der BO-Näherung erfordern, sind in der TD-MCDFT enthalten, da in der Herleitung keinerlei Näherung verwen-det wird. In der Praxis sind die zeitabhängigen Austauschkorrelationspotentiale, die ein Bestandteil der KS-Potentiale sind und die relevanten Vielkörpereffekte enthal-ten, nicht leicht zu berechnen. Diese müssen daher genähert werden. Der Erfolg der TD-MCDFT steht und fällt also mit der verwendeten Näherung.
TD-MCDFT wurde verwendet, um die gekoppelte Elektron-Kern-Dynamik beispiel-haft in einem eindimensionalen Wasserstoffmolekül-Ion zu beschreiben. In Abschnitt 2.3 wurde ein Algorithmus hergeleitet und diskutiert, um exakte KS-Potentiale für dieses System zu berechnen. Die komplexen KS-Orbitale wurden zunächst aus der numerisch berechneten quasi-exakten Dichte konstruiert. Durch Umkehrung des Split-Operator-Propagators wurden numerisch exakte KS-Potentiale berechnet. So-bald die genauen KS Potentiale bekannt sind, kann das Austauschkorrelationspo-tential einfach durch Subtraktion des HartreepoAustauschkorrelationspo-tentials und des externen PoAustauschkorrelationspo-tentials berechnet werden. Diese Methode ist nur für kleine Modellsysteme durchführbar, aber sie gibt einen Einblick in die Eigenschaften des Austauschkorrelationspotenti-als.
In den Abschnitten 2.4 und 2.5 wurde der Photodissoziationsprozess im Wasserstoff-molekül-Ion genutzt, um die KS-Potentiale in TD-MCDFT zu untersuchen. Zu-nächst wurde ein geeignet präparierter Zustand, der ohne externes Feld dissoziiert, betrachtet. Anschließend wurde das System in der Anwesenheit eines starken Feldes untersucht, welches Fragmentation verursacht. Es wurde festgestellt, dass die in die-sen Situationen resultierenden nuklearen KS-Potentiale eine stückweise Kombination von Grundzustandspotentialkurven (für kleine Kernabstände) und Potentialkurven angeregter Zustände (für große Kernabstände) sind.
Die Hartree-Näherung vernachlässigt die Korrelationsbeiträge. Das entsprechende „Hartree-only“-KS-Potential hat einen zu stark bindenden Charakter für große Kern-abstände. Dieser Mangel wirkt sich sowohl auf die statischen als auch auf die zeit-abhängigen KS-Berechnungen aus. Daher versagt die Hartree-Näherung sogar bei der qualitativen Beschreibung des Photodissoziationsprozesses. Wir schlussfolgern,
dass Korrelationseffekte eine wichtige Rolle bei der Beschreibung der Dynamik des Systems spielen.
Die Schwierigkeiten, die mit dem TD-MCDFT-Ansatz verbunden sind, können durch die Verwendung eines explizit korrelierten Ansatzes für die Vielteilchenwellenfunkti-on umgangen werden. In Kapitel 3 haben wir die zeitabhängige MultikVielteilchenwellenfunkti-onfiguratiVielteilchenwellenfunkti-ons- Multikonfigurations-Hartree(MCTDH)-Methode für ein System von gekoppelten Elektronen und Kernen eingeführt. In MCTDH-Methode wird die gesamte Wellenfunktion als Summe von Produkten aus Einteilchenfunktionen (SPFs) geschrieben, die selbst zeitabhängig sind. Eine kurze Einführung in die MCTDH-Methode wurde in Abschnitt 3.1 gege-ben und die Bewegungsgleichungen für das System wurden diskutiert. Wir hagege-ben zwei eindimensionale Modellsysteme beschrieben, nämlich das Wasserstoffmolekül-Ion mit einem elektronischen Freiheitsgrad und das Wasserstoffmolekül mit zwei Elektronen. Grundzustandseigenschaften wie Grundzustandsenergien, Kern- und Elektronendichten wurden unter Verwendung der Hartree-Methode (äquivalent zu Ein-SPF-MCTDH) und unter Verwendung der MCTDH-Methode mit zwei bis acht SPFs berechnet und in Abschnitt 3.2 mit exakten Rechnungen verglichen. Es wurde gezeigt, dass zwei bis vier SPFs genügen, um die Grundzustandseigenschaften des Wasserstoffmolekül-Ions zu beschreiben, während vier bis acht SPFs für das Wasser-stoffmolekül mit einem zusätzlichen Elektron benötigt wurden. Dies deutet darauf hin, dass sich die erforderliche Anzahl von SPFs mit der Anzahl der Freiheitsgrade erhöht.
In Abschnitt 3.3 wurde die Wechselwirkung des eindimensionalen Ions mit Laserpulsen verschiedener Wellenlängen untersucht. Das Wasserstoffmolekül-Ion wurde schon in der Vergangenheit ausgiebig erforscht. Deshalb war unser pri-märes Ziel nicht, die Physik dieses Modells zu studieren, sondern die Anwendbar-keit der MCTDH-Methode auf ein Modellsystem mit gekoppelter elektronischer und nuklearer Bewegung. Das Modellsystem, welches hier betrachtet wurde, ist recht einfach, aber sehr nützlich, da es eine vielfältige Dynamik aufweist und viele der wichtigsten Merkmale zeigt, die charakteristisch für das Verhalten von zweiatomi-gen Molekülen in starken Feldern sind. Das Scheitern der Hartree-Methode spiegelt sich in der zu schmalen Form der Grundzustandskerndichte wider. Das nukleare Po-tential, das von diesem Zustand erzeugt wird, ist stark bindend, und in Anwesenheit von moderat starken Feldern ist das Molekül nicht in der Lage, genügend Energie zu absorbieren, um diese bindende Kraft zu überwinden. Die Ergebnisse der MCTDH-Methode dagegen konvergieren mit zunehmender Anzahl der SPFs gegen die exakten Ergebnisse. Es wurde auch gezeigt, dass – abhängig von der Observablen – die An-zahl der Konfigurationen, die für numerisch konvergierte Ergebnisse erforderlich ist, sich mit der Spitzenintensität des Laserpulses erhöht, d.h. die Korrelation zwischen elektronischer und nuklearer Bewegung wird durch das Feld erhöht. Dennoch war die Anzahl der benötigten SPFs moderat.
Intensive Laserfelder können Ionisation oder Dissoziation eines Systems hervorru-fen. In einer exakten Rechnung wird die Wellenfunktion in eine vollständige Basis
viii
entwickelt, die eine genaue Beschreibung aller Details der sich entwickelnden Wel-lenfunktion ermöglicht. Im Falle von MCTDH kann die zugewiesene Anzahl der SPFs nicht eine vollständige Basis bilden, aber die Rechnung ergibt die bestmög-liche Beschreibung der Wellenfunktion für die fest vorgegebene Anzahl von SPFs. Mit zunehmender Anzahl der SPFs wird die Basis vollständiger und die Entwick-lung immer genauer. Die Verwendung der MCTDH-Methode bietet den Vorteil, dass die Genauigkeit der Berechnung systematisch verbessert werden kann, während die numerischen Anforderungen überschaubar im Vergleich zu exakten Rechnungen blei-ben.
Bei Betrachtung von vielteilchensystemen wie dem Wasserstoffmolekül, werden die exakten Berechnungen sehr zeitaufwändig und die Gittergröße ist begrenzt, so dass alternative Methoden benötigt werden. Abgesehen davon ist es von Interesse, die zusätzlichen Elektron-Elektron-Korrelationseffekte im Wasserstoffmolekül im Rah-men der MCTDH-Methode zu studieren. Diese Methode ermöglicht es, Dynamik auf großen Gittern zur Untersuchung von Ionisations- und Dissoziationsprozessen zu behandeln. Als Anwendungen haben wir in Abschnitt 3.4 den Zwei-Zentren-Interferenzeffekt bei der Erzeugung hoher Harmonischer (HHG) und den Loch-fraßeffekt bei der Ionisation des Wasserstoffmoleküls studiert. In Abschnitt 3.4.1 wurde zunächst die Zwei-Zentren-Interferenz mit dem Hartree- und dem MCTDH-Ansatz untersucht. Dies wurde mit den exakten Rechnungen für das eindimensionale Wasserstoffmolekül verglichen. Ein Minimum im HHG-Spektrum wurde für alle Re-chenmethoden festgestellt. Die HHG-Spektren aus den Hartree- und den Zwei-SPF-MCTDH-Rechnungen reproduzierten die Peaks an den richtigen Positionen, aber die Amplituden waren nicht in Übereinstimmung mit der exakten Rechnung. Die Zeit-Frequenz-Analyse der Emission der harmonischen Strahlung wurde unter Ver-wendung der Gabor-Transformation ausgewertet. Die Gabor-Analyse der exakten Rechnung zeigte, dass die Harmonischen in einem Bereich um die harmonische Ord-nung 30 unterdrückt wurden. Das Minmum wurde besonders klar durch Integration der Gabordichte über alle Emissionszeiten herausgestellt. Der Hartree-Ansatz zeigte die Unterdrückung der Harmonischen bei der 35. harmonischen Ordnung. Eine sol-che, recht gute qualitative Übereinstimmung war für die Hartree-Methode nicht un-bedingt erwartet worden, da diese keine sinnvollen Ergebnisse bei der Untersuchung des Wasserstoffmolekül-Ions geliefert hatte. Die Zwei-SPF-MCTDH-Rechnung sowie die MCTDH-Rechnung mit einer großen Anzahl von SPFs ergaben ein Minimum an der richtigen Position.
In Abschnitt 3.4.2 haben wir den Lochfraßeffekt in Wasserstoff- und Deuteriummole-külen für unterschiedliche Pumpwellenlängen vom ultravioletten bis zum infraroten Bereich untersucht und festgestellt, dass der Effekt bei kurzen Wellenlängen stär-ker ausgeprägt ist. Die Wechselwirkung dieser Moleküle mit ultrakurzen, intensiven Feldern führt zu schneller Ionisation und diese findet bevorzugt bei großen Kern-abständen statt. Eine Überlagerung des Grundzustandes und des ersten angeregten Vibrationszustandes oszilliert nach dem Pulsende weiter. Aus dieser periodischen
Schwingung konnten wir Schwingungsperioden extrahieren, die mit der Energiedif-ferenz zwischen den beiden niedrigsten Schwingungszuständen im Einklang waren. Diese Vibrationswellenpakete bleiben für extrem lange Zeit kohärent. Eine inter-essante Beobachtung wurde bei der Untersuchung des Lochfraßzustandes für unter-schiedliche Zeitverzögerungen in einem Pump-Probe-Schema gemacht. Wir fanden zusätzliche Oszillationen mit einer Periode, die kürzer als die gewöhnlichen Vibra-tionsperioden waren, vermutlich aufgrund der Bevölkerung höherer Schwingungsni-veaus.
In den TD-MCDFT- und den MCTDH-Rechnungen mussten wir absorbierende Randbedingungen anwenden, um Reflexionen von den Grenzen des Gitters zu ver-meiden. In Kapitel 4 haben wir die Auswirkungen von absorbierenden Rändern auf Mean-Field-Rechnungen untersucht, wobei die absorbierenden Ränder auf zwei unterschiedliche Arten behandelt wurden. Im Variationsansatz wurden die absor-bierenden Grenzen in Form eines komplexen Potentials direkt zum (Vielteilchen-)Hamiltonoperator addiert. Die Bewegungsgleichungen wurden dann mit dem Va-riationsprinzip hergeleitet. Im TDDFT-Ansatz werden die absorbierenden Grenzen stattdessen in die effektiven KS-Potentiale eingebaut, d.h. in die schon bestehenden Bewegungsgleichungen. Wir haben gezeigt, dass die Größe des Gitters fast keine Auswirkung auf die TDDFT-Ergebnisse hat, während die Ergebnisse für das Varia-tionsprinzip signifikant von der Gittergröße abhängen. Die effektiven Potentiale im Variationsansatz führten zu falschen Ergebnissen für kleine Gitter, da aufgrund der Absorption an den Gitterrändern eine effektive Erhöhung der Potentiale durch Re-normierung der Einteilchenwellenfunktion im Innenraum erfolgt. Die KS-Potentiale im TDDFT-Ansatz sind nicht von diesem Problem betroffen. So wurde gezeigt, dass sich das TDDFT-Konzept besser schlägt als die Variationsmethode, denn die TDDFT-Ergebnisse sind unempfindlich gegenüber der Gittergröße.
Abschließend möchten wir festhalten, dass eine korrekte Beschreibung der Elektron-Kern Kopplung, und nicht nur eine Beschreibung der Einzelbewegungen, notwendig ist um die Vorgänge in einem atomaren/molekularen System zu verstehen und um dessen Dynamik in Gegenwart eines externen Feldes zu beschreiben. Mit ihren inhä-renten Beschränkungen bei der Beschreibung der Elektron-Kern Kopplung versagt die Hartree Näherung bei der zuverlässigen Voraussage der Systemdynamik in An-oder Abwesenheit eines externen Feldes. Mit Hilfe der Vielteilchenmodellsysteme ha-ben wir zwei Ansätze vorgestellt, nämlich TD-MCDFT und MCTDH, die über die Hartree Näherung hinaus gehen und die Elektron-Kern Korrelationseffekte adäquat beschreiben.
Chapter 1
Introduction
Atoms, considered as indivisible a long time ago, are the basic building blocks of physical and biological systems and play a key role in the processes that occur in nature. In fact an atom [1] is divisible and is made up of a tightly packed nucleus in the center and a cloud of negatively charged electrons surrounding it. The nucleus is made up of the positively charged protons and neutrons which are neutral. The electrons are bound to the nucleus by the electromagnetic force. The hydrogen atom is the simplest of all atoms containing one negatively charged electron and a positively charged proton. The electric field that an electron in the first Bohr orbit of the hydrogen atom experiences is 5.14 × 1011 V/m. The atomic unit of intensity is defined using this value as 3.51 × 1016 W/cm2. Compared to this, the first operational laser developed by Maiman [2], after the theoretical description of lasers was developed by Schawlow and Townes [3], was rather weak. The attainable peak intensities in the 1960s were of the order of 109 W/cm2. Pulse intensities of the order of 1015 - 1018 W/cm2 were obtained during 1970-1990. Also the pulse duration was dramatically reduced from 100 picoseconds to 10 femtoseconds. Since then intensities of the order of 1021W/cm2 have been achieved and the pulse lengths attained are shorter than 1 femtosecond. As a result of this rapid development in the field of laser technology, nowadays, ultra-short high-intensity laser pulses are readily available.
Interactions between an atomic or a molecular system and a laser pulse give rise to interesting effects. Multiphoton processes, corresponding to the net absorption or emission of more than one photon in atomic or molecular transitions, are ob-served when atoms or molecules interact with the laser fields of sufficiently strong intensities. If the field strength of the laser pulse is much less than the atomic or molecular binding forces, the system dynamics arising from such field strengths can be explained using perturbation theory [4]. However, when the field strengths of the laser pulses are comparable to or larger than the typical atomic or molecular binding forces [5–7], the perturbation theory is no longer applicable. Due to the very strong field, the observed processes are highly nonlinear and a non-perturbative treatment
is required to address this issue. The interaction of atoms and molecules with such intense laser fields leads to a multitude of effects. In above-threshold ionization (ATI) [8], an electron escapes by absorbing more photons than necessary to over-come the ionization threshold. A particularly interesting case is the following: when an electron is accelerated and driven back to the core, processes such as excitation, further ionization, or high-order harmonic generation (HHG) [9] are observed. If the returning electron interacts with the parent ion inelastically, the process may lead to excitation or further ionization of the parent ion, potentially leading to fragmenta-tion in the case of molecules. Electrons may also scatter elastically from the parent ion. If the returning electron recombines with the parent ion by emitting a photon at higher-order multiples of the laser frequency, the process is known as HHG. HHG is important because it serves as a source of coherent extreme ultra-violet (XUV) radiation. HHG spectra from molecules provide information about the molecular structure [10] and have been used to retrieve electron orbitals [11–13].
For ultra-short high-intensity laser fields, the information about the multiphoton processes can be obtained by direct solution of the time-dependent Schrödinger equation (TDSE) [14]. The TDSE, although efficacious enough to describe all the relevant non-relativistic physics of atoms and molecules in strong laser fields, is too complex and time consuming to solve for all but very simple systems. For an ac-curate description of molecules in a strong field, the response of the electrons as well as the nuclei must be taken into account. Ideally, the dynamics of the different degrees of freedom should be treated fully quantum mechanically. The pioneering work of Born and Oppenheimer, also known as Born-Oppenheimer (BO) approx-imation, simplifies the problem by separating the nuclear and electronic motion [15–17]. However, the BO approximation cannot fully explain the complex interplay between the different degrees of freedom in situations where electrons are excited into high-lying bound or continuum states. The standard method of solving the full TDSE uses a representation of the wave function and Hamiltonian in an appropriate product basis. This approach works very well for small systems but the required computational resources grow exponentially with increasing number of degrees of freedom. Even simple molecular systems such as the H+
2 molecular ion pose a chal-lenge to theoretical description when they are subject to intense laser pulses because large grids or basis sets are required to capture the ionization and fragmentation dynamics in the non-perturbative regime. For H+
2, full non-BO quantum mechanical simulations incorporating the nuclear motion have been reported only for the case of frozen rotational degree of freedom. For the usual infrared fields, such simula-tions are limited to the case of the molecular axis aligned parallel to the laser field [18–20] while the case of arbitrary oriented three-dimensional H+
2 and H2 molecules have been treated in the fixed-nuclei approximation [21–26]. Recently the combined treatment of arbitrary orientation and nuclear motion has been reported [27, 28] for driving fields in the high-frequency (extreme ultraviolet) regime, which are com-putationally less demanding than the low-frequency fields. Another set of recent calculations [29, 30] includes the full rovibrational quantum mechanical motion of the nuclei, but uses a BO expansion including the three lowest electronic states.
3 As an alternative to the direct solution of the TDSE, approximate methods have been developed that maintain a fully quantum mechanical picture while alleviating the problem of exponentially growing computational requirements. Density func-tional theory (DFT) [31] is considered as a viable alternative. At the heart of static DFT lies the Hohenberg-Kohn (HK) [32] theorem. According to the HK theorem, to fully describe the stationary electronic system it is sufficient to know its ground-state density. From this quantity all observables can be obtained in principle. Density is a physical observable and it only depends on three spatial coordinates, in contrast to the N-particle many-body wave function, which is a complex function of 3N spa-tial coordinates. The HK theorem also describes a variational principle in terms of density, showing that the energy can be written as a density functional whose min-imum is attained at the exact density. A major breakthrough occurred when Kohn and Sham [33] proposed the use of an auxiliary system of non-interacting particles, known as the Kohn-Sham (KS) system, to obtain the density of the interacting par-ticles. The time-dependent analogue of the HK theorem is the Runge-Gross (RG) theorem [34]. The RG theorem shows that, for a given initial wave function, there is a unique mapping between the time-dependent external potential of a system and its time-dependent density. Again, a system of non-interacting particles can be used to determine the density, with the KS potential incorporating all the many-body effects of the system. It is the RG theorem that has been stated as the beginning of the modern-day time-dependent density functional theory (TDDFT) [35]. In the past two decades TDDFT has been used to describe various properties of atoms [36–38], clusters [39, 40], solids [41–43], semiconductor nanostructures [44–46], bio-logical systems [47–49], and fullerenes [50]. Attempts to use TDDFT in the optimal control theory are gradually gaining strong ground [51, 52]. TDDFT has been used to describe the dynamics of atomic/molecular [53–63], and cluster systems [64–66] in strong laser fields. Mostly, only the electrons have been treated fully quantum mechanically in the TDDFT applications mentioned above and the motion of nuclei has been neglected. However, the nuclear degrees of freedom play an essential role in a wide range of phenomena. Quantum transport [52, 67, 68] and superconductivity [69, 70] as well as ionization and dissociation dynamics of molecules in strong laser fields are examples where the coupling between the electrons and nuclei needs to be considered. Therefore a general framework is required that accounts for this coupling and treats the dynamics of those degrees of freedom fully quantum mechanically and on equal footing. To address this issue in the spirit of DFT, multicomponent DFT (MCDFT) has been proposed [71–73]. In Chapter 2, we first briefly discuss the time-dependent MCDFT (TD-MCDFT) and explain the multi-component Kohn-Sham (MCKS) scheme and the MCKS equations. The MCKS equations describe the system dynamics exactly, provided the Kohn-Sham (KS) potentials that govern the system are exact. In practice, one always has to work with approximations for these KS potentials. The KS potentials consist of three parts: (i) external poten-tials, (ii) Hartree potenpoten-tials, and (iii) exchange-correlation potentials. The trickiest part is to find an approximation for the exchange-correlation (xc) potentials. The xc potentials contain all the many-body effects of the system. To obtain insight into these xc potentials, we aim at finding exact potentials for small model systems.
If one has the knowledge of the exact KS potentials governing the system, then the xc potentials can be easily obtained. We describe an algorithm that constructs such exact KS potentials. With the help of a one-dimensional H+
2 molecular ion we describe photodissociation processes using a superposition of BO states or alterna-tively, by subjecting the model system to an external field. We show that the exact KS potentials reflect the dynamics and are piecewise similar to the BO potential energy surfaces. Neglecting the xc contributions from the KS potentials proves to be a fatal approximation as the Hartree potentials alone are not able describe the photodissociation dynamics, and we find that the system may not dissociate at all. The main purpose of constructing the exact KS potential is that its knowledge tells general properties, also for more complex systems. The exact potentials can also serve as a reference to check the quality of the multicomponent functionals to be developed in the future.
The time-dependent self-consistent field (TDSCF) approach [74] can also be con-sidered as a viable alternative. In the TDSCF approach, unlike DFT where the many-body wave function is represented as a functional of density, the many-body wave function is represented as a product of single-particle functions (SPFs) known as a Hartree product. This results in a set of coupled single-particle equations of motion for the SPFs while maintaining the quantum mechanical picture and abating the scaling problem. While the required computational effort is reduced significantly, the correlation between the degrees of freedom is lost, resulting in poor accuracy. To overcome the loss of correlation, the multiconfiguration TDSCF (MC-TDSCF) ap-proach has been proposed [75, 76]. The multiconfiguration time-dependent Hartree (MCTDH) method [77–79] is a variant of the MC-TDSCF approach. The MCTDH wave function is a sum of Hartree products. Here each Hartree product can be un-derstood as a configuration and the sum covers all possible configurations that arise from the set of single-particle functions (which are time-dependent themselves). Ap-plications of the MCTDH method are numerous. MCTDH has been used to calculate properties such as photodissociation and absorption spectra in various diatomic and polyatomic molecules [80–87]. Molecule-surface scattering [88–90], scattering cross-sections [91, 92] and nuclear dynamics during electron scattering processes [93, 94] have been described with the help of MCTDH. Computation of eigenenergies [95, 96] and eigenstates [97, 98] have also been performed. This indicates that MCTDH has been applied extensively to study the nuclear dynamics of multidimensional sys-tems. Its application has been extended to the case of fermionic systems [99–101], bosonic systems [102, 103] and mixtures of fermionic/bosonic systems [104, 105] as well. Similar to TDDFT, MCTDH is a contender in developing the solutions for op-timal control theory [106, 107]. In Chapter 3, we introduce the theory of MCTDH for a correlated system of an electron and nuclei and study the electron-nuclear correlation in molecules subject to strong external laser fields. In the presence of ultra-short intense laser pulses the electronic structure of the molecule is distorted to such an extent that it no more resembles its field-free structure. This leads to a strong interplay between the ionization and dissociation processes and the system dynamics is highly correlated. Using the example of a model H+
5 demonstrate that the MCTDH method can successfully describe the ground-state properties as well as multiphoton processes in the presence of strong fields with different wavelengths and intensities. Thus, it is a viable option to study larger systems. Furthermore, we study the two-center interference effect in HHG, i.e., the suppression of harmonics around a particular harmonic order, and the lochfraß effect which leads to the appearance of vibrational wave packets in the electronic ground state of H2 molecule upon strong-field ionization.
The purpose of the effective potential methods TDDFT and TDSCF is to over-come the scaling problem. In the numerical simulations describing the interactions of molecules with intense laser fields, inclusion of absorbing boundaries is nearly inevitable as most of the simulations are performed on finite sized grids. A large amount of required computational time may be saved if one works with small grids. However, numerical simulations on a small grid are often hampered by the spuri-ous reflections caused by the grid boundaries [108]. These reflections may interfere with the propagating wave packet in the inner region and in turn give unphysical results. Absorbing boundaries have been used to eliminate such reflections. Simi-lar to the direct solution of the TDSE, absorbing boundaries can be applied in the TDDFT and TDSCF approaches. The effects of including the absorbing boundaries in the above mentioned effective-field methods is studied in Chapter 4. We first derive and discuss the equations of motion for the model H+
2 molecular ion using the two approaches TDDFT and TDSCF. We then show that TDDFT and TDSCF, when including absorbing boundaries, differ for small and moderate grids. While the TDDFT approach gives consistent results independent of the grid size, the TD-SCF approach shows inconsistency. Effective-potential methods are widely used and therefore the present work has profound implications. The same deficiency may be found for the multiconfigurational methods with few configurations.
As we discussed earlier, even for the H+
2 molecular ion, solving the full three-dimensional TDSE is a daunting task. Model systems with reduced dimensions are often employed to understand the complex behavior of various systems. With the help of such model systems, one can reduce the complexity of the real system and solve the problem numerically exactly. In this work, we use the one-dimensional H+ 2 molecular ion and H2 molecule [109–115] to compare the MCDFT and the MCTDH calculations with exact calculations. Unless stated otherwise atomic units are used in the thesis.
Chapter 2
Multicomponent density functional
theory
In the past two decades DFT has been used extensively to describe atomic and molecular systems. The first step of DFT is to prove the HK theorem [32]. The HK theorem demonstrates that there exists a one-to-one correspondence between the external potentials and the density. From a given potential one can solve the Schrödinger equation and obtain the many-body wave function and the density. The HK states that from the density it is possible to obtain the external potential and hence the many-body wave function. It can be very hard to obtain the density of an interacting system. In modern DFT, one introduces an auxiliary system of non-interacting particles, with the same density as the system of non-interacting particles. This auxiliary system is governed by a local effective potential termed as Kohn-Sham (KS) potential [33].
Mostly the calculations done with DFT have only treated electrons quantum me-chanically while the nuclei have been considered as static. This means the electron-nuclear interaction only enters as a static contribution to the electronic KS potential. A classical treatment of nuclei is possible where the Newtonian equations for nuclei are solved at the instantaneous nuclear positions [116, 117]. When the atoms and molecules are subjected to strong external fields, a further improvement in the the-ory is required that treats the electron and nuclei on the same footing. To address this issue, multicomponent density functional theory (MCDFT) has been proposed for stationary systems [73, 119]. In the next section we review the MCDFT for the time-dependent problems, which has been introduced in [71–73].
2.1
Multicomponent density functional theory: An
overview
The Hamiltonian for the complete system of Ne electrons and Nn nuclei reads b H(t) = bTn(R) + bTe(r) + cWen(r, R) + cWnn(R) + cWee(r) + bUext,n(R, t) + bUext,e(r, t) (2.1.1) b Tn= Nn X ν=1 − ∇ 2 ν 2Mν ; bTe = Ne X j=1 −∇ 2 j 2m (2.1.2) c Wnn = 1 2 Nn X ν6=µ ZνZµ |Rν − Rµ| ; cWee = 1 2 Ne X j6=k 1 |rj − rk| (2.1.3) c Wen= − Ne X j=1 Nn X ν=1 Zν |rj− Rν| (2.1.4) b Uext,n = Nn X ν=1 Uext,n(Rν, t) ; bUext,e = Ne X j=1 Uext,e(rj, t) . (2.1.5) Eqs. (2.1.2) denote the kinetic-energy operators of the nuclei and electrons, re-spectively. Eqs. (2.1.3) and (2.1.4) represent the interparticle Coulomb interac-tions. Equation (2.1.5) represents all possible external potentials acting on the system. Here we note that the interaction between the electrons and nuclei are treated within Eq. (2.1.4) and do not contribute to the external potentials. All the electronic coordinates are labeled by {rj} ≡ r and the nuclear coordinates by {Rν} ≡ R. The above Hamiltonian provides a quantum mechanical description of all degrees of freedom and the evolution of the wave function is governed by the TDSE i∂tψ(t) = bH(t)ψ(t). The external potentials vanish for isolated systems. In this case, the Hamiltonian remains unchanged under translations, reflecting the conservation of total momentum.
2.1.1
Coordinate transformation
To describe such isolated systems, the translational and rotational symmetries must be accounted properly. To this end, two coordinate transformations are performed. First the laboratory coordinates of the electron are transformed to coordinates de-fined relative to the nuclear center of mass (CM). Then these relative coordinates are rotated into the body-fixed internal coordinates r′. The new electronic coordinates are written as
2.1.1 Coordinate transformation 9 where the term RCMN = PNn1
ν=1Mν
PNn
ν=1MνRν describes the nuclear center-of-mass. In Eq. (2.1.6) R defines a three-dimensional orthogonal matrix. It represents the Euler rotations [118]. The Euler angles (α, β, γ) are functions of the nuclear coor-dinates R and specify the orientation of the body-fixed coordinate frame. To keep the derivation flexible, the nuclear coordinates have been left unchanged
R′
ν = Rν. (2.1.7)
A detailed discussion and derivation of coordinate transformation is given in [35, 73, 119–121]. After the coordinate transformation the Hamiltonian reads
b
H(t) = bTn(R) + bTe(r′) + bTMPC(r′, R) + cWnn(R) + cWee(r′) + cWen(r′, R)
+ bUext,n(R, t) + bUext,e(r′, R, t) . (2.1.8)
In the transformed Hamiltonian of Eq. (2.1.8) bTn, bTe, cWee, cWnn, and bUext,n are left unchanged from their original form in Eqs. (2.1.2) and (2.1.3). bTMPC denotes mass-polarization and Coriolis contributions
b TMPC = Nn X ν=1 − 2M1 ν ∇Rν + Ne X j=1 ∂r′ j ∂Rν∇r ′ j !2 − bTn(R). (2.1.9) As the electronic coordinates are transformed to the CM-fixed coordinates, the mass-polarization terms appear. Including their diagonal elements into bTe/n changes the bare masses to the reduced masses. The transformation to body-fixed internal co-ordinates introduces additional electron-nuclear coupling terms known as Coriolis coupling. The electrons are now referred to a non-inertial coordinate frame which rotates with the nuclei, such that the electrons are subjected to a fictitious Coriolis force. The electron-nuclear Coulomb interaction term, after coordinate transforma-tion, reads c Wen(r′, R) = − Ne X j=1 Nn X ν=1 Zν |r′ j − R′ν| . (2.1.10) The term R′
ν = R(α, β, γ) (Rν−RCMN), in Eq. (2.1.10), is invariant under rotations and translation of the nuclear framework. The external potential acting on the electrons now depends on the electronic and nuclear coordinates and reads
b Uext,e(r′, R, t) = Ne X j=1 Uext,e(R−1 r ′ j + RCMN, t) := Ne X j=1 uext,e(r ′ j + R, t). (2.1.11)
Eq. (2.1.11) is no longer a one-body operator, but now acts as an effective interac-tion. The coordinate transformation described above reflects the internal symmetry of the system. The transformed Hamiltonian in Eq. (2.1.8) will now be used as the starting point in the formulation of MCDFT.
2.1.2
Defining the densities
The electronic and nuclear single-particle densities defined in terms of the inertial coordinates r and R are not useful as they reflect the symmetry of the corresponding external potentials. For example, the densities are constant for isolated systems and fail to describe the internal properties of the system.
The basic requirements for the electronic density ρ and nuclear density Γ are the following: 1) The basic electronic density should be a single-particle quantity. 2) The densities should be characteristic for the internal properties of the system. 3) The treatment of the nuclear degrees of freedom should allow for appropriate description of physical situations such as photodissociation. A set of densities which meets these requirements is given by [73, 119] Γ(R, t) = Z dr′|ψ(R, r′, t)|2 , (2.1.12) ρ(r′, t) = Ne Z dR Z dr′2· · · dr′Ne|ψ(R, {r ′ , r′2, · · · r′Ne}, t)| 2 , (2.1.13) where ψ(R, r′, t) is the full solution of the TDSE with the Hamiltonian of Eq. (2.1.8). Here we note that the densities in Eqs. (2.1.12) and (2.1.13) are defined as functions of the transformed coordinates. The electronic density in Eq. (2.1.13) gives the number density at position r′as measured from the certain orientation of the nuclear framework. This truly reflects the internal properties of the system. The nuclear degrees of freedom are described using the diagonal of the nuclear density matrix in Eq. (2.1.12).
2.1.3
The Runge-Gross theorem for multicomponent system
The Runge-Gross theorem [34] shows that in a system of many identical particles, for a given initial state, there is a one-to-one correspondence between time-dependent densities and time-dependent potentials. The theory is based on the TDSE i∂ψ(t) ∂t = bH(t)ψ(t). To set up the framework for TD-MCDFT, the RG theorem has to be extended to multicomponent systems.
The Hamiltonian of Eq. (2.1.8) is slightly modified to take the form b
H(t) = bTn(R) + bTe(r′) + bTM P C(r′, R) + cWee(r′) + cWen(r′, R)
2.1.3 The Runge-Gross theorem for multicomponent system 11 The potentials bVext, e and bVext, n are given by
b Vext, e(r, t) = Ne X j=1 vext, e(rj, t) , (2.1.15) b Vext, n(R, t) = vext, n(R, t) . (2.1.16)
The term bUext, e(r′, R, t) is treated as a fixed many-body term as it depends on both electronic and nuclear coordinates. The term bTMPC(r, R) is treated as an additional electron-nuclear interaction.
On the basis of the densities given in Equations (2.1.12) and (2.1.13), the Runge-Gross theorem which is the time-dependent analogue of the Hohenberg-Kohn theo-rem can be proved for multicomponent systems. Considering the time evolution of a molecule from an arbitrary initial state ψ(t0) = ψ0 under the influence of different external potentials, the Schrödinger equation defines a map which uniquely assigns the many-body wave function to the set of external potentials:
A : {vext, n(R, t), vext, e(r′, t)} −→ ψ(t). (2.1.17) Employing the densities of Equations (2.1.12) and (2.1.13), a second map is estab-lished between the many body wave function and the set of time-dependent densities,
B : ψ(t) −→ {Γ(R, t), ρ(r′
, t)} . (2.1.18)
The main issue here is the proof of the invertibility of the combined map C ≡ A ◦ B: C : {vext, n(R, t), vext, e(r′, t)} −→ {Γ(R, t), ρ(r′, t)} . (2.1.19) In [71], the Runge-Gross theorem has been extended to multicomponent systems. The statement is the following. Two sets of densitiesΓ′(R, t), ρ′(r′, t) andΓ(R, t), ρ(r′, t) , evolving from the same initial state ψ
0 under the influence of two sets of potentials v′
ext, n(R, t), v′ext, e(r′, t)
and vext, n(R, t), vext, e(r′, t)
differ infinites-imally after t0, provided at least one component of the potentials differs by more than a time-dependent function,
vext, n(R, t) 6= vext, n′ (R, t) + C(t) ∨
vext, e(r′, t) 6= vext, e′ (r ′
, t) + c(t) . (2.1.20)
Consequently, the rigorous one-to-one mapping between time-dependent densities and external potentials is established for a fixed initial state ψ0.
{vext, n(R, t), vext, e(r′, t)} 1−1
←→ {Γ(R, t), ρ(r′, t)} (2.1.21)
by solving the TDSE. The time-dependent density can be calculated from this wave function. If the same density can be reproduced by a system of non-interacting particles for some external potential then the density is said to be non-interacting v-representable. The non-interacting {ve}-representability has been proved only for electronic TDDFT [122]. Analogous to the static MCDFT [73], the non-interacting {Vn, ve}-representability is assumed to hold for the TD-MCKS scheme. Eq. (2.1.21) only assures the uniqueness and not the existence of the effective potentials. The one-to-one correspondence of Eq. (2.1.21) indicates that the full many-body wave function is determined by this set of densities up to within a time-dependent phase α(t),
Ψ(t) = e−iα(t)Ψ[Γ, ρ](t).˜ (2.1.22)
As a consequence of this, every observable O of the time-dependent system is a unique functional of the time-dependent densities
O[Γ, ρ](t) = h ˜Ψ[Γ, ρ](t)| bO| ˜Ψ[Γ, ρ](t)i , (2.1.23) and the ambiguity in the phase cancels out. The above functional implicitly depend on Ψ0, because the one-to-one correspondence can be established only for a fixed initial state. If this initial state is a non-degenerate ground state then the functional will depend on the densities alone, since Ψ0 is a unique functional of its densities {Γ0, ρ0}.
2.1.4
The time-dependent Kohn-Sham scheme for
multicom-ponent system
The multicomponent Kohn-Sham (MCKS) system consists of non-interacting elec-trons and nuclei. The Hamiltonian for such an auxiliary system, termed as Kohn-Sham Hamiltonian, is given by
b HKS(t) = bTn(R) + bTe(r′) + bVKS, n(R, t) + bVKS, e(r′, t). (2.1.24) Here, b VKS, n(R, t) = vKS, n(R, t), (2.1.25) b VKS, e(r′, t) = Ne X j=1 vKS, e(r′, t), (2.1.26)
and the potentials bVKS, n and bVKS, e are termed as Kohn-Sham potentials. Ac-cording to the multicomponent Runge-Gross theorem there is a set of potentials bVKS, n(R, t), bVKS, e(r′, t)
that reproduces a given set of densityΓKS, n(R, t), ρKS, e(r′, t)
. The potentials are therefore functionals of the densities Γ and ρ. This is due to the fact that the one-to-one correspondence between the time-dependent potential and densities can be established for any given interaction. The system considered here
2.1.4 The time-dependent Kohn-Sham scheme for multicomponent system 13 is noninteracting. This leads to the separation of the electronic and nuclear motion. If the initial KS wave function is chosen to be a product of a nuclear and electronic wave function, then the time-dependent KS wave function will also be in a product form.
ΦKS(R, r′, t) = χ(R, t)φ(r′, t) (2.1.27)
If χ and φ are suitably normalized, the corresponding densities are given by
Γ(R, t) = |χ(R, t)|2 , (2.1.28) ρ(r′, t) = Ne X j=1 |φj(r′, t)|2 , (2.1.29)
where χ and φj are the solutions of an Nn-particle nuclear and a single-particle electronic Schrödinger equation
i∂ ∂tχ(R, t) = − X α ∇2 α 2Mα + VKS, n[Γ, ρ](R, t) ! χ(R, t) , (2.1.30) i∂ ∂tφj(r ′ , t) = −∇ 2 2 + vKS, e[Γ, ρ](r ′ , t) φj(r′, t) . (2.1.31)
It is noteworthy that bVKS, n in the nuclear KS equation (2.1.30) is an Nn-body inter-action while bVKS, e is a one-body operator. Therefore, the electronic wave function is chosen as a Slater determinant built from the orbitals φj.
The key assumption in the MCKS scheme is that the local effective potentials {bVKS, n, b
VKS, e} exist such that the densities of the auxiliary system reproduce the exact densities of the fully interacting system. The effective potentials reproducing the exact densities can be decomposed as
VKS, n[Γ, ρ](R, t) = vext,n(R, t) + VHxc,n[Γ, ρ](R, t) , (2.1.32) vKS, e[Γ, ρ](r′, t) = vext,e(r′, t) + vHxc,e[Γ, ρ](r′, t) , (2.1.33) where vext,n and vext,e are the true external potentials as given in Eqs. (2.1.15) and (2.1.16). For each of the two potentials, the remaining term with subscript “Hxc” denotes the sum of Hartree, exchange and correlation potentials. The Hartree potentials account for the classical electrostatic interaction between the particles and written as [72] VH,n(R, t) = Z d3r′ Wen(r′, R) + uext,e(r′, R, t) ρ(r′, t), (2.1.34)
vH,e(r′, t) = Z d3R1· · · Z d3RNn Wen(r ′, R) + u ext,e(r′, R, t) Γ(R, t) + Z d3r′′ W ee(r′, r′′)ρ(r′′, t). (2.1.35)
Eqs. (2.1.10) and (2.1.11) have been used in above equations. The remaining xc contributions to the KS potentials contain all the many-body xc effects. The ex-pressions for Vxc,n and vxc,e as functional of the densities are defined as follows.
Vxc,n(R, t) = δAxc[Γ, ρ] δΓ(R) , (2.1.36) vxc,e(r′, t) = δAxc[Γ, ρ] δρ(r′) , (2.1.37)
where Axcis the xc part of the multicomponent action A = i lnhΨ0| bTce−i R
Cdt bH(t)|Ψ0i.
b
Tc denotes the time-ordering operator along the Keldysh time contour C. Ψ0 is the initial state of the system. For a detailed overview of the multicomponent action and the derivation of xc potentials, see [35]. Explicit functionals for the xc potentials are not easily obtained and have to be approximated.
Eqs. (2.1.28)–(2.1.33) are the essence of the time-dependent MCKS (TD-MCKS) scheme. This scheme implies that the dynamics of an interacting many-particle system can be calculated by propagating a set of coupled equations. We note that no approximations were made to derive the MCKS Eqs. (2.1.28)–(2.1.33) and the TD-MCKS equations include all many-body effects. The effective potentials in Eqs. (2.1.30) and (2.1.31) are time-dependent and adapt to the situation the system is subjected to. With this adaptability, the TD-MCKS potentials are able describe the non-adiabatic dynamics in the presence of strong fields that takes place over a broad range of potential energy surfaces (PESs), e.g. ionization, harmonic generation. This is unlike the time-dependent Born-Oppenheimer (BO) approximation, where the system dynamics is restricted to a fixed number of PESs. For the description of the BO approximation and PES please refer to Appendix A. The electronic properties do not depend parametrically on the nuclear configuration but the information is included through the functional dependence on Γ. Similarly any change in the electronic density is conveyed to the nuclear TD-MCKS potential through functional dependence on the electronic density ρ. The TD-MCKS Eqs. (2.1.30) and (2.1.31) are coupled and must be solved simultaneously. Even with all these simplifications the numerical solution may become quite demanding for large systems. Further simplification can be made by treating only the valence electrons dynamically while adding the core electrons to the nuclei via a suitable pseudopotential [123]. In the presence of a strong external field the nuclei cannot be considered as static and the nuclear motion (vibration and rotation) has to be taken in to account. The TD-MCKS scheme allows us to treat the electron and nuclear degrees of freedom on equal footing. However, the success of the proposed approach relies on the quality of approximations for the “Hxc” potentials.
2.2 The model system 15
2.2
The model system
In this section we set up a model system to understand the laser-matter interactions. Our objective is to study the dynamics of diatomic molecules subjected to intense ultrashort laser pulses. The particle dynamics is described non-relativistically and the dipole approximation is used. With this, the laser-matter interaction can be described using the TDSE i∂ψ(t)
∂t = H(t)ψ(t). The magnetic field component of the laser pulse is neglected. The external radiation is treated classically.
After performing the coordinate transformations as described earlier, the Hamilto-nian of the system transforms into
b H(t) = − ∇ 2 R 2µn + cWnn+ Ne X j=1 −∇ 2 rj 2 ! − 2M1 nuc Ne X i,j=1 ∇ri∇rj + bTC +cWee+ cWen− qnR− qe Ne X j=1 R−1rj ! E(t), (2.2.1)
where qn = Z1MM2−Znuc2M1 and qe = Z1+ZM2total+Mnuc are the reduced charges. Mnuc and Mtotal are the total nuclear mass and total mass respectively. E(t) is the electric field vector of the laser pulse. µn is the reduced mass of the nuclei and R = R1− R2 is the internuclear distance vector. Here R is the three-dimensional orthogonal matrix representing the Euler rotations [118] and bTC is the Coriolis coupling term. This Hamiltonian provides a basis to study the laser-matter interaction in diatomic molecules.
The strength of the laser fields we usually work with are strong enough so that the perturbation theory [124] treatment is not adequate. In the traditional ways, one can expand the full molecular wave function in terms of Born-Oppenheimer (BO) states (see appendix A) and solve the resulting time-dependent coupled equations for the nuclei. The nuclear dynamics is treated on a few lowest potential energy surfaces but a large number of electronic states is required as the field strength is strong and leads to a time consuming calculation. Ionization processes are expected to take place when working with strong fields. The Born-Oppenheimer approach is not capable of describing the ionization dynamics correctly. The highly excited and ionized electrons can be treated within the clamped nuclei approximation [125]. But this approach is constrained in a way that it leaves the nuclei fixed and cannot explain the dissociation dynamics. In [19], it was shown that the exact solution for the strong field dynamics is possible for H+
2 molecule. But even for this simple molecule, a full three-dimensional calculation is a daunting task.
2.2.1
One-dimensional H
+2molecular ion
We consider a H+2 molecular ion, perfectly aligned with the polarization axis of the external field while neglecting the Coriolis force. The nuclear motion is treated as one-dimensional. This model system is characterized by the Hamiltonian
b H(x, R, t) = bTn+ bTe+ Wen1D(x, R) + Wnn1D(R) + qexE(t) (2.2.2) with bTn = −2µ1n d2 dR2 and bTe = −2µ1 e d2
dx2. Here µn and µe denote the reduced mass
of nuclei and electron respectively. “x” denotes the electron coordinate and “R” denotes the internuclear distance. In our model the bare Coulomb potentials are replaced by soft-Coulomb potentials,
Wen1D(x, R) = −q 1 (x −R 2)2+ ae −q 1 (x + R 2)2+ ae , (2.2.3) Wnn1D(R) = √ 1 R2+ a n . (2.2.4)
Here ae = 1 and an = 0.03 are the soft-core parameters. qe = 2m2mnn+m+2e is the reduced charge and E(t) is the applied external field. The ground state is obtained by the split-operator method [126] and imaginary time propagation [127]. This ground-state wave function ψ(x, R) gives the exact densities Γ(R) and ρ(x) according to
Γ(R) = Z dx|ψ(x, R)|2 , (2.2.5) ρ(x) = Z dR|ψ(x, R)|2 . (2.2.6) The 1D H+
2 model system considered here allows to analyze the performance of various approximations methods. The computation effort is reduced significantly compared to the 3D H+
2. As the rotational degrees of freedom are not present, the phase space accessible by this model system is different from the real three-dimensional molecules. Therefore the model system provides only a qualitative description of the system dynamics.
2.3
Construction of the exact Kohn-Sham potential
The decomposition of the exact TDKS potential for nuclei and electron is given in Eqs. (2.1.32) and (2.1.33) respectively. The external potentials are given and the time-dependent Hartree potential contributions are easily obtainable from Eqs. (2.1.34) and (2.1.35). The time-dependent exchange-correlation potential terms
2.3 Construction of the exact Kohn-Sham potential 17 Vxc,n and vxc,e given in Eqs. (2.1.36) and (2.1.37), incorporating all the many-body effects, have to be approximated. If one finds a good approximation for the exchange-correlation potentials, the DFT calculation becomes accurate. In practice, one inserts the time-dependent densities into a known ground-state functional with the hope that it gives a reasonable description of the dynamics. In the presence of strong fields, where non-adiabatic processes can be important, such an approach may meet with astonishing success or a grave failure. In [62] Lein and Kümmel presented an approach to obtain the exact time-dependent KS and exchange-correlation po-tentials for the strong-field dynamics of a model He atom. Note that this approach is viable only for small model systems. The numerical computation is not feasible for large systems. We start our discussion with the time-dependent Schrödinger equation for the 1DH+
2 model system,
i∂ψ(x, R, t)
∂t = bH(t)ψ(x, R, t). (2.3.1)
In order to generate the exact KS potential two hurdles have to be overcome. First we have to generate the exact density of the system. Second, an algorithm to construct the KS potential from the density is needed.
By time-propagation of the initial state, we get the nuclear (electronic) density by integrating the absolute square of 2D wave function, |ψ(x, R, t)|2 over the electronic (nuclear) coordinates, Γ(R, t) = Z |ψ(x, R, t)|2dx, (2.3.2) ρ(x, t) = Z |ψ(x, R, t)|2dR. (2.3.3)
Similarly we evaluate the nuclear and electronic current from the following equations: jn(R, t) = 1 µn Re Z ψ∗(x, R, t) ∂ i∂Rψ(x, R, t)dx, (2.3.4) je(x, t) = 1 µe Re Z ψ∗(x, R, t) ∂ i∂xψ(x, R, t)dR. (2.3.5)
Equipped with density and current we can construct the KS orbitals as follows:
χ(R, t) =pΓ(R, t)eiαn(R,t), (2.3.6)
φ(x, t) =pρ(x, t)eiαe(x,t), (2.3.7)
where αn and αe are the phases of the respective orbital. From the Runge-Gross theorem we can say that the densities Γ, ρ of the interacting-particle system are identical to the densities ΓKS, ρKS of the non-interacting particles. Considering the equation of continuity
one can find that ∂je ∂x = ∂jKS, e ∂x , (2.3.9) ∂jn ∂R = ∂jKS, n ∂R . (2.3.10)
In 1D finite systems the current is zero at the large distances. This leads to the conclusion that the Kohn-Sham currents jKS,e, jKS,n are identical to the ordinary current je, jn which are needed to get Eqs. (2.3.11) and (2.3.12). This generalization is not valid in the 3D case.
The phases αn and αe are calculated up to a time-dependent constant by using the KS orbitals as follows. jn = 1 µn Re χ∗(R, t) ∂ i ∂ Rχ(R, t) jn = 1 µn RepΓ(R, t)e−iαn(R,t) ∂ i ∂ R p Γ(R, t)eiαn(R,t) jn = 1 µn RepΓ(R, t)e−iαn(R,t) p Γ(R, t)eiαn(R,t)∂ αn(R, t) ∂ R + e iαn(R,t)∂ p Γ(R, t) i ∂ R ! jn = 1 µn Γ(R, t)∂ αn(R, t) ∂ R
A rearrangement of above equation gives the equation for phase αn ∂αn(R, t)
∂R = µn
jn(R, t)
Γ(R, t). (2.3.11)
Similarly, the equation for phase αe is ∂αe(x, t)
∂x = µe
je(x, t)
ρ(x, t). (2.3.12)
After obtaining the KS orbitals, the KS potential for nuclei as well as electrons can be obtained by inversion of the split-operator time propagator [62]:
VKS,n(R, t) = − 1 2δtarcsin " Im e +i bTnδtχ(R, t + δt) e−i bTnδtχ(R, t − δt) !# + C(t), (2.3.13) vKS,e(x, t) = − 1 2δtarcsin " Im e +i bTeδtφ(x, t + δt) e−i bTeδtφ(x, t − δt) !# + c(t). (2.3.14) Here bTn and bTe are the single-particle kinetic energy operators for nuclei and elec-tron respectively. The operators e+i bTnδt and e+i bTeδt are applied by multiplication in
momentum space. Appropriate time-dependent constants C(t) and c(t) are added to the resulting potentials so that the potentials vanish in the asymptotic region
2.4 Dynamics on Born-Oppenheimer potential energy surface 19 (VS,n(r, t) → 0 for R → ∞ and VS,e(x, t) → 0 for |x| → ∞).
The nuclear correlation potential Vc,n and electronic correlation potential vc,e are defined by
VKS,n(R, t) = Wnn1D(R) + VH,n(R, t) + Vc,n(R, t), (2.3.15) vKS,e(x, t) = vH,e(x, t) + vc,e(x, t). (2.3.16) No exchange contribution is present as the model system consists of two distin-guishable degrees of freedom, namely the electron and the internuclear distance. Electron-electron Coulomb interaction is not present in the case of the H+
2 molecu-lar ion. In a general many-body system, an individual electron orbital is governed by the interaction of that electron with the nuclei as well as with all the other electrons. The interaction of the electron with the remaining electrons causes the presence of an exchange potential. The interaction of the electron with itself, known as self-interaction, also contributes to the Hartree and exchange potential. These self-interaction Hartree and self-exchange terms cancel each other, so they have been omitted in Eq. (2.3.16). The Hartree potentials in Eqs. (2.3.15) and (2.3.16) are therefore given by VH,n(R, t) = Z Wen1D(x, R)ρ(x, t)dx, (2.3.17) vH,e(x, t) = Z Wen1D(x, R)Γ(R, t)dR. (2.3.18)
Therefore the correlation potentials can be easily calculated from Eqs. (2.3.15) and (2.3.16) once the KS potentials are known.
2.4
Dynamics on Born-Oppenheimer potential
en-ergy surface
We first study the laser-field-free dynamics of 1D H+
2. To set up a suitable initial state we perform BO calculations for the model system discussed in Sec. 2.2.1. In Fig. 2.1(a) we plot the ground state BO potential energy surface (BOPES), the ground state KS potential, and the ground state Hartree-only KS potential VH, n(R) + Wnn1D. The corresponding nuclear densities are plotted in Fig. 2.1(b), respectively. The Hartree nuclear density plotted in Fig. 2.1(b) is obtained using a self-consistent Hartree calculation. We note that we have used the exact ground-state density to calculate the KS potential and Hartree-only KS potential in Fig. 2.1(a). The ground-state KS potential is obtained by inserting the ground-state nuclear density into the inverted KS Schrödinger equation,
VS,n(R) = 1 2µn 1 χ(R) d2χ(R) d R2 + const, (2.4.1)
1 2 3 4 5 6 Internuclear distance (a.u.)
-0.8 -0.6 -0.4 -0.2 0 0.2 Energy (a.u.) BO-PES KS potential Hartree-only KS potential 1 2 3 4
Internuclear distance (a.u.) 0 0.5 1 1.5 2 Γ (a.u.) BO density Exact density Hartree density (a) (b)
Figure 2.1: A comparison of the (a) ground state BO-PES, ground state KS po-tential, and the ground state Hartree-only KS potential and (b) BO ground state nuclear density, exact ground state nuclear density, and the Hartree nuclear density. where χ=√Γ. Excellent agreement is found when comparing the ground-state nu-clear KS potential with the ground state BOPES of the model system in Fig. 2.1(a) and the densities also agree very well in Fig. 2.1(b). The ground state Hartree-only KS potential, obtained from the exact density, follows the ground state BOPES and ground state KS potential for small internuclear distance but deviates for large internuclear distance. If this Hartree-only KS potential is used to obtain the ground-state nuclear density, the density will be similar to the Hartree density plotted in Fig. 2.1(b) with a narrow shape. Note that in the region where the density is close to zero, the inversion of nuclear KS potential leads to numerical instabilities. There-fore in such a region we cannot evaluate the nuclear KS potential reliably. This can be seen in Fig. 2.1(a) beyond the internuclear distances of 5.6 a.u.
We define three different initial states as given below in Eqs. (2.4.2), (2.4.3), and (2.4.4) for the real-time propagation study. We first obtain the nuclear ground-state wave function χBO, gr
n (R), which is calculated within the BO approximation. The initial state is then constructed by taking the product of χBO,gr
n (R) with an electronic state. The three cases considered here are
ψinitial, gr(x, R) = χBO, grn (R) · ΦBO, gre (x|R) , (2.4.2) ψinitial, ex(x, R) = χBO, grn (R) · ΦBO, exe (x|R) , (2.4.3)
ψinitial, ss(x, R) = χBO, grn (R) · α ΦBO, gre (x|R) + β ΦBO, exe (x|R)
, (2.4.4) where ΦBO,gr
e (x|R) and ΦBO,exe (x|R) are the electronic ground state and first excited state respectively. Here α=0.9 and β=0.4358 indicate the weights in the superposi-tion (α2+ β2 = 1). The initial state in Eq. (2.4.2) represents a bound state (within the BO approximation) that does not change in time unless perturbed by an external
2.4 Dynamics on Born-Oppenheimer potential energy surface 21
0 2 4 6 8
Internuclear distance (a.u.) -0.80 -0.60 -0.40 -0.20 0.00 Energy (a.u.)
Ground state BO-PES
0 2 4 6 8
Internuclear distance (a.u.) 1st excited state BO-PES
0 2 4 6 8
Internuclear distance (a.u.) BO nuclear ---- ground-state
(a) (b) (c)
Figure 2.2: Illustration of the three different initial states used for the time-dependent calculations for the 1D H+
2 molecular ion, see text.
field. For such a state the KS potential should be similar to the ground-state BOPES shown in Fig. 2.1(a). The initial state of Eq. (2.4.3) represents a dissociative state which leads to the wave packet motion on the excited-state BOPES. Therefore the corresponding KS potential should also be dissociative. Eq. (2.4.4) represents a superposition state of a part that remains bound in the ground state and a part that photodissociates. The different initial states mentioned in Eqs. (2.4.2), (2.4.3), and (2.4.4) are illustrated in Figs. 2.2(a), 2.2(b), and 2.2(c).
For the initial state Eq. (2.4.3), the nuclear density at different times is shown in Figs. 2.3(a), 2.3(b), 2.3(c), and 2.3(d) for the times 2 fs, 4 fs, 6 fs, and 8 fs. The corresponding nuclear KS potentials and the Hartree-only potentials are plotted in Figs. 2.3(e), 2.3(f), 2.3(g), and 2.3(h). The KS potential is similar to the first excited state BOPES while the Hartree-only KS potential fails to reproduce the repulsive KS potential. As said before, we cannot construct the exact KS potential reliably where the nuclear density is very small. This leads to the unphysical oscillation at small and large internuclear distances.
Eq. (2.4.4) describes an initial state where the electronic state is a superposition of ground and excited states (see Fig. 2.2(c)). This mimics the physical situation of photodissociation. Fig. 2.4(a) shows the nuclear density after 2 fs. The correspond-ing KS potential and Hartree-only KS potential are shown in Fig. 2.4(e). The major part of the wave packet remains bound in the ground electronic state while a small portion of the wave packet starts to move on the excited state BOPES. Therefore, the nuclear density splits up into two parts: (i) a part that remains bound in the ground state BOPES and (ii) a part that moves on the first excited state BOPES. This splitting of nuclear density can be seen at R = 3 a.u. in Fig. 2.4(a). The corre-sponding KS potential in Fig. 2.4(e) shows a bound potential in the range of 1-3 a.u.
0.00 0.50 1.00 1.50 2.00 Γ (a.u.) Nuclear density -0.80 -0.60 -0.40 -0.20 Energy (a.u.)
Excited state BOPES
0.00 0.50 1.00 1.50 2.00 Γ (a.u.) -0.70 -0.60 -0.50 -0.40 Energy (a.u.) 0.00 0.50 1.00 1.50 2.00 Γ (a.u.) -0.70 -0.65 -0.60 -0.55 Energy (a.u.) KS potential 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
0.00 0.50 1.00 1.50 2.00 Γ (a.u.) 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
-0.75 -0.70 -0.65 Energy (a.u.) Hartree-only KS potential (a) (b) (c) (d) (e) (g) (f) (h)
Figure 2.3: Nuclear density, KS potential, and Hartree-only KS potential at t=2 fs [panels (a) and (e)], t=4 fs [panels (b) and (f)], t=6 fs [panels (c) and (g)], and t=8 fs [panels (d) and (h)] obtained for the 1D H+
2 molecular ion starting from the initial state in Eq. (2.4.3).
and a dissociative potential in the range R = 3-5 a.u. The ground state and excited state BOPES are also plotted for comparison. The corresponding Hartree-only KS potential fails to reproduce this combination of bound and repulsive potential. At 4 fs, the two wave packets have fully separated as seen in Fig. 2.4(b). The corre-sponding KS potential in Fig. 2.4(f) now clearly resembles the ground state and excited state BOPES. Figs. 2.4(c) and 2.4(d) show the nuclear density at 6 fs and 8 fs respectively. Now the dissociation process becomes clearer with part of the den-sity remaining bound between 1-4 a.u. and the other part moving on the repulsive potential. The corresponding KS potentials for time 6 fs and 8 fs are shown in Figs. 2.4(f) and 2.4(g) respectively. Here we can clearly see that the KS potential is a piecewise combination of the ground state and excited state BOPES. The respective Hartree-only KS potentials do not reproduce this behavior.
2.5
Dynamics in presence of an external field
In this section we consider ionization of the 1D H+
2 molecular ion by a slowly varying external electric field, starting from the ground state. The interaction of molecule with the external field is given by the potential ve(x, t) and Equation (2.3.16) is replaced by,
2.5 Dynamics in presence of an external field 23 0.0 0.5 1.0 1.5 Γ (a.u.) Nuclear density -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.)
Ground state BOPES
0.0 0.5 1.0 1.5 Γ (a.u.) -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.)
Excited state BOPES
0.0 0.5 1.0 1.5 Γ (a.u.) -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.) KS potential 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
0.0 0.5 1.0 1.5 Γ (a.u.) 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
-1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.) Hartree-only KS potential (a) (b) (c) (d) (e) (g) (f) (h)
Figure 2.4: Nuclear density, KS potential, and Hartree-only potential for the su-perposition of ground-state and excited-state at t=2 fs [panels (a) and (e)], t=4 fs [panels (b) and (f)], t=6 fs [panels (c) and (g)], and t=8 fs [panels (d) and (h)] obtained for the 1D H+
2 molecular ion starting from the initial state in Eq. (2.4.4).
The direction of the field is chosen such that the electron can escape towards x = ∞. We subject the model system to an electric field E(t) that is increased slowly with time using a sin2 envelope is written as E(t) = E
0 sin2 π 2 t ton , with ton = 17 fs. The interaction with the field is truncated at a distance of xt= 50 a.u. from the nucleus to avoid numerical problems with the strongly accelerated electrons. The potential due to the external field can then be written as
ve(x, t) = qexE(t) for x < xt, (2.5.2)
= qextE(t) for x > xt, (2.5.3)
where qe=2+2m2mn+1n denotes the reduced charge.
In Figs. 2.5 and 2.6 we plot the nuclear density, KS potential, and Hartree-only KS potential for E0 = 0.01689 a.u. and E0 = 0.05340 a.u. respectively. We have plotted the nuclear density and potentials only up to 8 fs as for later times, the spread of the dissociated wave packet at large R makes the KS potentials completely illegible. The initial state is the exact ground-state wave function. Initially the model system approximately remains in the ground state due to the slow switching of the electric field. The nuclear ground state density is still concentrated around the equilibrium distance and the resulting KS potential is a binding potential. This can be seen in the plots of nuclear density (Figs. 2.5(a) and 2.6(a)) and the KS potential (Figs. 2.5(e) and 2.6(e)) for time 2 fs. The KS potential in Figs. 2.5(e) and 2.6(e) is in
10-15 10-12 10-9 10-6 10-3 100 Γ (a.u.) Nuclear density -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.)
Ground state BOPES
10-15 10-12 10-9 10-6 10-3 100 Γ (a.u.) -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.) KS potential 10-15 10-12 10-9 10-6 10-3 100 Γ (a.u.) -1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.) Hartree-only KS potential 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
10-15 10-12 10-9 10-6 10-3 100 Γ (a.u.) 1 2 3 4 5 6 7 8 9 10 11 12
Internuclear distance (a.u.)
-1.0 -0.8 -0.6 -0.4 -0.2 Energy (a.u.)
Nuclear Coulomb repulsion (a) (b) (c) (d) (e) (g) (f) (h)
Figure 2.5: Nuclear density, KS potential, and Hartree-only KS potential plotted for E0 = 0.01689 a.u. at t=2 fs [panels (a) and (e)], t=4 fs [panels (b) and (f)], t=6 fs [panels (c) and (g)], and t=8 fs [panels (d) and (h)] obtained for the 1D H+ 2 molecular ion. The initial state is the ground state.
agreement with the ground-state BOPES (plotted in black-dashed line) while the Hartree-only KS potential has an overbinding character. After ionization, the two bare proton are repelled away from each other due to the Coulomb repulsion force. This is visible in Figs. 2.5(b), 2.6(b), 2.5(c), 2.6(c), 2.5(d), and 2.6(d) at times 4 fs, 6 fs, and 8 fs where part of the nuclear density remains bound in the region of 1-5 a.u. and a small part of density moving away towards large internuclear distance. This phenomena, known as Coulomb explosion, is reflected in the KS potential at large R in Figs. 2.5(f), 2.6(f), 2.5(g), 2.6(g), and 2.5(h), 2.6(h) at times 4 fs, 6 fs, and 8 fs. From R = 5 a.u. onwards, the KS potential matches approximately with the nuclear Coulomb repulsion.
2.6
Discussion
In this chapter we have given a brief introduction to the MCDFT and TD-MCKS equations. Similar to TDDFT, TD-MCDFT also relies on the quality of approximation for the exchange-correlation (xc) potentials. In order to get insight into the exact KS potentials, in Section 2.3 , we discussed an algorithm that allows us to construct the exact time-dependent KS potentials for a small model system with coupled electronic and nuclear motion. With the help of the one-dimensional
