阿秒脈衝激發和遠紅外雷射驅動中稀有氣體原子的光子發射譜中的半脈衝週期震盪: 自作用修正的隨時變密度泛函理論計算

國立臺灣大學理學院應用物理所 碩士論文

Graduate Institute of Applied Physics College of Science

National Taiwan University Master Thesis

阿秒脈衝激發和遠紅外雷射驅動中稀有氣體原子的光子發射 譜中的半脈衝週期震盪:

自作用修正的隨時變密度泛涵理論計算

Subcycle Dynamics of Photon Emission Spectra of Rare Gases Atoms Excited by Attosecond Pulses and Driven by Near-Infrared

Laser Fields:

Self-Interaction-Free Time-Dependent Density-Functional-Theory Approach

周繼暉 Chi-Hui Chou

指導教授：朱時宜 博士 Advisor: Shih-I Chu, Ph.D.

中華民國 103 年 7 月

July, 2014

誌謝

首先要感謝指導教授朱時宜老師，老師是個有耐心和學術涵養豐富的教授，

感謝他細心的指導，能夠容忍並且修正我犯一些錯誤，讓我能夠提升自己，老師 儘管在強場物理已發展出精確的理論，仍然有熱誠於研究新的且更完整的理論，

那精神值得我們晚輩學習。

我也要感謝 Dr. John Heslar, 我的論文主要延伸至他的工作，感謝他在數值方 法還有 Fortran 上給我很多幫助，他平時幽默風趣。Dr. Ho, Dr. Tong, Dr. Telnov 等 訪問學者也給予我一些幫助。

實驗室的學長和博士後也給我很多幫助，能夠給我在學術和生涯一些建議的 佑航，幫我處理計算主機、latex 及繪圖軟體等問題的聖崙，研究上很認真的 Dr.

Lee，經常邀大家吃飯的曄琳，來自澳門和時常討論時事的焌梃，英文很好和在 CUDA 常和我討論的昶明，和我一起進來實驗室的書豪。還有幫我們處理瑣事的 助理秀翎和曉佩，幫實驗室丟垃圾的阿姨，感謝你們讓我能在碩士二年能夠順利。

除了實驗室外，我還要感謝碩一時一起打電動和討論功課的焌佑，已經結婚 且帶小孩的榮伸；以及室友們，身在美國深造的家銓，將要入職場的浩明，打理 環境的伯棟，熱愛調酒的柏聿，讓我在生活上能夠安穩的做研究。

最後感謝我的家人以及我教會的朋友們，你們在背後支持我，鼓勵我能夠堅 持自己的夢想，並且教導我學術研究以外的事。碩士結束只是個一個階段的開始，

希望我自己能夠完成我的夢想以及面對挑戰。

中文摘要

本文應用了自作用修正的隨時變密度泛涵理論和非微擾的數值計算方法，計 算阿秒脈衝激發和遠紅外雷射驅動中氫、氦和氖原子的光子發射譜，藉由改變阿 秒脈衝和遠紅外雷射的時間差，我們可以觀察到半雷射周期的震盪以及能階的改 變。這種現象已從吸收光譜實驗上觀察到，我們也計算出激發態的電子機率隨時 間差有著相同的週期震盪，我們從兩個光子吸收觀點來解釋這種現象。

關鍵字：阿秒脈衝、自作用修正時變密度泛涵理論、半雷射週期現象

Abstract

We present an ab initio method to study the sub-cycle dynamics of hydrogen, helium and neon atoms in near-infrared(NIR) laser fields subject to excitation by a single ex- treme ultraviolet attosecond pulse(SAP). We extended the self-interaction-free time- dependent density functional theory(TD-KLI-SIC) to describe multi-electron system and solve the time-dependent Kohn-Sham equations by time-dependent generalized pseudospectral(TDGPS) method. We calculated the photon emission spectra and population of several excited states as the function of the time delay between the NIR pulse and SAP. The phenomena can be explain by two-photon absorption.

keyword: attosecond, time-dependent density functional theory, subcycle dynamics

Contents

Contents ii

List of Figures iv

1 Introduction 1

1.1 Attosecond Physics . . . 1

1.2 Real-time observation with attosecond technology . . . 2

1.3 Sub-cycle AC stark shift . . . 3

1.4 Purpose of this work . . . 5

2 Theory and Method 6 2.1 Time-dependent Generalized Pseudospectral Method . . . 6

2.1.1 The eigenvalues problem . . . 6

2.1.2 Time propagation . . . 11

2.2 Time-Dependent Density Functional Theory . . . 15

2.2.1 Density Functional Theory and Kohn-Sham scheme . . . 15

2.2.2 Optimized Effective Potential method and Krieger-Li-Iafrate approximation . . . 17

2.2.3 KLI-SIC method . . . 18

2.2.4 TD-KLI-SIC method . . . 20

2.3 Implement of numerical methods on graphics processing unit . . . 22

2.3.1 GPU architecture . . . 23

2.3.2 Implementation on GPU . . . 26 2.3.3 Result . . . 27

3 Result and Discussion 28

3.1 Hydrogen . . . 28 3.2 Helium . . . 38 3.3 Neon . . . 44

4 Conclusions and Perspectives 50

Bibliography 51

List of Figures

1.1 Characteristic length and time scales for structure and dynamics.[1] . 2

1.2 Illustration of pump-probe (time-resolved) spectroscopy.[1] . . . 3

2.1 The basis function in correspond to domain.Reproduced from[2]. . . . 8

2.2 Data distribution of wave function. N_l and N_r is the total number of partial wave and the grid points of wave function. . . 12

2.3 Illustration of matrix-vector multiplication. . . 13

2.4 Illustration of the transformation in Eq. 2.29 . . . 14

2.5 Mask function for one dimensional spatial grid with absorbing bound- ary condition. Reproduced from[3]. . . 15

