卡氏座標下三維時變薛丁格方程式的平行化程式發展及其於雷射與分子間交互作用之應用研究

(1)

國

立

交

通

大

學

機械工程學系

博士論文

卡氏座標下三維時變薛丁格方程式的平行化程式

發展及其於雷射與分子間交互作用應用之研究

Parallel Solver for Three-dimensional Cartesian-grid Based

Time-Dependent Schrödinger Equation and Its Applications in

Laser-Molecule Interaction Study with Single-Action-electron Assumptio

n

研究生：李允民

指導教授：吳宗信博士

江進福博士

(2)

卡氏座標下三維時變薛丁格方程式的平行化程式發展及其於

雷射與分子間交互作用之應用研究

Parallel Solver for Three-dimensional Cartesian-grid Based

Time-Dependent Schrödinger Equation and Its Applications in

Laser-Molecule Interaction Study with Single-Action-electron Assumption

研究生：李允民 Student：

Yun-Min Lee

指導教授：吳宗信博士 Advisor：

Dr. Jong-Shinn Wu

江進福博士

Dr. Tsin-Fu Jiang

國立交通大學

機械工程學系

博士論文

A Thesis

Submitted to Department of Mechanical Engineering College of Engineering National Chiao Tung University

in partial Fulfillment of the Requirements

for the Degree of Doctor of Philosophy in Mechanical Engineering June 2009

Hsinchu, Taiwan

(3)

卡氏座標下三維時變薛丁格方程式的平行化程式發展及其於

雷射與分子間交互作用之應用研究

研究生：李允民 Student： Yun-Min Lee

指導教授：吳宗信 Advisor： Jong-Shinn Wu

江進福

Tsin-Fu Jiang

國立交通大學機械工程學系學系博士班 (熱流組)

中文摘要

我們以有限體積法來離散化一個在卡氏座標中的三維時變薛丁格方程式 (TDSE)，假設分子的原子核不移動且只有一電子有作用(SAE)，並用一個交錯時間法(stagger-time)來處理波方程式的實部與虛部在時間上演進，以平行計算的程式來分析此三維時變薛丁格方程式，此程式利用多層圖型切割法來作區域切割(domain decomposition)，以此分配不同處理器(processor)的計算區域。此程式以一個 H2+系統來作驗證，在沒有雷射(laser)的交互作用下，總電子機率和總能量的守恆。與其它二維薛丁格方程式的計算比較，在雷射的作用下離子化(ionization rates)的機率比較。此程式的平行化效率在使用 128 顆處理器下可以達到 75%。最後以 H2+系統，原子核距離 9 au，電子在不同雷射入射角度(χ= 0° 和 90°)下的機率隨時間的分佈，及原子核距離 2 au，和諧光譜(HHG)在不同雷射入射角度(χ= 0°, 30°, 60° and

(4)

90°)，以及不同雷射入射角度對 N2, O2及 CO2分子離子化機率的影響並和實驗資料

(5)

Parallel Solver for Three-dimensional Cartesian-grid Based

Time-Dependent Schrödinger Equation and Its Applications in

Laser-Molecule Interaction Study with Single-Active-electron Assumption

Student: Yun-Min Lee

Advisor: Dr. Jong-Shinn Wu

Dr. Tsin-Fu Jiang

Department of Mechanical Engineering

National Chiao Tung University

Abstract

A parallelized three-dimensional Cartesian-grid based time-dependent Schrödinger equation (TDSE) solver for molecules with single active electron assumption, assuming freezing the motion of nucleus is presented in this thesis. An explicit stagger-time algorithm is employed for time integration of the TDSE, in which the real and imaginary parts of the wave function are defined at alternative times, while a cell-centered finite-volume method is utilized for spatial discretization of the TDSE on Cartesian grids. The TDSE solver is then parallelized using domain decomposition method on distributed memory machines by applying a multi-level graph-partitioning technique. The solver is validated, using a H2+ molecule system, both by observing total

(6)

an aligned incident laser pulse. Parallel efficiency of this TDSE solver is presented and discussed, in which the parallel efficiency can be as high as 75% using 128 processors. Finally, examples of temporal evolution of probability distribution of laser incidence onto a H2+ molecule at inter-nuclei distance of 9 a.u. (χ= 0° and 90°) and spectral

intensities of harmonic generation at inter-nuclei distance of 2 a.u. (χ= 0°, 30°, 60° and 90°) and the angle effect of laser incidence on ionization rate of N2, O2 and CO2

(7)

誌謝

首先最要感謝的人是吳宗信老師，江進福老師，在我博士的求學中給我指導與幫助。學如逆水行舟，而老師則像是裝在舟上的電動馬達，帶我一路勇往直前奮力不倦，不僅於課業上給予指導，更教曉我做人做事的態度，如今輕舟已過萬重山，老師對我的照顧及付出，點滴在心，無限感激！另外要特別同學周欣芸和洪捷粲在課業及生活上的幫助，和實驗室學弟們的大力協助。最後特別感謝我的家人，如不是你們一直以來的包容及支持，我也不可能有今天的小小成就。當然還要感謝佩錡對我多年來的打氣與陪伴，妳就像是我手心的太陽，總在我忙碌無助時，適時給予我安定的力量，支持我一路走下去，認識妳是我最大的福氣，謝謝妳！最後，要感謝的人實在太多，如有被遺忘的朋友在此亦一併感謝，感謝大家於這多年來對我的照顧。在此也祝福所有還在學的朋友們都能夠在求學的路上綻放出各式各樣的光芒！ 2009 年 6 月 30 日

(8)

中文摘要

中文摘要... i

Abstract... iii

誌謝

誌謝... v

Table of Contents ... vi

List of Tables ... viii

List of Figures... ix

Nomenclature ... xvii

Chapter 1. Introduction ... 1

1.1 Background and Motivation ... 1

1.2 Literature Survey... 2

1.2.1 The hydrogen molecular ion ... 3

1.2.2 Multi-electron and multi-nuclear molecules... 6

1.3 Objectives ... 8

Chapter 2. Fundamentals of Quantum Physics ... 10

2.1 Time-dependent Schrödinger equation... 10

2.2 Time independent Schrödinger equation... 11

2.3 Expectation values (measurable physical quantities) ... 12

2.4 The Born-Oppenheimer approximation (BOA) ... 13

2.5 Single-Active-Electron assumption (SAE) ... 13

2.6 Laser field in the Time-dependent Schrödinger equation ... 16

Chapter 3. Numerical Methods ... 17

3.1 Three-dimensional Time independent Schrödinger Equation ... 17

3.2 Three-dimensional Time Dependent Schrödinger Equation ... 17

3.3 Effective potential of core electrons in molecule ... 18

3.4 Discretization of the Three-dimensional Time independent Schrödinger Equation ... 19

(9)

3.5 Discretization of the Three-dimensional Time Dependent Schrödinger Equation

... 23

3.5. Parallel Implementation of TDSE Solver ... 26

Chapter 4. Code Verification... 28

4.1. Ground state energy of H atom and H2+ ion molecule... 28

4.2. Energy Conservation of H2+ Without Laser Incidence ... 29

4.3. Ionization Rates of H2+ With Laser Incidence Along the H-H Axis... 29

Chapter 5. Results and applications ... 31

5.1 Parallel Performance of the 3D TDSE Solver... 32

5.2 Results of single electron and multi-nuclear molecular system... 33

5.2.1 Instantaneous Electron Probability Distribution of H2+ molecule ... 33

5.2.2 Spectra of Harmonic Generations of H2+ molecule ... 34

5.3 Multi-electron and multi-nuclear molecular system ... 35

5.3.1 The oriental effect of laser incidence on ionization rate of multi-electron and two-nuclear molecular system: N2... 36

5.3.2 The oriental effect of laser incidence on ionization rate of multi-electron and two-nuclear molecular system: O2... 37

5.3.3 The oriental effect of laser incidence on ionization rate of multi -electron and three-nuclear molecular system CO2... 39

Chapter 6. Conclusion ... 42

6.1 Summary of the major findings ... 42

6.2 Recommendation of future studies ... 42

(10)

List of Tables

Table 1. Ground state eigen-energy of H2+ molecule in different inter-nuclear

distance. The eigen-energy are calculated by J-D solver and compared to exact value listed in reference 76. These calculated eigen-energy are different from exact values [76] by less than 0.2%. ... 50 Table 2. Timing breakdown of a typical parallel simulation for 400 timesteps

(2,545,548 cells, time step size=0.01 a.u., laser intensity = 1014 W/cm2, and wave length =1064nm.) ... 51 Table 3. Simulation conditions for weaker laser incidence onto a H2+ molecule at

χ=0° and 90°... 52 Table 4. Simulation conditions for stronger laser incidence onto a H2

+

molecule at different angles of incidence (χ=0°, 30°, 60° and 90°)... 53 Table 5. Fitting parameters of Yukawa like soft-coulomb potential for N2 molecule. 54

Table 6. Fitting parameters of Yukawa like soft-coulomb potential for O2 molecule. 55

Table 7. Fitting parameters of Yukawa like soft-coulomb potential for CO2 molecule.

... 56 Table 8. Simulation conditions for laser incidence onto the N2, O2 and CO2 molecule

(11)

List of Figures

Figure 1. Sketch of the typical finite-element grid system projected in two-dimensional space. ... 58 Figure 2. Local number and coordinates of the finite-element grid system... 59 Figure 3. Sketch of the typical finite-volume grid system projected in

two-dimensional space. ... 60 Figure 4. (a)Typical grid system for 3D TDSE simulation (a slice through the

midplane). ... 61 Figure 4. (b) Typical slice of domain decomposition through midplane (16

processors). ... 62 Figure 4. (c) Typical surface domain decomposition (16 processors)... 63 Figure 5. Total electron probability and total energy variance of H2+ molecule

without laser incident, the internuclear distance is 9 au. ... 64 Figure 6. Applied electric field to the H2+ molecule along the H-H axis as a function

of time. Laser intensity = 1014 W/cm2, and wave length =1064nm... 65 Figure 7. Comparison of the ionization rates as a function of inter-nucleous distance,

obtained by the present parallelized 3D TDSE solver and previous 2D-axisymmetric TDSE solver for an aligned sub-femto-second linearly

(12)

polarized laser pulse interacting with a H2+ molecule (power

intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles). ... 66 Figure 8. Parallel efficiency of the present parallelized 3D TDSE solver as a function

of the number of processors. Case 1: 2.54M cell grid and case2: 14.8M cell grid. Speedup is normalized by 2 processor data... 67 Figure 9. Typical snapshots of the electron probability distribution over the

axisymmetric plane for a normally incident sub-femto-second linearly polarized laser pulse interacting when t=0 a.u. (0.00 cycle) with a H2+

molecule (R=9) (power intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles, angle of incidence=0°). ... 68 Figure 10. Typical snapshots of the electron probability distribution over the

axisymmetric plane for a normally incident sub-femto-second linearly polarized laser pulse interacting when t=810 a.u. (5.53 cycle) with a H2

(13)

(14)

polarized laser pulse interacting when t=0 a.u. (0.00 cycle) with a H2+

molecule (R=9) (power intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles, angle of incidence=90°) ... 75 Figure 17. Typical snapshots of the electron probability distribution over the

molecule (R=9) (power intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles, angle of incidence=90°) ... 76 Figure 18. Typical snapshots of the electron probability distribution over the

molecule (R=9) (power intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles, angle of incidence=90°). ... 77

(15)

Figure 19. Typical snapshots of the electron probability distribution over the axisymmetric plane for a normally incident sub-femto-second linearly polarized laser pulse interacting when t=2250 a.u. (15.36 cycle) with a H2+

molecule (R=9) (power intensity=1014W/cm2, wave length=1064nm, pulse duration=25 cycles, angle of incidence=90°). ... 79 Figure 21. Harmonic spectra of H2+ for different orientation angles: χ=0~90°. Laser

intensity = 5*1014 W/cm2, and wave length = 800nm. Internuclear distance =2.0 a.u. ... 80 Figure 22. Slice contour topology of Yukawa like soft-coulomb potential for N2

molecule on x=0 plane. The maximum value of the potential is about 4.16... 81 Figure 23. The 3D iso-surface contour of initial wave function for laser-N2 molecule

interaction. The black dots are the positions of the nuclear of N2 molecule.

(16)

Figure 24. The ionization yield to laser incidence angle in every 15∘for N2 molecule.

Laser intensity = 1.5*1014 W/cm2, and wave length = 820nm. Internuclear distance =2.075 a.u... 83 Figure 25. Ionization signal S(α) converted from measured ionization yield as a

function of the angleα between the polarization axes of the aligning and the ionizing beams for N2, and The peak laser intensities is1.5*1014 W/cm2. This

data is from reference [65]. Red solid line and orange dash dotted line are converted form experimental data by different method, green dotted line is from MO-ADK calculation... 84 Figure 26 Electron probability density distributions of N2 molecule under different