2.6 The effective potential rV_{ef f} with the LSDA and LSDA-KLI-SIC in neon and argon(left to right). . . 19

2.7 GPU architecture. . . 23

2.8 SM architecture. . . 24

2.9 GPU Memory. . . 25

2.10 Runtime of the TDGPS and TD-KLI-SIC on one GPU and 16-cores CPU with different atoms. (nvidia Kepler K20(GPU) and Intel(R) Xeon(R) CPU E5-2690 2.9GHz(16 cores CPU)) . . . 27

3.1 Illustration of SAP and NIR with time delay -5 fs. . . 29

3.2 Photon emission energy spectrum of the exicted states (2p[3s− 5s]

and 2p[3d− 4d] as a function of the time delay between the NIR pulse and SAP. The yellow color indicates the highest energy emitted. The color bars are represented by the log₁₀S(ω) of the spectral density in Eq. 3.7 . . . 31 3.3 Photon emission energy spectrum of the 1s2p excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 32 3.4 Photon emission energy spectrum of the 1s3p excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 32 3.5 Photon emission energy spectrum of the 1s4p excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 33 3.6 (left) Energy emitted near 1s2p transition (right) Energy emitted near

1s3p transition . . . 34 3.7 Illustration of Autler-Townes effect in hydrogen. When SAP comes

first, the electron can be exicted to np orbitals. The NIR is not weak for excited states, so it can make 2p3s and 2p3d transition happen. . 34 3.8 Population of several excited states as a function of the time delay

between the NIR pulse and SAP. The center frequency of NIR is 800 nm. . . 35 3.9 Illustration of two photon absorption. . . 36 3.10 Population of several excited states as a function of the time delay

between the NIR pulse and SAP.The center frequency of NIR is 800 nm. . . 37 3.11 Photon emission energy spectrum of the exicted states (1s[2p− 5p]

as a function of the time delay between the NIR pulse and SAP. The yellow color indicates the highest energy emitted. The color bars are represented by the log₁₀S(ω) of the spectral density in Eq. 3.7 . . . . 39 3.12 Photon emission energy spectrum of the 1s2p excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 40

3.13 Photon emission energy spectrum of the 1s3p excited state as a func- tion of the time delay between the NIR pulse and SAP. . . 40 3.14 Photon emission energy spectrum of the 1s4p excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 41 3.15 (left) Energy emitted near 1s2p transition (right) Energy emitted near

1s3p transition . . . 42 3.16 Population of several excited states as a function of the time delay

between the NIR pulse and SAP. . . 42 3.17 Photon emission energy spectrum of the excited states (2p[3s− 5s]

and 2p[3d− 4d] as a function of the time delay between the NIR pulse and SAP. The yellow color indicates the highest energy emitted. The color bars are represented by the log₁₀S(ω) of the spectral density in Eq. 3.7 . . . 45 3.18 Photon emission energy spectrum of the 2p3s excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 46 3.19 Photon emission energy spectrum of the 2p4s excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 46 3.20 Photon emission energy spectrum of the 2p3d excited state as a func-

tion of the time delay between the NIR pulse and SAP. . . 47 3.21 Photon emission energy spectrum of the 2p5s , 2p4d and 2p6s(down

to up) excited states as a function of the time delay between the NIR pulse and SAP. . . 47 3.22 Population of several excited states as a function of the time delay

between the NIR pulse and SAP. . . 49

Chapter 1 Introduction

In this chapter, we introduce the basic concept about attosecond physics and the purpose of our work.

1.1 Attosecond Physics

Attosecond(as,10⁻¹⁸ second) science and technology has been popular topic in recent year since the first generation of attoscecond pulse in experiments[4, 5]. This tech- nique allow us to study the electronic processes in atoms, molecules and surfaces at the attosecond timescales. There are two main progresses in attosecond science. The first one is that development of attosecond source. We want to generate more intense and shorter attosecond pulse. In our group, we focus on High-Order Harmonic gen- eration which can generate attosecond pulse. The second one is that application of attosecond pulse. There are many experiments on which they use attosecond pulse to discover ultra-fast the dynamic of electron[6, 7, 8, 9, 10, 11, 12].The complete discussion in attosecond physics can be found in [1]

1.2 Real-time observation with attosecond technol- ogy

In Fig. 1.1, this defines the characteristic time scale for motion of atoms in molecules to hundreds of femtosecond. The motion of individual electrons in semiconductor nanostructures, molecular orbitals, and the inner shells of atoms occurs from ten femtosecond to one attosecond. Motion of nuclei is even faster, on a zeptosecond time scale(10⁻²¹ second). The main application of attosecond pules is to observe ultrafast process in attosecond time scales.

Figure 1.1: Characteristic length and time scales for structure and dynamics.[1]

Real-time observation of ultrafast motion requires the ability to trigger and probe the process under scrutiny. Dynamical information is provided by an observable varying as a function of the delay between triggering and probing events in a pump- probe measurement in Fig. 1.2, where pump pulse triggers the process and probe pulse images the process at some delay time. This quantity varies on the time scale at which the motion occurs, affording the observer real-time access to the

process. If the observable can reveal the information on location of the moving particles, a series of freeze-frame pump-probe images allows retrieval of the ultrafast motion. Real-time observation and control of atomic motion relies on femtosecond laser techniques, using cycle-averaged quantities such as field amplitude and carrier frequency for triggering, probing and controlling dynamics. Attosecond technology allow us to improve the resolution of probing and control by orders of magnitude with sub-fs XUV pulse.

Figure 1.2: Illustration of pump-probe (time-resolved) spectroscopy.[1]

1.3 Sub-cycle AC stark shift

Stark shift is the shift of atomic energy level in a external field, which can be static(DC) or dynamic(AC). With the existing theory[13], when the monochromatic field with w_L far from atomic resonance, the AC stark shift is equivalent to the quadratic DC Stark shift. The energy shift is proportional to the cycle average of the laser field’s square:

M Ea =−α_a(ω_L)

2 hεL(t)²i = −α_a(ω_L)ε₀(t)²

4 (1.1)

where ε_L(t) = ε₀cos(ω_Lt) is the laser electric field. The polarizability α_a = P

k6=a

ω_ka|dka|² ω_ka² − ω²L

depends on both the dipole matrix elements d_ka coupling a with

other states k and the detuning of the laser frequency w_L. For highly excited states, we have

alim→∞α_a(ω_L) =− 1

ω_L² (1.2)

which leads to the well-known upshift of ionization threshold by the ponderomotive energy: U_p ≡ _4ω^ε202

L

.

Let|ni be the eigenstates of the time-independent Hamiltonian H0, such that H₀|ni = En⁰|ni(n = 1, 2, . . . ). (1.3) The energy of |ni under the external field ε(t) = ε0(t) cos(ω_Lt), based on second- order time-dependent perturbation theory, is given by

E_n(t) = E_n⁰ + ε(t)d_nn− iX

k6=n

Z t 0

dt⁰ε(t)ε(t⁰)e^{−iωkn(t−t0)}|dnk|² (1.4)

Supposed that laser pulse envelop is ε₀(t) = ε_pe^−|t|τp with duration τ_p, the integra- tion in Eq. 1.4 can be simplified for multicycle pulses:

δE_n(t) = 1

2ε₀(t)²X

k6=n

[ω_nk|dnk|² ω_kn² − ω²L

cos²(ω_Lt)− iω_L|dnk|² ω_kn² − ωL²

sin(2ω_Lt)] = 1

2ε₀(t)²[α_ncos²(ω_Lt)− iγnsin(2ω_Lt)].

(1.5)

where γ_n=P

k6=nωL|dkn|² ω²_kn−ω²L

specifies the subcycle changes in the population of n as it couples to k.

However, the average-cycle shift can be verified by using long-pulse( ns) but lacked time resolution. Recently, pump-probe measurements using probe laser pulses longer than the oscillation period of the Stark field revealed Stark shifts with time resolution 10 fs[14], but only the average-cycle shift can be measured.

To resolve the subcycle ac stark shift, we probe the atoms with single extreme ultraviolet attosecond pulse(SAP) with duration nearly 20 times smaller than the NIR Stark laser period. When the atom absorb an XUV photon, electrons can be moved to one of the excited states or to the continuum states because the wide

frequencies range of XUV pulse. When the NIR field act on the electron, the effects of the NIR field can be observed through the changes in the photo emission spectra and population of excited states. In most experiments, they use the transient absorption technique to observe the subcycle oscillation[15, 16, 17, 18]. But for our work, we calculate the photo emission spectra and observe the similar features.

1.4 Purpose of this work

We give the main targets of this work.

1.We perform 3D calculation of the rare gas atoms (H,He,Ne) in the NIR field subject to excitation by SAP. We have calculated the photon emission spectra and the population of excited states with respect to time delay between the SAP and NIR fields and then we found the subcycle oscillation.

2.We include that multi-electron correlation effect by using time dependent den- sity functioanl theory (TDDFT). Unlike other theory which only makes single-active- electron(SAE) model[6, 8, 9, 15, 17, 19] or time independent model potential[20].

Chapter 2

Theory and Method

In this chapter, we introduce the numerical method for solving PDE and the theory for many-electron system. In addition, we also give simple introduction about graphic processing unit(GPU), on which we can accelerate our calculations.

2.1 Time-dependent Generalized Pseudospectral Method

There are many numerical methods for solving PDE, finite difference(FD), finite element(FE), and spectral method[2, 21, 22]. GPS is under the category of spectral method. Spectral method is more accurate and efficient than FD and FE because it needs less grid points and converges faster. Moreover, the grid points of GPS are more denser near the origin, so we can describe Coulomb potential better than equal-spacing grid points.

2.1.1 The eigenvalues problem

The basic idea of spectral method is to approximate function f (x) by order N poly- nomial with orthogonal basis φ_j(x)

A(x)' XN

j=0

a_jφ_j(x) (2.1)

We need the coefficient a_j for defining the function, so it requires the values of function at N + 1 grid points x_j.

a_j = XN

i=0

ω_iA(x_i)φ_j(x_i) (2.2)

The approximation of function can be expressed by

A(x) = XN

j=0

A(xj)gj(x) (2.3)

where g_j(x) is the cardinal function given by

g_j(x) = XN

i=0

ω_iA(x_i)φ_j(x) (2.4)

With the definition of the cardinal function and the collocation points{xi}, we can obtain the differential operator matrix

d

dxA(x)|x=xi ' XN

j=0

A(x_j) d

dxg_j(x)|x=xj = XN

j=0

(D_N)_ijA(x_j) (2.5) The value of weight and the form of the cardinal function and the differential operator matrix are depend on basis function and grid points. How to choose the basis function and grid points is dependent on the problems, mainly the boundary condition and domain. We list the basis function in Fig 2.1

Figure 2.1: The basis function in correspond to domain.Reproduced from[2].

In this work, we choose the Legendre polynomials. And then we choose the grid points from different boundary condition[21].

Legendre-Gauss

x_j(j = 0, . . . , N )zeros ofP_{N +1};

ω_j = 2

(1− x²j)[P_{N +1}⁰ (x_j)]², j = 0, . . . , N.

Legendre-Gauss-Radau

x_j(j = 0, . . . , N )zeros ofP_N + P_{N +1}; ω₀ = 2

(N + 1)²

ω_j = 1− xj

(N + 1)²[P_N(x_j)]², j = 1, . . . , N.

Legendre-Gauss-Lobatto

x0 =−1, xN = 1, xj(j = 1, . . . , N − 1)zeros ofPN +1⁰ ;

ωj = 2

N (N + 1)[P_{N +1}(x_j)]², j = 0, . . . , N.

We use Legendre-Gauss-Lobatto grid points for the radical coordinate and Legendre- Gauss grid points for the angular coordinate.

For atomic or molecule calculation including Coulomb potential, the general prob- lem is the singularity at r = 0 and the long-range potential. We truncates the semi-infinite (0,∞) to the finite interval [rmin, r_max] to avoid the problem, but it needs a lot of grid points for convergence and accuracy and cost much time for time propagation.

The generalized pseudospectral method offers us the solution of above problem.

This method can be used in time-dependent problem. Follow the discussion, we use nonlinear mapping function[23, 24] to map the [−1, 1] to [rmin, r_max]

r(x) = r_m 1 + x

1− x + α (2.6)

where r_m and α = 2_r^rm

max is the mapping parameters. Here, we begin with the time- independent Schrödinger equation for Hydrogen atom.

ψ(r, θ, φ) = X

n,l,m

Cnlm

ϕ_nl(r)

r Ylm(θ, φ) (2.7)

Hˆ_l⁰ϕ_nl(r) = (−1

2∇²+ V_l(r))ϕ_nl(r) = E_nlϕ_nl(r) (2.8) where

V_nl(r) =−1

r +l(l + 1)

2r² (2.9)

Applying the nonlinear mapping for Eq. 2.6.

−1 2[ 1

r⁰(x)² d²

dx² − r⁰⁰(x) r⁰(x)³

d

dx]ϕ(r(x)) + V (r)ϕ(r(x)) = Eϕ(r(x)) (2.10) This equation is not symmetrical, so we choose ϕ(r(x)) =p

r⁰(x)f (x)

−1 2

1 r⁰(x)

d²

dx²f (x) + V_m(x)r⁰(x)f (x) + V (x)r⁰(x)f (x) = Er⁰(x)f (x) (2.11)

where

V_m(x) = 3(r⁰⁰)²− 2r⁰⁰⁰r⁰

8r⁰⁴ (2.12)

For special mapping Eq. 2.6, V_m(x) = 0. With a lot of effort, we finally translate the problem into the symmetric eigenvalue problem.

−1 2

1 r⁰(x)

d² dx²

1

r⁰(x)A(x) + V (x)A(x) = EA(x) (2.13) where

A(x) = r⁰(x)f (x) =p

r⁰(x)ϕ(x) (2.14)

We throw Eq. 2.3 into Eq. 2.13

−1 2

1 r⁰(x_i)

X

j

A(x_j) r⁰(x_j)

d²g_j(x) d²x

xi

+ V (xi)X

j

A(xj)gj(xi) = EX

j

A(xj)gj(xi) (2.15)

where[24]

g_j⁰⁰(xi) = d⁽²⁾_ij P_N(x_i)

P_N(x_j) (2.16)

and

d⁽²⁾_ij =− 2

(x_i− xj)²[i6= j, (ij) 6= (0N), (ij) 6= (N0)] (2.17)

d⁽²⁾_0N = d²_{N 0}= N (N + 1)− 2

4 (2.18)

d⁽²⁾_jj =−N (N + 1)

3(1− xj)²[j 6= 0, j 6= N] (2.19)

d⁽²⁾₀₀ = d²_{N N} = N (N + 1)[N (N + 1)− 2]

24 (2.20)

And the cardinal function has the property

g_j(x_i) = δ_ij (2.21)

and choose

A_i = s

2

(N + 1)(N + 2)

A(xi)

P_{N +1}(x_i) (2.22)

The whole equation is simplified to matrix form X

j

[(D₂)_ij + V (x_i)δ_ij]A_j = EA_j (2.23) The final wave function on the collection points is

ψ(r(x_i)) = A_i rp

ω_ir⁰(x_i) (2.24)

where the ω_i is the weight of the cardinal function. And then we can check the wave function normalized.

< ψ|ψ >=

Z A_iA^∗_i

r²ω_ir⁰(x_i)r²dr

=X

i

AiA^∗_i

r²ω_ir⁰(x_i)r²dr dxω_i

=X

i

AiA^∗_i

(2.25)

Finally, we solve the Schrödinger equation and get the eigenvalues and eigenfunc- tions of the system. Moreover, there are something we should be careful.The real wave function is in Eq. 2.24 although the summation of A_iA^∗_i is normalized by most eigensolver. It’s convenient to see Ai as the wave function numerically because some variables are eliminated by integral elements in Eq. 2.25.

2.1.2 Time propagation

The next step is to do time propagation by second-order split operator technique[25]

in spherical coordinates:

ϕ(r, t + δt)' exp(−i bH0δt/2) exp(−ibV (r, θ, t + δt)δt) exp(−i bH0δt/2)ϕ(r, t) + O(δt³) (2.26) The first term and third term is from time independent Hamiltonian. The the second term is the time dependent potential, namely the laser field.

Figure 2.2: Data distribution of wave function. N_l and N_r is the total number of partial wave and the grid points of wave function.

We define the initial wave function for hydrogen.

ψ(r_i, l = 0) = A(l = 0)_i

ψ(r_i, l = 1 . . . N_l− 1) = 0 (2.27)

where the i is the index of radical part of the wave function and l is the partial wave number from angular quantum number. N_l and N_r is the total number of partial wave and the grid points of wave function. In general, we set N_l = 32 and Nr = 256 for the hydrogen atom.

We define the exp(−i ˆH⁰δt/2) by the eigenvalues and the eigenfunctions of Hamilto- nian.

Sij(l) =X

n

< ri, l|n, l > exp(−iEnlδt/2) < n, l|rj, l > (2.28) where the eigenfunctions are < ri, l|n, l > and eigenvalues are Enl.

Figure 2.3: Illustration of matrix-vector multiplication.

We perform the operationexp(−iH⁰δt/2)ψ(r, 0)in Fig. 2.1.2, composed of matrix- vector multiplications for each l. The process can be paralleled by openmp or other parallel methods.

Before we do exp(−iV (r, θ, t)δt), we need to transform the representation(r, l) to (r, θ). The θ_j is the Gaussion-Legendre points we mention before.

< r_i, θ_j|ψ >=X

l

< r_i, θ_j|ri, l >< r_i, l|ψ > (2.29) The Eq. 2.29 can be illustrated in Fig. 2.1.2.

exp(−iV (ri, θ_j, t)δt)ψ(r_i, θ_j) (2.30) And we transform the representation (r, θ) back to (r, l) before we perform the oper- ation exp(−i ˆH⁰δt/2) again. To prevent reflection, the wave function are multiplied by mask function after each time step. We partition our finite spatial grid into an

Figure 2.4: Illustration of the transformation in Eq. 2.29

inner region, which is large enough to completely contain the finite system of interest , and a border region, where outgoing flux is to be absorbed.

Figure 2.5: Mask function for one dimensional spatial grid with absorbing boundary condition. Reproduced from[3].

Throughout all processes above, we get the next time step wave function and we can calculate the observables.

O(t) =< ψ(r, t)| ˆO|ψ(r, t) > (2.31) O(t) is the observables,like number of electron,electron density, dipole moment in length form and in acceleration form. Finally, we don’t need to save the time depen- dent wave function.

ψ(r, l) = ψ(r, l, δt) (2.32)

where we set δt as the initial wave function and then to get 2δt, . . . , t_{f inal}. The TDGPS method can be applied to many kinds of system.

2.2 Time-Dependent Density Functional Theory

2.2.1 Density Functional Theory and Kohn-Sham scheme

Density functional theory is the most popular method for many-electron system, atoms, molecules, and solid s. In 1964, Hohenberg and Sham develop the basic theorem of density functional theory[26]. We give simple concept about the theorem.

First theorem, we can represent the energy as a functional of the electron density

for given potential. We don’t need to use wave functionΨ(r₁, r₂, . . . , r_N), which has 3N variables in N-electron system, but use electron densityρ₍r), which has only 3 variables in N-electron system. That avoid the main computational difficulty. Second theorem, The energy functional is minimized by the ground state density. We can find the ground state potential by variantional principle. In Hohenberg-Kohn theorem, if we know the exact form of the universe functional, we can find the ground state by minimizing the functional. But we don’t know the exact functional and systematic way to find such functional. In 1965, Kohn and Sham develop systematic way to approximate the functional and find the ground state density.

Kohn and Sham develop another method[27] for the functional, by using non- interacting system as auxiliary system. Therefore, the ground-state wave function is the single slater-determinant.

Ψ = 1

√N !det[ϕ1ϕ2. . . ϕN] (2.33) And the density is

ρ(r) =X

σ Nσ

X

i=1

|ϕiσ(r)|² = ρ_α(r) + ρ_β(r) (2.34)

The energy functional is

E[ρ_α, ρ_β] = T_s[ρ] + J [ρ] + E_xc[ρ_α, ρ_β] + Z

v_ext(r)ρ(r)d³r (2.35) T_s[ρ] is the non-interacting kinetic energy functional, the J [ρ] is the Hartree en- ergy functional and the E_xc[ρ] is called exchange-correlation energy functional. The exchange energy is from the Pauli-expulsion and the correlation is from the ap- proximation of single slater-determinant. The complexity of system is inside the exchange-correlation functional. And finally we take the density derivative of the energy functional, we get the Schrödinger-like equation

Hˆ_KSϕ_iσ(r) = [−1

2∇²+ v_eff,σ(r)]ϕ_iσ(r) = ε_iσϕ_iσ(r), i = 1, 2, . . . , N_σ,

(2.36)

where v_{ef f,σ} is the effective KS potential and σ is the spin index. The effective potential is

v_{ef f,σ} = v_ext(r) + δJ [ρ]

δρ_σ(r) + δE_xc[ρ_α, ρ_β]

δρ_σ(r) (2.37)

where v_xc,σ(r) is the exchange-correlation potential vxc,σ(r) = δE_xc[ρ_α, ρ_β]

δρ_σ(r) (2.38)

The KS equations are solved self-consistently. One guesses the initial density at first and then solves the KS equation to get the new density from new orbitals until the convergence.

2.2.2 Optimized Effective Potential method and Krieger-Li- Iafrate approximation

The self-interaction is from the classical Coulomb repulsion. The effect of self- interaction should be cancelled by exchange-correlation functional, but for most of exchange-correlation functionals, the self-interaction correction is not consider. One of the most important error is the incorrect long-range tail of Kohn-Sham potential, which will affect the ionization energy. Therefore, the self-interaction correction is crucial for excited states.

Perdew and Zunger proposed the self-interaction correction(SIC)[28] by giving the approximate exchange-correlation energy functional E_xc[ρ_α, ρ_β],

E_xc^SIC[ρ_α, ρ_β] = E_xc[ρ_α, ρ_β]−X

σ Nσ

X

i=1

{J[ρiσ] + E_xc[ρ_iσ, 0]} (2.39) where ρ_iσ is the one-electron density of the ith KS spin orbital.

However, the SIC energy functional is explicit orbital-dependence, so for each electron orbital they have different potentials. That cause each orbital to be nonorthogonal and be complicated.

Another approach is the optimized effective potential method[29, 30] .In this approach, one solves the set of one-electron equations, similar to the KS equations in Eq. 2.36.

Hˆ_OEPϕ_iσ(r) = [−1

2∇²+ v^OEP_σ (r)]ϕ_iσ(r) = ε_iσϕ_iσ(r), i = 1, 2, . . . , Nσ

(2.40)

The optimized effective potential v_σ^OEP(r) is obtained by the orbitals{ϕiσ} in Eq. 2.40 which minimized the energy functional E[ϕ_iα, ϕ_jβ],

δE_xc[{ϕjσ}]

δv^OEP_σ (r) = 0 (2.41)

Eq. 2.41 can be written as, using chain rule for functional derivative, X

j

Z

d³r⁰δE^OEP[ϕ_iα, ϕ_i0β] δϕ_jσ(r⁰)

δϕ_jσ(r⁰)

v^OEP_σ (r) + c.c. = 0 (2.42) Eq. 2.42 leads to an integral equation that is complicated. Krieger, Li and Iafrate[31, 32, 33] make an approximate procedure to simplify the original OEP in- tegral equations into the set of linear equations. Although the KLI procedure can’t reach the exact exchange functional, it reduces the computational difficulty and the its result is pretty close to OEP method.

2.2.3 KLI-SIC method

The OEP method and KLI approximation uses the exchange part of the density functional contains a Hartree-Fock-like nonlocal functional.

E_x^exact[{ϕjσ}] = −1 2

X

σ Nσ

X

i,j=1

Z d³r

Z

d³r⁰ϕ^∗_iσ(r⁰)ϕ_jσ(r⁰)ϕ_iσ(r)ϕ_jσ(r)

|r − r⁰| (2.43)

Even though Eq. 2.43 provides more accurate exchange potential, it’s computa- tionally more expensive than the traditional DFT functional with only local func- tional. Therefore,we present the extension of KLI procedures to the SIC term[34, 35]

in Eq. 2.39. This new KLI-SIC procedure can speed up the static DFT calcula- tion and time dependent DFT calculation. This KLI-SIC procedure make the self- interaction-free effective potential orbital-independent. In other word, this avoid the problems with respect to nonorthogonal spin-orbitals. And the KLI-SIC procedure give the optimized effective potential with the correct long-range behavior(−1/r)in Fig. 2.6 and surprisingly improvement of ionization energy and excited states.

-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0

0 1 2 3 4 5 6 7 8

r*Veff(r) (a.u.)

r (a.u.)

LSDA-KLI-SIC LSDA

-18 -16 -14 -12 -10 -8 -6 -4 -2 0

0 1 2 3 4 5 6 7 8

r*Veff(r) (a.u.)

r (a.u.)

LSDA-KLI-SIC LSDA

Figure 2.6: The effective potential rV_{ef f} with the LSDA and LSDA-KLI-SIC in neon and argon(left to right).

Define the total energy functional with SIC to be E_SIC^OEP[{ϕiα, ϕ_jβ}] = E^OEP[{ϕiα, ϕ_jβ}] −X

σ

X

i

{J[ρiσ] + E_xc[ρ_iσ, 0]} (2.44)

where E^OEP[{ϕiα, ϕ_jβ}] is normal energy functional in Eq. 2.35.Following the OEP- KLI procedure, one finds that

v_SIC,σ^OEP (r) = vext(r) +

Z ρ(r⁰)

|r − r⁰|d³r⁰ +δE_xc[ρ_α, ρ_β]

δρ_σ(r) + vSIC,σ(r), (2.45) where

v_SIC,σ =X

i

ρiσ(r)

ρ_σ {viσ(r) + vⁱ_SIC,σ− viσ}, (2.46) v_iσ =−

Z ρiσ(r⁰)

|r − r⁰|d³r⁰− δExc[ρiσ, 0]

δρ_iσ(r) (2.47)

and

vⁱ_SIC,σ =< ϕ_iσ|vSIC,σ(r)|ϕiσ > (2.48) v_iσ =< ϕ_iσ|viσ(r)|ϕiσ > (2.49) The value of the v_SIC,σ is unknown, but we can solve it through linear equations

NXσ−1 i=1

(δ_ji,σ− Mji,σ)(vⁱ_SIC,σ− viσ) = v^s_jσ− vjσ (2.50)

where

M_ji,σ =

Z ρjσ(r)ρiσ

ρ_σ d³r (2.51)

and

v^s_iσ =< ϕiσ|

Nσ

X

j=1

ρ_jσ(r)v_jσ(r)

ρ_σ(r) |ϕiσ > (2.52) The highest occupied orbital dominates the potential at the long-range. We choose v^i=N_SIC,σ^σ = v_Nσ to make sure the potential has correct asymptotic behavior.

2.2.4 TD-KLI-SIC method

We extend the KLI-SIC method into the time dependent system. The basic theorem of time dependent density functional theory is from Runge and Gross[36], mainly the similar structure of HK theorem and KS scheme. We’ll give the main theorems of TDDFT.More detail proofs and discussions can be found in.

First theorem, there is a one-to-one correspondence between time dependent den- sity and time dependent potential for any fixed initial states. In general, potential, hamiltonian and wave function is the functional of the time dependent density. This is the basic existence theorem of TDDFT.

Second theorem, we define the action A of the many-body system as the functional of many-body wave function,

A[ψ] = Z t2

t1

dt < ψ(t)|i∂

∂t − ˆH(t)|ψ(t) > (2.53)

If the variation of the action is δA[ψ] =

Z _t2

t1

dt < δψ(t)|i∂

∂t− ˆH(t)|ψ(t) > + Z _t2

t1

dt < ψ(t)|i∂

∂t − ˆH(t)|δψ(t) >

(2.54) We do integration by parts in the second term:

δA[ψ] = Z _t2

t1

dt < δψ(t)|i∂

∂t− ˆH(t)|ψ(t) > + Z _t2

t1

dt < (i∂

∂t− ˆH(t))ψ(t)|δψ(t) > +i < ψ(t)|δψ(t) >^t2

t1

(2.55) The third term is zero because of the boundary condition on t₁ and t₂. The δA[φ] = 0 leads to the time dependent Schrödinger equation (i_∂t^∂ − ˆH(t))φ(t) = 0. With the first theorem, the wave function is the functional of density, so we define

A[ρ] = Z t2

t1

dt < ϕ[ρ](t)|i∂

∂t− ˆH(t)|ϕ[ρ](t) > (2.56) And we rewrite the Eq. 2.56 as

A[ρ] = A₀[ρ]− Z _t2

t1

dt Z

d³rρ(r, t)v(r, t) (2.57) The action A₀ is the universal functional from kinetic and electron-electron interac- tion term. The time dependent density can be solved from the variational principle.

δA[ρ]

δρ(r, t) = 0 (2.58)

Eq. 2.58 leads to the set of one-electron time dependent Schrödinger-like equation(TD Kohn-Sham equation).

(−1

2∇²+ v_{ef f,σ}(r, t))ϕ_iσ(r, t) = i∂

∂tϕ_iσ(r, t) (2.59) where the effective potential is

v_{ef f}(r, t) = v(r, t) +

Z ρ(r⁰, t)

|r − r⁰|dr⁰+ δA_xc[ρ]

δρ(r, t) (2.60)

The last term is the time dependent exchange-correlation potential.

The OEP method is also available in the time dependent system[37].

A_xc[ϕ]

v_{ef f,σ}^OEP[(r, t)] (2.61)

We present the TD-KLI-SIC method with adiabatic approximation[38]. The action is defined as

A^SIC_xc = Z _t2

t1

dtE_xc^SIC[ρ₀]|ρ0−→ρ(r,t) (2.62) The adiabatic approximation means that the system is only dependent on instant time. The derivative of the action A^SIC_xc leads to the time dependent potential,

v_xc,σ^SIC(r, t) = δE_xc[ρ_α, ρ_β]

δρ_σ(r, t) + vSIC,σ(r, t) (2.63) where

v_SIC,σ(r, t) =X

i

ρ_iσ(r, t)

ρ_σ(r, t){viσ(r, t) + vⁱ_SIC,σ(t)− viσ(t)}, (2.64) v_iσ(r, t) =−

Z ρ_iσ(r⁰, t)

|r − r⁰| d³r⁰−δE_xc[ρ_iσ, 0]

δρ_iσ(r, t) (2.65)

We solve the v_SIC,σ(r, t) with the same procedures as before.

Finally, we solve the TDKS Eq. 2.60 with TDGPS method. We define ˆH⁰ as Hˆ⁰ =−1

2∇²− Z

r + v_{σ,ef f}(r, 0) (2.66)

And

V (r, t) = v[ρ]ˆ _{σ,ef f}(r, t)− vσ,ef f(r, 0) (2.67) Therefore, we have to solve the static Kohn-Sham equation by the self-consistency.

And we follow the TDGPS method that we discuss before to propagate the wave function.

2.3 Implement of numerical methods on graphics processing unit

In this section, I’ll give simple introduction about GPU and how to implement TDGPS on GPU. I recommend this book for beginner[39].

2.3.1 GPU architecture

Fig. 2.8 shows the architecture of the CUDA-capable graphic processing unit(GPU), which is composed of many streaming multiprocessors(SMs). CUDA is the abbre- viation of compute unified device architecture developed by NVIDIA. CUDA is the hardware and software architecture for the programmers who can develop and ex- ecute the programs in C, C++, Fortran and other languages. The programmer organizes the threads in blocks and grids of blocks in the program, which is called kernel, compile and execute. The programmer don’t worry about how the GPU ex- ecutes the threads and only focus on how to organize the threads. When executing a kernel, the machine will distribute the threads to SMs in the blocks. And for each SMs, the threads in the same block will be distribute to stream processors(SPs) in a SM in Fig. 2.7. In the Tesla K20, there are 13 SMs and 192 SPs in single SM. In general, a group of 32 threads forms a warp to hide the latency, so the threads in a warp was executed by the same SM. If threads are not organized well, the perfor- mance of the GPU would be bad. Therefore, we need to know the architecture in order to exploit the the power from GPU.

Figure 2.7: GPU architecture.

Figure 2.8: SM architecture.

GPU has several memory spaces, mainly global, shared, register. The memory is crucial for speeding up the programs. The register memory is for each thread. It’s very small, of scale of ten bytes but the fastest. The shared memory is for each block, so the threads in the same block share the shared memory which is of the scale of kilo-bytes. The global memory is of the scale of giga-bytes and it can be accessed by any thread but it’s slower than shared memory and register memory. Even though the global memory is not faster, it’s still far faster than CPU memory. Among all the memory access, the slowest one is the communication between CPU and GPU, so we minimize data transferring between them.

Figure 2.9: GPU Memory.

2.3.2 Implementation on GPU

Algorithm 1 pseudocode for TDGPS on GPU Input: ψ(r, 0),exp(−iH⁰δt/2),Grid points(r_i, θ_j),δt Output: O(t)(Observables in each time step)

1: Allocate and Transfer GPU device memory for ψ(r, t),exp(−iH⁰δt/2),Grid points(r_i, θ_j)

2: for t = 0 to t = t_max do

3: Operation ψ₁(r, t) = exp(−iH⁰δt/2)ψ(r, 0)

4: Operation ψ₂(r, t) = exp(−iV (r, θ, t)δt)ψ1(r, t)

5: Operation ψ(r, δt) = exp(−iH⁰δt/2)ψ₂(r, t)

6: Operation < ψ₍r, δt)| ˆO|ψ(r, δt) >

7: ψ(r, 0) = ψ(r, δt)

8: end for

9: Transfer O(t) from GPU to CPU

Here we show the pseudocode for TDGPS on GPU. Because we solve Schrödinger equation at only one time, we can do that on CPU. Next step is to allocate and transfer GPU device memory for wave functions, operators and grid points. This is a basic and important programming technique. Just like preparing ingredients for the cooking. The crucial part is to perform time propagation on GPU. Fortunately, we use the CUBLAS library to perform matrix-matrix and matrix-vector multiplications and it really reduces the difficulty of programming. But we still have to program some parts which the CUBLAS can’t include. Finally, we transfer data from GPU to CPU and print the results. For TDDFT, we solve the static Kohn-Sham equations on CPU because we do that only one time. And in the each time propagation, we have to calculate the effective potential which is the functional of the time dependent density.

2.3.3 Result

Figure 2.10: Runtime of the TDGPS and TD-KLI-SIC on one GPU and 16-cores CPU with different atoms. (nvidia Kepler K20(GPU) and Intel(R) Xeon(R) CPU E5-2690 2.9GHz(16 cores CPU))

Here we show that we speed up our program on GPU. We use the intel mkl library and openmp on multi-core CPU and use the CUBLAS library and parallel the programs on GPU. In our workstation, there are four GPU cards and we can also execute all of these cards at the same time. But I want to say there are some limitations on GPU. There only 5 GB on each GPU card, so we can’t put the data more than 5 GB on single GPU card. Even for four cards, the maximum memory is still only 20 GB which is not large enough for some AMO problems with thousands grid points.

How to reduce the data on time propagation is the main question in our research.

Chapter 3

Result and Discussion

In this chapter, we show the result of the calculation and give explanations of the phenomena.

3.1 Hydrogen

Hydrogen atoms is the simplest system in all atoms. We don’t need to use any approximations about multi-electron effect(only one electron).

We solve the time-dependent SchrÃűdinger equation by TDGPS method.

i∂

∂tϕ(r, t) = [−1

2∇²− 1

r + v_ext(r, t)]ϕ(r, t) (3.1) The laser fields are polarized along z-axis:

v_ext(r, t) =−z[εX(t) + ε_L(t)] (3.2) The SAP field can be defined as follow:

ε_X(t) = F_Xexp(−2 ln(2)(t− td)²

τ_X² ) cos(ω_X(t− td)) (3.3) Here, F_X is the peak field strength of the SAP, τ_X = 140as is its full width at half maximum(FWHM), and ω_X = 13.6 eV is its central frequency(here we choose the laser frequency as the ionization energy of 1s orbital because we want to excite

atoms). The SAP peak intensity is 1 × 10¹⁰W/cm². The parameter t_d represents the time delay between the NIR and SAP; the negative time delay refers to the SAP arriving first. The NIR field has the form:

εL(t) = FLexp(−2 ln(2)t²

τ_L² ) cos(ωLt) (3.4)

Here, FLis the peak field of the NIR pulse, τL fs is its FWHM, and ωL is the central frequency of the NIR field(here we choose the laser wavelength as 800 nm [ω_L= 1.55 e.V] and 656 nm [ω_L = 1.89 e.V]). The NIR laser peak intensity is 1× 10¹²W/cm².

-0.006 -0.004 -0.002 0 0.002 0.004 0.006

-20 -15 -10 -5 0 5 10 15 20

-6e-05 0 6e-05

εNIR(t) (a.u.) εSAP(t) (a.u.)

Time (fs) Delay

SAP NIR

Figure 3.1: Illustration of SAP and NIR with time delay -5 fs.

After the time-propagation procedure, the dipole moment and the dipole accel- eration can be expressed as follow:

d(t) =hϕ(r, t)|r|ϕ(r, t)i (3.5)

a(t) =hϕ(r, t)|∇(1

r − vext(r, t))|ϕ(r, t)i (3.6) The spectral density of the radiation energy is given by the following expression:

S(ω) = 2 3πc³|

Z _∞

−∞

a(t) exp(−iωt)dt|² (3.7)

Here ω is the frequency of radiation, c is the velocity of light. S(ω) has the meaning of the energy emitted per unit frequency.

In the calculation, we use 128 radial and 32 angular grid points and the time step _ω ^2π

L1024 (nearly 0.1 a.u.). The maximum radius is 60 a.u. and we place absorber between 40 a.u. and 60 a.u. describe the ionization process. The time delay was varied in steps of 4td = 20as within the range of −20fs 6 td 6 20fs (2048 steps in total).

Figure 3.2: Photon emission energy spectrum of the exicted states (2p[3s− 5s] and 2p[3d − 4d] as a function of the time delay between the NIR pulse and SAP. The yellow color indicates the highest energy emitted. The color bars are represented by the log₁₀S(ω) of the spectral density in Eq. 3.7

In Fig. 3.2 we show the 3D plot of the photon emission spectrum as a function of t_d for the excited states 1snp(n ≤ 5). The higher excited states (1s4p and 1s5p) are shifted by the pondermotive potential U_p of the NIR field, where U_p = ε²_L

(2ω_L)²; for the field strength and frequency used, U_p = 0.17 eV.

Figure 3.3: Photon emission energy spectrum of the 1s2p excited state as a function of the time delay between the NIR pulse and SAP.

Figure 3.4: Photon emission energy spectrum of the 1s3p excited state as a function of the time delay between the NIR pulse and SAP.

Figure 3.5: Photon emission energy spectrum of the 1s4p excited state as a function of the time delay between the NIR pulse and SAP.

The density plots of the photon emission spectrum in Fig. 3.3-3.5 depict the transition from 1s2p, 1s3p and 1s4p as the function of t_d. We can observe the oscillation structure in the region where the NIR and SAP overlap. The period of the oscillation is 1.1 fs, which is half of the NIR laser optical cycle. This phenomenon was also observed in theoretical calculation where they use absorption spectra[12].

In Fig. 3.3 and Fig. 3.4 we observe the splitting of the lines near t_d∼ 10fs. The electron absorbs one XUV photon to np states and then absorbs more NIR photons to forbidden states, ns or nd. If we try the NIR with wavelength of 656nm, the transition is more obvious. The splitting has been known as result of Autler-Townes effect[40]. We can identify this splitting by the Hamiltonian without some of the excited states. For 1s2p transition, we choose the t_d = 10f s and remove 3s and 3d states in the Hamiltonian and for 1s3p transition we choose the t_d= 10f s and remove 2s states in the Hamiltonian in Fig. 3.2. The splitting disappears in both of them and make sure this splitting can be explained in terms of two-photon absorption and emission process. The SAP excites the ground state to 1snp states; then the NIR