laser incidence angle on x=0 plane at t=5 optical cycle (~13.78 fs). (a) initial (b)laser incidence angle of χ=0∘, 30∘, 60∘and 90∘... 85 Figure 27. The slice contour topology of Yukawa like soft-potential for O2 molecule

on x=0 plane. The maximum value of the potential is about 7.4. ... 86 Figure 28. The 3D iso-surface contour of initial wave function for laser-O2 molecule

interaction. The black dots are the positions of the nuclear of O2 molecule.

This is a πg type orbital. The orbital is symmetric to molecule center and

(17)

Figure 29. The ionization yield to laser incidence angle in every 15∘for O2 molecule.

Laser intensity = 1.3*1014 W/cm2, and wave length = 820nm. Internuclear distance =2.28 a.u... 88 Figure 30. Ionization signal S(α) converted from measured ionization yield as a

function of the angleα between the polarization axes of the aligning and the ionizing beams for O2, and The peak laser intensities is1.3*1014 W/cm2. This

data is from reference [65], Red solid line and orange dash dotted line are converted form experimental data by different method, green dotted line is from MO-ADK calculation... 89 Figure 31. 3D iso-surface contour of electron probability density distributions of O2

molecule under different laser incidence angle at t=10 optical cycle (~27.57 fs). (a) initial (b)laser incidence angle of χ=0∘, 30∘, 60∘and 90∘... 90 Figure 32. Slice contour topology of Yukawa like soft-coulomb potential for CO2

molecule on x=0 plane. The maximum value of the potential is about 5.66... 91 Figure 33. The 3D iso-surface contour of initial wave function for laser-CO2

molecule interaction. The black dots are the positions of the nuclear of CO2

molecule. This is a πg type orbital. The orbital is symmetric to molecule

(18)

Figure 34. The ionization yield to laser incidence angle in every 15∘for CO2

molecule. Laser intensity = 1.3*1014 W/cm2, and wave length = 820nm. Internuclear distance =2.28 a.u. ... 93 Figure 35. Ionization signal S(α) converted from measured ionization yield as a

function of the angleα between the polarization axes of the aligning and the ionizing beams for CO2, and The peak laser intensities is1.3*1014 W/cm2.

This data is from reference [65]. Red solid line and orange dash dotted line are converted form experimental data by different method, green dotted line is from MO-ADK calculation. ... 94 Figure 36. 3D iso-surface contour of electron probability density distributions of CO2

molecule under different laser incidence angle at t=10 optical cycle (~27.57 fs). (a) initial (b)laser incidence angle of χ=0∘, 30∘, 60∘and 90∘... 95

(19)

Nomenclature

r Position vector of electron.

R

Position vector of nucleus

H Hamiltonian operator of Schrödinger equation

( )

r

ψ

Wave function of Schrödinger equation 2

r

∇ Laplacian operator for electron 2

R

∇ Laplacian operator for nucleus

α Soft coulomb potential parameter

β Yukawa potential parameter eff

Z Effective charge of ion

E Eigen-energy

( )

N r

Shape function of finite element grid

( )

yukawa

V r

Yukawa like soft coulomb potential

Ω A small element

∂Ω Surface of a small element χ Laser incidence angle

(20)

Chapter 1. Introduction

1.1 Background and Motivation

In the past decade, the laser technology had been greatly improved. Technology of production of ultra-short laser pluses has been extensively developed and ultra-short high power laser pluses with durations about few optical cycles are now available to extend the science of atomic, molecular and optical physics [1]. In these cases, this had led to discovery of many new nonlinear nonperturbative optical phenomena and processes such as above threshold ionization (ATI), tunnelling ionization, and high-order harmonic generation (HHG). The frontier applications of these processes can provide a tool for molecule imaging [2-6], the light source in the region of X-ray to ultraviolet [7, 8]. From a theoretical viewpoint, such studies are extremely complex in the strong-field regime and have been of continuous interest for nearly two decades. In general, only results based on approximate theories such as the molecular strong-field approximation [9, 10] and tunnelling [11] models have been applied to calculate effects related to molecular orientation with respect to the light polarization vector. Such approximate theories are, however, often gauge dependent [10, 12] and limited in their applicability to describe complex processes. Despite several years of study, the interaction of intense laser pulses with atoms and molecules remains a very attractive

(21)

problem. Many problems in atomic-scale physics rely on the development of efficient numerical methods for the solution of the time-dependent Schrödinger equation (TDSE, introduced in chapter 2). Theoretical modelling of these phenomena is important to gain a fundamental understanding of atomic and optical processes at the microscopic level and to provide a scientific basis for the design and development of molecule imaging, light source in the region of X-ray to ultraviolet, which demands theoretical and large-scale computational modelling of all of these phenomena.

1.2 Literature Survey

There are numerous papers about laser-molecule interactions; we will arrange the papers according to the molecular system for systematic introduction.

First is the one electron and two-nucleus-centers molecule; the hydrogen molecular ion. Hydrogen molecular ion is the only and the simplest molecule of one electron molecules. Because, the hydrogen molecular ion can be solved theoretically with Born-Oppenheimer approximation (BOA, introduced in chapter 2), that is why hydrogen molecular ion is two fundamental and prototypical systems which can be used to understand and extend these fields of physics.

(22)

molecular systems are much difficult than one electron systems on the view point of computation. For practically implementation in theoretical analysis is always with approximations and modelling.

1.2.1 The hydrogen molecular ion

H2+ is the simplest molecule. It is usually used as a prototype of molecular

physics like the hydrogen atom in atomic physics. The quantum dynamics of H2+ and

other small molecules under the interaction of strong laser pulse has been a subject of continuous interest for more than two decades. Interesting phenomena such as bond-softening, bond-hardening, above threshold dissociation and Coulomb explosion were studied by Posthumus [13]. Especially, a recent experiment clearly showed that the ionization rates of H2+ at some characteristic bond lengths are greatly enhanced with

laser pulse of ∼1014

w/cm2 and duration of several hundreds femtoseconds. When the molecule is ionized, the two positively charged bare ions repel each other with strong Coulomb repulsion due to the small separation between them. All of a sudden the two ions fly apart rapidly [14]. This dramatic phenomenon was named either Coulomb explosion or charge resonant enhanced ionization.

From the computational viewpoint, the H2 +

quantum system is not simple at all. The time evolution of H2+ under laser field is in general 9 degrees of freedom in space

(23)

plus the time variable. This kind of time-dependent Schrödinger equation (TDSE) is unlikely to be solved even using the most advanced numerical scheme and computer, and the TDSE is very often be simplified by the Born-Oppenheimer approximation (BOA). This approximation will remove the dynamic effect of the two nuclei and reduce the dimension of the TDSE of H2+ from 9 to 3 degrees of freedom in space. All

papers as below mentioned are using with the BOA unless specified.

The H2+ has two fixed nuclei and a fast moving electron. With a linearly polarized

electric field along the molecular axis, the system is cylindrical symmetric in spatial coordinates in addition to the time variable. The time-dependent two-dimensional TDSE has been used often because it is computationally less demanding [e.g., 15-19]. However, to align the molecule totally oriented with the laser field is still very difficult (or impossible) experimentally and also the laser field may not be linearly polarized, which makes the time-dependent three-dimensional computation of TDSE strongly required in essence. Unfortunately, this kind of calculation is still very limited and still at its infancy due to the very high computational cost, even with the rapid advancement of computer technology. Until very recently, there have been three types of numerical methods used in solving the 3D TDSE for the study of laser-molecular interaction, which include: 1) First transformed the 3D TDSE from Cartesian coordinates into

(24)

technique with high-order Taylor’s series expansion for time propagation [20-21; Bandrauk’s group, 22-26] or Legendre pseudospectral discretization [Chu’s group, 20, 27]. 2) Solved the 3D TDSE in spherical coordinates using spherical harmonics basis expansion technique with split operator for time propagation [Madsen’s group, 28-32]. 3) Solved the 3D TDSE directly in Cartesian coordinates using parallel finite-difference method (FDM) [33] or finite-element method (FEM) [Collion’s group, 33-37; Bandrauk’s group, 38] with split operator or high-order Taylor’s series expansion for time propagation on distributed-memory machines. The above techniques have been used successfully to study the variation of ionization rate, higher-order harmonic generation (HHG) and electron probability distribution due to changes of internuclear distance, angle of laser incidence and laser intensity. For 3D problems, it is very time-consuming to apply methods of types 1 and 2 on a single-processor machine, however it is found difficult to parallelize the code on distributed-memory machines. In contrast, for FDM and FEM (type 3), it is relatively easy to parallelize the simulation code. However, it is well-known that with FDM it is rather clumsy to manipulate the grid distribution. On the other hand, it often requires matrix inversion using FEM [38], if non-orthogonal type basis functions are used, even with explicit time-marching scheme, which is very time-consuming and memory-demanding in practice. If orthogonal polynomials are used as the basis functions in FEM [33-37], then matrix

(25)

inversion becomes unnecessary. Even so, the number of non-zero entries in the matrix (or memory storage) resulting from the spatial discretization and number of machine operations per time step (or computational time) using FEM are both large, as compared to the cell-centered finite-volume method, which will be introduced later in Section 2. Thus, a general numerical method without the above-mentioned shortcomings, which can be used to solve the truly 3D TDSE for general diatomic molecule (and possibly beyond) under strong laser pulse, is still strongly desired in the physics community.

1.2.2 Multi-electron and multi-nuclear molecules

Even if under Born-Oppenheimer approximation, when a molecular system contains more than one electron, the numerical difficulty is dramatically increasing due to the non-linearity of Schrödinger equation and enormously computational cost. The theoretical modelling for the multi-electron molecule system is unavoidable.

There are some theoretical modelling analyses for a laser-molecule interaction system. The Ammosov-Delone-Krainov (ADK) theory [39] is based on the ionization rate of a hydrogen like atom in a static field, with modifications introduced. The ADK model can be used for atomic system and valid only for high laser intensity and low optical frequency [40]. The molecular-ADK (MO-ADK) theory [41] uses an ab-initio

(26)

applications on molecular system such as N2 and O2 [42]. For low laser intensity and

high optical frequency region, the lowest order perturbation theory (LOPT) can describe the ionization process. However, LOPT calculations are computationally demanding and thus only few systematic and correlated calculations are been performed [43]. Another popular model is strong field approximation (SFA). The traditional strong-field approximation is described as a transition from a field-free initial state to a final Volkov state ignoring thus the Coulomb interaction of the ejected electron with the remaining ion. Molecular effects have recently also been incorporated [44] into the model. And, there is some dispute about the use of length or velocity gauge and the applicability of a Coulomb-correction factor [9]. The above models have some intrinsically limitations. Such as, the ADK rate is optical frequency independent. The SFA does not contain any information on the possible influence of resonances.

The Density functional theory (DFT) is the most common theory for multi-electron molecule system. In analogy the DFT, Time-dependent DFT (TDDFT) can deal with this multi-electron molecule system under strong laser field, but the accuracy of TDDFT is dependent on the choice of the time-dependent exchange-correction potential because of the adiabatic assumption of exchange-correction potential [45]. Full orbital TDDFT computation is also very time-consuming as electron number is large. Due to above reasons, there are very few full orbital TDDFT studies about multi-electron atom or

(27)

molecule system under strong laser field. Chu[46, 47] performed some TDDFT simulations for small diatom molecule. Although Chu’s method is good but it can only be applied on diatom molecule. Bauer and Koval [48] provided a TDDFT solver for intense laser–atom interaction, but Bauer and Koval’s program is only for atoms.

1.3 Objectives

After studied these papers about laser-molecule interacting phenomena, we found that there were very few studies for molecule with more than two nuclei because of the special spatial coordinates were used or basis overlap integration were time consuming. The general purpose parallel three dimension TDSE solver with the capability to deal with multi-electron and multi-nuclear molecular system was rarely to find. The objectives of this thesis are as followed.

(1) To develop a parallelized 3D TDSE solver using Cartesian-grid based finite-volume method (FVM) and demonstrate and verify its capability by simulating the physics of H2+ under strong laser pulse. Then to

compare the data with others published data before.

(2) To make a model potential for multi-electron molecules. Then to couple with the developed parallel 3D TDSE solver and extended to treat

(28)

general polyatomic molecule with single-active-electron (SAE) under strong laser pulse.

(3) To make numerical analysis on multi-electron and multi-nuclear molecular system with laser incidence. To compare the data with others published data before.

This thesis begins with a basic introduction of quantum physics, derivation of model potential for multi-electron, then a detailed description of the finite-volume method used for discretizing the 3D TDSE. The development and implementation of parallelized TDSE solver is then outlined and simulations of a laser pulse incident along the axis of H2+ is then carried out to verify the accuracy of the codes. These simulations

are compared to the results of other studies applying symmetric TDSE in the literature. Results of simulated instantaneous electron probability and spectrum of harmonic generations are then presented, in which the laser is incident at different directions with respect to H2+ axis, to demonstrate the capability of the new TDSE solver. And we show

the applications of the new TDSE solver on multi-electron and multi-nuclear molecular system as: N2, O2 and CO2. Then make a comparison with other experimental and

(29)

Chapter 2. Fundamentals of Quantum Physics

Quantum theory grew out of an interplay of ground-breading experiments and radical theoretical proposals that were not based on accepted classical physics. Now it is an important theory for nano-scale and ultrafast physics studied. Here we make a quickly introduction to the basic equation (Schrödinger equation) in quantum physics, and the derivation of model potential Schrödinger equation.

2.1 Time-dependent Schrödinger equation

Schrödinger equation is fundamental equation of quantum physics; it is also called Schrödinger wave equation. Schrödinger equation provides the description of motion and interaction of particles at the small scales where the discrete nature of matter becomes important. In this scale, classical physics is fundamental break because of the continuous assumption of matter broken down.

Time dependent Schrödinger equation for n electrons and N nuclei in atomic unit can be expressed as:

(

)

(

) (

)

1 1 1 1 1 1 ,..., , ,..., , ,..., , ,..., , ,..., , ,..., , n N n N n N r r R R t i H r r R R t r r R R t t ψ ψ ∂ = ∂ (2-1)

(30)

i i i i

r =x i +y j+x k

is the position vector of i-th electron and R_j =x i_j +y j_j +x k_j

is the position vector of j-th nuclear. The HamiltonianH r

(

₁,..., ,r R_n ₁,...,R t_N,

)

is expressed as the sum of kinetic energy operator 2 2

1 1 1 1 2 i 2 j n N r R i= j= mi

−

∑

∇ −

∑

∇ for electrons and nuclei and coulomb potential of electron to electron

' ' 1 n i i≠ ri−ri

∑

, nuclear to nuclear ' ' ' N j j j j j j Z Z R R ≠ −

∑

and electron to nuclear

1 1 n N j i j _i _j Z r R = = − −

∑∑

. Multi-particle time dependent

Schrödinger equation is a non-linear equation and the dimension of multi-particle time dependent Schrödinger equation is equal to3

(

n+N

)

+ . Because of the high dimension 1 of the equation when

(

n+N

)

> , it is not easy to do numerical or theoretical analysis 1 without any assumption.

2.2 Time independent Schrödinger equation

If the Hamiltonian operator of Schrödinger equation is time independent (H r

(

₁,..., ,r R_n ₁,...,R t_N,

)

=H r

(

₁,..., ,r R_n ₁,...,R_N

)

). The Schrödinger equation can be solved by variables separation. The solution of the Time independent Schrödinger equation can be expressed as: ψ

(

r₁,..., ,r R_n ₁,...,R t_N,

) (

₌ϕ r₁,..., ,r R_n ₁,...,R_N

)

_⋅e−iE t⋅ _. And E in above equation is a real number. Then, we can rewrite eq(2-1) to:

(

1,..., ,n 1,..., N

)

(

1,..., ,n 1,..., N

) (

1,..., ,n 1,..., N

)

E⋅ϕ r r R R =H r r R R ⋅ϕ r r R R

(2-2) Above equation is frequently called stationary Schrödinger equation or time independent Schrödinger equation. The stationary Schrödinger equation is an

(31)

eigenvalue equation, and the corresponding eigen-energy and eigen-function are the system energy and wave function of the Schrödinger equation. And the wave function ( ϕ

(

r₁,..., ,r R_n ₁,...,R_N

)

) is often used as the initial wave function (ψ

(

r₁,..., ,r R_n ₁,...,R_N, 0

) (

=ϕ r₁,..., ,r R_n ₁,...,R_N

)

) for a time-dependent computation.

2.3 Expectation values (measurable physical quantities)

Wave function is not a measurable physical quantity. The measurable quantities can be obtained by defining the quantity operator f r

( )

. We consider a one-electron Schrödinger equation for example; the wave function is treated as a probability amplitude function. We can define a quantity operator f r

( )

, and the expectation value of the quantity operator can be obtained by integrating all coordinate space:

( ) ( ) ( )

3

F ψ r f r ψ r ∞ψ∗ r f r ψ r dr

−∞

= =

∫

⋅ ⋅ ⋅

. And the F is the measurable quantity. Some quantity operators are listed:

1. f r

( )

=1

, F is equal to total electron umber. 2. f r

( )

=H r

( )

, F is equal to total system energy 3.

( )

1 2

2

f r =− ∇

(32)

2.4 The Born-Oppenheimer approximation (BOA)

The calculations are often simplified with the Born-Oppenheimer approximation (BOA) where the nuclei are fixed during the light interaction, because the motions of nuclei are much slower than the electrons. It is known that the BOA is applicable when the optical period is smaller than the vibrational period. The equation (2-1) can be simplified to

(

)

(

) (

)

1 1 1 ,..., ,..., , ,..., , n n n r r t i H r r t r r t t ψ ψ ∂ = ∂ (2-3)

by neglecting the movement of nuclei.

2.5 Single-Active-Electron assumption (SAE)

It is not easy and practical to directly solve the Schrödinger equation of multi-electron molecule because the electron-electron interaction and high dimension of Schrödinger equation. In order to simplify the multi-electron molecule problem, we use the single-active-electron (SAE) assumption [49]. In the SAE approximation, the time-dependent Schrödinger equation for a single electron moving in the effective field generated by the nuclei and all the other electrons is solved numerically. SAE models where one reduces the dimensionality of the multi-electron problem by freezing the core electrons have proven to be very useful in cases where multiple electronic excitations are insignificant, and the SAE approximation is probably the most widely used

(33)

approach when studying phenomena such as single ionization, above-threshold

ionization (ATI) and high-harmonic generation (HHG) [50]. The key of SAE assumption is how to make a good effective potential, and the soft-coulomb potential method is the most common one.

The basic ideal of soft-coulomb potential is to remove the singularity of coulomb potential and represent the screened electrons of core electron. The soft-coulomb potential can be written as:

eff Z r R α − − + (2-4) eff

Z and α are a tuneable parameter. Z_eff and α are usually modified to fit the ground state energy to atom or molecule system and the asymptotic behavior of coulomb potential at r→ ∞

. The electron-nuclear coulomb potential and the electron-electron coulomb potential are replaced by the soft coulomb potential in the time dependent Schrödinger equation. Equation (2-4) is a simple and common approximation for soft coulomb potential, but sometimes the two tuneable parameters (Z_eff and

α

) can not satisfy both the ground state energy and the asymptotic behavior of coulomb potential at the same time. In order to overcome this problem, we further introduce the Yukawa potentials. Yukawa potentials had been studied for the bound states energy calculations [51]; it can provide a good model for screening effect of core

(34)

r core Z e r β − − ⋅ (2-5)

The Zcore in equation (2-5) is the effective charge of core nuclear. If we include the

Yukawa potentials to soft-coulomb potential, we can rewrite equation (2-4) to:

( )

_core r R yukawa Z V r e r R β α − − − = ⋅ − + (2-6a) and

( )

eff

(

1 r R

)

yukawa core Z V r Z e r R β α − − − = ⋅ + ⋅ − + (2-6b)

Now, there are four parameters (Z_eff, Z_core, α and β) for fitting the ground state energy and the asymptotic behavior of coulomb potential at r→ ∞

. Equation (2-6a) is used for screening effect of neutral atoms which interacts with electron in the mother molecule and the asymptotic behavior of coulomb potential at r→ ∞

will go to zero. And equation (2-6b) used for screening effect of ion with positive charge Z_eff in the

mother molecule and the asymptotic behavior of coulomb potential at r→ ∞

will

become to Zeff

r

−

. Theoretical analysis of one dimensional and two dimensional soft

coulomb potential models had been used for helium atom [52] and hydrogen molecule [53], and the three dimension model had been used for characteristic analysis of coulomb singularity in HHG [54].

(35)

Another type for the effective potential that is used for SAE model is the Hatree-Fork potential. The Hatree-Fork potential is calculated by time independent Hatree-Fork calculation. And the directly use it in the time dependent simulation. The potential is invariant during the time propagation. The Hatree-Fork potential for two electrons molecular system can be written as follow.

( )

2 3 ' ' ' i HF r V r dr r r φ = −

∫∫∫

(2-7)

To our knowledge, the SAE model with Hatree-Fork potential was used only for two electrons system. The Hatree-Fork potential had been used for H2 molecule [50, 55] and

for helium with basis expansion method [52].

2.6 Laser field in the Time-dependent Schrödinger equation

If we neglect the relativistic and the non-dipole effects in laser-molecule interactions, the laser field will therefore be described in dipole approximation by a spatial homogeneous vector electricA t

( )

or spatial homogeneous field E t

( )

with

( )

/

E t = −dA t dt . The magnetic component B t

( )

of the laser field, given by

( )

B t = ∇ ×A t

, will vanished by the dipole approximation. Then the operator of laser filed in the Time-dependent Schrödinger equation can be written as follow.

( )

,

E r t ⋅r

(36)

Chapter 3. Numerical Methods

3.1 Three-dimensional Time independent Schrödinger Equation

If we use the single-active-electron assumption (use the Yukawa like soft-coulomb potential), then the time independent Schrödinger Equation can be rewritten by combining the Yukawa like soft-coulomb potential equation (2-6), and expressed as follow:

( )

1 2

( ) ( )

2 yukawa E⋅ψ r =_− ∇ +V r _ψ r   (3-1) where r=xi +yj+xk

is the position vector of the electron, and 1 2 2

− ∇ is the energy operator. If the coefficient α and β in equation (2-6) are both equal to zero, then Yukawa like soft-coulomb potential equation will becomes a regular coulomb potential. And the parameter set (α=β = ) will be used in single electron molecule. 0

Equation (3-1) is an eigenvalue problem, and the solution of equation (3-1) is used as the initial wave function of a time-dependent evolution.

3.2 Three-dimensional Time Dependent Schrödinger Equation

We consider the linear time-dependent Schrödinger equation (TDSE) for a molecule with single-active-electron (SAE) under the incidence of a laser field. The TDSE can then be written as, with E r t

( )

,

(37)

(

)

(

)

2

(

)

(

)

, 1 , , , 2 r t i H r t V r t r t t ψ ψ ψ ∂   = = − ∇ +  ∂   (3-2)

and the Hamiltonian H is expressed as the sum of kinetic energy operator 1 2 2

− ∇ and potential energy operator V r t

(

,

)

. In addition, the potential energy operator under laser incidence can be generally expressed as

(

)

( )

1 1 , , N j _j V r t E r t r r R = = − + ⋅ −

∑

(3-3) where R_j

and N are the position vector of nucleus j and number of nuclei of the

molecule under consideration, respectively. And 1 1 N j= r−R_j

∑

is the electron-nuclear potential.

3.3 Effective potential of core electrons in molecule

In order to simplify multi-electron molecule Schrödinger equation, we assume that the distribution of core electron do not change during the laser field is active because of the core electron are tightly bound to nuclear in most molecules. We use SAE model and treat the effect of core electron and nuclei as effective potential. We model the effective potential by soft-coulomb potential. The α in equation (2-6a) and (2-6b) is determined by fitting the ground state energy of target molecule, and Z_eff and β are

(38)

potential, GAMESS[56] is a free package for ab-initio calculations; it can provide the molecular orbital information φ_i

( )

r and the index i is orbital number. We can

calculate the Hatree-Fork potential

( )

2 3 ' ' ' i HF r V r dr r r φ = −

∫∫∫

with this orbital

information.

Then we can modify the equation (3-3) to:

(

)

( )

1 , i i , core N r R i i _i _i Z V r t e E r t r r R β α − − = = − ⋅ + ⋅ − +

∑

(3-4a)

(

)

(

)

( )

1 , 1 i i , eff N r R core i i i _i _i Z V r t Z e E r t r r R β α − − = = − ⋅ + + ⋅ − +

∑

(3-4b)

for an Yukawa like soft-coulomb potential.

(

)

( ) ( )

1 1 , , N HF i _i V r t V r E r t r r R = = − + + ⋅ −

∑

(3-4c)

for Hatree-Fork potential.

It will be a general form for TDSE with single-active-electron assumption.

3.4 Discretization of the Three-dimensional Time independent Schrödinger Equation

For a generalized eigenvalue problem, the QZ [57] algorithm is the most popular one to solve it. The QZ algorithm code or program can be easily found in many

(39)

numerical mathematical libraries, such as LAPACK and IMSL. Unfortunately, the QZ algorithm is not efficient for large matrix problem. The computational cost is proportional to cube of matrix size. For a large matrix (>3000), another efficient algorithm is implemented in this problem. The Jocobi-Davidson algorithm [58] is an iterative algorithm, which solves the generalized eigenvalue problem iteratively. It was found that it is very efficient for large matrix problem for the first few selected eigenstates. The Jocobi-Davidson algorithm is only capable of dealing with the symmetric matrix problem. If we use the FVM to discretize equation (3-1), equation (3-1) in matrix form is not a symmetric matrix because of the potential term. In this thesis, we have applied finite-element-method to discretize equation (3-1) to form a symmetric matrix, which the Jocobi-Davidson algorithm can handle efficiently.

We use the finite-element method (FEM) to discretize equation (3-1) on structured non-uniform grid and the details are described as follows. A typical two-dimensional projected grid is shown in Figure 1. First, the wave function is approximated as

( )

i i

( )

i r c N r ψ

∑

(3-5) where the N r_i

( )

is the shape function or trial function for the ith grid point of

computational mesh, and c_i is the weightings for shape function N r_i

( )

. N r_i

( )

is equal to unity at grid point i and is equal to zero at any grid point j other than i,

(40)

( _i

( )

1 i N r =

∑

). N r_i

( )

is equal to 0, when r

is outside the corresponding element

which is under consideration.

There are eight grid points within a hexahedron. Numbering and coordinating are shown schematically in Figure 2. We employ the Lagrange polynomials which can satisfy the general properties of shape function as described earlier. The shape functions of a hexahedron can be written as follows:

( )

[

][

]

( )

[

][

]

( )

[

][

]

( )

[

][

]

( )

[

][

]

( )

[

][

]

( )

[

][

]

( )

[

][

]

1 1 1 1 2 0 1 1 3 0 0 1 4 1 0 1 5 1 1 0 6 0 1 0 7 0 0 0 8 1 0 0 1 1 1 1 1 1 1 1 N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV N r x x y y z z dV = − − − − = − − − = − − − − = − − − = − − − = − − − − = − − − = − − − − (3-6)

where dV is the volume of a hexahedron element. The parameters x₀, y₀, z₀, x₁,

1

y and z₁ are the local coordinates of a hexahedron element in Figure 1. Subscript of a shape functions denotes the local node number in a typical element as shown in Figure 2.

(41)

( )

1 2

( )

2 i i yukawa i i i i E⋅ c N r = − ∇ +_ V r _ c N r  

∑

(3-7)

Define the Galerkin weighted residual function as

( )

1 2

( )

( ; ) 2 i i i yukawa i i i i R r c =E⋅ c N r − − ∇ +_ V r _ c N r  

∑

(3-8)

Next, by applying the Galerkin weighted residual condition ( _{( ; )} _{( )} 3 ₀

i i

R r c N r dr =

∫

) on

equation (3-8) for each hexahedron element in the computational domain, we can obtain the following: .

( )

3

( )

1 2

( )

3 2 i i i j j i yukawa j j j j E N r c N r dr N r V r c N r dr Ω Ω   ⋅ = _− ∇ + _  

∑

∫

(3-9)

where ∂Ω is the boundary of integration volume _i Ω and _i N r_i

( )

=0

on ∂Ω . _i Further, by applying the divergence theorem, equation (3-9) is reduced to

( )

3 3 3 1 1 2 2 i i i i i j j j i j j i j j j r on j i yukawa j j j E N r c N r dr N r c N r N r c N r dr N r V r c N r dr Ω Ω ∂Ω Ω ⋅ = − ∇ + ∇ ∇ +

∑

∫

∑

∫

∑

∫

(3-10) since N r_i

( )

=0

on ∂Ω . The equation (3-10) then becomes: _i

( )

3 3 3 1 2 i i i i j j j i j j i yukawa j j j j E N r c N r dr N r c N r dr N r V r c N r dr Ω Ω Ω ⋅ = ∇ ∇ +

∑

∫

∑

∫

(3-11)

By applying equation (3-11) to each hexahedron in the computational domain, we can construct a generalized eigenvalue matrix equation into the form as Ax=λB x, which can then be solved numerically. In the present study, initial spatial distribution of wave function for the time-dependent Schrödinger equation is obtained by numerically

(42)

solving the generalized eigenvalue matrix equation as mentioned in the above using the Jocobi-Davidson algorithm [58].

3.5 Discretization of the Three-dimensional Time Dependent Schrödinger Equation

As introduced in chapter 1 that very few studies have focused on direct real-space discretization of the TDSE, except those using FDM [33] and FEM by Collin’s group [33-37]. It is generally agreed that for solving PDEs the memory storage and computational time by using FEM would be higher as compared to the cell-centered FVM for achieving the same solution accuracy, as mentioned in chapter 1. In this thesis, we solve the 3D TDSE directly on real-space coordinates using cell-centered FVM, which is much simpler in practical implementation and faster in simulation speed as compared to those using FEM.

By first dividing the volume of interest into several discrete cells and applying the standard finite-volume method [59] by taking volume integration to the TDSE, equation (3-2), in each discrete cell, we can obtain

(

)

3 1 2

(

)

(

)

3 , , , 2 i r t dr V r t r t dr t_Ωψ _Ω ψ ∂   = − ∇ +  ∂

∫

  (3-11)

where Ω represents the cell volume of interest. Next, by applying the divergence theorem, equation (3-11) is reduced to

(43)

(

)

3 1

(

)

(

) (

)

3 , , , , 2 i r t dr r t d s V r t r t dr t_Ωψ _∂Ω ψ _Ω ψ ∂ = − ∇ ⋅ + ∂

∫

(3-12)

where ∂Ω represents the cell surfaces of interest. We then apply cell-centered finite-volume scheme, in which the variable of wave function ψ is placed at the centroid of the cell. In the present thesis, we approximate the spatial part of equation (3-12) using Cartesian-grid based non-uniform hexahedral cells (or termed “non-conformal” mesh), which a typical sketch of mesh projected in two-dimensional plane is shown in Figure 3. Then, in each cell, equation (3-12) can be simply approximated as

(

)

(

)

(

) (

)

1 1 , , , , 2 s N c m i c c i i r t V r t s V r t r t V tψ Ω = ψ ψ Ω ∂ ∆ = − _∇ ⋅ ∆ _+ ∆ ∂

∑

(3-13)

where the subscript c and m represents the centroid of cell Ω and surface ∆ , s_i

respectively. In addition, Ns and ∆VΩ represents number of surfaces and volume of the cell under consideration. Note the gradient terms at the cell interface can be further approximated by central difference scheme using values of ψ at centroids of neighboring cells.

The time propagation term on the left-hand side of equation (3-13) is approximated using an explicit stagger-time algorithm following the idea presented by Visscher [60] for 1D time-dependent Schrödinger equation, in which the algorithm was shown to be 2nd-order accuracy in time for an uniform grid. The ideas are redescribed here for

(44)

(

)

(

)

R iI i H R iI t ∂ + = + ∂ (3-14)

where ψ( , )r t =R r t( , ) +iI r t( , ) , and R and I represents the real and imaginary part of the wave function ψ, respectively. By separating these two terms, equation (3-14) can be further reduced to R HI t ∂ = ∂ (3-15a) I HR t ∂ = − ∂ (3-15b)

Then, we can propagate the equation (3-15a) and equation (3-15b) alternatively in time using a leap-frog like explicit scheme, termed as “explicit stagger-time scheme”, as the following:

(

)

1 1 , , , , 0 ~ 2 2 R r t + ∆t=R r t − ∆t+ ∆ ⋅t HI r t t= T     (3-16a)

(

)

1 1 1 1 , , , , ~ 2 2 2 2 I r t + ∆t=I r t − ∆t− ∆ ⋅t HR r t t= − ∆t T − ∆t     (3-16b) where T is the total simulation time. Data are then synchronized by simple temporal interpolation of either R or I. In addition, absorbing type boundary conditions [61] are employed at the outer boundaries of computational domain, which is similar to previous studies in this regard. Note larger domain size is often necessary to delay the wave reflection from the numerical boundaries, which indeed is worthwhile of studying in the future. The present TDSE solver is designed to easily set up the arbitrary number of regions having different cell sizes centered around the molecule, where most refined cells are clustered in this region (as shown in Figure 3). In generating the

(45)

non-conformal grid, we always set the ratio of cell length between two neigboring cells at the interface refined regions as two, which further ensures numerical accuracy using this kind of mesh.

3.5. Parallel Implementation of TDSE Solver

The above algorithm is readily parallelized through decomposition of the physical domain into groups of cells which are then distributed among the parallel processors. Each processor executes the explicit stagger-time scheme in serial for all cells in its own domain. Parallel communication between processors is required when evaluating intrerfacial fluxes of wave function requiring cell-centered data to be transferred between processors. To achieve high parallel efficiency it is necessary to minimize the communication between processors while maintaining a balance between the computational load on each processor. In the present proposal, since we have adopted non-conformal mesh, it would be very difficult to have approximately equal number of cells in each processor simply using coordinate-based partitioning technique. Instead, we have used a publicly available mutli-level graph-partitioning library [62] for decomposing the computational domain. With this library, we can easily achieve the requirement of approximately having the same number of cells in each processor with any arbitrary non-conformal mesh. Figure 4a-4c, respectively, shows a typical

(46)

mid-plane for solving the TDSE, domain decomposition through the mid-plane, and surface domain decomposition for the 3D mesh. Most important of all, we would expect high parallel efficiency can be obtained since we have applied the explicit stagger-time scheme for time propagation having minimal communication load once the load is properly balanced among processors.

(47)

Chapter 4. Code Verification

4.1. Ground state energy of H atom and H2+ ion molecule

The hydrogen atom and hydrogen molecular ion are the simplest atomic system and molecular system, and the hydrogen atomic system has analytic solution. These are a good for verification casese. The exact eigen-energy for a hydrogen atom are 1₂

2n −

and n is the principle quantum number. The exact eigen-energy for a hydrogen molecular ion can be found in reference [63]. We first want to verify time independent Schrödinger equation solver (Jocobi-Davidson algorithm) by comparing the exact and calculated eigen-energy. The exact ground state energy (n=1) of hydrogen atom is -0.5, we make a ground state energy calculation for H atom, and use 343000 mesh elements. The calculated eigen-energy is -0.4994098 a.u. and compares to exact value -0.5. The difference between the two values is 0.1180%.

For the H2+ ion molecule, we make a ground state energy calculation in different

inter-nuclear distance from 0.2 to 3 a.u. The results are listed in table 1. These calculated eigen-energys are different from exact values [63] by less than 0.2%.

These results are quite good as an initial state for a time-dependent Schrödinger equation problem. For higher accuracy requirement, it can be reached by increasing the

(48)

solver for the initial wave function calculation with model potential that conducted in chapter 3.

4.2. Energy Conservation of H2 +

Without Laser Incidence

As a first verification of the implementation of the parallel TDSE solver, we have monitored the time evolution of total energy conservation of the H2+ molecule without

laser interaction. Simulation conditions include: 2.5 million hexahedral cells, 192 a.u. of length in x and y directions and 288 a.u. length in z direction (molecular axis), time-step size of 0.01 a.u. and inter-nucleus distance (R) of 9 a.u. Number of processors is kept as 16, unless otherwise specified. We calculated total energy as a function of time, in which the distribution of ground-state wave function is used as the initial condition. As figure 5 shows that total electron probability (1.0) and total energy (-0.6216 eV) are both nearly conserved with very small variance of 0.001% and 0.04%, respectively, even after long-time propagation (35 fs), which demonstrates the parallel TDSE solver is implemented correctly at least without considering the externally applied electric field.

4.3. Ionization Rates of H2+ With Laser Incidence Along the H-H Axis

As a further verification, we apply an intense laser pulse field to the H2+ molecule

(49)

predicted previously using axis-symmetric TDSE solvers. Applied laser pulse is the same as that in [64], as shown in Figure 6, but related data are repeated here for completeness. Envelope of the laser pulse is E f t₀

( )

cos

( )

ωt and f t

( )

is defined as follows,

( )

1 1 1 1 2 1 2 1 2 1 1 1 cos , 0 2 1, 1 1 cos , 2 2 0, t t t f t t t other π τ τ τ τ τ π τ τ τ τ τ     − ≤ ≤           ≤ ≤ +  =     − + ≤ ≤ +           ,τ₁=5,τ₂=15 (4-1)

Simulation conditions include: laser intensity of 1014 W/cm2, wavelength of 1064 nm, 5.5~7.5 million hexahedral cells, 224~254 a.u. of length in z direction, and time-step size of 0.01 a.u.. Figure 7 illustrates the comparison of ionization rates as a function of inter-nuclei distance among different studies. We employed equation (6) of reference [64] to calculate the ionization rate as a function of time, and then obtained the presented ionization rate by averaging over time during which the evelope levels at highest amplitude (5-20 optical cycles). Results show that our results agree reasonably well with previous simulation data using axisymmetric TDSE codes. The well-known phenomenon, which has been observed experimentally, such as “Coulomb explosion” at R=9 a.u. is clearly reproduced by our simulation as others’ axisymmetric code. Briefly

卡氏座標下三維時變薛丁格方程式的平行化程式發展及其於雷射與分子間交互作用之應用研究

國

立

交

通

大

學

機械工程學系

博士論文