Multielectron effects in strong field processes in moleculesbyYuqing Xia
B.A., University of Science and Technology of China, 2010
M.S., University of Colorado at Boulder, 2013

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Department of Physics
2016This thesis entitled:
Multielectron effects in strong field processes in molecules
written by Yuqing Xia
has been approved for the Department of Physics

Prof. Agnieszka Jaroń-Becker

Prof. Andreas Becker

Date

The final copy of this thesis has been examined by the signatories, and we find that both thecontent and the form meet acceptable presentation standards of scholarly work in the above
mentioned discipline.

iii
Xia, Yuqing (University of Colorado at Boulder)Multielectron effects in strong field processes in molecules
Thesis directed by Prof. Agnieszka Jaroń-Becker

Laser technology has experienced a rapid evolution in available intensities, frequencies, andpulse durations over the last three decades. Many new laser induced phenomena in atoms havebeen discovered, such as multiphoton ionization, above-threshold ionization, high-order harmonic
generation etc. For the interaction with atoms, usually only one electron in the outermost shellis assumed to be active (called single-active-electron approximation) while all other electrons are
considered to remain frozen in their initial states. Due to the extra degrees of freedom (vibration
and rotation) and the more complex structures, the interaction of molecules with intense laser pulses
reveals many new features. Recent experiments have indicated that electrons from inner valenceorbitals of molecules can have significant contributions to ionization and high harmonic generation.
Theoretical analysis of these processes in molecules faces the challenge to extend previous theories
developed for the atomic case by including the multielectron character of the molecular target.
In this thesis we systematically investigate multielectron effects in the interaction of moleculeswith intense laser light. To this end, we apply time-dependent density-functional theory to solve themultielectron Schrödinger equation and analyze highly nonlinear processes such as high harmonic
generation, laser-induced ionization and nonadiabatic electron localization. Based on the results ofour numerical simulations we predict a new feature in the harmonic spectra of molecules, namely
the occurrence of fractional harmonics in the form of Mollow sidebands. Such additional peaks inthe spectra appear due to a field-induced resonant coupling of an inner valence orbital with theoutermost orbital in a molecule. Furthermore, we show that the theoretical explanation of recent
experimental data for the ellipticity of high harmonics in N2 and CO2 require the systematic consideration of all inner valence shells as well as the proper alignment distribution in the experiment.
We also show that the coupling of molecular orbitals in the field can lead to an enhancement of

iv(inner-shell) ionization, potentially leading to a population inversion in the ion, as well as nonadiabatic electron dynamics, where the electron can be trapped at one side of the molecule over several
field cycles. Finally, we present the development of a new intense-field theory based on the Floquet
theorem with complex Gaussian basis sets and show results of first applications for ionization of
simple systems.

Dedication

To my love and family

vi

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisor Prof. Agnieszka JarońBecker for her excellent guidance, caring, patience, and providing me with an excellent atmosphere
during my Ph.D study and research. She taught me a lot how to question thoughts and express
ideas. Her patience and broad knowledge in Chemistry and Physics helped me to overcome manycrisis situation in our project and finish this dissertation. She is also very warmhearted and always
be ready to help us no matter what problem we met in both research and daily life.
Next, I would like to thank Prof. Andreas Becker, who is like a co-advisor for me in thegroup. He is a very nice and easygoing person but with rigorous attitude towards science. I am
deeply grateful to him for the discussion and always get new ideas because of his keen intuition in
Physics.
I would also like to thank my previous and current group members: Dr. Carlos HernandezGarcia, Daniel Weflen, Hongcheng Ni, Jing Su, Andrew Spott, Michelle Miller, Cory Goldsmith,
Benjamin Miller, Ran Reiff, Erez Shani, Quynh Nguyen, as well as Felipe Cajiao. We not only
share our ideas on projects but also any interesting things in life. They helped create a pleasantworking environment and make the research full of fun.
The Janus supercomputer and the clusters at JILA have given me continuous computational
support during these years. Without these computation resources, I would not be able to finish
my research. So I am grateful to all the people working for the Janus supercomputer and the staff
members of JILA computing group.

vii
Finally, and most importantly, I would also like to thank my family and my friends for the
support they provided me. In particular, I must acknowledge my wife. Thanks for her encouragement when I was frustrated and thanks to her to let me have a such lovely son. All in all, without
their love, encouragement I would not have accomplished this thesis.

viiiContents

Chapter
1 Introduction

1

2 Time-dependent density-functional theory

6

2.1

Multielectron system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

The Hohenberg-Kohn theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3

Kohn-Sham equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4

Time-dependent density-functional theory . . . . . . . . . . . . . . . . . . . . . . . . 13

2.5

Exchange correlation functionals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5.1

LDA and LB94 functionals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.6

The pseudopotential approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.7

Implementation of TDDFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.7.1

Comparison of exchange-correlation functionals . . . . . . . . . . . . . . . . . 22

3 Multielectron effects in high-order harmonic generation from molecules
3.1

25Classical and semi-classical description of HHG . . . . . . . . . . . . . . . . . . . . . 26
3.1.1

Three-step model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.1.2

Strong field approximation (SFA) . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.1.3

Saddle point approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.2

High-order harmonic generation from molecules . . . . . . . . . . . . . . . . . . . . . 32

3.3

Fractional harmonics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

ix
3.3.1
3.4

3.5

Open shell molecules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Ellipticity of high-order harmonics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.4.1

Ellipticity from one-electron system H+
2 . . . . . . . . . . . . . . . . . . . . . 56

3.4.2

Complex ellipticity pattern due to two-electron effects . . . . . . . . . . . . . 58

3.4.3

Alignment angle average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.4.4

Comparison with experimental data for N2 . . . . . . . . . . . . . . . . . . . 60

3.4.5

Inner shell contributions to ellipticity of harmonics . . . . . . . . . . . . . . . 66

3.4.6

Ellipticity of harmonics from CO2 . . . . . . . . . . . . . . . . . . . . . . . . 69

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4 Multielectron effects in strong-field ionization from molecules

75

4.1

Role of ionization potential and orbital symmetry . . . . . . . . . . . . . . . . . . . . 76

4.2

Resonance enhanced ionization of open-shell molecules . . . . . . . . . . . . . . . . . 78

4.3

Multiorbital contributions and enhanced ionization from inner-shell orbitals of N2 . . 80
4.3.1

Total ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

4.3.2

Ionization of N2 in filamentation experiments . . . . . . . . . . . . . . . . . . 90

4.4

Ionization from other molecules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

4.5

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

5 Nonadiabatic electron localization

96

5.1

Coupling induced nonadiabatic electron dynamics . . . . . . . . . . . . . . . . . . . . 98

5.2

Dependence of nonadiabatic behavior on laser intensity

5.3

Nonadiabatic localization in two-color laser field . . . . . . . . . . . . . . . . . . . . 102

5.4

Electron localization in other molecules . . . . . . . . . . . . . . . . . . . . . . . . . 104

5.5

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

6 Complex Gaussian basis for non-Hermitian Floquet formalism
6.1. . . . . . . . . . . . . . . . 101

108

Floquet theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

x
6.2

Complex basis sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

6.3

Gaussian basis integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.3.1

Fundamental integral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

6.3.2

Recursion relation for general primitive Gaussian basis integrals. . . . . . . 117

6.4

Technical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6.5

Optimization for complex Gaussian basis . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.5.1

6.6

6.7

Optimization algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.6.1

One photon ionization of H atom . . . . . . . . . . . . . . . . . . . . . . . . . 126

6.6.2

Two photon ionization from H atom . . . . . . . . . . . . . . . . . . . . . . . 130

6.6.3

Two photon ionization from H+
2. . . . . . . . . . . . . . . . . . . . . . . . . 131

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

7 Summary

134

Bibliography

136

xi

TablesTable
3.1

Orbital energies (eigenvalue) of N+
2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

6.1

Results for one-photon ionization cross sections compared to those from the literature. F is root mean square of laser electric field in atomic units. All the calculations
above are performed for n = 0, 1 except for F = 0.05, in which n is −1, 0, 1, 2 due tothe strong coupling at high intensity. The superscripts of the reference values: (a)stand for [243] and (b) for [233] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

6.2

The exponents of complex Gaussian basis sets for different laser frequencies. Thesebasis sets are tested for the range of laser intensity from 7 × 108 W/cm to 1.8 ×
1014 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

6.3

Complex Gaussian exponents for two photon ionization . . . . . . . . . . . . . . . . 133

xii

FiguresFigure2.1

Density plots for H atom in a 800 nm laser field at intensity of 1 × 1014 W/cm2 withmask function (left) and imaginary absorbing potential (right). . . . . . . . . . . . . 21

2.2

Comparison of population from valence orbitals of N2 between LB94 (left) and LDA(right) functionals. Laser field is at wavelength of 400 nm and intensity of 5 ×
1013 W/cm2 . Only the relevant KS orbitals which possess an important response to
the laser field are shown with their labels. . . . . . . . . . . . . . . . . . . . . . . . . 22

2.3

Comparison of high harmonic spectra of N2 for LDA and LB94 functionals. Laser
field is at wavelength of 400 nm and intensity of 1 × 1014 W/cm2 with polarization
parallel to the molecular axis.

2.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Projection between coupled orbitals of N2 in a 400 nm laser field at intensity of
1 × 1014 W/cm2 for LDA and LB94 functionals. . . . . . . . . . . . . . . . . . . . . . 23

3.1

Three step model for high-order harmonic generation (left) and typical harmonic
spectra from noble gas atom (right), adapted from [144]. . . . . . . . . . . . . . . . . 27

3.2High harmonic spectra minimum due to the quantum interference from the wavefunction at the two nuclei in H+
2 . The laser field is at the wavelength of 800 nm andintensity of 3 × 1014 W/cm2 with pulse length of 20 cycles. Molecular orientation
angle is 30 degrees with respect to polarization of the laser field. . . . . . . . . . . . 33

xiii
3.3

The high harmonic spectra of N2 from its two outermost states HOMO (3σg ) (bluedash) and HOMO-1 (1πu ) (green dots). The total spectra from these two orbitals arealso plotted (red solid line). The wavelength of the laser field is 800 nm at intensityof 2 × 1014 W/cm2 with polarization perpendicular to the molecular axis. . . . . . . 34

3.4

Energy level diagram (left) and orbitals (right) of the valence shell of N+
2 . . . . . . . 36

3.5Time dependent dipole for HOMO-2 of N+
2 for the case when the molecular axis
is perpendicular (a) and parallel (b) to the laser electric field at wavelength of 400
nm and intensity of 2 × 1014 W/cm2 . The amplitude of the laser field (on arbitraryscale) is shown by dashed line for comparison. The Fourier transform spectrum of
the residual dipole after the laser pulse for parallel case (b) is plotted in (c). . . . . . 37

3.6High harmonic spectra of N+
2 , aligned parallel (a) and perpendicular (b) to thepolarization axis of a laser pulse at 400 nm and 2 × 1014 W/cm2 . The insets show
an enlarged view of part of the spectrum, which clearly exhibit the first and second
order Mollow sidebands displaced by the Rabi frequency Ωr and 2Ωr from the odd
harmonics for the aligned molecules. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.7

Illustration of Mollow triplets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.8

Comparison of the distance of the sideband positions, i.e. first and second order
fractional harmonics, with respect to the corresponding odd harmonic (averaged
over harmonics bellow 11th) and the Rabi transition frequency for different laser
field intensities. Laser wavelength of 400 nm and orientation of molecule parallel tothe laser field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.9

HHG spectra from N+
2 for laser field at wavelength of 600 nm (left) with intensity of
1×1014 W/cm2 and at wavelength of 800 nm (right) with intensity of 2×1014 W/cm2 .
Molecular axis is oriented parallel to the polarization of laser field. Pulse length is
20 fs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

xiv
3.10 High harmonic spectrum of N+
2 aligned perpendicular to the polarization axis of 60
fs laser pulse at wavelength of 800 nm and intensity of 2 × 1014 W/cm2 . The inset
provides an enlarged view of the part of the spectrum.

. . . . . . . . . . . . . . . . 43

3.11 Time dependence of the projection between the coupled orbitals of N+
2 during theinteraction with the laser pulse for the case of parallel orientation (left) and perpendicular orientation (right). Laser field parameters: wavelength of 800 nm andintensity of 2 × 1014 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.12 Comparison of harmonic spectra for N+
2 calculated within modified SFA and TDDFT.
Laser parameters: λ=400 nm and I = 1 × 1014 W/cm2 . Polarization parallel to themolecular axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.13 Harmonic spectra from HOMO (a) and HOMO-2 (b) of spin down component and
HOMO-2 (c) of spin up component and population of coupled orbitals (d) (ionizationfrom other orbitals is negligible). Laser field is at 400 nm and 2 × 1014 W/cm2 withpolarization parallel to the molecular axis.

. . . . . . . . . . . . . . . . . . . . . . . 46

3.14 Difference between the full time-dependent electron density of N+
2 , aligned along(left) and perpendicular (right) to the polarization of laser, and the initial electron
density as a function of time and the position along the molecular axis. The density
is integrated over the spatial dimensions transversal to the molecular axis. Laserparameters as in Figure 3.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.15 Wavelet analysis of high harmonic generation in N+
2 aligned perpendicular to thepolarization direction of laser field. Laser parameters are the same as in Fig. 3.6. . . 50
3.16 Wavelet analysis of high harmonic generation in N+
2 aligned parallel to the polarization direction of the laser field. Laser parameters are the same as in Fig. 3.6. . . . . 51
3.17 Orbital image of CO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.18 Fractional harmonics from CO+
2 with laser field at wavelength of 350 nm and intensityof 1 × 1014 W/cm2 at parallel (left) and perpendicular orientation (right). . . . . . . 53
3.19 Orbital images of C2 H4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

xv
3.20 Fractional harmonics from C2 H+
4 with laser field at wavelength of 400 nm and intensity of 1 × 1014 W/cm2 with polarization parallel (left) and perpendicular (right) to
the C-C bond.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.21 Orbital images of NO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.22 Fractional harmonics from NO2 with laser field at wavelength of 400 nm and intensityof 2 × 1014 W/cm2 with polarization direction parallel to the oxygen-oxygen axis.

. 55

3.23 The amplitudes of parallel and perpendicular components (left) and phase differenceδ (right) of 57th harmonic order of H+
2 as a function of the alignment angle in a laser
field at wavelength of 800 nm and intensity 3 × 1014 W/cm2 . . . . . . . . . . . . . . 58
3.24 The amplitude ratio r (left), phase difference δ (middle), and the ellipticity (right) of
high-order harmonics from H+
2 in a laser field at wavelength of 800 nm and intensity
3 × 1014 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.25 The same as Figure 3.24 but for H2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.26 Illustration diagram for angle average procedure. . . . . . . . . . . . . . . . . . . . . 61
3.27 Ellipticity of high-order harmonics as a function of the harmonic orders for threemolecular alignment angles at λ = 800 nm and peak intensity of Io = 1×1014 W/cm2 .
Comparison of results in which no average (left panel) and a proper average (middle
panel) over the distribution of the alignment in the molecular ensemble has beentaken into account and experimental results (at 2 × 1014 W/cm2 )(right), adaptedfrom [165]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.28 Ellipticity from N2 with 800 nm laser field at intensity of 8 × 1013 W/cm2 (left) and1 × 1014 W/cm2 (right), adapted from [154]. . . . . . . . . . . . . . . . . . . . . . . . 64
3.29 Ellipticity of high-order harmonics as a function of the molecular alignment angle andthe order of the harmonics at λ = 800 nm and peak intensity I0 = 7.5×1013 W/cm2 .
Comparison of results in which no average (left panel) and a proper average (rightpanel) over the distribution of the alignment in the molecular ensemble is made. . . 65

xvi
3.30 Mechanism studied in [154]. Electron is ionized from HOMO and recombines to
HOMO-1. The hole at HOMO-1 is created by the 800 nm laser pulse which resonantly
couples HOMO-1 and HOMO. Figure is adapted from [154]. . . . . . . . . . . . . . . 66
3.31 Ellipticity of harmonics as function of harmonic order (vertical axis) and angle between molecular axis and laser polarization direction (horizontal axis). Comparison of results from the full calculation (upper left panel) and contributions from
HOMO, HOMO-1 and HOMO-2 (upper right panel), HOMO and HOMO-1 (lower
right panel) and HOMO only (lower left panel). Laser parameters: wavelength of
800 nm, peak intensity of 1 × 1014 W/cm2 and pulse duration of 20 fs. . . . . . . . . 67
3.32 Same as Fig. 3.31 but for a laser wavelength of 1000 nm.

. . . . . . . . . . . . . . . 68

3.33 Intensity ratio of perpendicular to parallel component. Comparison between experimental result (a) adapted from [165], and our simulation result (b). . . . . . . . . . 69
3.34 Illustration of high harmonic ellipticity orientation angle (a), experimental results ofellipticity orientation angle from CO2 (b), adapted from [165], and our simulation
results (c).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

3.35 Ellipticity of high harmonics as function of alignment angle for CO2 for a laser fieldat λ = 800 nm and I = 1.5 × 1014 W/cm2 . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.36 Ellipticity of harmonics as a function of alignment angle for CO2 . Starting with theresults from HOMO only (upper left), contributions from inner valence orbitals areadded step by step in the other panels (as indicated). A proper average over the
distribution of the alignment in the molecular ensemble is made. Applied laser field
is at wavelength of 800 nm and intensity of 1.5 × 1014 W/cm2 with 30 fs pulse length. 72
3.37 Projection of coupled inner valence orbitals HOMO-3 (4σg ) to HOMO-1 (3σu ) (a)and HOMO-5 (3σg ) to HOMO-4 (2σu ) (b). Laser wavelength of 800 nm and intensityof 2 × 1014 W/cm2 with orientation angle of 20o are considered. . . . . . . . . . . . . 73

xvii
4.1

Ionization from KS orbitals of N2 for parallel (a) and perpendicular (b) alignment
with respect to external laser field polarization. The laser field parameters: wavelength of 600 nm and intensity of 5 × 1013 W/cm2 . . . . . . . . . . . . . . . . . . . . 77

4.2+
Ionization of KS orbitals of openshell molecules: NO2 (a), CO+
2 (b), and C2 H4 (c).

The laser parameters and coupling orbitals are illustrated in the captions below eachfigure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3

The ionization probability of valence KS orbitals of N2 . Only the relevant KS orbitals
which possess an important response to the laser field are shown with their labels.
Left: Laser polarization is parallel to the molecular axis at 400 nm and intensity of 2×
1014 W/cm2 , inducing a coupling of HOMO (3σg ) and HOMO-2 (2σu ). The curves
for 1πux and 1πuy overlap because they are degenerate. Right: Laser polarizationis perpendicular to the molecular axis at wavelength of 800 nm and intensity of
2 × 1014 W/cm2 , inducing a coupling between HOMO (3σg ) and HOMO-1 (1πu ). . . 81

4.4

Comparison of ionization from N2 between LDA (left) and LB94 (right) functionals.
Laser field has a wavelength of 400 nm and intensity of 5 × 1013 W/cm2 with polarization parallel to the molecular axis of N2 . Only the relevant KS orbitals which
exhibit an significant response to the laser field are shown with their labels. . . . . . 82

4.5

Comparison of ionization from N2 between sin2 (left) and trapezoidal (right) envelope. Wavelength is 400 nm and intensity is 8 × 1013 W/cm2 with polarization
parallel to the molecular axis of N2 . Only the relevant KS orbitals which possess animportant response to the laser field are shown with their labels. . . . . . . . . . . . 83

4.6

Ionization of KS orbitals of N2 with moving nuclei. LB94 functional and trapezoidal
envelope are used in simulation. The laser field is set at wavelength of 400 nm andintensity of 5 × 1013 W/cm2 with polarization parallel to the molecular axis. . . . . . 84

xviii
4.7

Ionization yields of states of N2 at the time of zero electric field. Only the relevant
states which possess an important response to the laser field are shown with theirlabels. Left: Laser field is parallel to the molecular axis at 400 nm and intensity of
2 × 1014 W/cm2 inducing a coupling of HOMO (3σg ) and HOMO-2 (2σu ). Right:
Laser field is perpendicular to the molecular axis at 800 nm and intensity of 2 ×
1014 W/cm2 inducing a coupling of HOMO (3σg ) and HOMO-1 (1πu ).

4.8. . . . . . . 85

Upper pannels: Ratio of ionization yields between HOMO-2 and HOMO (left) and
HOMO-1 and HOMO (right) of N2 as a function of laser wavelength at intensityof 5 × 1013 W/cm2 and a pulse length of 14 fs. Lower panels: Ionization yields for
HOMO and inner valence orbitals as a function of intensity at wavelength of 400 nm(left) and 800 nm (right). Panels on left are for parallel alignment, while panels onright are for perpendicular alignment. . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.9

Instantaneous eigenvalues of Kohn-Sham Hamiltonian for N2 valence orbitals for 400
nm laser field at parallel orientation (left) with intensity of 5 × 1013 W/cm2 and 800
at perpendicular orientation (right) with intensity of 5 × 1013 W/cm2 . Only therelevant KS orbitals, which exhibit significant response to the laser field, are shownwith their labels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.10 The projection of propagated orbital 2σu (t) onto the coupled initial orbital 3σg (t =0). Laser field is at wavelength of 400 nm and intensity of 5×1013 W/cm2 for parallelorientation of molecular axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.11 Total ionization yield from N2 as a function of wavelength. The laser electric field isaligned along the molecular axis. Intensity is equal to 5×1013 W/cm2 . . . . . . . . . 89
4.12 Measured fluorescence spectra obtained from laser filaments in air and normalized
to the same interaction volume for λ= 800 nm (black) and λ= 400 nm (gray).
Fluorescence from the first band obtained for filaments at 400 nm is about 6.4 times
larger compared to that at 800 nm. Adapted from [184].

. . . . . . . . . . . . . . . 90

xix
4.13 The spectra captured with (blue) and without (green) the pump pulse. Laser parameters Adapted from [186]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.14 Ionization yield of N2 orbitals with two-color laser field: 400 nm with intensity of1012 W/cm2 and 800 nm with intensity of 1014 W/cm2 . Both laser field are inparallel to the molecular axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.15 Ionization of valence orbitals of CO2 at 350 nm laser field and intensity of 1 ×1014 W/cm2 in parallel to the molecular axis. Only the relevant KS orbitals which
possess an important response to the laser field are shown with their labels. . . . . . 93
4.16 Ionization of valence orbitals of C2 H4 (left) and projection between the coupled states(right). The laser field is at wavelength of 400 nm and intensity of 1 × 1014 W/cm2
with polarization perpendicular direction to the molecular axis but in parallel tothe plane where the molecule lies. Only the relevant KS orbitals which possess animportant response to the laser field are shown with their labels. . . . . . . . . . . . 94
4.17 Ionization of KS orbitals of C2 H6 (left) and the projection of coupled states. Laser
field parameters: 200 nm and 1×1014 W/cm2 , electric field is parallel to the carboncarbon molecular axis. Only the relevant KS orbitals which exhibit significant response to the laser field are shown with their labels. . . . . . . . . . . . . . . . . . . 95
5.1

Ionization and fragmention patterns represented by the mass spectra of hexatriene,
decatetraene and β-carotene for the interaction with laser pulse at 800 nm and 1450
nm. Adapted from [40]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

5.2

Multiburst ionization from H+
2 at internuclear distance of R=7.0 a.u. along themolecular axis in a strong laser field (shown in d) using (a) 1D model and (b, c) 2D
model. Adapted from [203] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

5.3

The electron density of HOMO-2 of N+
2 along molecular axis for different laser wavelength: 400 nm (left), 600 nm (middle), and 800 nm (right). The white curve is theelectron population transfered from HOMO-2 to HOMO.

. . . . . . . . . . . . . . . 100

xx
5.4

13 W/cm2
Electron density of HOMO-2 from N+
2 at different laser intensity: 1 × 10(a), 5 × 1013 W/cm2 (b), 1 × 1014 W/cm2 (c), and 3 × 1014 W/cm2 (d). Laser
wavelength is fixed at 400 nm.
5.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

The electron density plot from HOMO-2 of N+
2 in a two color laser field set up: 400
nm laser field at I = 1 × 1012 W/cm2 and 800 nm laser field at I = 1 × 1014 W/cm2 103

5.6

The projection from HOMO-2 to HOMO of N+
2 for two-color laser field: 400 nm
laser field at I = 1 × 1012 W/cm2 and 800 nm laser field at I = 1 × 1014 W/cm2 . . . 104

5.7

Orbital density plot of the coupled orbitals for different molecules. . . . . . . . . . . 106

6.1

Effect of complex coordinate rotation. Bound states and thresholds are invariant. As
continuum states rotate, complex resonance eigenvalues may be exposed. Adaptedfrom [227]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

6.2

Diagram of the Floquet program representing the relationship between classes. . . . 120

6.3

Distribution of exponents α of Gaussian basis sets (red vertical line) after optimization for the H atom in a laser field with ω = 0.6 (left), ω = 1.0 (middle), ω = 2.0√
2(right). The solid blue curve is the normalized coefficient α−3/2 e−s /4α for s = 2Ek
and the other two dashed blue curves are for electron energy of 2Ek and Ek .

6.4

. . . . 126

Quasienergy trajectory of H atom in a laser field at intensity of 4.387 × 1013 W/cm2
with frequency ω = 0.6. Results in the top two panels are for 3-complex Gaussianbasis sets while those in the bottom row are for 5-Gaussian basis sets. Figures on
the left show the imaginary part of quasienergy as a function of the rotation angle θand panels on the right show the trajectory of the quasienergy. A stationary point
is found with ε = −0.498 − 0.00109i a.u.. . . . . . . . . . . . . . . . . . . . . . . . . 127

6.5

Trajectory of quasienergy of H for different laser parameters. F is the mean squareroot of laser electric field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

xxi
6.6

Trajectory of quasienergy (a-c) for two photon ionization from H atom in a laser
field with ω = 0.3 and F = 0.025 and the value of exponents of complex Gaussianbasis after optimization are indicated by the vertical lines in (d). . . . . . . . . . . . 131

6.7

13 W/cm2 .
Cross section of two-photon ionization from H+
2 at laser intensity I=1.76×10

The solid curve is adapted from [211, 244] and the red dot is our calculations. . . . . 132

Chapter 1

Introduction

Since the invention of the laser in 1960s, it has been used extensively to study the responseof quantum systems such as atoms, molecules, and solids to the light field, which has extended
our knowledge about the mechanism of electron response to coherent radiation. At low lightintensity the light-matter interaction is dominated by the process of one- or two-photon absorption
or emission, which can be sufficiently analyzed with the lowest-order perturbation theory [1]. The
capability of producing high intense light fields over a wide frequency range in the form of short laser
pulses opened a new domain in the study of light-matter interaction. At intensities of 1013 W/cm2 ,the strength of the interaction between electron and light electric field can compete with that of
the interaction between electron and Coulomb field. Therefore, the interaction between quantum
systems and light field can no longer be treated as a perturbation to the system [2]. Due to theintense laser field, higher order processes, corresponding to the net absorption and emission of
more than one photon, called multiphoton processes [3], start to play a significant role. Novel
phenomena such as above-threshold ionization (ATI) [4], high-order harmonic generation (HHG)[5, 6], laser-induced correlated electron emission [7, 8] have been observed as a consequence of
nonlinear nonperturbative effects.
For the atomic system, the interaction with laser field has been studied extensively and isconsidered to be well understood. Several ab-initial theoretical techniques as well as systematic
approximation methods have been developed and the results agree well with experimental observations. The strong-field ionization from atoms can be modeled successfully by the strong-field2
approximation (SFA, also known as Keldysh-Faisal-Reiss (KFR) theory [9, 10, 11]) or by quasistatic tunneling models [9, 12, 13]. In these models, it is assumed that only one of the boundelectrons responds to the external field and all others are frozen in their initial states, which iscalled single-active-electron (SAE) approximation. Such assumption is justified by the fact of large
spacing between electronic energy levels in atoms. Another assumption behind these models isthat the active electron follows the oscillation of the electric field adiabatically and thus most of
strong-field phenomena are interpreted with quasi-static pictures.
Two important fundamental processes in the interaction between atomic systems and strong
field, namely ionization and high harmonic generation are studied very well. For ionization at
relatively low intensity of a light field, an electron can only be ionized if the energy of the absorbed
photon is larger than the ionization potential. When the intensity of the light field increases, multiphoton ionization becomes probable, in which the electron can absorb more than one photon to accumulate enough energy to reach the continuum. At sufficiently high intensities (I > 1011 W/cm2 ),it was found that the ejected electron can absorb photons in excess of the minimum number requiredfor ionization, which is called above-threshold ionization (ATI), as discovered by Agostini et al. [4]in 1979. If the laser field is sufficiently strong and the frequency is low enough, ionization can beinterpreted by using a quasi-static model. In this case, the bound electron experiences an effective
potential barrier formed by adding to the atomic potential the potential of the instantaneous laser
electric field. The ionization rate can be approximated by the static-limit tunneling rate as givenin [14]. The separation between multiphoton ionization and tunneling ionization can be associatedwith the Keldysh parameter γK [9] given by
s
γK =

Ip
,
2Up

(1.1)

where Ip is the ionization potential and Up is the average energy of a free electron in a laser field.
The multiphoton ionization is usually more likely to happen for γK > 1 while tunneling ionization
dominates for γK < 1.
Another important process in strong field is high-order harmonic generation (HHG). The

3interaction between atoms and sufficiently intense laser field can lead to the emission of radiation
at high-multiples, or harmonics, of the angular frequency ω of the laser field. Since the discoveryof harmonic generation in 1961 by Franken et al. [15], the study of high-order harmonic generation
has attracted considerable interest because it extends the range of laser light sources to shorter
wavelengths and thus can provide a source of very bright, short-pulse, high-frequency coherent
radiation, an attractive low-cost alternative to large scale facilities like synchrotrons and freeelectron lasers. Nowadays, the spectra of HHG can be extended to keV regime [16, 17]. The
ultrabroad bandwidth of the spectra makes it possible to produce attosecond pulses, which canin principle be used to detect electron dynamics directly at its natural time scale. This opened a
new research area: attosecond science [18, 19, 20]. Much of earlier work in experiment as well as
in theory in the nonperturbative regime of interaction in atomic systems with radiation have beensummarized in [21, 22, 23, 3, 24].
When moving from atoms to study molecules, because of the extra degrees of freedom (rotation and vibration), the multi-center nuclear frame, and a more complex electronic energy levelstructure, new phenomena such as Coulomb explosion, bond softening, vibration trapping, optically triggered explosion, above-threshold dissociation, and quantum interference were discovered[25, 26, 27, 28, 29, 30, 31]. Most of the previous theories developed for studying atoms needmodification to analyze the interaction of molecules with intense laser light. For example, electron
wavepackets emitted from different centers in molecules can interfere destructively leading to a suppression of molecular ionization [30, 31] and a minimum in the spectra of high harmonics [28, 29].
On the other hand, the single-active-electron (SAE) model can break down in the interaction between molecules and intense laser fields because the valence orbitals of molecules usually haveionization energies close to each other and thus nonperturbative processes may involve electrons
from more than one orbitals [32, 33, 34, 35, 36, 37]. Recently, the alignment [38] of molecules madeit possible to control the intersection angle between molecular axis and the polarization of laser
field and thus made it possible to study angle-dependent properties. The studies of ionization and
HHG dependence on the alignment angle reveal that such processes strongly depend on the sym-

4
metry of orbitals. As a result, molecular imaging techniques based on HHG were developed, which
has became an important topic, since Itatani et al. [39] have shown that HHG may be exploited
to tomographically image the highest occupied molecular orbital (HOMO) for diatomic molecules.
Besides this, quasi-static picture can also break down in the interaction between molecules andintense laser field, leading to nonadiabatic effects. Compared to atomic systems, the electronic
states in molecules are usually much closer to each other and thus a nonadiabatic behavior is easier
to induce and observe in molecules because of the coupling of states. The nonadiabatic behavior isalso easier to observe in large molecules because the time for an electron inside a large molecule totraverse is longer which is in favor of nonadiabatic response to external laser field [40, 41, 42].
Because of these complex phenomena in molecules, it is important to study the responseof molecules in intense laser field to help us understand the mechanisms behind these phenomena
better. In this thesis, we analyze multielectron effects in high-order harmonic generation, ionizationand nonadiabatic electron localization. Since the solution of multielectron Schrödinger equation for
the interaction between a molecule and a laser field at nonperturbative intensities is numericallynot achievable, we use the time-dependent density-functional theory. In particular, we study thenonlinear processes for scenarios where the laser field couples states inside the molecule. The thesis
is organized as follows:
In Chapter 2 we introduce the time-dependent density-functional theory (TDDFT). TDDFT
transforms the problem of solving multielectron Schrödinger equation to the problem of solving
single-particle Kohn-Sham equations, which can significantly reduce the computation time. Related
concepts like Hohenberg-Kohn theorem, Kohn-Sham equation, exchange-correlation functional, etc.
are discussed in this Chapter.
Chapter 3 is devoted to analyze multielectron effects in high-order harmonic spectra frommolecules. The basic three-step model and analytical strong-field-approximation (SFA) method
are introduced. Specifically, we first investigate the harmonic spectra from open-shell moleculesinteracting with a laser field, which resonantly couples two states in the molecule. Two new
effects are identified and analyzed. Rabi flopping induced by the laser field leads to fractional

5
harmonics in the form of Mollow sidebands over the whole spectrum of harmonics and modificationof the traditional three-step model for the HHG process. Furthermore, we analyze the ellipticityof harmonics from N2 and CO2 . Ellipticity patterns obtained from our simulations are in good
agreement with experimental results. By examining contributions from each valence orbital wesurprisingly find significant contributions of inner valence orbitals, which have not been reported
previously.
In view of the strong modification of the high-order harmonic generation process due to theresonant coupling, we analyze its influence on the ionization process in Chapter 4. Specifically, we
study the coupling of two occupied valence orbitals for molecules like N2 , CO2 , etc. and observe
significant enhanced ionization from inner valence orbitals. Consequently, the approximation thatinner valence electrons are treated as frozen cannot be applied in such cases. We further showthat our results can explain recent observations of population inversion in molecules in recent
filamentation experiments.
The coupling between states in strong field can further lead to nonadiabatic electron dynamics
in laser driven molecules. In Chapter 5, we study the mechanism behind the nonadiabatic dynamics
resulting from coupled states by analyzing the time-dependent changes in the electron density in thelaser field. Two competing processes which together lead to the nonadiabatic effects are analyzed
with respect to the different type of coupled orbitals as well as the laser intensity on the nonadiabatic
processes.
In Chapter 6, we introduce an intense-field method based on basis sets. To this end, weapply the Floquet theorem to time-dependent Schrödinger equation (TDSE) to transform the timedependent problem into a time-independent problem. To make the Floquet method capable forapplications to large molecular systems, we analyze the time-independent problem using Gaussianbasis sets. Complex Gaussian basis sets are used to represent the continuum states. An optimization
procedure for the exponents of the basis sets is developed. First results for one- and two-photon
cross sections are in excellent agreement with those of other calculations.
Finally, we present a summary of the thesis in Chapter 7.

Chapter 2

Time-dependent density-functional theory

The interaction with high intensity laser field has lead to the discovery of many new phenomena in atoms and molecules, which puts a demand for more advanced theoretical approaches
to help us understand the mechanism behind it. A real challenge in theoretical description is thatthe interaction of electrons with the time-dependent electric field cannot be described within perturbation theory because of the high field strength that is comparable to Coulomb field. On theother hand, direct numerical integration of Schrödinger equation is impossible for a system withmore than two electrons with current available high power computers. By now, most numericalsimulations have been performed within single-active-electron (SAE) approximation, where electrons are assumed to be independent and only one electron is considered to be ’active’ while all
the remaining electrons are frozen in their initial states. On the other hand, considerable progress
towards numerical solution of the six-dimensional Schrödinger problem for double ionization has
been achieved with state-of-the-art computations [43, 44, 45] in which two electrons can be ’active’.
Another theoretical approach for strong field problems is Intense-field Many-body S-matrix
Theory (IMST). The earliest forms of S-matrix theory that accounted for the non-perturbative
intense-field matter interaction, are analogs of the well-known Keldysh-Faisal-Reiss (KFR) model,
also called Strong-Field Approximation (SFA), in which the interaction with laser field is emphasized [10, 9, 11]. The IMST which has been developed to extend this model has been used to study
multiphoton ionization from outer and inner valence shells of diatomic and polyatomic molecules[46, 47, 30, 48, 49, 50]. However, e.g. the ionization saturation intensity for polyatomic molecules

7
calculated based on IMST only qualitatively agrees with experimental results. This suggests thationization from these molecules might show a breakdown of single-active-electron approximation.
To investigate such multielectron systems, quantum dynamical version of the configuration interaction (CI) and multiconfiguration self-consistent field (MCSCF) methods of electronic structure
calculations have been developed. These interesting approaches, time-dependent CI with single excitations (TD-CIS) [51, 52, 53] and time-dependent Hartree-Fock (TDHF) [54] as well as multiconfiguration time-dependent Hatree-Fock (MCTDHF) [55, 56, 57, 58], are able to treat the correlated
motion of many electrons in nonequilibrium situations. However, they were only applied to smallor artificial model molecules. Another promising method is time-dependent density-functional theory (TDDFT) [59, 60, 61, 62, 63, 64], which in principle has the ability to describe strong fieldprocesses for polyatomic molecules. Otobe and Yabana [59] have studied ionization of benzene
and acetylene with TDDFT for the first time. The method has also been applied by Castro et al.
[60, 65, 66, 67, 68] and Chu et al. [69, 70, 71, 72] for smaller molecules such as hydrogen molecule,nitrogen molecule, sodium dimer, and helium trimer.
The TDDFT method does not only provide a way to study multielectron system such aslarge molecules but also gives access to inner valence orbitals in order to investigate their roles ina specific mechanism. In the following sections, the basic concepts within TDDFT are going to bereviewed. More detailed information is available at [73, 74, 75, 76, 77, 78].

2.1

Multielectron system
In quantum mechanics, a multielectron system with N electrons can be described fully by

its wavefunction Ψ(t, r1 , r2 , · · · , rN ) which satisfies the multielectron Schrödinger equation in timedependent form



N 
N
2 ∇2
X
X
∂
~
i
i~ Ψ(t, r1 , r2 , · · · , rN ) = −+ v(ri ) + v laser (t, ri ) +U (ri , rj ) Ψ(t, r1 , r2 , · · · , rN ),
∂t
2m
i=1

i<j=1

(2.1)

8
or time-independent form



N 
N
2 ∇2
X
X
~
i

−
+ v(ri ) +U (ri , rj ) Ψ(r1 , r2 , · · · , rN ) = EΨ(r1 , r2 , · · · , rN ),
2m
i=1

(2.2)

i<j=1

where U (ri , rj ) is the electron-electron interaction which has the form
Û =

N
X

N
X

U (ri , rj ) =

i<j=1

i<j=1

e2
|ri − rj |(2.3)

for a Coulomb system, where we defined an operator Û for the electron-electron interaction. Similarly, we can define a kinetic energy operator for the electrons
T̂ =

N
X

~2 ∇2i2m(2.4)

X −Qk e|ri − Rk |

(2.5)

−

i=1

and nuclear-electron interaction energy
V̂ =

N
X
i=1

v(ri ) =i,kwhere Qk and Rk are the nuclear charge and nuclear position of k th nuclei, respectively. Theinteraction between electron and external laser field is included in the term v laser (t, ri ). In the timedependent multielectron Schrödinger equation (2.1), we ignored the motion of the nuclei and its
interaction with the laser field, since the time-scale of the motion for nuclei is very large compared tothe motion of electrons. Thus, the whole wavefunction that describes both nuclei and electrons canbe decomposed into two parts: one for nuclei and one for electron, which leads to Eq. (2.1) satisfied
by the electron wavefunction. Such approximation is called Born-Oppenheimer approximation [79].
Although the computation power has been improved dramatically over the last decades, it is still
impossible to solve the multielectron Schrödinger equation (2.1) directly for more than two-electron
systems in realistic time. An alternative approximation method is TDDFT, an extension of densityfunctional theory (DFT), which provides rather precise results within realizable computing time.
Before going to TDDFT, we first look at the Hohenberg-Kohn theorem which establishes thefoundation of DFT.

9

2.2

The Hohenberg-Kohn theorem
DFT is based on the Hohenberg-Kohn (HK) theorem [80] which states that the ground stateof a system is uniquely determined by its ground state one-electron densityˆˆˆ
no (r) = N dr2 dr3 · · · drN Ψ∗o (r, r2 , r3 · · · rN )Ψo (r, r2 , r3 · · · rN ),(2.6)

where the subscript o refers to ground state. This means Ψo (r, r2 , r3 · · · rN ) is a functional of thedensity and thus all ground-state observables are functionals of the one-electron density no (r).
Hohenberg and Kohn showed in their paper [80] that if Ψo was not determined by the groundstate density no , then there is a contradiction to the variational principle. The HK theorem has
been scrutinized carefully over the years, for example, M. Levy [81] and E. Lieb [82] proposed
independently constrained-search proof of the theorem. A strong form of the Hohenberg-Kohn
theorem was discovered and given by [83, 84, 85]ˆdr∆n(r)∆v(r) < 0,(2.7)

which means we cannot keep either the change of external potential ∆v(r) or the change of groundstate density ∆n(r) to be zero when the other is changed. So, the HK theorem gives the following
conceptual relationshipno (r) → v(r) → Ψo (r1 , r2 , · · · , rN ),(2.8)

which says that the ground density uniquely determines the potential which uniquely determines
the ground state of the system. Based on the above information, the ground state wavefunction isa functional of the ground state densityΨo (r1 , r2 , · · · , rN ) = Ψ[no (r)](2.9)

and thus any observable of the ground state is also a functional of no (r)O[no (r)] = hΨ[no (r)]|Ô|Ψ[no (r)]i.

(2.10)

One of the important observables is total energy
E[no ] = T [no ] + U [no ] + V [no ].

(2.11)

10

2.3

Kohn-Sham equation
Since the system is uniquely determined by the electron density, it is sufficient to determine

the correct electron density to calculate any interesting physical quantities. The electron density canbe obtained by solving the multielectron Schrödinger equation. However, if it would be possible
to obtain it from an auxiliary single-particle Schrödinger equation, then the computation cost
would dramatically decrease. Fortunately, such auxiliary equation exists and it is called KohnSham equation. Before the derivation of the Kohn-Sham equation, let us first consider a system ofnoninteracting electrons in external potential vs (r).
Ĥ = T̂ + V̂sˆV̂s =drn(r)vs (r)(2.12)The ground state wavefunction Ψo of the N-particle system can be approximated by a Slater
determinant φ1 (r1 ) · · ·1..
Ψo = √ det .
N!

φ1 (rN ) · · ·φN (r1 ) 


..

.


φN (rN )(2.13)

where orbital φi (ri ) is the solution of single-particle Schrödinger equation for the ith particle−~2 2∇ φi (ri ) + vs (ri )φi (ri ) = εi φi (ri ).
2m(2.14)For simplicity, the subscript i presents both spatial and spin quantum numbers. The ground state
density of the N-noninteracting-particle system is given bynos (r) =NX|φi (r)|2 ,(2.15)

i=1

in which, subscript s is added to differentiate it from the true interacting system. The total energyof the system is
Eos =

N ˆ
X
i=1

drφ∗i (r)(−~2 ∇2 )φi (r) +2m

ˆ
drvs (r)nos (r).

(2.16)

11
According to HK theorem, we already know Eos and the second term on the right of Eq. (2.16)
is a functional of density, then the term for kinetic energy should also be an functional of density.
This implies that the orbital is a functional of the ground state energy
φi = φi [nos ].

(2.17)

Now, assume that for a real interacting system with external potential vext (r), its ground state
density no (r) can be obtained from the ground state density of some noninteracting system with asingle-particle external potential vs (r) (related but different from vext (r)). In this way, it transforms
the problem from solving a multielectron Schrödinger equation to several non-interacting singleparticle Schrödinger equations, and thus it significantly reduces the computation time. The question
whether such a corresponding single-particle potential vs (r) exists for any arbitrary density is called
v-representability problem and has been early recognized [86]. This is a very important question
because HK theorem only guarantees that there cannot be more than two potentials for each density
but it does not guarantee that such a potential exists for a given density. Fortunately, it is known
that for discretized system the density is ensemble v-representable. A density is said to be ensemble
v-representable, if there is a mix of degenerated ground states for some potential which has thatdensity [87, 88, 89].
For an interacting system, we rewrite the exact energy of the system
E[no ] = T [no ] + U [no ] + V [no ]= Ts [no ] + UH [no ] + Exc [no ] + V [no ],(2.18)

where Ts is the total kinetic energy for non-interacting system defined as
Ts [no ] =N ˆX

drφ∗i (r)i(−~2 ∇2 )φi (r),
2m(2.19)no (r)no (r0 ).
|r − r0 |(2.20)and the Hartree energy UH is defined ase2
UH [no ] =
2

ˆ

ˆ
dr

dr0

12
The exchange-correlation energy Exc contains the difference Tc = T − Ts , which is called correlationenergy, and the difference U − UH . The exchange-correlation energy is usually decomposed into
two parts Exc = Ex + Ec , where Ex , the exchange energy, is due to the Pauli principle and Ec ,the correlation energy, is due to correlation (Tc is only a part of Ec ). The exchange energy can be
expressed in terms of the single-particle orbitals
q2 X
Ex = −
2

ˆ

ˆ
dr

dr0

jk

φ∗j (r)φ∗k (r0 )φj (r0 )φk (r)|r − r0 |(2.21)

which has no general exact expression in terms of density. The derivative of exchange-correlationenergy with respect to density is called exchange-correlation potential
vxc [no ] =

δExc
δno

(2.22)and the exact form of this potential is unknown. By now, there exist hundreds of exchangecorrelation functionals for different systems [90].
Based on the variational principle with respect to electron density, we have0=

δE[no ]=
δno (r)=δTs [no ] δV [no ] δUH [no ] δExc [no ]+++δno (r)δno (r)δno (r)δno (r)δTs [no ]+ vext (r) + vH (r) + vxc (r),δno (r)(2.23)

δV
H
is the external potential in which the electrons move and the term vH = δU
where vext (r) = δnδnois the Hartree potential. Similarly, for a system of non-interacting particles moving in a potentialvs (r),0 ==δEs [no ]δTs [no ] δVs [no ]
=
+
δn(r)
δno (r)δno (r)δTs [no ]+ vs (r).
δno (r)(2.24)

Comparing Eq. (2.23) and Eq. (2.24), if we set
vs (r) = vext (r) + vH (r) + vxc (r),(2.25)

then, both minimizations should yield the same density n(r) = ns (r). Therefore, the density of theinteracting many-particle system in a potential vext (r) can be obtained by solving the equations

13of a non-interacting system in potential vs (r), if we know the exchange-correlation potential. The
orbitals of the auxiliary non-interacting system satisfy the Kohn-Sham (KS) equation
 2 2

~ ∇
−
+ vs (r) φi (r) = εi φi (r)2m(2.26)and the density of the original system is obtained via
n(r) = ns (r) =NXfi |φi (r)|2 ,(2.27)

i=1

where fi is the electron population in ith orbital.
Since the Hartree potential vH and the exchange-correlation potential vxc depend on thedensity n and, thus, depend on the solution φi ; while the solution φi is determined by the potential
vs , the KS equation (2.26) must be solved self-consistently. It usually starts with an initial guess
for n(r), followed by calculating the potential vs , and then obtaining the solution φi , which will
yield a new density n(r). With the new density, we can start the process again until the result
converges.
Once the converged ground state no is obtained, we can calculate any interesting quantity
like the total energy
Eo =

N
X
i

e2
εi −
2

or the total dipole of the system

ˆ
dr

dr0

no (r)no (r0 )|r − r0 |(2.28)

ˆ
d=

2.4

ˆ

drno (r)r(2.29)

Time-dependent density-functional theory
So far, we have considered a static external potential only. For a system subject to a time-dependent external potential, e.g. that of a laser field, the system is described by the time-dependent
Schrödinger equationi~∂Ψ(r, t) = (T̂ + Û + V̂ )Ψ(r, t)∂t(2.30)The kinetic energy T̂ is the same as in Eq. (2.4) and the electron-electron interaction Û is definedin Eq. (2.3). The external potential V̂ not only includes the nuclear potential V̂ne as defined in Eq.

14(2.5) but also includes the external electric field which is typically of the form
vLaser = f (t)E sin(ωt)

NXri · α,(2.31)

i=1

where α, ω, and E are polarization, frequency, and electric field amplitude of the laser. The function
f (t) is called pulse envelope and describes the shape of the laser pulse in time. For Eq. (2.31),
there are two assumptions applied:• The laser field is treated classically, which is valid as long as the density of photons islarge and the individual quantum nature of photons can be ignored. For all simulations
presented in this thesis, this is the case.
• The dipole approximation is applied, in which we assume that the wavelength of the laser
field is much larger than both the size of the system and the path that is traveled byelectron in one period of the laser field. Also this assumption is valid for all simulations
presented in this thesis.
Similar to the static case, although the time-dependent Schrödinger equation provides anaccurate description of physics, it is, however, impossible to solve it directly for multielectron
system. It is, therefore, natural to ask whether the Hohenberg-Kohn theorem holds for the timedependent case as well. Fortunately, Runge and Gross [61] have proven that if two potentials v(r, t)
and v 0 (r, t) differ by more than a purely time-dependent function c(t), then the time-dependent
density produced with these two potentials cannot be the same:v(r, t) 6= v 0 (r, t) + c(t) ⇒ n(r, t) 6= n0 (r, t),(2.32)

which implies that there exists a one-to-one mapping between the density and the external timedependent potential. Again, such density is not obtained directly but with the auxiliary timedependent Kohn-Sham equation
 2 2

∂
~ ∇
i~ φi (r, t) = −+ vks (r, t) φi (r, t),∂t
2m(2.33)

15
where the Kohn-Sham potential is defined as in Eq. (2.25) except that the external potential
includes both the nuclear potential in Eq. (2.5) and the laser field vlaser ,vext = vnu + vlaser .

(2.34)The time-dependent density is then obtained from the time-dependent Kohn-Sham orbitalsn(r, t) =NX

fi |φi (r, t)|2(2.35)

i=1

2.5

Exchange correlation functionals
As shown above, the Kohn-Sham equation is simple and similar to a single-particle Schrödinger

equation. However, it contains a term νxc (r) which takes account of the multielectron effects, although its explicit formula does not exist. To make DFT and TDDFT a practical theory, different
approaches are proposed to design density functionals [91]:• Local density approximation (LDA). This group of functionals has the form
ˆ
LDA
Exc
[n] =

exc [n(r)]dr,(2.36)

where exc (n) is function of electron density only. The analytic formula is derived throughanalysis of the exchange-correlation energy in a system of uniform electron gas and isapplied to the case of non-uniform densities directly or with empirical modifications [92,93, 94, 78, 95].
• Density-gradient expansion (DGE). LDA is exact for a uniform electron gas and quite
accurate for solids but it is less satisfactory for atoms and molecules. The natural stepbeyond LDA is a formal expansion of Exc in gradients of the density, which was initially
proposed by Hohenberg and Kohn [15]. The general form of such group of functionals is
given by
DGE
Exc
[n] =

ˆ h
i
(1)
(2)
2
e(0)
xc (n) + exc (n)∇n + exc (n)|∇n| + · · · dr.

(2.37)

16
The derivation of the coefficients ekxc (n) is very mathematically involved [96, 97, 98, 99].
Since it is based on Taylor expansion, these functionals can only be applied for systems
with very slowly-varying densities.
• Constraint satisfaction. This approach is to apply some constraints to the property of theexchange-correlation potential. DGE gives only modest improvements over LDA and theDGE functionals have the wrong analytic behavior of the exchange-correlation potential.
The failure of the second-order DGE (divergence of the exchange-correlation potential,
overcorrection of the correlation energy, etc) were analyzed in detail in [100, 101, 102].
Accordingly, functionals with desired analytic properties; such as the asymptotic behaviorof exc and νxc , upper and lower energy bounds, density scaling transformation, etc. [103,104, 102] were designed, which have the general form,ˆExc [n] = exc (n, ∇n, ∇2 n, τ, · · · )dr.

(2.38)

To eliminate the asymptotic divergence of vxc (r), various modifications of the second-order
DGE have been proposed. Most of them fall into the category of generalized gradientapproximations (GGA) [103, 105, 106].
• Modeling the exchange-correlation hole. This type of functionals are designed based on
the exchange-correlation hole model, forming one of the largest and most diverse groups.
The general analytic form is the same as in Eq. (2.38). The exchange-correlation may beapproximated by a truncated Taylor series expansion in ρ, ∇ρ, ∇2 ρ or modeled from a
problem that can be analytically solved [102, 107, 104].
• Empirical fits. Such group of functionals are devised by fits to experimental results foratoms and molecules with analytic forms of Exc [ρ] which might be borrowed from otherfunctionals or simply postulated without derivation [108, 109]. The general form of thistype of functionals isˆEex [ρ] =

X
k

Ck

2
e(k)
xc (ρ, ∇ρ, ∇ ρ, τ, · · · ; ak , bk , · · · )dr,

(2.39)

17
where Ck ,ak ,bk ,· · · are fitting parameters.
• Mixing exact and approximation exchange. These functionals are also called hybrid functionals and have the formˆ
hybridExc[ρ] = exactTTaex (r) + b(eDF(r) + eDF(r)) dr,xc(2.40)

where the mixing coefficients a and b either depend on r or are kept constant [110, 111].
There exist hundreds of density functionals so far designed for the calculations of different
systems, from atoms to molecules, from liquids to solids. There is no universal functional, so
appropriate functional must be selected for a specific problem. For our simulations, we have done
systematic studies to compare results for different functionals and found that LDA yields goodresults and is also very fast compared to application of other functionals. For more accurate
results, for example in the case of ionization from atoms and molecules, it is necessary to chosefunctionals with correct asymptotic tail, 1/r, such as LB94 [112, 113] functionals. The following
subsection is a brief introduction to LDA and LB94 functionals.

2.5.1

LDA and LB94 functionals

As mentioned above, the local density approximation (LDA) has the form
ˆ
LDA
Exc
[ρ] =

ρ(r)νxc [ρ]dr,(2.41)

where the exchange-correlation potential νxc (ρ) is often separated into two parts: exchange potentialνx (ρ) and correlation potential νc (ρ). They can be obtained from the theory of a uniform electrongas.
The details about the derivation of the exchange energy for a uniform electron gas are analyzed by Bloch [92] and Dirac [93]. The exact LDA exchange formula for a uniform electron gas
isˆ
ExLDA [ρ] = −Cxρ4/3 (r)dr,(2.42)

18
where Cx = 43

3 1/3
. The correlation energy is more complicated than the exchange energy becauseπthe exact analytic form of the correlation energy can only be derived for two limiting cases. The
first one is the high-density (weak correlation) limit of a spin-compensated uniform electron gas
νc = AGB ln rs + B + rs (C ln rs + D),(2.43)

where rs is defined as

rs =


3 1/3
ρ
 1,
4

(2.44)

which is interpreted as the radius of a sphere that contains the charge of one electron. GellMann and Brueckner [114] evaluated the coefficients AGB and B, while Carr and Maradubin [115]evaluated C and D. For the low-density (strong correlation) case, Nozières and Pines [116] and Carr[117] obtained the correlation energy in the form of1
νc (rs ) =
2

Uo
U1
U2
+ 3/2 + 2 + · · ·
rs
rs
rs

!
, rs  1,(2.45)

where Uk are known constants. The exact form for intermediate densities is unknown and the value
was obtained through Monte Carlo simulations for uniform electron gas by Ceperley and Alder[118]. In order to connect the low density and high density form as well as to reproduce the value
for the intermediate density range, Perdew and Zunger [119] proposed the following fitting formula,


γ


if rs ≥ 11/2
νc (rs ) = 1 + β1 rs + β2 rs,(2.46)A ln rs + B + Crs ln rs + Drs if rs < 1
where γ, β1 , β2 , A, B, C, and D are parameters. There also exist some other fitting formulas
proposed by Vosko et al. [120] and Perdew et al. [106]. In our simulations, Eq. (2.46) is used andwe obtain high accuracy results in comparison with experiments.
One shortage of LDA is the fast decay of atomic potential which does not have correct asymptotic behavior, namely 1/r. This is however important to obtain accurate results for ionization.
To improve it, Leeuwen and Baerends [112] proposed a model of exchange potential that has thecorrect 1/r potential tail:LBνxc(α, β; r) = ανxLDA (r) + νcLDA (r) −βx2σ (r)ρ1/3 (r),1 + 3βx(r) ln{x(r) + [x2 (r) + 1]1/2 }(2.47)

19
where the last term in Eq. (2.47) is the gradient correction with x(r) = |∇ρ(r)|/[ρ(r)]4/3 , which
ensures the proper asymptotic behavior. The parameters α = 1 and β = 0.05 were obtained by a
fit to the exact exchange-correlation potential for the beryllium atom.

2.6

The pseudopotential approximation
For large systems, like large molecules and solids, the concept of pseudopotential (PP) is

usually applied so that only the wavefunction for valence electrons is explicitly calculated. The
idea is that only valence electrons participate in chemical bonding in molecules and solids, while
the inner electrons remain as an atomic-like configuration and their orbitals do not change too much,
if the environment of the atom is changed. Such technique leaves only the valence electron density
to be determined in the self-consistent cycle for the system which greatly reduces the computationtime.
The effective potential in the original KS equation is vs [n] = vext [n] + vH [n] + vxc [n], where
the density is the full electron density n(r) of the system. Using the pseudopotential, the electron
density is replaced with the density of valence electrons nv (r). The effect of core electrons is taken
P P . Therefore,
into account by replacing the external potential vext with a pseudopotential vext

PP
vsP P [nv ] = vext
+ vH [nv ] + vxc [nv ].

(2.48)

P P can be obtained by applying the following steps:
The pseudopotential vext(1) Determine an effective pseudopotential vsP P so that the single-particle wavefunction fromvsP P agrees with the valence wavefunction obtained from the all-electron KS equations
outside a cut-off radius rc for the same atom. In other words, the valence orbital densities
obtained through these two approaches are the same.
(2) Then the pseudopotential for the atom is equal toPP= vsP P − vH [nv ] − vxc [nv ].
vext(2.49)

20
There are many methods to generate the effective pseudopotential vsP P . The most popularmethods were proposed by Troullier and Martines [121], Bachelet, Hammann and Schlüter[74], Kleinman and Bylander [122], and Vanderbilt [123]. Nowadays one can generate thepseudopotential with these methods through the open source program atomic pseudopotential engine (APE) [73].

2.7

Implementation of TDDFT
To solve the Kohn-Sham equation (2.33), the wavefunction is discretized in both space andtime with uniform spacing step ∆x and time step ∆t, which converts the Kohn-Sham equation intoa matrix equation. The initial wavefunction can be obtained by solving the eigenvalue problem selfconsistently. An initial guess is needed to start the self-consistent process. For molecules, the initial
Kohn-Sham orbital is usually obtained as a linear combination of atomic orbitals. Such atomic
orbitals usually implemented are pseudo-orbitals that are used to generate the pseudopotential[124].
Once the initial ground state orbitals are obtained, the wavefunction of each orbital is propagated forward with the method called enforced time-reversal symmetry (ETRS),
φn+1 (r) = exp(−i∆tHn+1 /2)exp(−i∆tHn /2)φn (r).

(2.50)

It depends on the term Hn+1 which can only be known exactly via the solution φn+1 . The way to
avoid such problem is to estimate the exact solution by performing φ∗n+1 = exp(−i∆tHn )φn andthen set Hn+1 = Hn+1 [φ∗n+1 ]In general, in our simulations, the propagation time step and spatial grid spacing are set to0.03 a.u. and 0.3 a.u., respectively. There is a relationship for the propagation time step ∆t andspatial spacing ∆x which must be satisfied in order to get stable results
2

(∆x)2
> π.
∆t

(2.51)

The condition is based on the following argument. The maximum energy determined by the grid

2
1
π 2
π
spacing is kmax
2 = 2 ∆x , while the maximum energy determined by the time step is Emax = ∆t ,

21
2

and it requires Emax > kmax
to keep the propagation to be stable which yields the relationship of2
Eq. (2.51).
When a laser field is applied, an atom or molecule can be ionized and some part of the
wavefunction then propagate toward the boundary which will lead to reflection and unphysicaleffects. Therefore, such ionized wavepacket from the system should be absorbed once it reaches
the boundary of the simulation box. There are two types of absorbing boundaries. One procedure
is to multiply the wavefunction at each time step with a so-called mask function. The maskfunction is a function with value 1 for the region outside the absorbing region, decaying from 1
to 0 inside the absorbing boundary. The other method is to add an imaginary potential in theabsorbing boundary region. The wavefunction inside the region of imaginary potential is thendamped and thus absorbed. To compare the effect of these two boundaries, we have performed
simulations for 1-D H atom in a 800 nm laser field with intensity of 1 × 1014 W/cm2 using these
two boundaries. In both cases, the simulation box size and absorbing boundary width are set to
92 a.u and 13 a.u., respectively. Results are presented in Fig. 2.1. For mask absorbing boundary(left), the wavefunction shows a strong interference pattern due to the interference between outgoing wavepacket and reflected wavepacket from boundary while for imaginary absorbing potential(right), there is no interference pattern because the out-going wavepacket is well absorbed without
reflection.

Figure 2.1: Density plots for H atom in a 800 nm laser field at intensity of 1 × 1014 W/cm2 withmask function (left) and imaginary absorbing potential (right).

22
The pseudopotential we used is the Martine-Troullier pseudopotential [121]. The simulation
with TDDFT is carried out in two stages: first, we perform a DFT calculation to obtain the
ground state, and then we start TDDFT calculation by propagating each Kohn-Sham orbital in
the laser field. Although self-consistent calculations should be carried out for each time step, these
are actually not performed at each time step due to the computational costs. The results are
however checked by comparing with those for smaller time step. In each time step, information likeeigenvalue, laser induced dipole, population, projection of each orbital are calculated and stored
for analysis.

2.7.1

Comparison of exchange-correlation functionals

In this thesis, the physical quantities we are interested in include spectra of high harmonic
generation and ionization yield. So, we mainly compared results for this two properties obtainedwith LDA and LB94 exchange-correlation functionals. As an example, we first compare the ionization yield from N2 at 400 nm for these two functionals in Fig. 2.2. The LB94 functional has acorrect tail, 1/r, and thus the ionization is not as strong as that for the LDA functional. However,
the relative ionization yields from valence orbitals are nearly the same for both functionals. So
unless we are interested in comparing ionization yield with experimental results, LDA is already

2.000
1.995
1.990
1.985
1.9800

2σu
3σg

4

8 12
Time [fs]16

Population

Population

good enough to give us insights in a physical mechanism. The high harmonic generation spectra

2.000
1.995
1.990
1.985
1.9800

2σu
3σg

4

8 12
Time [fs]16

Figure 2.2: Comparison of population from valence orbitals of N2 between LB94 (left) and LDA(right) functionals. Laser field is at wavelength of 400 nm and intensity of 5 × 1013 W/cm2 . Only
the relevant KS orbitals which possess an important response to the laser field are shown with theirlabels.

23
obtained with LDA and LB94 functional are also compared for the same laser parameters, as shown

HHG yield (arb. units)in Fig. 2.3. Although the ionization yields from the two functionals have different values, the am-

10-3

LB94
LDA

10-5
10-7
10-9 02

4 6 8 10 12
Harmonic orderFigure 2.3: Comparison of high harmonic spectra of N2 for LDA and LB94 functionals. Laser
field is at wavelength of 400 nm and intensity of 1 × 1014 W/cm2 with polarization parallel to themolecular axis.

plitudes of each high harmonic order are nearly the same. As we will discuss later, at wavelengthof 400 nm, the laser field will induce a coupling between HOMO-2 (2σu ) and HOMO (3σg ) for
parallel orientation. We therefore also looked at the projection between the two coupled orbitals
to see if the two different functionals give us the same Rabi flopping behavior. The comparison
is presented in Fig. 2.4. For both functionals, we see an oscillation of projections with nearly­®
| 3σg (t =0)|2σu (t) |2

1.0
0.8
0.6
0.4
0.2
0.00

LB94
LDA

5

10 15
Time [fs]20

Figure 2.4: Projection between coupled orbitals of N2 in a 400 nm laser field at intensity of
1 × 1014 W/cm2 for LDA and LB94 functionals.

the same period. The small difference in the period comes from a slightly different eigenvalue oforbitals calculated with the two functionals. Based on these comparisons, we conclude that LDA

24
in general , is sufficient for the calculation presented in this thesis, since we are mainly interested in
identifying multielectron effects in laser induced process. If it turns out to be necessary to computeusing LB94 functional, we will show the comparison in detail between these two functionals in thecorresponding Chapter.

Chapter 3

Multielectron effects in high-order harmonic generation from molecules

High-order harmonic generation (HHG) is one of the highly nonlinear nonperturbative processes in intense laser field science. It was first observed for neon gas by McPherson et al. [5] in1987. Typically, HHG is generated by focusing an ultrashort, high-intense, near-IR laser pulse into
a noble gas at intensities of about 1014 W/cm2 . The intensity spectra of high harmonics emitted
by an atom reveal some characteristic features. Perturbation theory predicts that the intensitiesof the first few harmonics rapidly decreases with the increasing order of harmonics [125]. However,
the experimental observations and calculations with non-perturbation theories show that the typical spectrum of the harmonic intensities has a fast initial decay followed by a region with fairly
constant harmonic intensities, which is called ’plateau’. At the end of the plateau is the so-called
cutoff, beyond which the harmonic intensity drops dramatically.
Over the last three decades, HHG has progressed dramatically, becoming a very active research topic [126, 35, 127]. One application of HHG is to serve as table-top XUV light source.
By now, the plateau of HHG can be extended to keV regime [16, 17]. The ultrabroad bandwidth of the spectra makes it possible to produce attosecond pulses, which can in principle beused to detect electron dynamics directly at its natural time scale. This opened a new researcharea: attosecond science [128, 18, 19, 20]. HHG itself contains information about the electronic
structure of the target and therefore it is extensively used to study the dynamics in atoms andmolecules[129, 130, 131, 132]. As a spectroscopic tool, HHG is successfully applied to image staticand dynamic properties of molecules and solid state materials [133, 134, 135, 136, 137, 32]. Molec-

26
ular imaging with HHG became an important topic, since Itatani et al. [39] have shown that HHG
may be exploited to tomographically image the highest occupied molecular orbital (HOMO) for
diatomic molecules. Concerning bigger molecules, there have also been performed experiments on
HHG for butane, propane and and on allene and acetylene [138]. Going to larger molecules was
motivated by the possibilities of imaging of molecules within HHG. One of the problems, that canbe met in pursuing this goal is related to the fact that in larger molecules one can have interferences between electron wavepackets that can be traced to different orbitals. Several experiments[33, 32, 34, 139] indeed suggest that electrons from more than one orbital contribute during suchprocess.
In the non-perturbative intensity regime it is a great challenge to do theoretical analysis
for interaction with multielectron systems. Without the machinery of perturbation theory one isforced to compromise. One of the most popular approaches assumes the Single-Active-Electron(SAE) approximation. Methods used within this approximation range from tunneling theories, theso called strong field approximation (also called Lewenstein approximation) to numerical solutionsof Schrödinger equation with effective single-electron potential modeling the rest of the molecule.
The approaches like tunnel theory and strong field approximation were very successful for atoms.
Clearly, these approaches do not include multielectron effects or electron correlation.
In this Chapter, we are going to analyze multielectron effects in the HHG spectra frommolecules using TDDFT framework. In particular, we consider the spectra at driving laser frequency matching the energy difference of orbitals in molecules.

3.1

Classical and semi-classical description of HHG

3.1.1

Three-step model

HHG can be understood using the semiclassical three-step model [140, 141] as shown in Fig.
3.1 (left). In the first step, an intense laser field tilts the Coulomb potential and an electron cantunnel through the barrier into the continuum, then the free electron is accelerated by the laser field.

27
When the laser field reverses direction, the electron can return to the parent ion, and recombine
emitting excess energy in the form of high-order harmonic photon. Such process occurs every
half-cycle and it produces a pulse train of radiation. A simple semiclassical theory based picture
captures not only HHG properties but also other interesting phenomena, such as non-sequential
ionization [142] and above threshold-ionization (ATI) [143].

Figure 3.1: Three step model for high-order harmonic generation (left) and typical harmonic spectrafrom noble gas atom (right), adapted from [144].

The generated harmonic spectra have particular characteristic features as shown in Fig.
3.1(right). The harmonic intensity decays quickly as a function of harmonic order followed by aregion with fluctuating but fairly constant intensity, which is called a plateau. The plateau ends
with the so called cutoff, where the harmonics intensities drop several orders of magnitude. The
cutoff can be estimated based on the classical formula [140, 141]
Ecutof f = Ip + 3.17Up
2

(3.1)

2

e E
where Ip is the ionization potential and Up = 4mω
2 is the ponderomotive potential, which is definedoas the average kinetic energy of a free electron in a linearly polarized electric field with amplitude
E and frequency ωo .
Since HHG is a quantum process, its description is contained in the solution of the timedependent Schrödinger equation (TDSE). In general, the solution of TDSE is however very time
consuming. Therefore, models, in particular those using the strong field approximation (SFA), are28
popular and discussed in the following section.

3.1.2

Strong field approximation (SFA)

The strong field approximation is based on the following assumptions: (1) only one electronin the outermost shell is involved in HHG (SAE approximation); (2) only the ground state andcontinuum states are involved in the process; (3) ground-state depletion due to ionization is ignored,and (4) Coulomb interaction between electron and parent ion has no influence on the wavepacket
in the continuum. Following these assumptions, we can derive the SFA formula for HHG (similar
to that given in Lewenstein et al. [141]).
For an electromagnetic field applied to one electron system, the corresponding time-dependent
Schrödinger equation is given byi∂t Ψ(r, t) = [Ho + V + Hint ] Ψ(r, t),(3.2)

where Ho is the free particle Hamiltonian, V is the Coulomb potential, and Hint is the interaction
between the particle and laser field which has the form, Hint = eE · r in length gauge and Hint =
2

e
e
−m
p · A + 2m
A2 in velocity gauge. The exact solution has the form,ˆ ∞Ψ(t, r) = Ψo (t, r) +G(t, t0 )Hint (t0 , r)Ψo (t0 , r)dt0 ,

(3.3)

−∞

where Ψo (t, r) is the initial (ground) state. G(t, t0 ) is the Green’s function describing the propagationfrom t0 to t. Inserting Eq. (3.3) back into Eq. (3.2) shows that the Green’s function satisfies theequation,i∂t G(t − t0 ) = (Ho + V + Hint )G(t − t0 ) + δ(t − t0 ).

(3.4)

The exact solution depends on the potential V . Let us now approximate the solution by ignoring
the effect of Coulomb potential, V , in a strong field. The solution corresponding to the propagationof a free particle in the presence of laser field has the form,
G(t, t0 ) = −iΘ(t − t0 )Go (t − t0 ),

(3.5)

29
where the step function Θ(t − t0 ) arises due to the function δ(t − t0 ) in Eq. (3.4). Then Go satisfies
the equation,i∂t Go = (Ho + Hint )Go .

(3.6)

This equation can be solved in velocity gauge,1i∂t Ψ(p, t) = (p − eA)2 Ψ(p, t),2
withiΨ(p, t) = C exp(−2

ˆ t(p − eA(t00 ))2 dt00 ),

(3.7)

(3.8)

−∞

and can be transferred back to the space domain giving,iΨVp (r, t) = C exp{ip · r −2

ˆ t(p − eA(t00 ))2 dt00 },

(3.9)

−∞

which is called Volkov state [145]. The coefficient C can be obtained through normalization condition,hΨ(r, t)|Ψ(r, t)i = δ(p − p0 ).

(3.10)The solution in velocity gauge is then given by,ΨVp (r, t) =



ˆ
1
i t
00 2 00(p − eA(t )) dt .
exp ip · r −2 −∞(2π)3/2(3.11)The solution in length gauge can be obtained through the gauge transformation:
A = A + ∇χ

(3.12)∂
χ∂t(3.13)

φ = φ0 −Ψ(r, t) = Ψ0 (r, t) exp [iχ(r, t)] .

(3.14)In order to have φ = eE · r and A = 0, we set χ = −eA · r. Thus, the solution in length gauge isequal to,ΨVp (r, t) =
=



ˆ

1
i t
00 2 00
exp
i
(p
−
eA)
·
r
−
p
−
eA(t
)
dt
2 −∞
(2π)3/2


ˆ t

1
i
00 2 00
exp −p − eA(t ) dt |p − eAi,
2 −∞(2π)3/2(3.15)(3.16)

30
where|p − eAi = exp [i(p − eA) · r] .

(3.17)

Thus, the approximated Green’s function then isˆGV (t, t0 ) = −iΘ(t − t0 )d3 p|ΨVp (tihΨVp (t0 )|.

(3.18)

Inserting the Green’s function back into Eq. (3.3), we getˆ ∞Ψ(t, r) = Ψo (t, r) +GV (t, t0 )Hint (t0 , r)Ψo (t0 , r)dt0 .

(3.19)−∞The time-dependent laser induced dipole of the system can then be evaluated asˆ ∞d(t) = hΨ|(er)|Ψi = hΨo |(er)|Ψo i +dt0 hΨo (t)|(er(−iΘ(t − t0 )G(t, t0 )(V ) ))Hint |Ψo (t0 )i + c.c.,−∞(3.20)
where the terms containing the second and higher orders of Hint are ignored. The first term isconstant or zero while the second term is the interesting part and gives information about theinduced oscillating dipole. We denote it asˆ ∞dt0 hΨo (t)|(er(−iΘ(t − t0 )GV (t, t0 )))Hint |Ψo (t0 )i

I=
−∞

ˆ

2

=e i

ˆ t
dp

ˆ t

−∞

dt0 hΨo (t)|r|ΨVp (t)iE(t0 ) · hΨVp (t0 )|r|Ψo (t0 )iˆdpD∗ (p − eA(t), t)E(t0 ) · D(p − eA(t0 ), t0 )
 ˆ t



00 1
00 2
× exp −idtp − eA(t ) + Ip
2
t0
ˆ t
ˆ

=e2 i
dt0 dpD∗ (p − eA(t), t) × exp iS(p, t, t0 ) × E(t0 ) · D(p − eA(t0 ), t0 ),
2

=e i

0

dt

−∞

(3.21)−∞

where the action S is defined as
ˆ t



1
00 2
p − eA(t ) + Ip ,2

(3.22)

D(p − eA(t0 ), t0 ) = hp − eA(t0 )|r|Ψo (r, t0 )i.

(3.23)

0

S(p, t, t ) =

00

dt
t0



and D is defined as

31
In Eq. (3.23), we recognize that hp − eA(t0 )|r|Ψo (r, t0 )i represents the dipole moment for thetransition from state Ψo to the continuum at time t0 , which is usually referred to as ionization
dipole moment, while D∗ = hΨo (r, t)|r|p − eA(t)i is referred to as recombination dipole moment.
The action defined in Eq. (3.22) is the phase accumulated by the wavepacket in the continuum.
Based on the interpretation here, Eq. (3.21) represents the three-step model: an electron is ionizedat time t0 into the continuum, then propagates in the laser field and accumulates the phase S during
the time period from t’ to t, before it finally recombines back to the ground state Ψo at time t.

3.1.3

Saddle point approximation

The direct integration over p in Eq. (3.21) is difficult, not only because it is an integration
over three dimensions but also because strong oscillations from the time dependent phase term partat large p require more precise integration method for accurate results. Saddle point approximation
is often used for such oscillating integrals. It can be motivated in simple terms by noting that thehighly oscillating contributions will cancel each other and only the stationary point has an influenceon the result of the integral. The stationary condition for this integral is equal toˆ t∇p S = ∇p1dt00 ( (p − eA)2 + Ip ) = 0,2t0which yields that´t
pst = e t−τA(t00 )dt00
τ,(3.24)(3.25)

where τ = t − t0 . Using this approximation, andˆ ∞
2π 3/2 ∗I 'e idτD (pst − eA(t), t) × E(t − τ ) · D(pst − eA(t − τ ), t − τ )×
iτ
0
 

ˆ
1 2
1 t
2 00
00
exp −i − pst τ +A (t )dt + Ip τ
.
2
2 t−τ
2



(3.26)The high harmonic spectra then can be obtained as Fourier transform of the second derivativeof the time-dependent dipole in Eq. (3.26) and we will compare the results obtained in this way
within SFA to our TDDFT simulations later.

32

3.2

High-order harmonic generation from molecules
High-order harmonics are currently routinely generated from atoms. HHG can, however,also be generated from molecules and the three-step model has also been used to interpret thebasic physics of the process in molecules. The study of HHG from molecules can provide a deeper
understanding of the process by investigating the contribution from each step of the three-stepmodel. For example, recent experimental [146] and theoretical [30] studies reveal that ionization
rates depend on the molecular symmetry. If the other two steps do not counteract the effects
in step 1, then the HHG yield will be increased when ionization is enhanced, which is confirmedby the experimental observation of high-order harmonics from aligned N2 [39]. For step 2 and3, one can expect quantum interferences, where the recombination amplitudes of the wavepacket
at different nuclei can constructively or destructively interfere with each other as pointed out inrecent experimental [147] and theoretical studies [30, 46]. Destructive interference in the HHG
from diatomic molecules is well studied within the two-center interference model [148, 149, 150].
Our TDDFT simulations also reproduce such destructive interference phenomena, as it can be seenin the Fig. 3.2. The results represent the HHG spectrum from H+
2 in an intense laser field at
wavelength of 800 nm and intensity of 3 × 1014 W/cm2 . The orientation angle, which is the angle
between molecular axis and laser polarization vector, is equal to 30o . The two-center interference
model predicts a destructive interference minimum at 31th harmonic order, as can be observed inour results as well.
Another difference between atoms and molecules is that molecules are not isotropic systems.
HHG spectra from an ensemble of randomly aligned molecules is often similar to those from atoms[151, 152, 153]. However, for aligned molecules, achieved through addition of a relatively weakpulse creating a rotational wave packet, HHG depends strongly on the orientation angle [147].
The alignment techniques make it also possible to use HHG as spectroscopic tool and to image
orbitals. The analysis, however, often relies on approximations such as semiclassical treatment of
the electron-field interaction, the single-active-electron approximation and the three-step model.

33

Figure 3.2: High harmonic spectra minimum due to the quantum interference from the wavefunction
at the two nuclei in H+
2 . The laser field is at the wavelength of 800 nm and intensity of 3 ×
2
14
10 W/cm with pulse length of 20 cycles. Molecular orientation angle is 30 degrees with respect
to polarization of the laser field.

34
One of the example is the proposed method of tomographic imaging of molecular orbitals [39]. It isbased on the single-active-electron approximation, simplified strong-field approximation in which
only recombination matrix element has been extended to molecules and the neglect of the Coulomb
effects. One of the attempts to extend this method by Patchkovskii et al. [32] shows that theconsideration of fully anti-symmetrized multielectron wavefunctions and electron relaxation leads
to a molecular orbital tomography image, that is a mix of Dyson orbital and exchange-correlation
contribution from inner shells. This illustrates how extension beyond the simplest form of SFA canbecome rather cumbersome. Although possible, the extension still is limited by SFA in a form of
separation into three steps and artificial isolation of the electron dynamics in these steps.
In recent years, it has been shown that, in particular in the case of molecules, the generated
harmonic spectra often incorporate more features than predicted by the basic three-step SFA model.
For example, the relevance of multielectron contributions for the interpretation of experimental
data on the ellipticity of high-order harmonics from nitrogen molecules has been demonstrated[154, 155]. This results from the fact that in molecules often there are several orbitals energetically

HHG yield (arb. units)close together in the neutral molecule, as well as the cation. Therefore, the emission of electrons

10-4
10-5
10-6
10-7
10-8
10-9 0

Total

HOMO

HOMO-1

5 10 15 20 25 30
Harmonic orderFigure 3.3: The high harmonic spectra of N2 from its two outermost states HOMO (3σg ) (bluedash) and HOMO-1 (1πu ) (green dots). The total spectra from these two orbitals are also plotted(red solid line). The wavelength of the laser field is 800 nm at intensity of 2 × 1014 W/cm2 withpolarization perpendicular to the molecular axis.

35from different orbitals or the coupling between different orbitals during the interaction with thelaser field, become more likely than in the case of atoms. For example, TDDFT results for HHG
spectra from HOMO and HOMO-1 of N2 are plotted in Fig. 3.3 and compared with the total
harmonic spectra. In this case, contributions from both HOMO and HOMO-1 are comparable in
the total spectrum.
In the rest of this Chapter we identify several multielectron effects in HHG from molecules.
Specifically, we first predict and analyze the presence of fractional harmonics in the HHG spectrafrom molecules. The fractional harmonics show up as sidebands to the usual integer odd harmonics.
They appear when the driving laser wavelength is tuned to the transition between one of the innervalence orbitals and a hole in the highest occupied molecular orbital (HOMO) and is, hence, mostprominent in the interaction of open-shell molecular ions with intense laser radiation. We also
analyze the ellipticity of HHG from N2 , CO2 and compare our results with experimental results,
yielding very good agreement after averaging over the alignment angle which accounts for imperfectalignment present and reported in experiments. The contributions from each valence orbital to theharmonic ellipticity are analyzed. Specifically, the results suggest that for N2 rather strongly boundelectrons in the HOMO-2 orbital have significant contribution to the total HHG spectra and forCO2 nearly all valence orbitals contribute significantly to HHG process and influence the harmonics
ellipticity. In our simulations, all valence orbital electrons are considered active, while the effectof all the other electrons is represented by the pseudopotential. Physical quantities like ionization,high harmonic spectrum, ellipticity, and both total and partial contributions from each orbitals areanalyzed and compared. From such comparison, we can analyze the behavior of each orbital andidentify the significance of their contributions to the processes.

36

3.3

Fractional harmonics1

3.3.1

Open shell molecules

In this section, we study HHG from open shell molecules in their ground states and equilib+
rium geometry such as N+
2 . Initial wavefunction for N2 is obtained by removing one electron fromthe HOMO of N2 and then by performing optimization to find the equilibrium internuclear distanceand corresponding wavefunction. Generating HHG from the ion instead of the neutral molecule canbe a way to extend the plateau cutoff since the ionization potential is larger. Experiments [157, 158]
on HHG from ions using capillary discharged plasma, have shown that the cutoff is indeed greatly
extended. Another scenario, where HHG from ions plays a role, is to study filamentation whereplasma is naturally generated [159, 160].
Fig. 3.4 shows the orbitals and energy level diagram of N+
2 . The HOMO and HOMO-2 hasan energy difference corresponding to wavelength of about 400 nm, while HOMO and HOMO-2
has an energy difference corresponding to 800 nm. Thus, we can study HHG from N+
2 at typicallaser wavelengths, which can couple resonantly two orbitals in the molecular ion.

Figure 3.4: Energy level diagram (left) and orbitals (right) of the valence shell of N+
2.

1

Part of the results presented in this section are in [156].

37
3.3.1.1

Nonadiabatic behavior of laser induced dipole

The high harmonic generation from a system can be calculated through the Fourier transformof the laser induced dipole. Within the TDDFT formalism, the calculation of laser induced dipole
is equal todtot =

Xdi ,(3.27)iwhere di is the dipole from ith orbital,ˆdi = hφi (t)|r|φi (t)i =rni (r, t)dr,(3.28)

where φi (t) and ni (t) are the wavefunction and corresponding partial density for the ith orbital,
respectively.
HOMO and HOMO-2 in N+
2 have an energy gap corresponding to about 400 nm and basedon the symmetry of the orbital they can only be coupled via one-photon transition, when the laser
electric field is parallel to the molecular axis. In Fig. 3.5, we present the results of simulations

(a)

(b)

(c)

Figure 3.5: Time dependent dipole for HOMO-2 of N+
2 for the case when the molecular axis isperpendicular (a) and parallel (b) to the laser electric field at wavelength of 400 nm and intensityof 2 × 1014 W/cm2 . The amplitude of the laser field (on arbitrary scale) is shown by dashed line forcomparison. The Fourier transform spectrum of the residual dipole after the laser pulse for parallelcase (b) is plotted in (c).

with a laser field set to the wavelength of 400 nm and intensity of 2 × 1014 W/cm2 with trapezoidal
envelope and pulse length of 16 fs. The laser induced dipole behaves adiabatically (i.e. following
the external laser field) in (a) for perpendicular orientation but does not follow the laser field for theparallel orientation (b). The latter nonadiabatic response is due to the coupling between HOMO

38
and HOMO-2. This can be seen from the analysis of the dipole after the laser pulse has finished.
In our simulations, we have propagated the wavefunction 3 fs after the end of the laser pulse. In
Fig. 3.5, the residual dipole after laser pulse is shown: for adiabatic case it is very small (a),while it is prominent for the parallel setup (b). To confirm that such residual dipole is due to thenon-zero population transfer between HOMO-2 and HOMO after the laser pulse, Fourier transform
is applied to the residual dipole and a dominant peak at 0.11 a.u. (corresponding to the energy
difference between HOMO and HOMO-2) in the spectrum is found as shown in Fig. 3.5(c). We
will now analyze the signatures of the nonadiabatic response in the HHG spectra.

3.3.1.2

Fractional harmonics from N+
2

Because of the relationship between HHG spectrum and the oscillating dipole,
P (ω) =

ω4d(ω) · d∗ (ω).
12πεo c3

(3.29)

one can expect additional features in the spectrum for the coupling between orbitals. In Fig. 3.6, thespectrum for the perpendicularly oriented molecular ion at 400 nm shows the traditional spectrum
consisting of odd harmonics only (panel a), while in the spectrum for the parallel aligned molecule
there appear first and second order sidebands for each of the odd harmonics in the spectrum (panelb).
The sidebands do not show up at the integer multiples but at fractionals of the fundamental
laser frequency and are, hence, termed “fractional” harmonics. For the present laser parameters
they are related to the transition between HOMO-2 (2σu ) to HOMO (3σg ) in this open shellmolecule. The intensity of the fractional harmonics is surprisingly high, which can be attributed
to a strong transition and large transition dipole matrix element describing the coupling between
these two molecular states.
The mechanism responsible for the sidebands is analogous to that for Mollow sidebands inquantum optics [161], and can be intuitively explained with the so-called dressed states. As shownin Fig. 3.7, if two states are coupled resonantly by an external laser field, they become degenerate

39

(a)

(b)

Figure 3.6: High harmonic spectra of N+
2 , aligned parallel (a) and perpendicular (b) to the polarization axis of a laser pulse at 400 nm and 2 × 1014 W/cm2 . The insets show an enlarged view ofpart of the spectrum, which clearly exhibit the first and second order Mollow sidebands displacedby the Rabi frequency Ωr and 2Ωr from the odd harmonics for the aligned molecules.

40
and create two dressed states |+, ni and |−, ni, where n is the photon number, separated by anenergy of Rabi frequency. The possible transitions between these dressed states are illustrated inFig. 3.7. There are three different possible emission energy values ωo , ωo + Ωr , and ωo − Ωr , whereΩr is Rabi frequency. Although in the illustration only first order sidebands are presented, higher
order sidebands have also been observed in simulations and experiments [162].

Figure 3.7: Illustration of Mollow triplets.

In the insets on the right of Fig. 3.6, we present enlarged views of parts of the spectra. These
clearly show the relation between the separation of the first and second order sideband peaks in
the HHG spectrum and the Rabi transition frequency Ωr .

Ωr =

µE
.
~

(3.30)

We studied the splitting of fractional harmonics more systematically by varying the laser field
intensity, which influences the Rabi transition frequency directly and the results of simulations areshown in Fig. 3.8. One can see that the first order splitting and the calculated Rabi transition
frequency are in almost perfect agreement with each other, while the second order splitting is givenby the double of the Rabi transition frequency. This confirms that the fractional harmonics arestrongly related to the Rabi transition frequency.
The results in Fig. 3.6 involve the transition between the same type of orbitals as in thepreviously studied case of H+
2 [163], where a 1064 nm laser field is applied to couple the so-calledFrequency [ωlaser]

41

1.2 1st order sideband
frequency
1.0 Rabi
2rd order sideband
0.8 2xRabi frequency
0.6
0.4
0.2
0.0 0.04 0.05 0.06 0.07 0.08
E [a.u.]Figure 3.8: Comparison of the distance of the sideband positions, i.e. first and second order
fractional harmonics, with respect to the corresponding odd harmonic (averaged over harmonics
bellow 11th) and the Rabi transition frequency for different laser field intensities. Laser wavelengthof 400 nm and orientation of molecule parallel to the laser field.

05

10 15 20
Harmonic order25

HHG yield (arb. units)

HHG yield (arb. units)

42

10-5
10-6
10-7
10-8
10-9
10-10
10-11

10-5
10-6
10-7
10-8
10-9
10-10
10-11

0

10 20 30 40
Harmonic order

50

Figure 3.9: HHG spectra from N+
2 for laser field at wavelength of 600 nm (left) with intensity of
1 × 1014 W/cm2 and at wavelength of 800 nm (right) with intensity of 2 × 1014 W/cm2 . Molecularaxis is oriented parallel to the polarization of laser field. Pulse length is 20 fs.

Charge-Resonance (CR) states HOMO (1σg ) and Lowest Unoccupied Molecular Orbital (LUMO)(1σu ) at stretched internuclear distances, in the regime of Charge Resonance Enhanced Ionization(CREI). In contrast, in our studies the fractional harmonics effect appears for the interaction of
the molecular ion with the field at equilibrium distance and the coupling between valence orbitals,
showing that the phenomenon is more general than considered before.
We have also performed simulations for the HHG spectra for laser wavelengths of 600 nm
and 800 nm for the molecule aligned parallel to laser polarization vector. Because the wavelengthis shifted away from the resonant wavelength 400 nm, we are entering near resonance (600 nm) and
off-resonance (800 nm) cases. As can be seen in Fig. 3.9 (left), for 600 nm the fractional harmonics
are still present in the spectra but are not as intense as those for 400 nm. But at 800 nm, thespectrum (Fig. 3.9 (right)) consists almost only of odd harmonics. Such simulations confirm thatthe fractional harmonics result from coupling between the orbitals. Since the fractional harmonics
can still be seen in the near-resonant case (600 nm), the fractional harmonics should be observablein experiments, even if the laser frequency does not match exactly the energy difference between
the orbitals.
Mollow sidebands in high harmonic generation have so far been predicted for the coupling
between orbitals of σ-symmetry only. However, in general, the phenomenon appears for the coupling

43
between molecular orbitals of different type of symmetry as well. Since the symmetry of the orbitals
governs the alignment dependence of the transition matrix element, the strengths of the respective
fractional harmonics depends on the alignment between the molecular axis with respect to the

HHG yield (arb. units)polarization vector.

10-5
10-6
10-7
10-8
10-9
10-10
10-11

10-8
10-9
10-10
10-11

0

10

21

20 30 40
Harmonic order

23

25

50

27

29

60

0.30
0.25
0.20
0.15
0.10
0.05
0.00
0

­

5

3σg |2σu®

­

3σg |1πu1015Time [fs]®

Projection

Projection

Figure 3.10: High harmonic spectrum of N+
2 aligned perpendicular to the polarization axis of 60
fs laser pulse at wavelength of 800 nm and intensity of 2 × 1014 W/cm2 . The inset provides anenlarged view of the part of the spectrum.

20

1.0
0.8
0.6
0.4
0.2
0.0
0

­

5

3σg |2σu®

­

3σg |1πu1015Time [fs]®

20

Figure 3.11: Time dependence of the projection between the coupled orbitals of N+
2 during the interaction with the laser pulse for the case of parallel orientation (left) and perpendicular orientation(right). Laser field parameters: wavelength of 800 nm and intensity of 2 × 1014 W/cm2 .

For example, fractional harmonics due to the coupling between orbitals of πu - (HOMO-1)and σg -symmetry (HOMO) in N+
2 are induced for the interaction with a 800 nm laser pulse andperpendicular alignment of the molecule. The respective results of our calculations for the HHG
spectrum in Fig. 3.10 clearly reveal the occurrence of the Mollow sidebands at these parameters.

44
For comparison the results for parallel orientation with the same parameters are shown in Fig.
3.9(right). One can still see small splitting which, however, comes from the weak coupling between
HOMO and HOMO-2. The coupled orbitals can be identified by projecting one orbital to the otherduring the simulation and results are shown in Fig. 3.11. For the parallel case, the projection from
HOMO-1 (1πu ) to HOMO (3σg ) is zero while there is almost 25% population transfer between
HOMO-2 (2σu ) and HOMO. For the perpendicular case, the HOMO and HOMO-1 coupling is verystrong leading to 80% population transfer between these two orbitals while the projection from
HOMO-2 to HOMO in this case is nearly zero. One can conclude that the fractional harmonics
for the perpendicular case arise due to the coupling between HOMO and HOMO-1 and the minor
fractional harmonics for the parallel case come from weak coupling between HOMO-2 and HOMO.

3.3.1.3

Modified SFA model

Since we now understand that the fractional harmonics arise from the Rabi flopping inducedby the laser field, we seek for a model to incorporate such mechanism into the strong-field approximation formalism. Suppose, the electron oscillates between the two states |ψ1 i and ψ2 i with Rabi
frequency Ωr . Then the corresponding state in the SFA model is|Φi = C1 (t)|ψ1 i + C2 (t)|ψ2 i,(3.31)

where C1 (t) = cos( Ω2r t) and C2 (t) = −i sin( Ω2r t), which are obtained by modeling the Rabi floppingin a two-level system. By substituting Eq. (3.31) into Eq. (3.26) we obtain
I = I11 + I22 + I12 + I21 ,
where

ˆ ∞
Iij =

dt0 hψi (t)Ci (t)(−iqrΘ(t − t0 ))G(t, t0 )Hint Cj (t0 )|ψj (t0 )i.

(3.32)(3.33)−∞

The four terms in Eq. (3.32) correspond to four HHG channels. The first two terms I11 and I22
represent the case, where the initial state and the final state are the same, while the last two terms,
I12 and I21 , represent the scenario in which the electron is ionized from one state and recombines into

45
the other. This model captures the resonant nature of the process and it illustrates the necessary
modification of the ordinary HHG mechanism. The comparison of results from the modified SFA
and the TDDFT simulations is shown in Fig. 3.12. As one can see, the modified SFA does produce
the fractional harmonics at exactly the same positions as the results from TDDFT, which further
confirms that the fractional harmonics are due to the coupling effect.

Figure 3.12: Comparison of harmonic spectra for N+
2 calculated within modified SFA and TDDFT.
Laser parameters: λ=400 nm and I = 1 × 1014 W/cm2 . Polarization parallel to the molecular axis.

3.3.1.4

Enhanced harmonic generation

Our TDDFT calculations also allow us to examine the harmonic spectra contributions fromeach orbital in N+
2 , separately. The results are shown in Fig. 3.13 for the parallel alignment case.
Interestingly, the strength of high harmonics from HOMO-2 (spin up component) is much larger
than those from spin down components (N+2 is created in our calculation by removing one spinup electron from HOMO). This can be explained by looking at the ionization from each orbital.
In Fig. 3.13(d), it can be seen that the ionization yield from spin-up HOMO-2 is much larger
than ionization from other orbitals. The resonance between HOMO-2 and HOMO enhances theionization, which is one part of the explanation. But, the ionization yield from spin-up HOMO-2

46

(a)

(b)

(c)

(d)

Figure 3.13: Harmonic spectra from HOMO (a) and HOMO-2 (b) of spin down component and
HOMO-2 (c) of spin up component and population of coupled orbitals (d) (ionization from other
orbitals is negligible). Laser field is at 400 nm and 2 × 1014 W/cm2 with polarization parallel tothe molecular axis.

47should not be larger than the ionization yield from spin-down HOMO, if just a resonance process
plays a role. The reason why the ionization yield for spin-up HOMO-2 is even larger than for
spin-down HOMO can be found from the analysis of the orbital energies as shown in Table 3.1.
For the HOMO the orbital eigenvalue is approximately equal to the vertical ionization potential.
Table 3.1: Orbital energies (eigenvalue) of N+
2
orbital
HOMO-3(2σg )HOMO-2(2σu )HOMO-1(1πu )HOMO(3σg )spin down-1.488381-0.939705-0.849458-0.848556

spin up-1.465770-0.883796-0.831411-0.779313(unoccupied)

As can be seen from the table, the spin-up HOMO’s ionization potential is lower than that of thespin-down HOMO. Therefore, the spin-up HOMO-2 coupled to the spin-up HOMO has a dominant
ionization yield in this case, resulting in a dominant HHG yield.
Please note that in our numerical simulations the spin-up and spin-down orbitals are not
restricted to be the same. Such method is called unrestricted open shell TDDFT [164]. The
advantage of unrestricted calculations is that they can be performed very efficiently. However, thewavefunction calculated in this way may no longer be the eigenfunction of the total spin operatorhS 2 i, and if that happens an error called ’spin contamination’ is introduced into the calculations.
The spin contamination can lead to lowering of the computed total energy because of the extra’freedom’, but more often it leads to a higher total energy since higher energy states are being mixed
into the solution. Therefore it is important to check, whether such relative energy shift between
spin-up and spin-down components is due to the spin contamination. The spin contamination
can be checked, for example, via the calculation of hS 2 i. According to [164], the total spin of
spin-polarized DFT can be calculated as
Nβ ˆ
Nα X
X
hS i = Ms (Ms + 1)~ + ~ Nβ − ~| φαi (r)φβj (r)d3 r|2 .
2

2

2

2

(3.34)

i=1 j=1
2
For N+
2 , the value hS i should be 3/4. We checked the calculations for different functionals and therelative error is less than 1% and thus we conclude that spin contamination does not have influence

48
on our calculations. So we can trust the energies presented in the Table 3.1.
As we mentioned before, the cutoff of HHG from ions is extended because of high ionizationthpotential. The cutoff for HHG from N+
2 in our simulations for the 400 nm is around 11 harmonic,

corresponding to Ip = 1.0 a.u., which means the cutoff is actually corresponding to the ionization
potential of spin-up HOMO-2. Without resonance, the dominant HHG is from spin-down HOMO
and thus for the cut-off, relevant ionization potential is that of spin-down HOMO. Therefore, theresonance further extends the cutoff of HHG. Because of the enhanced ionization due to resonance,
it seems such a resonance is a way to enhance HHG yield and extend the cutoff simultaneously.

3.3.1.5

Nonadiabatic electron density localization

The strong coupling between the orbitals does not only leave its footprints in the high harmonic spectra, but does strongly influence the dynamics of the electron density. More specifically,
the emergence of fractional harmonics comes along with profound variation of the electron dynamics inside the molecule as well as during its excursion in the continuum throughout the harmonic
generation process. First, the molecular orbital coupling leads to a laser induced nonadiabatic
electron dynamics in the molecule, where the electron instead of directly following the changes in
the oscillating electric field experiences lags and dwells at one of the nuclei for longer than a half
period of the laser field cycle. The time scale of the corresponding localization of the electron isgoverned, as the location of the fractional harmonics in the spectrum, by the transition frequency,
i.e. Rabi frequency, and can therefore in principle be controlled via the laser parameters. This nonadiabatic dynamics can be visualized in our calculations via the difference of the time-dependent
total electronic density of the propagated molecular wavepacket and the initial electron density
ρ(z, t) − ρ(z, 0), where,ρ(z, t) =Xρi,σ (z, t)(3.35)

dxdy|φiσ (r, t)|2 .

(3.36)

i,σ

with

ˆ
ρi,σ (z, t) =

49
The electron density ρ(z, t) is obtained by integrating over the two spatial dimensions transversal
to the molecular axis as shown in Eq. (3.36). We present the resulting density for N+
2 as a functionof the electron position along the internuclear axis and as a function of time in Fig. 3.14. The laser

Figure 3.14: Difference between the full time-dependent electron density of N+
2 , aligned along (left)
and perpendicular (right) to the polarization of laser, and the initial electron density as a functionof time and the position along the molecular axis. The density is integrated over the spatial
dimensions transversal to the molecular axis. Laser parameters as in Figure 3.6.

parameters and the orientation of the molecular axis (along the electric field vector) are the same
as for the results in Fig. 3.6. For the purpose of comparison the electric field is included at the
top. For the sake of comparison, we present changes in the density difference for the case where no
coupling is present (right) and one can observe typical adiabatic dynamics.
The result clearly shows that in the Rabi flopping regime the electron density does not
swap adiabatically from one side of the molecular ion to the other every half cycle of the electricfield (left). Instead dynamical localization islands can be observed along the internuclear axis inparticular. However, also the electron density localized around nuclei shows nonadiabatic behavior.
For adiabatic case (right) one can typically observe changes in electron density primarily around thenuclei. The frequency of the associated nonadiabatic oscillation is estimated to be approximately
equal to the Rabi frequency of the transition between the two states of σ-symmetry. We will further
explore the nonadiabatic dynamics in Chapter 5 of this thesis.

50
3.3.1.6

Modification of the time-frequency analysis of HHG

Since the Rabi flopping induces additional dynamics during the HHG process, we would liketo investigate the HHG process in time domain in more details. Time-frequency analysis (also called
wavelet analysis) is such a tool which makes it possible to examine how high harmonic generation
changes with time. The time-frequency spectra is obtained through formula Eq. (3.37)
ˆd(ω, τ ) = d(t) exp(−(τ − t)2 /2σ 2 ) exp(−iωτ )dt.

(3.37)

In time-frequency analysis, instead of performing fast Fourier transform (FFT) directly to the laser
induced dipole, a Gaussian type function exp(−(τ − t)2 /2σ 2 ) is added, which constraints the FFT
only to be applied to the dipole in the time window defined by the Gaussian function. The resultof the time-frequency analysis for nonresonant case is presented in Fig. 3.15. From the Figure

Figure 3.15: Wavelet analysis of high harmonic generation in N+
2 aligned perpendicular to thepolarization direction of laser field. Laser parameters are the same as in Fig. 3.6.

we can see that there is a HHG burst during every half laser cycle. For each harmonic order,
there are two recombination times, which are related to two types of trajectories recognized fromclassical calculations. In typical time frequency analysis the pattern has ”lambda” type shape,where the left ’arm’ is corresponding to short trajectories, while the right one is corresponding tolong trajectories, as indicated in Fig. 3.15(right)
The nonadiabatic electron dynamics in the molecule leads to modifications of this traditional
semiclassical picture of HHG, as can be seen from the time frequency analysis in Fig. 3.16. In

51

Figure 3.16: Wavelet analysis of high harmonic generation in N+
2 aligned parallel to the polarization
direction of the laser field. Laser parameters are the same as in Fig. 3.6.

52contrast to the typical wavelet analysis for high harmonic generation, which exhibits the return ofthe electron and harmonic emission via short and long trajectories, the generation of the frequencies
in the harmonic spectrum as a function of time in this case is more complex. Obviously, anyinterpretation in terms of classical paths becomes cumbersome. Instead, the time-frequency pattern
is modulated by the nonadiabatic dynamics and temporal variation of the electron localization inside
the molecule and, not surprisingly, shows an oscillating variation of the pattern with the time scale
related to the transition frequency between HOMO and HOMO-2 in the nitrogen molecular ion,
which is shown overimposed in the lower part of spectrum in Fig. 3.16(left).

3.3.1.7

Fractional harmonics from other molecules

Present section is devoted to illustration that our findings are not restricted to the nitrogenmolecular ion. Indeed, from the results of our calculations we observe the appearance of fractional+
harmonics in other larger open shell molecules as well, like CO+
2 , C2 H4 , and NO2 .

CO+
2
2
2
2
4
2
3
The valence shell electron configuration of CO+
2 is (3σg ) (2σu ) (4σg ) (1πu ) (3σu ) (1πg ) .

To induce coupling, we find that the 1πu and 1πg can be coupled at parallel orientation with awavelength around 350 nm, which is available in experiments. The coupled orbitals are shown in
Fig. 3.17. The harmonic spectra for both coupling and non-coupling case are shown in Fig. 3.18.

(a) 1πu(b) 3σu(c) 1πgFigure 3.17: Orbital image of CO2 .

The coupling between orbitals is very strong, leading to appearance of intense fractional harmonics

53(left panel) for parallel case, while for the perpendicular case the harmonic spectrum is regular and

10-4

HHG yield (arb. unit)

HHG yield (arb. unit)

presents odd harmonics only (right panel).

10-6
10-8
10-10
10-12 1

5
9
13
Harmonic order

17

10-4
10-6
10-8
10-10
10-12 1

5
9
13
Harmonic order

17

Figure 3.18: Fractional harmonics from CO+
2 with laser field at wavelength of 350 nm and intensityof 1 × 1014 W/cm2 at parallel (left) and perpendicular orientation (right).

C2 H+
4While the previously studied cases are for the coupling between inner valence shell orbitals
and HOMO, in C2 H+
4 two inner valence orbitals can be coupled. One electron is removed from 1b3g ,which leaves the electron configuration as (2ag )2 (2bu )2 (1b2u )2 (3ag )2 (1b3g )1 (1b3u )2 . With a 400 nm
laser field, a coupling can be induced between 1b2u and 1b3g . The orbital images of C2 H+
4 of thecoupled orbitals are shown in Fig. 3.19. We analyze the generated harmonic spectrum with driven(a) 1b2u(b) 1b3g

Figure 3.19: Orbital images of C2 H4 .

laser field at 400 nm with polarization parallel to the axis along carbon-carbon bond. Fractional
harmonics are present again due to the coupling, but since the coupling is not as strong as in CO+
2
and N+
2 , only first order splitting fractional harmonics are observed in this case, shown in Fig.

54
3.20(left). For the case where laser polarization is perpendicular to the carbon-carbon bond, only

10-4
10-6
10-8
10-10
10-12 1

5
9
13
Harmonic order

HHG yield (arb. units)

HHG yield (arb. units)

odd harmonics are present.

17

10-4
10-6
10-8
10-10
10-12 1

5
9
13
Harmonic order

17

Figure 3.20: Fractional harmonics from C2 H+
4 with laser field at wavelength of 400 nm and intensity2
14
of 1 × 10 W/cm with polarization parallel (left) and perpendicular (right) to the C-C bond.

NO2
Open shell molecules can be created by removing an electron from neutral molecules, whichis, however, difficult for experimentalists to obtain. We, therefore, selected also a neutral open shellmolecule NO2 in our studies, which has the electron configuration as (5a1 )2 (1b1 )2 (1a2 )2 (4b2 )2 (6a1 )2 ,where the orbital 6a1 is partially occupied. The valence orbitals’ images are presented in Fig. 3.21.
Considering the nonlinear geometry, it is interesting to see if the fractional harmonics are stillpresent. A laser at 400 nm with intensity of 1014 W/cm2 is applied to couple the orbital 4b2 andthe orbital 6a1 with polarization parallel to the oxygen-oxygen bond. The coupling is, however,not very strong. Instead of separated peaks in the spectrum, we therefore get a broadening of eachodd harmonics, as shown in Fig. 3.22.

3.4

Ellipticity of high-order harmonics 2
Since the high-order harmonic generation has recently become an important XUV source[17], the ellipticity of HHG is particularly interesting for experimentalists because the control of
ellipticity for soft-X ray is very difficult compared to lasers at wavelengths within the visible light2

Parts of the results in this section are published in Y. Xia and A. Jaron-Becker, ”Multielectron contributions in
elliptically polarized high-order harmonic emission from nitrogen molecules”, Opt. Lett. 39, 1461 (2014).

55(a) 4b2(b) 6a1

Figure 3.21: Orbital images of NO2 .

Figure 3.22: Fractional harmonics from NO2 with laser field at wavelength of 400 nm and intensityof 2 × 1014 W/cm2 with polarization direction parallel to the oxygen-oxygen axis.

56
range. It would be appealing, if specific elliptically polarized laser light can be generated from thesource directly. Since atoms are spherically symmetric systems, the generated HHG from noble gas
atoms by linearly polarized laser field is also linearly polarized. HHG from molecules, however, isnaturally elliptically polarized because of the non-spherical symmetric orbital. Recent experiments
have shown that elliptically polarized HHG can be obtained from molecules N2 and CO2 [165, 166].
In this section we present results of numerical calculations for the polarization and ellipticity
of high-order harmonics from molecules like H+
2 , H2 , N2 , and CO2 . We show that for simple
one-electron system, the ellipticity has a large value for harmonics at angles where destructively
interference occurs. The ellipticity pattern for larger molecules is more complex and multielectron
effects start to play a role. For multielectron systems, our calculated ellipticity of high-orderharmonics from N2 and CO2 are in good agreement with the data of recent experiments [165, 154]
and reveal the importance of multielectron contributions for N2 , where at least three different
orbitals (HOMO, HOMO-1, and HOMO-2) have influence on the ellipticity and for CO2 , where
nearly all valence electron orbitals play a role. In particular, we show that the contributions frominner-valence orbitals are most significant at wavelengths at which there is a resonance between
two orbitals. Furthermore, we confirm that it is essential to take proper account of the alignmentof the internuclear axis with respect to the polarization direction in the molecular ensemble.

3.4.1

Ellipticity from one-electron system H+
2

The high harmonic spectrum d(ω) is obtained through fast Fourier transform of the laser
induced dipole d(t), see Eq. (3.27). For linear molecules, the laser induced dipole usually has twocomponents: the parallel and perpendicular components categorized by the direction with respect
to the direction of the laser electric field. The ellipticity of a given harmonic is determined by[167, 150]s=p
1 + 2r2 cos(2δ) + r4p,
1 + r2 + 1 + 2r2 cos(2δ) + r4
1 + r2 −(3.38)

57
where r is the amplitude ratio of perpendicular component to parallel componentr = |d⊥ (ω, α)| / dk (ω, α) ,(3.39)and δ is the relative phase between the two components,δ = arg [d⊥ (ω, α)] − arg dk (ω, α) .

(3.40)

From Eq. (3.38), the maximum ellipticity occurs when the amplitude ratio r=1 and the relative
phase δ = 90o , which means the strength of perpendicular component d⊥ has to be comparable
to the parallel component dk . This however is usually not satisfied because the electron oscillates
strongly along the laser field and thus produces strong parallel and weak perpendicular components.
We first study the ellipticity of HHG from the simplest one-electron system H+
2 . In Figure
3.23 with a laser field at wavelength of 800 nm and intensity of 3 × 1014 W/cm2 , which is also
investigated by Chu et al. [150], we compare the amplitude and phase difference of parallel andperpendicular components of 57th harmonic. With the change of alignment angle of molecularaxis to the laser polarization axis, the parallel component has a minimum around 50o while theperpendicular component has a maximum value around the same angle and comparable to theparallel component, which indicates r ≈ 1 around 50o . On the other hand, the phase differencechanges dramatically from −0.6π to 0.4π. Based on our analysis, we therefore expect maximum
ellipticity around 50o while close to zero elsewhere.
The angle of 50o for this harmonic actually corresponds to the occurrence of an interference
minimum for parallel component at this harmonic according to the two center interference model
where minimum occurs when1
d cos θ = (n + )λ.
2(3.41)

Here d is the internuclear distance and λ is wavelength of ionized electron wavepacket. If the
wavefunction has different sign on the two center, then 1/2 should be dropped. According to theanalysis, we would expect the maximum ellipticity will occur around the interference minimum of
the individual harmonics in HHG of H+
2 , as shown in Figure 3.24. It can be seen that for each1.0

d (57)d (57)

10

30
50
70
Alignment angle [degree]

Phase difference [π]Amplitude [arb. units]

58

0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000

57

0.5
0.0
0.5
1.0 10

90

30
50
70
Alignment angle [degree]90

Figure 3.23: The amplitudes of parallel and perpendicular components (left) and phase differenceδ (right) of 57th harmonic order of H+
2 as a function of the alignment angle in a laser field at
wavelength of 800 nm and intensity 3 × 1014 W/cm2 .

harmonic the amplitude ratio is close to 1 and the phase difference also changes dramatically near
the interference minimum.
1.0

1.5
1.0
0.5
0.0

10

30
50
70
Alignment angle [degree]

90

1.0

27
29
31
33

0.5
0.0
0.5
1.0 10

27
29
31
33

0.8
Ellipticity

27
29
31
33

Phase difference [π]Amplitude ratio r

2.5
2.0

0.6
0.4
0.2

30
50
70
Alignment angle [degree]

90

0.0

10

30
50
70
Alignment angle [degree]90

Figure 3.24: The amplitude ratio r (left), phase difference δ (middle), and the ellipticity (right) of
14 W/cm2 .
high-order harmonics from H+
2 in a laser field at wavelength of 800 nm and intensity 3×10

3.4.2

Complex ellipticity pattern due to two-electron effects

We have shown that the maximum ellipticity of harmonics for a simple one-electron system
arises around the interference minimum for H+
2 . We will now consider the influence of two-electron
effects in H2 . Considering that the two electrons are in the ground state 1σg , we do not expect too
much difference for the amplitude and phase difference pattern. We have performed a simulation forH2 with the same laser parameters as for H+
2 . The results in Fig. 3.25 do share some similar pattern
for the amplitude ratio: maximum amplitude ratio and phase change around the interference
position which appears at the angle of 40o . However, the ratio here is not as smooth as that for

59
+
H+
2 . Concerning the ellipticity, while the results for H2 show clearly a strong ellipticity aroundthe interference minimum, the ellipticity pattern for H2 is more complex although the ellipticity

1.0

27
29
31
33

10

30
50
70
Alignment angle [degree]

90

1.0

27
29
31
33

0.5
0.0
0.5
1.0 10

27
29
31
33

0.8
Ellipticity

16
14
12
10
8
6
4
2
0

Phase difference [π]Amplitude ratio r

maximum is still located near the position of interference minimum.

0.6
0.4
0.2

30
50
70
Alignment angle [degree]

90

0.0

10

30
50
70
Alignment angle [degree]90

Figure 3.25: The same as Figure 3.24 but for H2 .

Such comparison of the results for these simplest molecules indicates that the ellipticity canbe strongly influenced by multielectron effects. In case of the two-electron molecule, the ellipticity
pattern is already considerably more complex than that for the one-electron molecule. For moleculeswith more than two electrons, we therefore expect even more complicated features since effects likeinterference from orbitals with different symmetry and coupling between the orbitals may occur.

3.4.3

Alignment angle average

In the following sections, we are going to compare our simulation results with experimental
data. In experiment, molecules are often aligned using a laser pulse which initiates a rotational
wavepacket. The molecules in the ensemble are then aligned at certain time according to therotational period. The alignment in the experiment is however usually not perfect. The distributionof alignment is measured via hcos2 θi. We account for this distribution by carrying out simulations
for the angles in the distribution in order to do angle average and compare with observations inexperiments.
To obtain the ellipticity at different alignment angle, the time-dependent Kohn-Sham equationi∂φk (r, t) = −∇2 φk (r, t) + (VKS [ρ(r, t)] + U (r, t)) φk (r, t)∂t(3.42)

60
is solved for each alignment angle. The parallel and perpendicular components of dipole momentdk (θ, t)) and d⊥ (θ, t) are obtained via Eq. (3.28). To take account of the distribution of alignment
angles in molecular ensemble in the experiment, it is then straightforward to average the dipole
moments using the reported alignment distributions [167, 165].
Let the probe laser electric field be along the z-direction (x=0,y=0,z=1) as shown in the
diagram in Fig. 3.26 and the pump laser field lies in the y-z plane with a polar angle α and thus
the electric field direction for the pump laser field is (0, sin α, cos α). After applying the pump laser
field, the linear diatomic molecule is oriented in the direction (sin θm cos φm , sin θm sin φm , cos θm ).
The intersection angle between the pump laser field and the molecular axis is denoted as θ. Thenthe orientation distribution of molecules with respect to the pump laser field is only a function ofθ, denoted as ρ(θ). The angle θ can be obtained viacos θ = sin α sin θm sin φm + cos α cos θm .

(3.43)

Thus the alignment distribution for the molecular axis at orientation angle (θm ,φm ) with respect
to the pump probe angle α isρ(α, θm , φm ) = ρ(θ(α, θm , φm )),(3.44)and the harmonic spectrum after angle average is given byd¯⊥/k (t, α) =

ˆ π ˆ 2πd⊥/k (t, θm , φm )ρ(α, θm , φm ) sin θm dθm dφm .
0(3.45)

0

For the calculation of d⊥/k (t, θm , φm ), in case of a diatomic molecule, we only perform simulations
for several polar angles θm at φm = 0. For other non-zero φm , the component d⊥/k (t, θm , φm ) canbe obtained via a transformation of the calculation results at φm = 0.

3.4.4

Comparison with experimental data for N2

Measurements of the polarization and ellipticity of high-order harmonics generated in nitrogenmolecules have been reported recently. Zhou et al. [165] carried out an experiment on molecule
N2 to study the ellipticity of the generated harmonics, in which a 30 fs probe pulse at intensity of61

Figure 3.26: Illustration diagram for angle average procedure.

62
2×1014 W/cm2 was focused on the molecule N2 which was aligned by a pump pulse. The ellipticity
dependence as a function of the pump-probe angle was studied by varying the time delay between
pump and probe pulse. In the experiment, the alignment parameter hcos2 θi was measured to be
0.65, as compared to 0.3 for random alignment. The measured ellipticity, shown in Fig. 3.27(c), isvery large as compared to that observed in a previous experiment [168].
To explain the experimental observations, we have performed simulations with similar laserparameters and the results are angle averaged using the measured distribution. We compared resultsfor different intensities and achieved a good agreement with experiment, which was performed at a
relatively low intensity 1 × 1014 W/cm2 . In fact, it is difficult to measure the intensity at the focal
spot accurately. Therefore, it is not surprising that results from theory and experiments agree at
different laser intensities. We attribute the sensitivity of the ellipticity as a function of intensity
to the multielectron character within the present theoretical method. It has been shown by Son
et al. [150], and confirmed by us in test calculations, that for one and two-electron molecules theellipticity varies more smoothly.
Fig. 3.27 compares our simulation results with angle average (middle) and without (left) as
well as experimental results (right). The angle averaged results are in good agreement with theexperimental observations showing a maximum ellipticity of about 0.4 for alignment angles around40o − 60o . As far as we know, this is the first time an agreement between results from ab-initio
calculations with experimental data for ellipticity has been found. In the following subsection wewill further analyze the contributions from different orbitals.
We note that the ellipticity value is pretty large before angle average. For example, for analignment at 40o , the ellipticity of the 19th harmonic is close to 0.8, while after angle average,
the maximum ellipticity is only 0.45 for the 21th harmonics. On the other hand, the ellipticity
structure is smoother and has similar trend for different harmonics after the average, which meanslots of detailed information get lost. Since obtaining and manipulating the ellipticity of HHG is an
important goal, our simulations suggest that improving the alignment technique is one way to gethighly elliptically polarized harmonics. Furthermore, an improved alignment technique would make

1
0.8
0.6
0.4
0.2
0
11

40o
50o
60o

15 19 23 27
Harmonic Order

Ellipticity

Ellipticity

63

0.5
0.4
0.3
0.2
0.1
0
11

40o
50o
60o

15 19 23 27
Harmonic Order

Figure 3.27: Ellipticity of high-order harmonics as a function of the harmonic orders for threemolecular alignment angles at λ = 800 nm and peak intensity of Io = 1 × 1014 W/cm2 . Comparisonof results in which no average (left panel) and a proper average (middle panel) over the distributionof the alignment in the molecular ensemble has been taken into account and experimental results(at 2 × 1014 W/cm2 )(right), adapted from [165].

64
it possible to access the rich information in the ellipticity, which might provide another approachto study electronic structure in molecules.
Another experiment on ellipticity in HHG from N2 is reported by Mairesse et al. [154], asshown in Fig. 3.28. Our TDDFT simulation results with and without angle average are presented inFig. 3.29. Again, we see that TDDFT simulation are in excellent agreement with the experimental
data. In this case, intensity 7.5 × 1013 W/cm2 used in our simulation is just slightly lower than
the reported experimental intensity of 8 × 1013 W/cm2 . Mairesse et al. [154] analyzed harmonic

Figure 3.28: Ellipticity from N2 with 800 nm laser field at intensity of 8 × 1013 W/cm2 (left) and1 × 1014 W/cm2 (right), adapted from [154].

generation from different HHG channels due to the dynamics in the ion. As shown in Fig. 3.30,
they considered the situation when an electron, liberated from HOMO, returns to the parent ion,and may find the hole in HOMO-1 orbital, which is created by the resonant coupling of HOMO and
HOMO-1 by the 800 nm laser field for large alignment angle between molecular axis and probe laser.
In their calculations, ionic ground state X (hole in HOMO), A (hole in HOMO-1), and B (a hole in
HOMO-2) are considered. In the paper, the authors studied the relative phase between ionizationfrom HOMO and HOMO-1 and showed that a π phase difference is necessary to reproduce theexperimental results. While in the paper the laser-induced subcyle coupling between ionic states ismodeled explicitly and decoupled from the HHG process, it is interesting to study such effects withTDDFT which capture the full multielectron dynamics. In the following subsection, we are going

65

Figure 3.29: Ellipticity of high-order harmonics as a function of the molecular alignment angle andthe order of the harmonics at λ = 800 nm and peak intensity I0 = 7.5 × 1013 W/cm2 . Comparisonof results in which no average (left panel) and a proper average (right panel) over the distributionof the alignment in the molecular ensemble is made.

66
to present our simulation results which reveal the contribution from each valence shell orbital andtheir influence on the ellipticity of HHG.

Figure 3.30: Mechanism studied in [154]. Electron is ionized from HOMO and recombines to
HOMO-1. The hole at HOMO-1 is created by the 800 nm laser pulse which resonantly couples
HOMO-1 and HOMO. Figure is adapted from [154].

3.4.5

Inner shell contributions to ellipticity of harmonics

In our full TDDFT calculations the response of all electrons in the nitrogen molecule is considered. Due to the ab-initio character of TDDFT calculations and in contrast to other theoretical
approaches, as e.g. the strong-field approximation used in [154], the recombination and ionization
events cannot be separated in the analysis and the laser driven dynamics is not restricted to part ofthe pulse. In order to analyze the role of transitions to individual orbitals we have further calculated
also the responses from a subset of occupied orbitals in the molecule. In Fig. 3.31 we compare
the results of such calculations with those of the full calculations for the same laser parameters as
in Fig. 3.27. In each panel the ellipticity of the harmonics is plotted as function of the harmonicorder and the molecular alignment angle. In order to illustrate the importance of the inner-valenceorbitals no ensemble average of the molecular alignment was performed in this set of calculations.
The role of inner-shell contributions becomes obvious from the comparison of the results ofthe full calculation (upper left panel) and those obtained from the HOMO orbital only (lower left

67

Figure 3.31: Ellipticity of harmonics as function of harmonic order (vertical axis) and angle between
molecular axis and laser polarization direction (horizontal axis). Comparison of results from thefull calculation (upper left panel) and contributions from HOMO, HOMO-1 and HOMO-2 (upper
right panel), HOMO and HOMO-1 (lower right panel) and HOMO only (lower left panel). Laserparameters: wavelength of 800 nm, peak intensity of 1 × 1014 W/cm2 and pulse duration of 20 fs.

68
panel). Neither for the low-order nor for the high-order harmonics the HOMO contribution agrees
well with the full results. A significant improvement is achieved when both contributions from the
HOMO and the HOMO-1 are included in the calculations (lower right panel). However, certain
features in the results of this partial calculation, such as the strong maximum for the 21st harmonicat an orientation angle of about 30o , are not present in the full calculation. Further inclusion of
the contributions from the HOMO-2 orbital (upper right panel) finally leads to a more satisfactory
agreement with the full results. Therefore, we may point out that we observe a strong influence ofmultielectron effects and contributions from, at least, three orbitals in the nitrogen molecule over
the whole HHG spectrum. This applies, in particular, also for the low-order harmonics which could
not be investigated in the previous theoretical studies based on strong- field approximation.

Figure 3.32: Same as Fig. 3.31 but for a laser wavelength of 1000 nm.

69
In our calculations we have further noticed that a certain inner-shell channel becomes particularly important at driving laser wavelengths, for which there is a resonance between the inner-shellorbital and the HOMO. This can be seen from the results presented in Fig. 3.32 which wereobtained for a laser wavelength of 1000 nm. At this wavelength the (unperturbed) HOMO and
HOMO-2 levels in the neutral nitrogen molecule are coupled via a two-photon resonance transition.
The results show that the effect of the HOMO-2 contributions on the ellipticity of the harmonics
is even stronger than at 800 nm (Figure 3.31).

3.4.6
3.4.6.1

Ellipticity of harmonics from CO2
Comparison with experimental data for CO2

We have also investigated ellipticity of HHG from CO2 , another extensively studied molecule.
Similarly, we have done a simulation for CO2 with 800 nm laser field at intensity of 1.5×1014 W/cm2
with orientation angle from 0o to 90o at a step of 10o . To compare with experimental results[165], angle average is applied as well. The harmonic intensity ratio of perpendicular to parallel
components is calculated and presented in Fig. 3.33(b). The intensity ratio for harmonic orderfrom 17th to 23rd is in the range from 0.02-0.07 which agrees with experimental data in Fig. 3.33(a).
The orientation angle φ of the elliptically polarized high harmonics, illustrated in Fig. 3.34(a), is(b)

Figure 3.33: Intensity ratio of perpendicular to parallel component. Comparison between experimental result (a) adapted from [165], and our simulation result (b).

70also calculated and shown in Fig. 3.34(c). The orientation angle is in opposite direction to theprobe-pump angle and has a minimum at 23rd harmonic order. Such features are in agreement withdata of experimental measurements shown in Fig. 3.34(b). Experimentally observed ellipticity from(c)Figure 3.34: Illustration of high harmonic ellipticity orientation angle (a), experimental results ofellipticity orientation angle from CO2 (b), adapted from [165], and our simulation results (c).

CO2 is relatively small [165]. We also obtained the ellipticity with ensemble angle averaged results,
shown in Fig. 3.35. The ellipticity is weak except for some strong ellipticity around the cutoff.
Such weak ellipticity is partly due to the ensemble angle average effect as we have seen for N2 .
Another reason is the role of multiorbital contributions, as we will discuss in the next subsection.

Figure 3.35: Ellipticity of high harmonics as function of alignment angle for CO2 for a laser fieldat λ = 800 nm and I = 1.5 × 1014 W/cm2 .

71
3.4.6.2

Inner valence shell orbital contributions to ellipticity of harmonics

We study the contribution from each inner valence by looking at the response from HOMO aswell as adding those from deeper orbitals. The results are shown in Fig. 3.36. No angle average isperformed here. The ellipticity of high harmonics from CO2 is influenced nearly by all the valenceshell orbitals. The ellipticity of HHG from the HOMO (1πg ) is strong as shown in Fig. 3.36(a)and we also find that based on the two-point interference model, the strong ellipticity happens
around the position of destructive interference. However, after adding contribution from HOMO-1
to HOMO, the ellipticity becomes weaker, as presented in Fig. 3.36(b). As we add contributions
from inner orbitals one by one, the ellipticity is gradually getting smaller. Therefore, in this casethe multiorbital contributions gradually lead to a weak ellipticity while for N2 we did not observe
similar effect.
The reason for the strong contribution from inner valence orbitals is that HHG depends
on the ionization which is influenced not only by the ionization potential but also the alignment
angle. Although the HOMO has relatively small ionization potential and thus is easier to ionize,
its ionization rate is suppressed at 0o and 90o because of quantum destructive interference, andhas a maximum ionization rate at 45o . The HOMO-1 and HOMO-2 have very similar ionization
potential but the HOMO-1, which is σu orbital, is easier to ionize at 0o while the HOMO-2, whichis πu orbital, is easier to ionize at 90o (ionization dependence on symmetry will be discussed innext Chapter in detail). So, the HOMO-1 and HOMO-2 have strong contributions to the HHG atangle close to 0o and 90o , respectively. As for the HOMO-3 (σg ), the contributions to HHG comesfrom the strong coupling between HOMO-3 and HOMO-1 at 800 nm close to 0o . Looking at the
projection between these two orbitals in the laser field, shown in Fig. 3.37(a), there is a significant
coupling between these two orbitals (even at off resonant frequency), which means that HOMO-3
plays a role as important as HOMO-1 and thus influences the ellipticity strongly. Finally, HOMO-4
and HOMO-5 contribute to ellipticity slightly. At 800 nm, the two orbitals are however coupledsignificantly and there is 40% population transfer between these two orbitals, as seen from the right72(a) (1πg )4(b) (3σu )2 (1πg )4(c) (1πu )4 (3σu )2 (1πg )4(d) (4σg )2 (1πu )4 (3σu )2 (1πg )4(e) (2σu )2 (4σg )2 (1πu )4 (3σu )2 (1πg )4(f) (3σg )2 (2σu )2 (4σg )2 (1πu )4 (3σu )2 (1πg )4Figure 3.36: Ellipticity of harmonics as a function of alignment angle for CO2 . Starting with theresults from HOMO only (upper left), contributions from inner valence orbitals are added step bystep in the other panels (as indicated). A proper average over the distribution of the alignment in
the molecular ensemble is made. Applied laser field is at wavelength of 800 nm and intensity of
1.5 × 1014 W/cm2 with 30 fs pulse length.

73

Figure 3.37: Projection of coupled inner valence orbitals HOMO-3 (4σg ) to HOMO-1 (3σu ) (a) and
HOMO-5 (3σg ) to HOMO-4 (2σu ) (b). Laser wavelength of 800 nm and intensity of 2×1014 W/cm2
with orientation angle of 20o are considered.

panel in Fig. 3.31. Based on our simulations, we can therefore conclude that the dynamics fromand between different orbitals during the HHG in CO2 is even more complex as compared to N2 .

3.5

Summary
In this Chapter, we first have shown that fractional harmonics in the form of Mollow side-

bands occur in high-order harmonic spectra of open shell molecules at equilibrium geometries. The
fractional harmonics appear due to a one-photon coupling (Rabi flopping) between one of the innervalence orbitals and the outer open shell induced by the driving laser. The phenomenon is foundand illustrated for parallel as well as perpendicular transitions, di-, tri- and polyatomic molecules,
coupling of different type of orbitals, and different orientations of the molecule. Thus, it is shownthat Mollow sidebands are much more general than the previously studied case of stretched H+
2in the regime of charge resonance enhanced ionization. Moreover, the time-dependent electron
density and the time-frequency analysis have shown that the appearance of fractional harmonics
comes along with a nonadiabatic electron dynamics and localization in the molecule as well as abreakdown of the traditional semiclassical picture of high-order harmonic generation.
We also analyzed the ellipticity of high harmonic spectra starting from simple one-electron
molecule H+
2 to multielectron molecules N2 and CO2 . The results obtained within the timedependent density-functional theory reveal the importance of multielectron contributions from

74
different orbitals for the interpretation of recent experimental data for the ellipticity of higherorder harmonics. It is shown that the full multielectron results are in excellent agreement with theexperimental data, if contributions from inner valence orbitals over the whole HHG spectrum aretaken into account. We have further confirmed that it is essential to take proper account of the
distribution of the molecular axis with respect to the polarization direction to find agreement withthe observations.

Chapter 4

Multielectron effects in strong-field ionization from molecules

The ionization of atoms and molecules is another fundamental strong-field process. It provides
information about the system from which electron is ionized. For example, streaking and RABBITT
are experimental techniques in which measurement of ionized electrons are used to gain dynamical
information through pump-probe experiments [169, 170].
Compared to atoms, for molecules the quasistatic picture of tunneling ionization needs to bemodified to incorporate the multi-center nature of atomic potentials in molecules. There are indeed several fundamental differences between the ionization of atoms and molecules. Earlier studies[171, 172, 173] revealed experimentally that the ionization signal of molecules and noble gas atoms
with comparable ionization potentials behave similarly, for example, the ionization yields of thepair N2 and Ar as a function of intensity are similar. Later it was found that O2 , which has a similar
ionization energy as Xe, violates such rule because of the destructive interference in the electron
emission from anti-bonding πg (HOMO) orbital in O2 [30]. Therefore, the ionization of molecules
is more complex than ionization of atoms, but contains also more information about the electronic
structure. Another fundamental difference is charge-resonance-enhanced ionization (CREI) discovered for H+
2 , which occurs at a critical internuclear distance Rc , greater than equilibrium distanceRe [174, 175].
One common approximation in theory of strong-field ionization is to assume that only one
active electron responds to the external laser field, which is called, as mentioned before, singleactive-electron (SAE) approximation. Recent experimental results on ionization and dissociation

76
of Br2 show however that during the dissociation process the valence orbitals rearrange. Ionizationfrom πu orbital is observed at early dissociation stage, while ionization from σg orbital as well asother valence orbitals is observed at large internuclear distances [176]. Such phenomena, orbital
rearrangement during dissociation process and changes in the ionization yields from orbitals, arealso observed in simulations for stretched C2 H2 [177]. TDDFT simulations show that inner valenceorbitals can contribute considerably to total ionization at high intensity [63, 178]. Furthermore,in experiments on filamentation in N2 and CO2 [179, 180], it is found that there is an unusualhigh ionization yield from the inner valence orbitals compared to the ionization of HOMO, whichleads to a population inversion in the ion. Thus, we need to understand the multielectron effects inionization process. However, the theoretical treatment remains difficult due to the complexity of
direct simulation for multielectron systems performed by solving multielectron Schrödinger equation. So by now, there are still not many theoretical studies on multielectron effects in strong-field
ionization.
In this Chapter, we are going to investigate the ionization of molecules using TDDFT simulations. This gives us the opportunity to analyze the response from each Kohn-Sham orbital ofa multielectron molecule in view of the role of the orbital symmetry and the dependence on thelaser field parameters. Specifically, we will study the influence of resonant coupling between valenceorbitals on the strong-field ionization.

4.1

Role of ionization potential and orbital symmetry
Within TDDFT approach the Kohn-Sham orbitals are propagated in the laser field. Thus,

the population of the ith orbital can be obtained by integrating the corresponding wavefunction in
the simulation boxˆnks
i (t) =|φi (r, t)|2 dr,

(4.1)

box

where φi is the wavefunction of the ith Kohn-Sham (KS) orbital. We also denote the dynamical
th KS orbital explicitly with superscript ks, which is necessary to separate
population nks
i (t) of i

77
it from state population in resonant case as we will show later. The corresponding ionization
probability of ith KS orbital is equal toPiks (t) = 1 − nksi (t).

(4.2)

The ionization yield from different orbitals can strongly depend on orbital symmetry andionization potential. Ionization from inner valence orbitals is usually at least an order of magnitude
smaller than ionization from HOMO, because of the larger ionization potential as confirmed by our
simulation results in Fig. 4.1. In the simulations, a 600 nm laser field at intensity 5 × 1013 W/cm2
is applied to N2 with polarization parallel (left) and perpendicular (right) to the molecular axis. In

(a)

(b)

Figure 4.1: Ionization from KS orbitals of N2 for parallel (a) and perpendicular (b) alignment withrespect to external laser field polarization. The laser field parameters: wavelength of 600 nm andintensity of 5 × 1013 W/cm2 .

both cases, the HOMO (3σg ) is the strongest ionized orbital as one would expect due to its lowest
ionization potential. The ionizations from other orbitals are negligible in this case. However, we
can still see that the second strongest ionized orbital is different for the two cases: it is HOMO-2(2σu ) for parallel orientation while for the perpendicular alignment it is HOMO-1 (1πu ). Although
HOMO-1 has slightly lower ionization potential compared to HOMO-2, its ionization is muchsmaller than HOMO-2 in parallel case. This is because the ionization not only depends on theionization potential but also strongly on the orbital symmetry and molecular orientation. The
symmetry of HOMO-1 influences its ionization in two aspects. HOMO-1 is a πu orbital (see Fig.
3.4) and it possesses a nodal plane along the molecular axis. If the polarization of laser field is78
parallel to the molecular axis, the ionized wavepackets from parts of the orbital above and below
the nodal plane interfere destructively with each other due to the opposite phase and thus suppress
the ionization, while this is not the case for perpendicular alignment. The ionization of HOMO-2(σu orbital) is also influenced in this way but destructive interference occurs in the perpendicularcase. The second effect comes from the density distribution, namely, how much it is extended
from the center. The σ orbitals are more extended along molecular axis, while π-orbitals are more
extended in perpendicular direction. Because of these factors, HOMO-2 is easier to ionize forparallel molecular orientation, while the HOMO-1 is easier to ionize for perpendicular alignment.
We can also see such effects by looking at the ionization from HOMO. In Fig. 4.1, the two figures
are generated for the same parameters except the orientation angle. The ionization probability of
HOMO in parallel case is obviously larger than that for perpendicular orientation.

4.2

Resonance enhanced ionization of open-shell molecules
We have previously studied the resonance induced fractional harmonics in open-shell molecules.

We have shown that for N+
2 with one spin-up electron removed from HOMO, the resonant laser field
at wavelength of 400 nm strongly couples HOMO-2 and HOMO and can enhance the ionizationfrom HOMO-2 strongly, shown in Fig. 3.13(d). For other open-shell molecules, in the same coupling
regime where we have seen fractional harmonics, we can also observe the enhanced ionization. The
results are shown in Fig. 4.2. For NO2 in Fig. 4.2(a) where the HOMO (|6a1 ↑i) is not occupied,the 400 nm laser field induces a coupling between HOMO-1 (|4b2 ↑i) and HOMO (|6a1 ↑i) andthus |4b2 ↑i has the relative strongest ionization compared to that of other orbitals. For CO+
2 in
Fig. 4.2(b), where an electron in |1πg ↑i is removed, a 350 nm laser field induces a coupling between
HOMO and HOMO-1 in parallel direction to the molecular axis. The ionization of valence orbitals
are more complex. There are three different types of orbitals strongly involved in the ionization
process. For the spin-down components, the ionization behavior is like what we expect based on
the ionization potentials and symmetries of valence orbitals. The HOMO-1 (|3σu ↓i) is strongly
ionized at parallel orientation compared to the other two types of orbitals: HOMO (|1πg ↓i) and79(a) NO2 with |6a1 ↑i unoccupied. λ=400 nm,
2
I = 2 × 1014 W/cm , |4b2 ↑i and |6a1 ↑i arecoupled.

(b) CO+
2 with |1πg ↑i partially occupied. Laser
2parameters: λ=350 nm, I = 1 × 1014 W/cm .
|1πu ↑i and |1πg ↑i are coupled. For the πg
and πu orbital, ionization is normalized by initial number of electrons.

(c) C2 H+
4 with |1b3g ↑i unoccupied. Laser pa2
rameters: λ = 400 nm and I = 1 × 1014 W/cm .
|1b3g ↑i and |1b2u ↑i are coupled.
+
Figure 4.2: Ionization of KS orbitals of openshell molecules: NO2 (a), CO+
2 (b), and C2 H4 (c).
The laser parameters and coupling orbitals are illustrated in the captions below each figure.

80
HOMO-2 (|1πu ↓i). For the spin-up component, the HOMO (|1πg ↑i) is only partially occupied, thecoupling between HOMO-2 (|1πu ↑i) and HOMO (|1πg ↑i) leads to a strong ionization of HOMO-2(|1πu ↑i). For C2 H+
4 in Fig. 4.2(c), where HOMO-1 (|1b3g ↑i) is unoccupied, the ionization fromHOMO (|1b3u i) is still the strongest ionized orbital due to its lowest ionization potential. A 400 nm
laser field induces a coupling between HOMO-1 (|1b3g ↑i) and HOMO-3 (|1b2u ↑i), which enhances
the ionization of HOMO-3 (|1b2u ↑i). But its ionization yield is still relatively smaller than that of
HOMO.
Such resonance enhanced multiphoton ionization (REMPI) is a very important process formolecules and it is often used in experiments as a spectroscopy tool for mass spectrometry [181].
Besides REMPI set up for open-shell molecules, we are also interested in studying enhanced ionization from inner-shell orbitals for closed-shell molecules. In the next section we therefore present
results for ionization of N2 at equilibrium internuclear distance.

4.3

Multiorbital contributions and enhanced ionization from inner-shellorbitals of N2
In section 4.1, we have seen a typical case in which the ionization probability of N2 is domi-

nantly determined by the orbital with the lowest ionization potential. In such case, it is reasonable
to assume that only one active electron or one orbital responds to the external field. However, wehave seen before for HHG that the resonant coupling between orbitals can lead to new effects. In
particular, we will now analyze if there is any effect on ionization when the laser frequency is tuned
to the energy difference for two occupied orbitals in neutral molecules.
The energy differences between HOMO-2 and HOMO, HOMO-1 and HOMO of N2 correspondto wavelengths of about 400 nm and 800 nm, respectively. These are typical laser wavelengths usedin experiments nowadays. We have carried out simulations at these wavelengths and studied theionization of N2 with 400 nm laser field with polarization parallel to molecular axis and 800 nm
laser field with polarization perpendicular to the molecular axis. The time-dependent ionizationsof valence KS orbitals are presented in Fig. 4.3. Compared to the case with 600 nm laser field81
shown in Fig. 4.1, we notice a significantly enhanced ionization from inner orbital HOMO-2 (2σu )

0.200

2σu
1πux
1πuy
3σg

0.150
0.100
0.050
0.00004

8 12Time [fs]16

Ionization probability

Ionization probability

for parallel orientation and HOMO-1 (1πu ) for perpendicular alignment.

0.080
0.060
0.040

2σu
1πux
1πuy
3σg

0.020
0.0000 4 8 12 16 20 24 28
Time [fs]Figure 4.3: The ionization probability of valence KS orbitals of N2 . Only the relevant KS orbitals
which possess an important response to the laser field are shown with their labels. Left: Laser
polarization is parallel to the molecular axis at 400 nm and intensity of 2 × 1014 W/cm2 , inducing acoupling of HOMO (3σg ) and HOMO-2 (2σu ). The curves for 1πux and 1πuy overlap because they
are degenerate. Right: Laser polarization is perpendicular to the molecular axis at wavelength of
800 nm and intensity of 2 × 1014 W/cm2 , inducing a coupling between HOMO (3σg ) and HOMO-1(1πu ).

The results shown above are calculated with LDA functional, and field having trapezoidal
envelope with fixed nuclei during the propagation of the molecular wavepacket. The LDA functionaldoes, however, not have a correct 1/r tail in the potential. Thus, we would like to clarify that suchenhanced ionization from inner shell orbitals is not caused by the incorrect asymptotic potential
tail. Since LB94 functional [112] has a correct 1/r tail, we have also performed simulations with thisfunctional for the same laser parameters at wavelength of 400 nm and intensity of 5 × 1013 W/cm2
and polarization parallel to the molecular axis of N2 . The results are presented in Fig. 4.4. As
expected, because of the correct 1/r tail, the ionization calculated with LB94 functional (right) issmaller than the result from LDA calculation (left). But the ionization behavior and the relative
orbital contributions in both cases are very similar. The ionization of coupled orbitals (2σu and
3σg ) are dominant in both cases and the relative ionization yield is also similar. Therefore, the 1/r
tail has no influence on the enhanced ionization from inner orbitals. In our calculations, trapezoidal
envelope of the laser pulse is used, because we need a pulse of relatively short duration to reduce

82

Figure 4.4: Comparison of ionization from N2 between LDA (left) and LB94 (right) functionals.
Laser field has a wavelength of 400 nm and intensity of 5 × 1013 W/cm2 with polarization parallelto the molecular axis of N2 . Only the relevant KS orbitals which exhibit an significant response to
the laser field are shown with their labels.

83
the simulation time. However, trapezoidal envelope has a point of discontinuity in the derivative,

0.040
0.030
0.020
0.010
0.0000 4 8 12 16 20 24 28
Time [fs]

2σu
1πux
1πuy
3σg

Ionization probability

Ionization probability

which is not physical and might produce some unphysical results. We have checked this carefully

0.040
0.030
0.020
0.010
0.0000

2σu
1πux
1πuy
3σg

4

8 12 16 20
Time [fs]Figure 4.5: Comparison of ionization from N2 between sin2 (left) and trapezoidal (right) envelope.
Wavelength is 400 nm and intensity is 8 × 1013 W/cm2 with polarization parallel to the molecularaxis of N2 . Only the relevant KS orbitals which possess an important response to the laser field
are shown with their labels.

by comparing the result between sin2 and trapezoidal envelopes in Fig. 4.5. Compared withtrapezoidal envelope (right), the relative ionization from HOMO-2 to HOMO is still enhanced, but
not as strong as for trapezoidal envelopes, which is understandable because the amplitude of laser
field for sin2 envelope changes with time which influences the coupling.
The results shown above are obtained in simulations with fixed nuclei during the propagationof the molecular wavepacket. Since the molecule vibrates and the corresponding movement of thenuclei plays a role in the electron dynamics, we would also like to see if the results change forsimulations without fixed nuclei. If the effect disappears then such enhanced ionization is notmeasurable in an experiment. To this end, the nuclei are propagated on a mean potential surface
corresponding to an electronic superposition state [182]. The results is shown in Fig. 4.6, theenhanced ionization from HOMO-2 is still present for simulations with moving nuclei. The relative
ionization of HOMO-2 to HOMO is not as strong as that for fixed nuclei because the energy gap
between the two states can change slightly with the internuclear distance.
The way we calculate the ionization from Kohn-Sham orbital based on Eq. (4.2) need tobe carefully examined when a coupling is present. Without any coupling, the wavefunction of theIonization probability

84

0.012
0.010
0.008
0.006
0.004
0.002
0.000
0 4 8 1216202428
Time [fs]

2σu
1πux
1πuy
3σg

Figure 4.6: Ionization of KS orbitals of N2 with moving nuclei. LB94 functional and trapezoidal
envelope are used in simulation. The laser field is set at wavelength of 400 nm and intensity of
5 × 1013 W/cm2 with polarization parallel to the molecular axis.

85
KS orbital remains in the initial state during the propagation. In contrast, when a coupling ispresent, the wavefunction will periodically evolve to the other coupled state. Therefore, when we
calculate the population in a state, integration over the KS orbital will include the population ofboth coupled states. A better way to calculate the population in one of the coupled states is toproject the KS orbital to that state
nst
i =|

N
X

hφi (t = 0)|φj (t)i|2 ,

(4.3)

j

where we add superscript st for the population to separate it from the population of KS orbital.
Actually, for coupled states, we only need to sum over the coupled KS orbitals. For non-coupledks
state, nst
i will give the same value as ni . On the other hand, the KS orbital φi (t) at the time ofnonzero electric field can be strongly deformed by the electric field and thus we can project the KS
orbitals φi (t) at the time of zero-electric field only. Then the ionization can be obtained via
Pist (t) = 1 − nst
i (t).

(4.4)

We have calculated the state ionization according to Eq. (4.4). The results for the same situation

0.20

2σu
1πux
1πuy
3σg

0.15
0.10
0.05
0.000

4

8 12
Time [fs]16

Ionization probability

Ionization probability

as in Fig. 4.3 is presented in Fig. 4.7. Compared to Fig. 4.3, we get similar ionization yields

0.08

2σu
1πux
1πuy
3σg

0.06
0.04
0.02
0.000

4

8 12 16 20
Time [fs]Figure 4.7: Ionization yields of states of N2 at the time of zero electric field. Only the relevant
states which possess an important response to the laser field are shown with their labels. Left:
Laser field is parallel to the molecular axis at 400 nm and intensity of 2 × 1014 W/cm2 inducing acoupling of HOMO (3σg ) and HOMO-2 (2σu ). Right: Laser field is perpendicular to the molecularaxis at 800 nm and intensity of 2 × 1014 W/cm2 inducing a coupling of HOMO (3σg ) and HOMO-1(1πu ).

86
for each state, which means such correction is not significant in our simulations. From Fig. 4.7,
we see that the strongest ionized state for N2 is HOMO (3σg ) and HOMO-2 (2σu ) for parallelorientation (left), while for perpendicular orientation the strongest ionized state are HOMO (3σg )and HOMO-1 (1πu ) (right), which confirms the fact that there is an enhanced ionization from innervalence states when a coupling is induced by laser field.
Such an enhanced ionization from inner valence orbitals arises due to the laser induced
coupling. To support this conclusion, we have done a systematic analysis by scanning the laser
wavelength from 200 nm to 1000 nm and show in Fig. 4.8 the ionization ratio between the ionization
yield for the two coupled orbitals in both cases (panels left for parallel alignment and panels rightfor perpendicular alignment). From the upper two panels in Fig. 4.8, one can see that the ionizationfrom inner valence orbital is enhanced over a relatively wide wavelength range with a peak ratio
close to 1 at the optimal wavelength for the coupling. We also studied the ionization ratio at
coupling wavelength as function of laser intensity. As we can see from the lower panels in Fig. 4.8,
the two coupled orbitals have nearly the same ionization yield and increase simultaneously with theincrease of laser intensity, which also indicates the coupling instead of the intensity indeed plays arole here.
In order to gain further understanding, we study the instantaneous eigenvalues of valenceorbitals in the laser field. In our simulations, after each time-propagation step, the eigenvalueof each orbital is calculated and the results for the two cases in N2 are shown in Fig. 4.9. The
eigenvalue of orbitals that do not couple with any other states oscillate in phase with the oscillationsof the electric field of the laser. On the other hand, the eigenvalues of all orbitals have a decreasing
trend, which can be understood by the fact that with increase of ionization in the system, it
becomes harder to ionize the system further and thus the eigenvalues become smaller. For thetwo coupled orbitals, their eigenvalues oscillate with opposite phases and around the same mean
average value with the same oscillation amplitude. Since the ionization potential is proportional tothe eigenvalue of each orbital, at every transition cycle, the two orbitals are equally ionized. Suchoscillating ionization of each orbital can also be seen from Fig. 4.3. Based on this argument, such

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
200 300 400 500 600 700 800 900 1000
Wavelength [nm]Ionization ratio HOMO-1/HOMO

Ionization ratio HOMO-2/HOMO

87

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
200

400

600
800
Wavelength [nm]

1000

0.38
0.40
0.42
0.44
0.46
0.48
0.50
0.520 2 4 6 8 10 12 14 16
Time [fs]

2σu
1πux
1πuy
3σg

Eigenvalue

Eigenvalue

Figure 4.8: Upper pannels: Ratio of ionization yields between HOMO-2 and HOMO (left) and
HOMO-1 and HOMO (right) of N2 as a function of laser wavelength at intensity of 5 × 1013 W/cm2
and a pulse length of 14 fs. Lower panels: Ionization yields for HOMO and inner valence orbitals
as a function of intensity at wavelength of 400 nm (left) and 800 nm (right). Panels on left are forparallel alignment, while panels on right are for perpendicular alignment.

0.38
0.40
0.42
0.44
0.46
0.48
0.500

2σu
1πux
1πuy
3σg

5 10 15 20 25
Time [fs]Figure 4.9: Instantaneous eigenvalues of Kohn-Sham Hamiltonian for N2 valence orbitals for 400 nm
laser field at parallel orientation (left) with intensity of 5 × 1013 W/cm2 and 800 at perpendicularorientation (right) with intensity of 5 × 1013 W/cm2 . Only the relevant KS orbitals, which exhibit
significant response to the laser field, are shown with their labels.

88mechanism can be detected by measuring the number of excited ions (N+2 with a hole in HOMO-2)
as a function of time in a pump-probe setup. At times longer than the oscillation period, the twoorbitals on average have the same ionization yield.
The reason why the coupled orbitals’ eigenvalues oscillate around each other can be analyzed
by looking at the projection of the KS orbitals to the other initial states |hφj (t = 0)|φi (t)i|2 . From­®
| 3σg (t =0)|2σu (t) |2

1.0
0.8
0.6
0.4
0.2
0.00 2 4 6 8 10 12 14Time [fs]Figure 4.10: The projection of propagated orbital 2σu (t) onto the coupled initial orbital 3σg (t = 0).
Laser field is at wavelength of 400 nm and intensity of 5 × 1013 W/cm2 for parallel orientation ofmolecular axis.

the projection in Fig. 4.10, we can see that the propagated KS orbital 2σu (t) periodically evolves
to the other state. Therefore, in the laser field, the coupled KS orbital becomes a mixture of thesetwo coupled states.
For the states in the laser field, it is often convenient to use photon-dressed states within the
Floquet picture. Floquet states or laser dressed states have been successfully used to describe theelectron response to a strong laser field [3, 183] and are capable to provide an intuitive picture formany physical processes. In the Floquet picture, each state |φi has an infinite number of copies thatare coupled to n number of photons, denoted as |φ, ni. For the resonance case in our simulation,
the two coupled dressed states |3σg , ni and |2σu , n + 1i become degenerate in the laser field andform two new photon-dressed states which are mixed states of the original states. According tothis picture, ionization from the new photon-dressed state indicates that the electron is ionizedsimultaneously from the two original coupled states. In our simulations, the ionization yields fromthe two coupled states are indeed close to each other but not exactly. This deviation is due to the89finite pulse duration, since the Floquet picture is strictly applicable for fields of constant amplitudeonly. An approach to include Floquet picture in the TDDFT formalism will be discussed in Chapter
6.

4.3.1

Total ionization

Since in the case of resonant coupling the ionization is strongly enhanced from an inner valence
state, a general question is whether the total ionization is also enhanced and if the consequencesof the coupling can be measured in the experiments. This is difficult to check because it is notpossible to switch off the coupling effects during the interaction of the laser pulse with the molecule
within our simulations. If we force to keep all the inner shell orbitals except the HOMO frozen,
then the system is not physical because the HOMO is coupled to HOMO-2 and can result in an
unphysical population larger than 2. Instead, we look at the change of the total ionization yield
as a function of laser wavelengths as shown in Fig. 4.11. In the corresponding calculation, a 14 fs
laser pulse is applied with trapezoidal envelope which has 10 fs constant plateau and 2 fs switch on
and off regions. The total ionization decreases smoothly with increasing wavelength. But at the

Figure 4.11: Total ionization yield from N2 as a function of wavelength. The laser electric field isaligned along the molecular axis. Intensity is equal to 5×1013 W/cm2 .

resonance wavelength we observe indeed a slightly enhanced ionization.

90
4.3.2

Ionization of N2 in filamentation experiments

Enhanced ionization from HOMO-2 of N2 in interaction with 400 nm laser field has been
observed in experiment [184]. In the experiment, a 400 nm laser field with power of 2.3 mJ
propagated in air to generate a filament, and the fluorescence from the filament is measured. As
a comparison, the fluorescence from filament generated by 800 nm laser field with the same powerhas been measured. The fluorescence spectra normalized to the same interaction volume are shownin Fig. 4.12. A strong peak from the first negative band, corresponding to the transition from 3σg

Figure 4.12: Measured fluorescence spectra obtained from laser filaments in air and normalized tothe same interaction volume for λ= 800 nm (black) and λ= 400 nm (gray). Fluorescence from thefirst band obtained for filaments at 400 nm is about 6.4 times larger compared to that at 800 nm.
Adapted from [184].

to 2σu in N+
2 , is observed for 400 nm which is about 6.4 times of the signal obtained with the 800
nm laser field. This experiment result suggests an efficient inner-shell ionization for the 400 nm
laser field. Such observations agree qualitatively with our simulation results.
The enhanced ionization from an inner valence shell also gives us a deeper insight into themechanism behind recent observation of abnormal population inversion between ground and excitedionic states in filamentation experiments on N2 [185, 180, 186, 179]. In the experiment [186], a 800
nm laser field with a power of 3.8 mJ is focused into a chamber full of N2 . A filament is createdalong the laser propagation direction, in which N2 is strongly ionized. Then, a weak probe pulse
at 400 nm with power of 0.1 µJ is applied and collimated with the 800 nm laser, and the forward

91
fluorescence spectra is measured as shown in Fig. 4.13. The interesting result is that after applying

Figure 4.13: The spectra captured with (blue) and without (green) the pump pulse. Laser parameters Adapted from [186].

the 800 nm laser field, the weak 400 nm laser field acting as a seed pulse leads to a strong stimulated
emission at 390 nm which correspongds to the transition from HOMO (3σg ) and HOMO-2 (2σu ).
The other peak around 428 nm corresponds to the transition between the same electronic states
but different vibrational states. The stimulated emission indicates that there is a strong ionizationfrom HOMO-2 and thus population inversion was created for N+
2 after applying the 800 nm laser
field. Such results are hard to understand because the HOMO-2 is deeply bound. Different models
were proposed to explain such abnormal behavior, however, the debates about the mechanism isongoing [187]. Our simulation results show that for the coupling wavelength at about 400 nm,
HOMO-2 can indeed be strongly ionized, which helps us understand better how the deeply boundorbital HOMO-2 can participate in the ionization process.
Since in the experiment the ionization from HOMO-2 is due to the 800 nm laser field which
is different to our previous simulations, we have performed another simulation for a two-color laser
pulse interacting with N2 (a weak 400 nm laser field at 1 × 1012 W/cm2 and a strong 800 nm laser
field at intensity of 1 × 1014 W/cm2 ) to help us confirm the mechanism. We would like to point out
here that the weak 400 nm laser field in our simulation is not to simulate the same effect as in theexperiment. The 400 nm laser field in the experiment only acts as a seed to initiate the stimulation

92
process. The 400 nm laser field in our simulation is considered to arise from the filament whichbroadens the spectrum of the pulse. Thus, we model the contribution of 400 nm in the filament
with a second weak laser field at 400 nm in the simulations. Our simulation results are presented
in Fig. 4.14. Although the 400 nm laser field is very weak, there is still a strong ionization from

Figure 4.14: Ionization yield of N2 orbitals with two-color laser field: 400 nm with intensity of1012 W/cm2 and 800 nm with intensity of 1014 W/cm2 . Both laser field are in parallel to themolecular axis.

inner shell orbital. Therefore, our simulations suggest that the spectral broadening of the laser
pulse in the filament might play an important role in inducing coupling between the orbitals andan enhanced ionization.

4.4

Ionization from other molecules
Based on our further simulation results, the enhanced ionization from inner shell orbitals is avery general phenomenon that is not only observed for N2 , but also for other molecules such as CO2 ,C2 H4 , and C2 H6 . For CO2 , which has an electronic configuration (1σg )2 (1σu )2 (2σg )2 (3σg )2 (2σu )2(4σg )2 (1πu )4 (3σu )2 (1πg )4 , we have found that a coupling between 1πu and 1πg can be inducedwith a 350 nm laser field at parallel orientation. The ionization of the strongly adapted states
are shown in Fig. 4.15. For CO2 , the dominant ionization at this wavelength is not the HOMO

93

­
®
| πg (t =0)|πu (t) |2

1.0
0.8
0.6
0.4
0.2
0.0
0 2 4 6 8 10 12 14Time [fs]Figure 4.15: Ionization of valence orbitals of CO2 at 350 nm laser field and intensity of
1 × 1014 W/cm2 in parallel to the molecular axis. Only the relevant KS orbitals which possess
an important response to the laser field are shown with their labels.

94(1πg ) but HOMO-1 (3σu ). This is because at parallel orientation, HOMO-1 (3σu ) is easier toionize than HOMO (1πg ) and HOMO-2 (1πu ). The ionization from HOMO-2 (1πu ) and HOMO
(1πg ) is suppressed due to the symmetry of the orbitals. We note that for CO2 , there also exist
filamentation experiments [180], in which inversion population is observed, corresponding to an
enhanced ionization from HOMO-1, in agreement with our results.
We also studied other polyatomic molecules like C2 H4 . A 400 nm laser field is applied
parallel to the carbon-carbon bond, which couples HOMO-3 (1b2u ) and HOMO-1 (1b3g ), as canbe seen from the projection in Fig. 4.16. HOMO (1b3u ) is still the dominantly ionized orbital.
However, we see that the ionization from HOMO-3 is enhanced and becomes comparable to that of
HOMO-1. Compared to ionization from HOMO, their ionization yields are relatively small, due tothe symmetry of both HOMO-3 and HOMO-1, which suppresses the ionization along the carboncarbon bond axis. Thus in this case, although a coupling is induced and enhances the ionization,
the effect remains relatively small. For C2 H6 , a laser field with wavelength of 200 nm and intensity­®
| 1b3g(t =0)|1b2u(t) |2

1.0
0.8
0.6
0.4
0.2
0.00

5

10
15
Time [fs]20

Figure 4.16: Ionization of valence orbitals of C2 H4 (left) and projection between the coupled states(right). The laser field is at wavelength of 400 nm and intensity of 1×1014 W/cm2 with polarization
perpendicular direction to the molecular axis but in parallel to the plane where the molecule lies.
Only the relevant KS orbitals which possess an important response to the laser field are shown with
their labels.

of 1 × 1014 W/cm2 is applied to induce coupling between HOMO-1 (3a1g ) and HOMO-3 (2a2u )with polarization parallel to the carbon-carbon bond. As shown from Fig. 4.17(left), the top threestrongly ionized orbitals are HOMO (1eu ), HOMO-1 (2a2u ), and HOMO-3 (2a1g ). The ionization

95
from HOMO is actually slightly suppressed due to its symmetry. The ionization of HOMO-3 isstrongly enhanced and approaches ionization rate of HOMO-1 as a result of the coupling.

­®
| 3a1g(t =0)|2a2u(t) |2

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.00

5

10 15
Time [fs]20

Figure 4.17: Ionization of KS orbitals of C2 H6 (left) and the projection of coupled states. Laser fieldparameters: 200 nm and 1 × 1014 W/cm2 , electric field is parallel to the carbon-carbon molecularaxis. Only the relevant KS orbitals which exhibit significant response to the laser field are shownwith their labels.

4.5

Summary
In this Chapter, we have shown that the ionization yield strongly depends on orbital symmetry

and ionization potential. Ionization from inner valence state is usually orders of magnitude smaller
than ionization from HOMO. However, this does not hold true when a laser field induces a strong
coupling between inner valence orbitals and HOMO. In that case, strong ionization from innervalence orbitals can be observed. Such phenomenon is rather general and can be observed for
di- and polyatomic molecules. Results of our simulations provide better understanding how innervalence orbitals participate in the ionization process. The results are useful to help us analyze
observations in filamentation experiments. As ionization is strongly related to dynamical electron
localization, we will study electron localization in the coupling regime in the next Chapter.

Chapter 5
Nonadiabatic electron localization 1

Since the invention of laser in 1960s, research on light-matter interaction in quantum systemsuch as atoms, molecules and solids extends our knowledge about the mechanism of electron response to the light field. As a consequence, new nonlinear nonperturbative phenomena like ATI,HHG were discovered, as discussed in previous Chapters. Usually, it is assumed that electrons
follow the change in the oscillating electric field of the laser adiabatically. This means, for example, that the electron density will swap periodically from one side to the other in an atom ora molecule. On the other hand, more recently, dynamic electron localization and nonadiabatic
behavior has been seen for such high intense laser fields. Laser induced nonadiabatic dynamics
can influence the excitation and ionization behavior of molecular systems. Therefore, nonadiabatic
electronic effects must be taken into account when attempting to image or control molecular dynamics through ionization, HHG or fragmentation. The nonadiabatic effects are also prominent
in processes such as vision [189], intersystem crossing [190, 191, 192], and many photochemical
reactions [193, 194, 195]. In order to achieve the goal of coherent control of such processes, it isimportant to study and understand nonadiabatic effects.
Adiabatic electron response is usually assumed in atoms, especially in rare gas atoms, due
to the large spacing of electronic energy states. Compared to atomic systems, the electronic statesin molecules are usually much closer to each other and thus a nonadiabatic behavior is easier toinduce and observe in molecules. The nonadiabatic behavior is also easier to observe in large
1Part of the results in this Chapter are presented in the review [188].

97
molecules because the time for an electron inside a large molecule to traverse is longer, in favor of
inducing nonadiabatic response to external laser field [40, 41, 42]. Such nonadiabatic response inlarger, conjugated molecules, like linear conjugated systems and aromatic molecules were studied
by Lezius et al. [40, 41] and Markevitch et al. [196, 197], respectively. Lezius et al. examined themass spectra of three linearly fully conjugated all-trans-hydrocarbons of increasing length. They
showed that only two peaks are found for the smallest molecules indicating stable production of
the first two charge states of ions, resulting in little fragmentation. On the other hand for large
molecules the mass spectrum shows extensive fragmentation of molecules, indicating nonadiabatic
electron dynamics. However, the spectrum also strongly depends on the wavelength, as shown in
Fig. 5.1. Stronger fragmentation at wavelength of 800 nm is observed than longer wavelength of

Figure 5.1: Ionization and fragmention patterns represented by the mass spectra of hexatriene,
decatetraene and β-carotene for the interaction with laser pulse at 800 nm and 1450 nm. Adaptedfrom [40].

1450 nm. It was pointed out by Lezius et al. that at 800 nm, the duration of the field oscillation
was insufficient for the electron to traverse the large molecules adiabatically, while the 1450 nm
field can afford the electron wavepacket with sufficient response time. Intramolecular nonadiabatic
behavior is also found for small conjugated molecules [198, 54, 199, 200, 201]. For example, the98interaction of the anthracene 1+ cation with a strong electric field of laser leads to nonadiabatic
instantaneous dipoles and such nonadiabatic behavior becomes stronger at high intensities [200].
Even for simple one-electron system H+
2 , strong nonadiabatic electron localization can beobserved. Zuo and Bandrauk [174, 163] investigated the coupling between charge-resonant states
1σg and 1σu for different internuclear distance. At large internuclear distances, these two states
become nearly degenerate and strongly couple to the laser field and can lead to CREI at somecritical internuclear distance [202, 174, 163]. What is important is that CREI coincides with theregime where nonadiabatic behavior is present. Previous studies have shown that the degeneracyof the relevant states can lead to the result that population is preferentially trapped in the uphillwell [175]. The importance of such nonadiabatic behavior becomes obvious when one examines the
temporal dynamics of the electron in the CREI regime. The time resolved studies of ionizationof H+
2 and electron behavior demonstrate that the ionization was not only enhanced but that thetiming of peak ionization could be changed as well as shown in Fig. 5.2. Unlike ordinary case, thelocalized electron does not follow the laser field adiabatically from one nucleus to the other. The
electron dynamics instead exhibits two strong localization events on each side for every half laser
cycle as shown in Fig. 5.2 (c). This led to the discovery of the so-called multiburst ionization [203].

The study of the response of electron to laser field can ultimately provides a convenient tool
for manipulation of the intramolecular electron dynamics and control of the molecular dynamics or
potentially provides new tools for control of chemical reactions. In this Chapter, we are going tostudy the electron dynamics in the coupling regime where enhanced ionization as well as fractional
harmonics are observed as discussed in previous Chapters.

5.1

Coupling induced nonadiabatic electron dynamics
As we have discussed, fractional harmonics can be generated for N+
2 with 400 nm laser field

for parallel orientation due to the additional Rabi transition induced electron dynamics. At the
coupling regime, the laser-induced instantaneous dipole behaves nonadiabatically, as it has been

99

Figure 5.2: Multiburst ionization from H+
2 at internuclear distance of R=7.0 a.u. along the molecular axis in a strong laser field (shown in d) using (a) 1D model and (b, c) 2D model. Adaptedfrom [203]100
shown in Fig. 3.5(b). Since the coupled orbitals are 2σg and 3σu , which can lead to electron
localization, we would also expect strong nonadiabatic effects in this case. In order to visualize theelectron localization directly, we analyzed the density of the HOMO-2 (spin-up) along the molecularaxis by integrating over the other two dimensionsˆρ2σu ↑ (z, t) = dxdy|φ2σu ↑ (r, t)|2 .

(5.1)

The density plots of the coupled orbital HOMO-2 of N+
2 along molecular axis are presented in

Figure 5.3: The electron density of HOMO-2 of N+
2 along molecular axis for different laser wavelength: 400 nm (left), 600 nm (middle), and 800 nm (right). The white curve is the electron
population transfered from HOMO-2 to HOMO.

Fig. 5.3 for different laser wavelengths: 400 nm (left), 600 nm (middle), and 800 nm (right). The
electron localization at resonant case has a strong nonadiabatic behavior but becomes adiabatic,
when the laser frequency is off-resonant (800 nm). For the near resonant wavelength of 600 nm, theelectron density behaves more adiabatically but still presents some nonadiabatic effect since it does
not oscillate symmetrically. This confirms that the nonadiabatic behavior of the time-dependent
electron density comes from the coupling between states 2σu and 3σg . Such nonadiabatic behavior
can be understood by decomposing it into two competing processes. The 400 nm laser field first
couples the two states, which leads to nonadiabatic localization as part of the electron dynamics andconsequently the localized electron oscillates along the molecular axis at Rabi transition frequencyΩr . Second, the laser field itself can drive the electron periodically along the molecular axis withthe laser frequency. The second process is indeed the adiabatic response of the electron to theelectric field of the laser pulse. In Fig. 5.3(left), we also plotted the projection from HOMO-2 to
HOMO which has a frequency of Ωr . In the nonadiabatic case, the electron is basically oscillating

101
from one side of the molecule to the other side with the Rabi frequency. On the other hand, there
is also a modulation in the electron density corresponding to the peaks of the laser electric field.
Thus, the electron behavior can be understood as a superposition of such two processes. In thenext subsection we discuss the effect of the laser intensity on the nonadiabatic electron dynamics.

5.2

Dependence of nonadiabatic behavior on laser intensity
As we mentioned at the outset of this Chapter, nonadiabatic behavior is mainly found forinteraction at high laser intensities. Therefore, it is interesting to analyze how the electron behavior
changes with intensity. To this end, we have done simulations for N+
2 at 400 nm laser field withparallel orientation, which will induce the coupling between HOMO-2 and HOMO. Four different
intensities are investigated: 1×1013 W/cm2 , 5×1013 W/cm2 , 1×1014 W/cm2 , and 3×1014 W/cm2 .
The results are presented in Fig. 5.4. The electron behavior is very different at different laser

13 W/cm2 (a),Figure 5.4: Electron density of HOMO-2 from N+
2 at different laser intensity: 1 × 10
2
2
2
13
14
14
5 × 10 W/cm (b), 1 × 10 W/cm (c), and 3 × 10 W/cm (d). Laser wavelength is fixed at
400 nm.

102
intensities. At relatively low intensities of 1 × 1013 W/cm2 , the electron is almost following thelaser field adiabatically. However, when laser intensity increases, nonadiabatic behavior begins
to show up. As we mentioned above, the interaction with the laser field leads to two competing
processes: electron oscillates at Rabi frequency due to the coupling effect and electron oscillates
following the electric field. Such two processes depend strongly on the electric field amplitude.
At relative low intensity 1 × 1013 W/cm2 , the electron’s behavior is dominated by the second
process and thus behaves more adiabatically. While at stronger laser field, such two processes arecomparable to each other which results in a strong nonadiabatic response to the laser field andmore complex dynamics.
Our simulation results indicate that the coupling between different states may not necessarily
always result in nonadiabatic behavior. At low intensity, the electron can still behave adiabatically,following the laser field. The two competing processes govern the behavior of the electron inside
a molecule and the nonadiabatic response is easier to see at strong intense laser field. Throughthe control of the Rabi frequency via the electric field strength, one can control the nonadiabatic
dynamics. All of the results that have been shown so far are obtained for one laser pulse. In thenext section we present results for interaction of the molecule with bichromatic laser pulse with
one of the frequencies tuned to induce the coupling leading to nonadiabatic dynamics.

5.3

Nonadiabatic localization in two-color laser field
While the two competing processes in previous section are two aspects of the dynamics driven

by the same laser pulse, it is interesting to separate such two processes to control them separately. In
this section, we therefore show results of simulations with two-color laser field: a 400 nm (resonant)weak laser field at intensity of 1 × 1012 W/cm2 and a 800 nm (nonresonant) intense laser field at
intensity of 1×1014 W/cm2 . The simulations were designed to use the 400 nm laser pulse to control
the coupling effect, while using the 800 nm laser field to control the second process. Of course, thetwo effects still cannot be totally separated because the resonant laser field still has influence on
both processes. However, we hope such two-color laser field setup can enable us to control it more

103or less because of the extra control parameters.
We again study the system N+
2 and the laser polarization of both fields are parallel to themolecular axis. The results are presented in Fig. 5.5. Since the intensity of the 400 nm laser field

Figure 5.5: The electron density plot from HOMO-2 of N+
2 in a two color laser field set up: 400
nm laser field at I = 1 × 1012 W/cm2 and 800 nm laser field at I = 1 × 1014 W/cm2

is two orders of magnitude smaller than the intensity of the 800 nm field, the second process ismainly influenced by the 800 nm laser field. A strong nonadiabatic response can be still seen inthis case, as compared to the previous results for a single laser field at wavelength of 400 nm andlowest intensity of 1 × 1013 W/cm2 , where the electron behaved more adiabatically. Therefore,
relative strength of the two competing processes indeed strongly depends on the laser wavelength.
What is more interesting is that we have found that during a relative long period (about 20 fs),
the electron is on average located more frequently on one side (Fig. 5.5(left)) or on the otherside (Fig. 5.5(right)), during the interaction with a 90 fs laser pulse. Such periodic localizationpreference from one side of the molecule to the other side can be analyzed via the projection of thetwo coupled orbitals, as for example shown in Fig. 5.6. The projection has two different oscillating
frequencies: the high-frequency oscillating component is caused by the 800 nm field and the slow
oscillating component is caused by the 400 nm laser field and corresponds to the Rabi frequency.
Because of the low intensity, the Rabi frequency is very small in the present case. Therefore, theelectron moves very slowly from one side to the other side due to the low intensity of the 400 nm
laser pulse, which explains why the electron can locate more frequently on one side for such a long­®

3σg |2σu |2

1.0
0.8
0.6
0.4
0.2
0.0
0 10 20 30 40 50 60 70 80 90
Time [fs]104

Figure 5.6: The projection from HOMO-2 to HOMO of N+
2 for two-color laser field: 400 nm laser
field at I = 1 × 1012 W/cm2 and 800 nm laser field at I = 1 × 1014 W/cm2 .

time period.
Control of electron localization in H+
2 at large internuclear distance has been reported previous
in [204, 205]. The coupling of the nearly degenerated 1σg and 1σu states can lead to permanent
trapping of the electron on one side. The control of such permanent trapping is very appealing since
it can influence strongly related processes such as ionization, dissociation as well as other chemical
reactions. Although our simulations with bichromatic field are still far from showing permanent
trapping of the electron on one side, it however helps to understand better the mechanism behind
it and opens up the possibility to study the conditions of permanent trapping in multielectron
molecules.

5.4

Electron localization in other molecules
According to the discussion in previous sections, the electron dynamics can be strongly in-fluenced by the coupling of the orbitals. For N+
2 , the coupling can be induced between for σg - andσu -type orbitals. In Chapter 3 and Chapter 4, we have shown fractional harmonics and enhanced
ionization for other molecules as well. For some of these molecules the efffect are due to the couplingof orbitals with other symmetries.
For C2 H+
2 , a 760 nm laser field is used to couple 2σu and 3σg , which has the same typeof orbital coupling as N+
2 . The density plot of 2σu is shown in Fig. 5.7(a), we see a strong

105+
nonadiabatic dynamics of electron localization but is different compared to N+
2 . For C2 H4 , one cancouple orbitals 1b2u and 1b3g at 400 nm. The orbital density for 1b2u is shown in Fig. 5.7(b) for thelaser intensity of 1×1014 W/cm2 . The electron response is almost adiabatic with small nonadiabatic
signatures. It is surprising when compared with N+
2 at the same resonant wavelength and intensity,where nonadiabatic behavior is already clearly visible. This indicates that the dynamic localization
can also depend on the orbital type. For CO+
2 , where 1πu and 1πg orbitals are coupled for alaser field of wavelength of 350 nm and intensity of 1 × 1014 W/cm2 , the electron density (in Fig.
5.7(c)) almost follows the laser field adiabatically. We also studied the density plot of the coupledorbital 4b2 and 6a1 of NO2 , shown in Fig. 5.7(d). The non-adiabatic effect is weak compared toN+
2 at the same laser intensity. Our simulations are currently limited to smaller molecules dueto the computational cost. However, there are no physical arguments which would prohibit the
dynamical localization for any open-shell molecule. We have chosen those molecules which due totheir electronic properties are more easily accessible for experimental studies. As we pointed out,
the effects can strongly depend on the laser intensity and the type of coupled orbital. The choiceof systems was also related to visualization of the dynamics, which is not easy for nonsymmetricmolecules or σ-π coupling of the orbitals.

5.5

Summary
In this Chapter, we systematically investigated the nonadiabatic electron localization due tothe coupling of 2σu and 3σg in N+
2 . The nonadiabatic response of electron is determined by twocompeting processes: coupling induced oscillations with Rabi frequency and laser induced oscillations with laser frequency. Such competing processes can strongly depend on the laser intensity
and thus by changing intensity one can change from adiabatic to nonadiabatic regime for electron
dynamics. Moreover, we have explored the possibility of controlling such nonadiabatic behavior
with two-color laser pulse and have shown that in such simulations the electron can reside on one
nuclei more frequently for relatively long time. The electron localization is also examined for othertype of couplings in the regions where we have found strong fractional harmonics and enhanced

106(a) Orbital density of 2σu of C2 H+
2 . Laser field:2
λ = 760 nm and I = 1 × 1014 W/cm(b) Orbital density of 1b2u of C2 H+
4 . laser field: λ =2
400 nm and I = 1 × 1014 W/cm(c) Orbital density of 2σu of CO+
2 . Laser field: λ =2
350 nm and I = 1 × 1014 W/cm(d) Orbital density of 4b2 of NO2 . Laser field: λ =2
400 nm and I = 1 × 1014 W/cm

Figure 5.7: Orbital density plot of the coupled orbitals for different molecules.

107
ionization. For several open-shell molecules nonadiabatic response has been found and illustrated.

Chapter 6

Complex Gaussian basis for non-Hermitian Floquet formalism

A lot of important results in strong field physics of atoms were obtained with theories basedon the Floquet theorem. This type of approach was initially developed to obtain exact solutions of
the time-dependent Schrödinger equation which describes the interaction between a quantum system and a continuous intense light source which is strictly periodic in time. The wavefuction of thesystem describing the photo-induced effects is represented in an extended Hilbert space spannedby photon dressed states, which are constructed as a direct product of system and field states.
This approach allows to transform the original time-dependent problem into a time-independent
formulation yielding quasi-stationary eigenstates which can be expressed in harmonics of the carrier
frequency of the field. Instead of dealing with time propagation one has the advantage of applicationof fast diagonalization algorithms. The Floquet ansatz was successfully applied to non-perturbative
studies of the bound-bound transitions in atoms and molecules as well as laser-assisted collisions[206, 183, 207, 208]. Moreover, the non-Hermitian version of the Floquet method provides a description of the decaying states via complex scaling [209]. It can therefore be used in the context
of resonance as well as ionization, for example it has been applied to study multiphoton ionizationof neutral atoms [210] and molecules (e.g. in [211, 212]), multiphoton detachment of negatively
charged ions, ac-Stark shifts, laser-induced chemical bond softening and hardening in a molecule[213], as well as multiphoton dissociation [183, 214], Coulomb explosion, and high harmonic generation [215, 216]. Floquet-based methods not only give results in agreement with experiments, they
allow, due to the inherent, intuitive picture of laser-dressed states, for straightforward interpreta-109tion of the results, e.g. determination of the states responsible for transient resonances in highharmonic spectrum [207] as well as identification of key mechanisms. In this way it can help in thedevelopment of simplified models.
Formally, Floquet method is limited to time-periodic problems such as interaction of systems
with cw-lasers. In recent years, this limitation of the Floquet theory has been overcome and Floquet
theory has been adapted to the treatment of laser fields whose amplitudes are varying in time. For
example, reasonable results were obtained already for relatively short laser pulses composed of 5
optical cycles. Moreover, it has been shown that with slight modifications one can reformulate theproblem to describe the interaction with ultrashort laser pulses [217].
The way to make such approach applicable to larger polyatomic molecules is to combine
the time-dependent density-functional theory (TDDFT) with the Floquet method [206, 183, 218]
The first applications of this method both in linear response regime [219] and in nonperturbative
regime [59, 60, 62, 63, 64] are promising. However, in the nonperturbative regime only grid-based
calculations were performed and it is very interesting to explore possibilities of the Floquet-based
density-functional theory (DFT) method using the basis-set approach. The Floquet-dressed statesmight give us better physical picture to understand the multiorbital contributions from inner valenceorbitals to the high harmonic generation as well as ionization.
In this thesis, we apply Floquet theorem to solve the time-dependent Schrödinger equation(TDSE) with Gaussian basis sets. Complex basis set is used to enable the description of boundfree and free-free transitions. We designed an optimization process which yield a very accurate
calculation of ionization rate. First applications regarding one- and two-photon ionization will bepresented.

6.1

Floquet theorem
A system subject to a time-dependent external potential, is described by the time-dependent

Schrödinger equationi~∂Ψ(r, t) = Ĥ(r, t)Ψ(r, t),
∂t

(6.1)

110
where Ĥ = T̂ + Û + V̂ . According to the Floquet theorem [220], if the external potential V isperiodic in time V (t + T ) = V (t), then the (ground state) solution of Schrödinger equation (6.1)
can be written as superposition of so-called quasienergy states given byΨλ (r, t) = e−iελ t Φλ (r, t),(6.2)

Φλ (r, t + T ) = Φλ (r, t),

(6.3)

where

and ελ is called quasienergy. Substituting Ψλ (r, t) into Eq. (6.1), we obtain an eigenvalue equation
for the quasienergy:Hˆ (r, t)Φλ (r, t) = λ Φλ (r, t),

(6.4)

∂
Hˆ (r, t) = Ĥ(r, t) − i .
∂t

(6.5)

where

Because of the periodicity of Φλ (6.3), we have the following transformation,0ελ = ελ + nω0

φλ (r, t) = exp(inωt)φλ (r, t),

(6.6)
(6.7)

which means that Floquet states are equivalent, if the quasienergies only differ by nω where n is anarbitrary integer number. Such property is the same as for Bloch eigenstates in solid state physics.
Since Φλ (r, t + T ) = Φλ (r, t) is periodic in time, its temporal part can be expanded into
Fourier seriesΦ(r, t) =∞X

e−inωt fn (r).

(6.8)

n=−∞

The integer n represents the number of photons absorbed (for positive value) or emitted (fornegative value). The space part fn (r) can be either represented by the value on a numerical gridor spanned by basis sets, such as Sturmian functions [221], Slater functions [211], and Gaussianfunctions [211, 222]. Considering the great success of real Gaussian basis sets for representation of

111bound states, especially for large molecular systems, we have chosen Gaussian basis sets to span
the spatial part of the wavefunction. If we denote the basis set functions as χj (r), we can write
Φλ (r, t) =

X

λ
Cnj
exp(−inωt)χj (r).

(6.9)

nj

Substituting it into Eq. (6.4), we obtain the following matrix equation
Xh

i
X
λ
λ
hχi |Ĥ (m−n) |χj i − mωδnm Sij Cnj= ελ
Sij Cmj,nj(6.10)jwhere
1Ĥ (r) =Tˆ TnĤ(r, t) exp(inωt)dt,(6.11)oand Sij = hχi |χj i is the overlap matrix. The integral hχi |Ĥ (m−n) |χj i can be decomposed ashχi |Ĥ (m−n) |χj i = hχi | −

X
∇2
−Zk
1
|χj iδnm +hχi ||χj iδnm + E · hχi |r|χj i(δn,m+1 + δn,m−1 )2m|r − Rk |2
k(6.12)
for one electron molecular systems where Rk denotes the nuclei’s positions and E is the externallaser electric field. The three terms on the right hand side are called kinetic-energy integral,
nuclear-attraction integral and dipole integral, respectively. We will show how to calculate themfor Gaussian basis sets with recursive relationship in section 6.3. When applied to the timedependent Kohn-Sham equation (2.33) to solve for the case of a multielectron system, the matrix
elements include two additional parts: electron-electron repulsion integral hχi χn |νH |χj χm i andexchange-correlation integral hχi |νxc |χj i. Eq. (6.10) can be efficiently solved using PETSc or
SLEPc numerical libraries [223, 224] which are based on methods to solve eigenvalue problems forlarge matrices.

6.2

Complex basis sets
Physical phenomena can be either studied via Hermitian formalism or non-Hermitian for-malism. Although Hermitian formalism of quantum mechanics can handle most of problems inphysics, non-Hermitian formalism can provide a straightforward explanation to some phenomena

112
which might be very hard or impossible to be explained within the framework of Hermitian formalism such as exceptional points, where two or more degenerate resonance states have coalesced,multiphoton ionization of atoms or multiphoton dissociation of molecules, where bound-free or
free-free transitions are present.
There are two popular ways to obtain a non-Hermitian Hamiltonian. One is to apply complex
coordinate rotation [225, 226], r → r exp(iα) (complex rotation is usually applied to radius r in
spherical coordinate), to the Schrödinger equation (6.1)
i∂Ψ(r exp(iα), t)= Ĥ(r exp(iα), t)Ψ(r exp(iα), t),∂t(6.13)

where the Hamiltonian operator Ĥ(r exp(iα), t) becomes non-Hermitian. We can again, following
Floquet theorem, write the wavefunction Ψ(r exp(iα), t) asΨλ (r exp(iα), t) = exp(−iελ t)Φλ (r exp(iα), t).

(6.14)

Substituting it back into the Schrödinger equation, we getHˆ (r exp(iα), t)Φλ (r exp(iα), t) = ελ Φλ (r exp(iα), t).

(6.15)

ελ is now a complex quasienergy and the imaginary part of quasienergy Im(ε) is related to theionization rate Γ of that state via Γ = −2Im(ελ ).
The influence of the transformation r → r exp(iα) on the spectrum of a many-body Hamiltonian [227] is that both bound state eigenvalues and the complex eigenvalue associated with resonant
states are independent of the rotation angle α while the continuum states are rotated by angle of
2α into the lower half complex energy plane, as shown in Fig. 6.1. If α is large enough, then thecomplex resonance eigenstates will be exposed and thus we can study the resonance states throughthe complex coordinate rotation.
The other way to construct non-Hermitian Hamiltonian is to use complex basis set which isequivalent to complex coordinate rotation [228] according to the following transformation
˚
f (r, θ, φ)H(reiα , θ, φ)g(r, θ, φ)r2 sin(θ)drdθdφ
˚
−3iα
=e
f (re−iα , θ, φ)H(r, θ, φ)g(re−iα , θ, φ)r2 sin(θ)drdθdφ,(6.16)

113

Figure 6.1: Effect of complex coordinate rotation. Bound states and thresholds are invariant. As
continuum states rotate, complex resonance eigenvalues may be exposed. Adapted from [227].

where f (r, θ, φ) and g(r, θ, φ) are real L2 basis functions. In the above equation, no complex
conjugate is applied to left basis set f (r, θ, φ) compared to Hermitian case, which is called cproduct [229] and denoted as (f |H|g). The complex basis is more appealing compared to thecomplex-coordinate rotation technique. The rotation of coordinate to complex plane is hard tointerpret in the case of molecules while there is no such interpretation problem for the complex
basis. Since only continuum is necessary to be rotated into the lower half complex energy plane,
only representation of the continuum by the complex basis is necessary and the basis for the boundstates can still be kept real. So previous results of Gaussian basis for the description of boundstates can be used directly without any modification. On the other hand, different complex basis
can have different rotation angle and thus provides even some flexibility.

6.3

Gaussian basis integral

6.3.1

Fundamental integral

As mentioned above, the spatial basis set functions χj (r) can be chosen as Sturmian functions,
Slater functions or Gaussian functions. The matrix element evaluation for Sturmian and Slater
functions is computationally demanding and thus limits their extension to larger molecular systems.
Gaussian basis sets on the other hand have several properties, as illustrated in this section, whichmake it possible to extend calculations to large systems.

114
The primitive Gaussian function we consider has the form2l

φa (r) = N xlax yay zalz e−αk |ra | ,(6.17)

where ra = r−A, which is centered at position A with complex exponential αk and the momentum
vector (lx , ly , lz ). We use subscript a here to denote the parameters (position A, momentum
vector, and exponential αk ) for a Gaussian basis. In the following, we will just use the subscript a
to represent a Gaussian basis. The normalization coefficient of Gaussian basis is(2α)(lx +ly +lz )/2+3/4
N=q.
Γ(lx + 12 )Γ(ly + 12 )Γ(lz + 12 )(6.18)

From the fact that the Hamiltonian contains only one- and two-electron operators, it is straightforward to show [230] that most of the matrix elements can be written in terms of integrals of thegeneral form¨
φa (r1 )φb (r1 )f (|r1 − r2 |)φc (r2 )φd (r2 )dr1 dr2 .

(ab|cd) =(6.19)For the Gaussian basis, we first consider a more fundamental integral (integral for s-state) [231]
¨
I=

2

2

2

2

e−α|r1 −A| e−β|r1 −B| f (|r1 − r2 |)e−γ|r2 −C| e−δ|r2 −D| dr1 dr2 ,(6.20)

where we ignored the normalization coefficient for convenience. All the matrix integrals in Eq.
(6.12) can be written in this general form. For example, the overlap integral,
ˆ
S=I

overlap=22e−α|r−A| e−γ|r−C| dr,(6.21)

can be obtained by setting
β = δ=0
f (|r1 − r2 |) = δ(|r1 − r2 |),(6.22)

in Eq. (6.20). And the nuclear-attraction integral,
ˆ
INE =

2

e−α|r−A| e−β|r−B|

2

−Z
dr,
|r − C|(6.23)

115
can be obtained from Eq. (6.20) using
δ=0γ=∞
γ 3/2 −γ|r−C|2e= δ(r − C)γ→∞ π 3/2
Z
f (|r1 − r2 |) = −.
|r − C|lim(6.24)For the electron-electron repulsion integral, f (|r1 − r2 |), we set
f (|r1 − r2 |) =1.
|r1 − r2 |(6.25)The dipole integral can be transformed into the overlap matrix with higher angular momentum
Gaussian basis sets, as we will show later.
To evaluate integral Eq. (6.20), we first notice that the product of two Gaussian functions isstill a Gaussian function2e−α|r1 −A| e−β|r1 −B|

2

2= GAB e−ξ|r1 −P|

2

2= GCD e−η|r2 −Q| ,e−γ|r2 −C| e−δ|r2 −D|

2(6.26)

where
GAB = eαβ
− α+β|A−B|2,ξ = α + β,P =αA + βB,
α+βGCD = eγδ− γ+δ|C−D|2,

η = γ + δ,Q =γC + δD.
γ+δ

(6.27)

Because of this property, the four center problem (A, B, C, D) in integral (6.20) can be simplified
to a two-center problem (P, Q)¨
I = GAB GCD

2

2

e−ξ|r1 −P| f (|r1 − r2 |)e−η|r2 −Q| dr1 dr2 .

(6.28)

116
Replacing the three factors of the integrand in Eq. (6.28),
ˆ  3/2π2−ξ|r1 −P|2−3e
= (2π)
e−k1 /4ξ eik1 ·(r1 −P) dk1 ,ξˆf (|r1 − r2 |) = (2π)−3 f˜(k2 )eik2 · (r1 − r2 )dk2 ,(6.29)ˆ  3/2π2e−k3 /4ξ eik3 ·(r1 −Q) dk3 ,η(6.31)e−η|r2 −Q|2−3

= (2π)(6.30)

substituting (6.29)–(6.31) into (6.28) and integrating over r1 , r2 , k1 , and k3 , Eq. (6.20) becomes a
simple integralˆ
k2 (ξ+η)GAB GCD− 2 4ξηI=ee−ik2 ·R f˜(k2 )dk2
3/2
8(ξη)
ˆ
πGAB GCD ∞2=u sin ue−u /4T f˜(u/R)du,2(ξη)3/2 R3 0

(6.32)

whereT= θ 2 R2

θ2 = ξη/(ξ + η)
R = Q − P.

(6.33)For the overlap integral, we obtainI

overlap


=

π
α+γ

3/2

2

e−αγ|A−C| /(α+γ) .

The kinetic-energy integral is connected to the overlap integral byˆ
1 2 −γ|r−C|2
kinetic−α|r−A|2I= e− ∇r edr2ˆ122= − ∇2C e−α|r−A| e−γ|r−C| dr
2
1= − ∇2C I overlap
2

3/2
2
3αγ
π
=
e−αγ|A−C| /(α+γ) .
α+γ α+γ(6.34)(6.35)

1
With f (|r1 − r2 |) = |r1 −r, the result for the electron-electron repulsion integral is2|

I EE = GAB GCD
ˆ 1
F0EE (T ) =02

2π 5/2F EE (T )ξη(ξ + η)1/2 0

e−T u du.

(6.36)(6.37)

117
This one dimensional integral (6.37) can be evaluated very efficiently using the GSL integration
library [232]. For the nuclear attraction integral in Eq. (6.24), we immediately obtain
2π
I N E = GAB 3/2 F0EE (T ).
ξ
6.3.2(6.38)

Recursion relation for general primitive Gaussian basis integrals

All the integrals presented so far are for primitive Gaussian basis functions with zero angularmomentum. For Gaussian basis functions with higher angular momentum, their integrals can beobtained from those for the primitive s-Gaussian basis function using recursion relationship. Before
we the introduce the actual recursion relationship, it is useful to look at a simplified notation [231].
In Eq. (6.19), a general two-electron integral over the two electron operator f (|r1 − r2 |) is given.
Such an equation defines an inner product between two functions(ab| = φa (r1 )φb (r1 )(6.39)|cd) = φc (r2 )φd (r2 ).

(6.40)

using the Dirac notation, where (ab| and |cd) are referred as ”bra” and ”ket”, respectively [231],and (ab|cd) is referred as ”braket”. It is helpful to generalize Eq. (6.39) and Eq. (6.40) to derive
the recursion relation. Though the generalized form is more complicated, it should be kept in mind
that bras and kets here are simply functions of the positions of electron 1 and 2, respectively. The
generalized primitive bra is given by2

2(abp| = DA DB e−α(r−A) e−β(r−B) Πi=x,y,z (i − Ai )ai (i − Bi )bi ξ pi /2 Hpi [ξ 1/2 (i − Pi )],(6.41)

where Hn is the nth Hermite polynomial defined as
H0 (x) = 1
H1 (x) = 2x
Hn+1 (x) = 2xHn (x) − 2nHn−1 (x).

(6.42)The three symbols abp in the bra represent the primitive Gaussian basis χa (r) centered at position
A, Gaussian basis χb (r) centered at position B, and the Hermite function centered at position P.

118
Similar general form can be given for ket. Then, the general primitive integral we are concernedwith can be written as(ab|cd) = (ab0|cd0)¨
22=e−α(r1 −A) e−β(r1 −B) Πi=x1 ,y1 ,z1 (i − Ai )ai (i − Bi )bi ·22

f (|r1 − r2 |) · e−α(r2 −C) e−β(r2 −D) Πi=x2 ,y2 ,z2 (i − Ci )ci (i − Di )di dr1 dr2 ,(6.43)

where we dropped the coefficients DA , DB , DC , and DD , which can be easily added back later.
Applying the properties of Hermite polynomials in Eq. (6.42), there exists a three-term relationship((a + 1i )bp| = pi (ab(p − 1i )| + (Pi − Ai ) (abp| + (2ξ)−1 (ab(p + 1i )| ,(6.44)

where 1i is a unit index vector with a non-zero element 1 in direction i (i can be x, y, or z). Thus,((a + 1i )bp| means that the angular momentum in Eq. (6.41) at direction i becomes ai + 1, whileothers remain unchanged. Similar three-term relationships also exist for the ket. Such three-term
relationship for a, b, c, and d can be applied recursively until all of the angular moments become
zero. Then, the integral (6.43) is finally converted into the evaluation of a serials of¨

22e−α(r1 −A) e−α(r1 −B) Πi=x1 ,y1 ,z1 ξ pi /2 Hpi [ξ 1/2 (i − Pi )]·(00p|00q) =2

2

f (|r1 − r2 |) · e−α(r1 −C) e−α(r1 −D) Πi=x2 ,y2 ,z2 ξ qi /2 Hqi [ξ 1/2 (i − Qi )]dr1 dr2¨
2=GAB GCDe−ξ|r1 −P| Πi=x1 ,y1 ,z1 ξ pi /2 Hpi [ξ 1/2 (i − Pi )]f (|r1 − r2 |)2e−η|r2 −Q| Πi=x2 ,y2 ,z2 ξ qi /2 Hqi [ξ 1/2 (i − Qi )]dr1 dr2 .

(6.45)

Notice that Hermite polynomials satisfy2e−ξ(r−p) Πi=x,y,z ξ ri /2 Hri (ξ 1/2 (i − Pi )) = (−1)r Πi=x,y,z= Πi=x,y,z

dri −ξ(i−Pi )2e
dri i

dri −ξ(i−Pi )2
e
.
dPiri

(6.46)

so Eq. (6.45) can be simplified as(00p|00q) = Πi=x,y,z

dpi dqiGAB GCD
dpi Pi dqi Qi

¨

2

2

e−ξ|r1 −P| f (|r1 − r2 |)e−η|r2 −Q| dr1 dr2 ,(6.47)

119
where the integral part is exactly the fundamental integral we considered before. So, we canimmediately get the final form of the integral as
dpi dqi GAB GCD(00p|00q) =Πi=x,y,z p
d i Pi dqi Qi 8(ξη)3/2
dsi GAB GCD=(−1) Πi=x,y,z s
d i Ri 8(ξη)3/2qˆ

k2 (ξ+η)e− 2 4ξηˆ

e−ik2 ·R f (k2 )dk2

k2 (ξ+η)e− 2 4ξηe−ik2 ·R f (k2 )dk2 ,(6.48)(6.49)

where we combined the two derivatives using the fact that R = P − Q and si = pi + qi . To simplify
the notation, we denote the last equation (6.49) as (S). For the overlap integral, we haveoverlap

(S)


π 3/2
2
Hrx (θRx )Hry (θRy )Hrz (θRz )e−θR θr .
=ξ+η


(6.50)

To avoid the derivative evaluation in Eq. (6.49), we further define(S)overlap
=
m




π 3/2 r 2 m2
θ (2θ ) Hrx (θRx )Hry (θRy )Hrz (θRz )e−θR .
ξ+η

(6.51)

asThen we can obtain a two-term recursion relationship for (S)overlapm(S)overlap= Ri (S − 1i )overlap− (si − 1)(S − 2i )overlap
m
m+1
m+1 .

(6.52)

Similar two-term recursion relationships can be found for the other type of integrals as well. So
far, we can summarize the algorithm using the calculation of Gaussian basis set as(ab|cd) = (ab0|cd0) ⇒ (00p|00q) ⇒ (S) ⇒ (0)m .

(6.53)

Instead of implementing the recursion function directly to evaluate each recursion relationship, it
is implemented with a non-recursive function to avoid extra time caused during calling recursivefunction and stack overflow. We can further summarize the strategy to calculate the Gaussian basis
integrals as follows:(1) Determine the largest m according to the momentum index of integral (ab|cd).
(2) Form a table with elements [0]i for i=0, ..., m.
(3) Calculate (S) via recursion relationship (6.52).

120(4) Form table (00p|00q) = (−1)q (S) for p=0,1,...p, and q=0,1,...,q.
(5) Calculate (ab|cd) via recursion relationship (6.44).
The evaluation of Gaussian basis set integral in this form is very efficient and is implemented as aclass GaussianBasisEval in our program (see Fig. 6.2).

6.4

Technical details
The main program structure is shown in Fig. 6.2. The IO class is a class that can help

Figure 6.2: Diagram of the Floquet program representing the relationship between classes.

other classes to parse the input and output without worrying about race condition in parallel
computation. The GaussianBasisEval class is used to calculate the Gaussian integral in Eq. (6.20)
while the ContractGaussianBasis class which is built on the top of GaussianBasisEval is used to121calculate the matrix element with contracted Gaussian basis sets. The main task of Hamiltonianclass is to calculate the Floquet matrix elements in Eq. (6.4). Parallelization with MPI is applied
to speedup the matrix element calculation process. After the matrix is created, EigenSolver classwill check the matrix property and choose an appropriate eigensolver to solve the matrix withparallel package SLEPc [224]. Everything except the optimization class is written in C++. The
optimization module is implemented with Python since many optimization packages are available
for Python. The optimization process will generate an input file and read the output from theprogram step by step to find the optimal Gaussian basis sets.
Based on our test calculations, there are several things to be taken into consideration:(1) Because of the periodic property of Floquet states in Eq. (6.3), we can no longer selectFloquet states based on the order of their energies. One way to avoid such problem is toconfine the solution ε ∈ (−ω/2, ω/2] but it still does not clarify from which states each
Floquet state is formed, which is an obstacle for the automatic optimization process. One
possible solution is to first calculate all bound states of interest without external laser field
in which case there is no periodic quasienergy problem and then project Floquet states
Φλ (r, t) to the bound state ψb (r) of interested via Eq. (6.54) and select the Floquet statewith highest projection value
ˆ T

ˆ
dt

d3 r|ψb (r)Φf (t, r)|2

(6.54)o(2) We also find that the imaginary part of quasienergy without external field using rotatedcomplex Gaussian basis sets has a non-zero value. Such phenomena arise from the incompleteness of basis sets and its effect can be reduced by increasing the basis sets as shownin [221, 233]. However, for a Gaussian basis set, it is not easy to improve the accuracyof bound states and the best ground state energy for hydrogen atom we got is -0.4998.
Because of this problem, we subtract the imaginary part of eigenvalue ε in laser field by theimaginary part of eigenvalue εo with zero-external field. Although the results for the cases

122in this thesis are accurate, the procedure should be dealt with carefully for other situations,especially for processes at high intensity.

6.5

Optimization for complex Gaussian basis
The complex coordinate rotation is quite straight forward in numerical calculations. Chu andTelnov have applied it to study multiphoton ionization, dissociation, etc. extensively for one- andtwo-color laser field [183]. However, so far it is only applied for problems in 1 or 2 dimensions. For
the calculation of large molecules with grid method in 3D, it is even impossible with current supercomputer because the matrix to solve is very huge. For example, for simulation box with 200 grid
points in each dimension and a 10 photon block, the matrix to be solved is 8 × 107 by 8 × 107 , whichis even impossible to store directly in the super computer. Especially, the computational costs will
be huge during self-consistent calculation after applying it to DFT. With complex basis set, it ispossible to get high precision results with relatively small basis. For example, a basis set with 21
Gaussian functions is applied to calculate the photoionization cross sections from N2 [222] and it
is even possible to obtain accurate results for photoionization cross sections with only one or twocomplex Gaussian function as shown in [234].
With complex coordinate rotation technique, eigenvalues of resonances can be found since
they do not change if the rotation angle reaches a critical value. The complex Gaussian basis is
usually obtained by rotating their exponent αe−iθ directly. Since complex basis sets are obtainedby rotating real Gaussian basis sets, according to the equivalence between complex coordinate
rotation and complex basis sets in Eq. (6.16), the eigenvalue of resonance states can be stationary
after some critical rotation angle. For a simple system, with large enough basis, it is possible by
optimizing θ to find the stationary angle, at which the quasienergy does not change with rotation
angle anymore. However, for complex systems, usually both α and θ have to be optimized. An
efficient optimization process for the exponent α is very important for the application of suchtechnique. For the optimization object function, the natural choice is the quasienergy, which has
been used in [235, 236, 237, 238]; while optimization on polarizability is applied in [222, 239, 234].

123
As for the optimization for the Gaussian exponent α, there is still no common choice yet, which
therefore is the most difficult part. A geometrical series for Gaussian exponent α was chosen in[238, 233]. Sometimes, the value of α is chosen such that a complex eigenvalue with real partsin the expected energy range is guaranteed [222]. And recently, Matsuzaki et al. proposed to fit
Gaussian functions to Slater function, where the exponent of the Slater function is determined by
the energy of ionized electron [240].
During past decades, various real Gaussian basis sets are proposed to describe the boundstates and yield relatively high precision for stationary calculations. However, there is still nosystematic efficient way to find precise Gaussian basis sets describing continuum states in a laser
field. This is because firstly, the basis set parameters strongly depend on the laser parameters,
especially the laser frequency. Secondly, unlike in the stationary case, in which basis parameters
for bound states can be optimized via minimization of bound state energy, in the laser field,∂εFloquet states only satisfy a stationary condition which means ∂α= 0 where ε is quasienergy andiαi are parameters to be optimized. This makes the optimization process more difficult than justminimizing the state energy. Thirdly, the optimization problem is a non-convex problem whichmakes the process rather difficult.

6.5.1

Optimization algorithm

In our calculations, we take real contracted Gaussian basis sets to describe the bound states
and add several additional complex basis sets to describe the continuum wavepackets. Since theGaussian basis sets are not complete, the stationary property usually only holds for a small rangeof angles, typically around 30o . So, one way to optimize the exponents of Gaussian basis set is
to calculate the eigenvalue for each angle from 0o to 80o degree and study the trajectory to see if
a stationary point is found. If no stationary point is present, then one changes the exponents ofGaussian basis set and repeats the above processes until a stationary point is found.
Such optimization is not very efficient, because we have to compute the eigenvalues at eachangle from 0o to 80o to see if the stationary property is present with the trial exponents of Gaussian

124basis set. In such process it is difficult to implement an auto-optimization process since we cannot
easily formulate a objective function for such method. Alternatively, we can assume that thestationary points always occurs around 30o which is a good assumption because the critical angle
reported in most papers is lower than 30o . In this way we only need to perform calculations around
30o without a full calculation.
To develop an objective function for optimization, we require that the first derivative of
the eigenvalue with respect to rotation angle is zero at 30o , because the eigenvalue at stationary
point (assumed to be 30o ) is independent of the rotation angle. According to our calculations, thereal part of eigenvalue strongly depends on the real Gaussian basis sets which are already known
and barely change for the complex Gaussian basis sets. Therefore, we only require the imaginary
eigenvalue Im(ε) to be stationary.
The stationary property is not guaranteed by the condition for the first derivative only. As
more higher derivatives are zero as well, as better the conclusion that the stationary conditionis achieved. To save computation time, we choose to optimize the first and second derivatives ofimaginary eigenvalue with respect to the rotation angle only. In our program, the first derivative
is calculated with the spline interpolation function. During the interpolation process, we set theboundary condition such that the second derivative at the ends of interpolated points are zero. We
have chosen our objective function:
"

 2
2 #
∂Im(ε) 2∂ Im(ε)
+C
∂θ
∂θ2

,

(6.55)

θ=30owhere we added C to adjust the relative weight between the first and second derivatives. Both firstand second derivative are squared to convert the stationary problem to a minimization problem.
The optimization process will stop when the value of Eq. (6.55) is less than a critical value.
However, the ionization rate γ = −2Im(ε) is strongly influenced by the laser intensity and can varyby several orders of magnitude. Thus, the optimized value of Eq. (6.55) can also vary dramaticallyand the critical value to stop optimization process need to be changed accordingly for different

125systems and laser intensities. Instead of doing this, it is better to optimize the relative value:
"

 2
2 #
1
∂Im(ε) 2∂ Im(ε)+C,(6.56)∂θ∂θ2(Im(ε))2oθ=30which does not depend strongly on the laser intensity and we do not have to change the criterion
for optimization for different intensities.
The optimization algorithm applied is a sequential penalty derivative-free method for nonlinear constrained optimization [241], which is provided in Python package pyOpt [242]. The
optimization problem is a highly non-convex problem and thus the final optimal result depends on
the initial values of exponents of Gaussian basis sets. To overcome such an initial value problem
we followed the idea of R. Matsuzaki et al. [240], who used Slater basis sets in their calculation.
The radial part of the Slater function has the form as
R(r) = N rn−1 e−sr .

(6.57)

With complex exponent s, the behavior of the Slater function is more similar to a plane wave eikr
than a Gaussian function because the oscillating frequency of Slater function only depends on s
while for Gaussian function the oscillating frequency increases with r. Since the ionized wavepacket
can be approximated by a plane wave eikr , we can choose Slater function to represent such a plane√
wave by setting s = k = 2E. We also want to find a similar connection between the energyof the ionized wavepacket and exponents of Gaussian basis set. Actually, there exists a integral
transformation between Slater function and Gaussian functionse−sr = √
2 π

ˆ ∞

2

2

α−3/2 e−s /4α e−αr dα,(6.58)

0which can be interpreted such that a Slater function e−sr can be written as a sum of Gaussian2functions e−αr with coefficient C = √s
2

(π)

2

α3/2 e−s /4α . The coefficient C constraints the range of

exponents α, and thus it tells us for certain value s, where the value of the exponents {αi } of theGaussian basis sets should be. Therefore, instead of fitting Slater function with Gaussian basis set
as in [240], we assign the initial value of {αi } according to C. To prove that this strategy is useful,
we have done several calculations and the optimized value of {αi } is shown as vertical line in Fig.

126
6.3. The solid curve line represents the normalized coefficient C for different ionization energy Ek
while the other two dashed curves are those for electron energy of 2Ek and Ek /2, respectively.
Calculations are performed for H atom in laser field of frequency ω=0.6 (left), 1.0 (middle), and
2.0 (right). From the calculation, the exponents {αi } of Gaussian basis sets do fall in the range
1.0
0.8
0.6
0.4
0.2
0.010-3 10-2 10-1 100 101
alpha

1.0
0.8
0.6
0.4
0.2
0.0 10-2

10-1

100
alpha

101

102

1.0
0.8
0.6
0.4
0.2
0.0
10-2

10-1

100
alpha

101

Figure 6.3: Distribution of exponents α of Gaussian basis sets (red vertical line) after optimization
for the H atom in a laser field with ω = 0.6 (left), ω = 1.0 √(middle), ω = 2.0 (right). The solid2
blue curve is the normalized coefficient α−3/2 e−s /4α for s = 2Ek and the other two dashed blue
curves are for electron energy of 2Ek and Ek .

defined by the coefficient C for several very different laser frequency. We, therefore, choose values
uniformly distributed in the range defined by C as the initial optimization value.

6.6

Results and discussion

6.6.1

One photon ionization of H atom

We first calculated the one-photon ionization of the hydrogen atom at intensity 4.387 ×
1013 W/cm2 with laser frequency ω = 0.6 a.u. The real Gaussian basis set used is 6-31++G∗∗ ,which represents the ground s-state well. For one-photon ionization, no higher bound states are
required to get accurate results. Complex Gaussian basis sets including both s and p states areadded to represent the continuum states. P states are necessary since the final angular momentumis 1 after one photon transition from s state. The complex s state is also necessary to get good
stationary result based on our calculations. s and p Gaussian basis sets share the same exponents
α. For the test, we first tried complex Gaussian basis sets with 3 and 5 states. The initial value isassigned according to the coefficient C with Ek = 0.1 (energy of ionized electron is Ek = ω − Ip ).
The number of Floquet block n is also tested for convergence. For this one-photon transition and127
laser intensity 4.387 × 1013 W/cm2 , n=0, 1 is sufficient to achieve a converged result. For higher
intensity, more Floquet blocks are needed as shown later.
After optimization, we calculated the eigenvalues with the optimized exponents {αi } for therotation angle from 0o to 80o . The results are shown in Fig. 6.4, where the top two panels are
0.0000

Eimag [a.u.]−0.0002−0.0004−0.0006−0.0008−0.0010−0.00120

10 20 30 40 50 60 70 80
Angle [degree]

0.0000

Eimag [a.u.]−0.0002−0.0004−0.0006−0.0008−0.0010−0.00120

10 20 30 40 50 60 70 80
Angle [degree]

Figure 6.4: Quasienergy trajectory of H atom in a laser field at intensity of 4.387 × 1013 W/cm2
with frequency ω = 0.6. Results in the top two panels are for 3-complex Gaussian basis sets while
those in the bottom row are for 5-Gaussian basis sets. Figures on the left show the imaginary partof quasienergy as a function of the rotation angle θ and panels on the right show the trajectory ofthe quasienergy. A stationary point is found with ε = −0.498 − 0.00109i a.u..

results obtained with 3 complex Gaussian basis sets and the bottom two panels are results obtainedwith 5 complex Gaussian basis sets. The Figures on the left show the imaginary part of quasienergy
as function of rotation angle. In both cases, we get a relatively flat curve around 30o which means
both first and second derivatives are close to zero as expected. On the other hand, although we only
minimize the objective function at 30o , the imaginary part of quasienergy is actually unchanged

128
ω = 0.6,F = 1 × 10−4
ω = 0.6, F = 7.08 × 10−4
ω = 0.6, F = 0.05
ω = 1.0, F = 7.08 × 10−4
ω = 2.0, F = 1 × 10−4

our results
0.13540.13560.13990.03380.0046

reference values
0.1378(a)0.1377(b)0.1384(a)0.0333(b)0.0044(a)Table 6.1: Results for one-photon ionization cross sections compared to those from the literature. F
is root mean square of laser electric field in atomic units. All the calculations above are performed
for n = 0, 1 except for F = 0.05, in which n is −1, 0, 1, 2 due to the strong coupling at high intensity.
The superscripts of the reference values: (a) stand for [243] and (b) for [233]for a much wider range of rotation angle. The two panels on the right show the trajectory of thequasienergy. In both cases, there is an angle, around which the trajectory of the quasienergy does
not change too much and is trapped, which is identified as stationary point, as indicated in theFig. 6.4(right). The stationary point for the 3 complex Gaussian basis sets is around 70o while forthe 5 complex Gaussian basis sets it is around 30o . From the trajectory plot, the trajectory of 5complex Gaussian basis sets has a better stationary quality than the other case. Therefore, in thefollowing calculation, we will use 5 Gaussian basis sets.
To compare the results with other papers [243, 233], we calculated the one-photon crosssections via the imaginary part of the quasienergy. For a M photon process, the cross section canbe calculated via
σM =

~Γω
,
IM

(6.59)

where Γω is the ionization rate and I is the laser intensity. The one-photon cross sections at
different intensities and frequencies are calculated and illustrated in the Table 6.1. Our results
obtained with only 5 complex Gaussian basis sets (each set includes both s and p components) arein good agreement with those of other papers. The exponents {αi } for different frequencies is given
in Table 6.2. We find that the optimized complex Gaussian basis sets do not depend much on theintensity of the laser field but only strongly depend on the frequency. For example, the results forthree different intensities but at the frequency ω = 0.6 are obtained with the same complex Gaussianbasis sets. for a broad range of intensities from 7 × 108 W/cm2 to 1.8 × 1014 W/cm2 , keeping the

129
ω = 0.6
5.623413
0.199059
0.030132
0.028864
0.006410

ω = 1.0
10.00000
0.630957
0.199526
0.125893
0.031623

ω = 2.0
7.617912
0.841466
0.621906
0.158489
0.100000

Table 6.2: The exponents of complex Gaussian basis sets for different laser frequencies. These basis
sets are tested for the range of laser intensity from 7 × 108 W/cm to 1.8 × 1014 W/cm2 .

excellent agreement with reference values. This is due to the fact that it is the laser frequency that
determines the energy of ionized outgoing wavepackets. The trajectories of eigenvalues for all of
this calculations are shown in Fig. 6.5. For all the cases, the stationary property around stationary

0 1e 7
1
2
3
4
5
6
7
08 07 06 05 04 03 02 01 00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.000.00
Ereal [a.u.] 4.989e 1

0.81e 7
1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4 5
5
5
5
0
0
0
0.52 0.52 0.51 0.51 0.50 0.50 0.49
Ereal [a.u.]

Eimag [a.u.]

Eimag [a.u.](a) F = 1 × 10−4 a.u., ω = 0.6a.u. (b) F = 7.08 × 10−4 a.u., ω =0.6a.u.

Eimag [a.u.]

0.0 1e 8
0.2
0.4
0.6
0.8
1.0
1.2
1.4
08 07 06 05 04 03 02 01 00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.000.00
Ereal [a.u.] 4.989e 1

Eimag [a.u.]

Eimag [a.u.]

point is very strong, demonstrating the validity of our calculations.
0.00000.00050.00100.00150.00200.00250.00300.0035 5
9 90 85 80 75 70 65 60 55
0.49 0.49 0.49 0.49 0.49 0.49 0.49 0.49 0.49
Ereal [a.u.](c) F = 0.05a.u., ω = 0.6a.u.

0.51e 10
0.0
0.5
1.0
1.5
2.0
2.5
2
8
0
6
4
0.50 0.50 0.50 0.50 0.50
Ereal [a.u.](d) F = 7.08 × 10−4 a.u., ω = (e) F = 1 × 10−4 a.u., ω = 2.0a.u.
1.0a.u.

Figure 6.5: Trajectory of quasienergy of H for different laser parameters. F is the mean square rootof laser electric field.

The ionization rates obtained using this method are very accurate and the computation timeis short. The ionization rate from numerical simulations is often obtained by applying a long laser

130
field with constant amplitude and then plot the time dependent population to obtain the slope as afunction of time which is equal to the ionization rate. With the Floquet method, the optimization
procedure for the one-photon ionization is very fast and typically optimization can be achieved
within 10 minutes. Furthermore, since the continuum states depend on the energy of the ionizedwavepacket, the optimization has to be done once only for a study as a function of laser intensity.
The calculation of the ionization rate itself takes less than 1 second, which is very fast comparedto numerical propagation method. On the other hand, Floquet method, also gives direct access to
the Stark shift which is usually difficult to obtain with other methods.

6.6.2

Two photon ionization from H atom

While for the one-photon ionization, it is very easy to get the stationary value, the caseof two photon ionization is more complicated. The complex Gaussian basis sets can be obtained
through optimization process in the same way as for one photon ionization. The maximum angularmomentum l should be 2 instead of 1 for this case. The bound states can be still represented with
6-31++G∗∗ Gaussian basis set. However, now further bound states are required since the twophoton ionization can go via intermediate states which are needed to be represented. Therefore,
for the two-photon ionization case, we have to add more Gaussian basis sets to represent possible
excited states and intermediate states.
For H atom, currently there are no real Gaussian basis sets for higher excited states available.
Instead of just using five Gaussian basis sets as in the one photon case, we added 15 complex
Gaussian basis sets including both s and d states during the optimization. The initial value foreach basis is still assigned evenly in the range determined by coefficient C.
We calculated the two-photon ionization of H atom at a laser frequency ω = 0.3 and intensity
4.3 × 1013 W/cm2 . After optimization, we calculated the eigenvalue for different rotation angle andthe trajectory is plotted in Fig. 6.6. Compared to one photon ionization, the imaginary part ofquasienergy is only flat for the angle between 30o and 60o . The stationary point can be found fromthe trajectory plot in Fig. 6.6(c). The optimized value of α is indicated by the vertical line, as131shown in Fig. 6.6(d).

(a)

(b)

(c)

(d)

Figure 6.6: Trajectory of quasienergy (a-c) for two photon ionization from H atom in a laser fieldwith ω = 0.3 and F = 0.025 and the value of exponents of complex Gaussian basis after optimization
are indicated by the vertical lines in (d).

Compared to the case of atom, for molecules, although the number of available Gaussianbasis sets for each atom does not change, the total number of Gaussian basis sets increases with thenumber of atoms and thus these Gaussian basis sets can represent more states than in the atomcase. Therefore, we do not need provide many complex Gaussian basis sets in the molecular caseas we will illustrate it in the next subsection.

6.6.3

Two photon ionization from H+
2We have also applied our method to the simplest molecule H+
2 . We have used the same real
Gaussian basis 6-31++G∗∗ , which represents the 1s, 2s, 2p, 3s and 3p state. For the continuum, 5

132
complex Gaussian basis sets including both s and p states centered on each atom are added. The
complex Gaussian basis sets on the two H atom are the same.
To compare with other papers, we have performed calculations for laser intensity I=1.76 ×
1013 W/cm2 at different frequencies ω. After optimization, we obtain the stationary quasienergy
and calculate the cross section via Eq. (6.59) with M = 2. The results are plotted (red dots) andcompared to reference values (blue curve) in Fig. 6.7. Our results are in good agreement withthe results reported in other papers [233, 244] over a broad range of frequencies even for the 2σu

Cross section [cm4 w−1 ]resonance corresponding to the first peak in the figure. Compared to the two-photon calculation

10-31
10-32
10-33
10-34
10-35
10-36
10-370.5 0.6 0.7 0.8 0.9 1.0 1.1
Photon energy [a.u.]

13 W/cm2 .
Figure 6.7: Cross section of two-photon ionization from H+
2 at laser intensity I=1.76 × 10
The solid curve is adapted from [211, 244] and the red dot is our calculations.

for H atom, we used the same real Gaussian basis set but here we only need 5 complex Gaussianbasis sets because of more states can be represented by the Gaussian basis sets. This simple test
implies that this Floquet method with Gaussian basis sets may work even for large molecules. The
optimal values of exponents {αi } of complex Gaussian basis sets in our calculation for H+
2 are given
in Table 6.3.

133
Table 6.3: Complex Gaussian exponents for two photon ionization
ω=0.6
7.02329e-1
4.09294e-2
2.34079e-2
8.35401e-4
7.91926e-4

6.7

ω=0.656.45221e-1
4.09330e-2
2.23400e-2
1.84340e-3
2.08791e-4

ω=0.7
6.45993e-1
4.09468e-2
2.33864e-2
5.20973e-3
5.67170e-5

ω=751.53936e1
6.023692.30414e-2
4.54144e-3
1.18183e-4

ω=0.832.55236e1
1.391903
1.000000
7.71358e-2
7.29241e-4

ω=0.8452.46937e1
1.072119.95687e-1
8.08699e-2
7.29315e-4

ω=0.9
7.33636
2.09922e-1
2.98224e-1
8.73335e-2
9.47368e-4

Summary
We have shown that the application of a combination of real and complex Gaussian basis setsin the Floquet method can help efficiently compute accurate quasienergies over a broad range of laser
intensities and frequencies. A very general optimization procedure to optimize complex Gaussianbasis sets has been developed. Although the optimization problem is a non-convex problem andthe optimization results depend on the initial value of the exponents of Gaussian basis sets, wefind that it is still possible to get accurate results with initial values based on the energy of ionizedwavepackets. The cross sections for one-photon and two-photon ionization for both H and H+
2 arein agreement with those reported in literature. For multiphoton ionization, more basis sets will beneeded for the excited and intermediate states. However, with the same available real Gaussianbasis sets, more states can be represented for molecules than for atoms by the same Gaussian basis
sets. Although we only applied it to one-electron case, it is straightforward to move to multielectron
system by applying such technique in the Kohn-Sham equations directly in the future.

Chapter 7

Summary

In this thesis, we have systematically studied multielectron effects in high-order harmonic
generation, ionization, and nonadibatic electron localization in the interaction between molecules
and an intense laser field by applying time-dependent density-functional theory to solve the multielectron Schrödinger equation. The results of our numerical simulations and theoretical analysis
showed that, in particular, the field-induced coupling between different molecular orbitals can havea strong influence on these fundamental processes.
We first analyzed multielectron effects in the high-order harmonic generation spectra in Chapter 3. We found that the laser induced coupling can lead to fractional harmonics in the form of
Mollow sidebands over the whole spectra of harmonics and modification of the traditional threestep model for the HHG process. Such phenomenon has been analyzed in detail for the case ofN+
2 but is found to appear in general for other open-shell molecules as well. We further analyzed
experimental data for the ellipticity of harmonics from N2 and CO2 . Our multielectron simulation results are in good agreement with experimental results. By analyzing the response from all
valence orbitals to the ellipticity pattern, we have found that contributions from inner valence orbitals play a significant role in the HHG process for molecules, which indicates a breakdown of thesingle-active-electron model.
In Chapter 4 we studied the coupling effects on ionization which is another important strongfield process. Our results show that the laser-induced coupling between inner valence orbital and
HOMO can significantly enhance the ionization from inner shell orbitals. These results do not135
only provide another indication of the importance of multielectron effects in strong-field interaction
with molecules, but they may also explain the population inversion in ions, recently observed infilamentation experiments.
Since electron localization in molecules can strongly modify the ionization dynamics, wehave further studied the effects of the field-induced coupling of states on the electron dynamics inmolecules. Analysis of results presented in Chapter 5 show that two competing processes, in which
the expected oscillation of the electron density with the change of external field is superimposed bya dynamics due to the Rabi flopping of the coupled states, lead a strong nonadiabatic dynamics.
The mechanism and the control of the two competing processes via the parameters of the field are
analyzed.
In Chapter 6, we have introduced a new intense-field method using the Floquet theorem
and complex Gaussian basis sets. Floquet theorem is applied to transform the time-dependent
Schrödinger equation into a set of time-independent equations. Complex Gaussian basis sets are
employed to represent the continuum states and an optimization procedure has been developed toobtain the exponents of complex Gaussian basis sets. The application of this method to one- andtwo-photon ionization yield good agreement with results in literature and the possibility to extend
the method to larger molecules has been discussed.

Bibliography[1] S. V. Popruzhenko. Keldysh theory of strong field ionization: history, applications, difficulties
and perspectives. J. Phys. B: At. Mol. Opt. Phys., 47:204001, 2014.
[2] F. H. M. Faisal. Theory of multiphoton processes. Plenum Press, 1987.
[3] C. J. Joachain, N. J. Kylstra, and R. M. Potvliege. Atoms in intense laser fields. Cambridge
University Press, 2014.
[4] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman. Free-free transitions
following six-photon ionization of xenon atoms. Phys. Rev. Lett., 42:1127–1130, 1979.
[5] A. McPherson, G. Gibson, U. Jara, H. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K.
Rhodes. Studies of multiphoton production of vacuum-ultraviolet radiation in rare gase. J.
Opt. Soc. Am. B, 4:595, 1987.
[6] X. F. Li, A. L’Huillier, M. Ferray, L. A. Lompré, and G. Mainfray. Multiple-harmonic
generation in rare gases at high laser intensity. Phys. Rev. A, 39:5751–5761, 1989.
[7] B. Walker, B. Sheehy, L. F. DiMauro, P. Agostini, K. J. Schafer, and K. C. Kulander. Precision
measurement of strong field double ionization of helium. Phys. Rev. Lett., 73:1227–1230, 1994.
[8] A. Talebpour, C.-Y. Chien, Y. Liang, S. Larochelle, and S. L. Chin. Non-sequential ionizationof Xe and Kr in an intense femtosecond Ti:sapphire laser pulse. J. Phys. B: At. Mol. Opt.
Phys., 30:1721, 1997.
[9] L. Keldysh. Ionization in the field of a strong electromagnetic wave. Sov. Phys. JETP,
20:1307, 1965.
[10] F. H. M. Faisal. Multiple absorption of laser photons by atoms. J. Phys. B, 6:L89, 1973.
[11] H. R. Reiss. Effect of an intense electromagnetic field on weakly bound system. Phys. Rev.
A, 22:1786, 1980.
[12] A. M. Perelomov, S. V. Popov, and M. V. Terent’ev. Ionization of atoms in an alternating
electric field. Sov. Phys. JETP, 23:924, 1966.
[13] M. V. Ammosov, N. B. Delone, and V. P. Krainov. Tunnel ionization of complex atoms andatomic ions in an alternating electromagnetic field. Sov. Phys. JETP, 64:1191, 1986.
[14] L. D. Landau and E. M. Lifshitz. Quantum Mechanics. Addison Wesley, 1958.

137[15] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich. Generation of optical harmonics.
Phys. Rev. Lett., 7:118–119, 1961.
[16] T. Popmintchev, M. C. Chen, D. Popmintchev, S. Alisauskas, G. Andriukaitis, T. Balciunas, A. Pugzlys, A. Baltuska, M. Murnane, and H. Kapteyn. Bright coherent attosecondto-zeptosecond kiloelectronvolt X-ray supercontinua. In Laser Applications to Photonic
Applications, page PDPC12, 2011.
[17] T. Popmintchev, M. C. Chen, D. Popmintchev, P. Arpin, S. Brown, S. Alisaukas, G. Andriukaitis, T. Balciunas, O. Mücke, A. Pugzlys, A. Baltuska, B. Shim, S. E. Schrauth,A. Gaeta, C. Hernandez-Garcia, L. Plaja, A. Becker, A. Jaroń-Becker, M. Murnane, and
H. Kapteyn. Bright coherent ultrahigh harmonics in the keV X-Ray regime from mid-infrared
femtosecond lasers. Science, 336:12878, 2012.
[18] I. P. Christov, M. M. Murnane, and H. C. Kapteyn. High-harmonic generation of attosecond
pulses in the ”single-cycle” regime. Phys. Rev. Lett., 78:1251, 1997.
[19] M. Hentschel, R. Kienberger, Ch. Spielmann, G. A. Reider, N. Milosevic, T. Brabec,
P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz. Attosecond metrology. Nature,
414:509, 2001.
[20] P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, Ph. Balcou, H. G. Muller, and
P. Agostini. Observation of a train of attosecond pulses from high harmonic generation.
Science, 292, 2001.
[21] K. Midorikawa. High-order harmonic generation and attosecond science. Jpn. J. Appl. Phys.,
50:090001, 2011.
[22] T. Brabec and F. Krausz. Intense few-cycle laser fields: Frontiers of nonlinear optics. Mod.
Phys., 72:545, 2000.
[23] P. Corkum and F. Krausz. Attosecond science. Nat. Phys., 3:381, 2007.
[24] N. B. Delone and V. P. Krainov. Multiphoton processes in atoms. Springer, 2000.
[25] Z. Vager, R. Naaman, and E. P. Kanter. Coulomb explosion imaging of small molecules.
Science, 244:426–431, 1989.
[26] P. H. Bucksbaum, A. Zavriyev, H. Muller, and D. W. Schumacher. Softening of the H2+molecular bond in intense laser fields. Phys. Rev. Lett., 64:1883, 1990.
[27] A. Giusti-Suzor, X. He, O. Atabek, and F. H. Frederic. Above-threshold dissociation of H2+
in intense laser fields. Phys. Rev. Lett., 64:515, 1990.
[28] A. Etches, M. B. Gaarde, and L. B. Madsen. Two-center minima in harmonic spectra fromaligned polar molecules. Phys. Rev. A, 84:023418, 2011.
[29] M. Lein, N. Hay, R. Velotta, J. P. Marangos, and P. L. Knight. Interference effects in highorder harmonic generation with molecules. Phys. Rev. A, 66:023805, 2002.
[30] J. Muth-Böhm, A. Becker, and F. H. M. Faisal. Suppressed molecular ionization for a classof diatomics in intense femtosecond laser fields. Phys. Rev. Lett., 85:2280, 2000.

138[31] A. Jaroń-Becker, A. Becker, and F. H. M. Faisal. Saturated ionization of fullerenes in intenselaser fields. Phys. Rev. Lett., 96:143006, 2006.
[32] S. Patchkovskii, Z. Zhao, T. Brabec, and D. M. Villeneuve. High harmonic generation andmolecular orbital tomography in multielectron systems: Beyond the single active electron
approximation. Phys. Rev. Lett., 97:123003, 2006.
[33] B. K. McFarland, J. P. Farrel, P. H. Bucksbaum, and M. Guhr. High harmonic generationfrom multiple orbitals in N2 . Science, 21:1232, 2008.
[34] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. Corkum, andM. Yu. Ivanov. High harmonic interferometry of multi-electron dynamics in molecules.
Nature, 460:972, 2009.
[35] B. B. Augstein and C. Figueira DE Morisson Faria. High-order harmonic generation indiatomic molecules: Quantum interference, nodal structures and multiple orbitals. Mod.
Phys. Lett. B, 26:1130002, 2012.
[36] D. Dundas. Multielectron effects in high harmonic generation in N2 and benzene: Simulation
using a non-adiabatic quantum molecular dynamics approach for laser-molecule interactions.
J. Chem. Phys, 136, 2012.
[37] D. Dundas and J. M. Rost. Molecular effects in the ionization of N2 , O2 , and F2 by intenselaser fields. Phys. Rev. A, 71:013421, 2005.
[38] M. Leibscher, I. Sh. Averbukh, and H. Rabitz. Molecular alignment by train of short laser
pulses. Phys. Rev. Lett., 90:213001, 2003.
[39] J. Itatani, J. Levesque, D. Zeidler, H. Niikura, H. Ppin, J. C. Kieffer, P. B. Corkum, andD. M. Villeneuve. Tomographic imaging of molecular orbitals. Nature, 432:867–871, 2004.
[40] M. Lezius, V. Blanchet, D. M. Rayner, D. M. Villeneuve, A. Stolow, and M. Yu. Ivanov.
Nonadiabatic multielectron dynamics in strong field molecular ionization. Phys. Rev. Lett.,
86:51, 2001.
[41] M. Lezius, V. Blanchet, M. Yu. Ivanov, and A. Stolow. Polyatomic molecules in strong laser
field: Nonadiabatic multielectron dynamics. J. Chem. Phys., 117:1575, 2002.
[42] A. Stolow. The three pillars of photo-initiated quantum molecular dynamics. Farad. Discuss.,
117:1575, 2013.
[43] J. Parker, K. T. Taylor, C. W. Clark, and S. Blodgett-Fort. Intense-field multiphoton ionization of a two-electron atom. J. Phys. B: At. Mol. Opt. Phys., 29, 1996.
[44] H. G. Muller. Low-frequency above-threshold ionization of a model two-electron atom. Opt.
Exp., 8, 2001.
[45] J. S. Parker, B. J. S. Doherty, K J. Meharg, and K. T. Taylor. Time delay between singly
and doubly ionization wavepackets in laser-driven helium. J. Phys. B: At. Mol. Opt. Phys.,
36, 2003.
[46] A. Jaroń-Becker, A. Becker, and F. H. M. Faisal. Ionization of N2 , O2 , and linear carbon
clusters in a strong laser pulse. Phys. Rev. A, 69, 2004.

139[47] A. Becker and F. H. M. Faisal. S-matrix theory of two-electron momentum distribution
produced by double ionization in intense laser fields. Opt. Exp., 8, 2001.
[48] J. Muth-Böhm, A. Becker, S. L. Chin, and F. H. M. Faisal. S-matrix theory of ionization ofpolyatomic molecules in an intense laser pulse. Chem. Phys. Lett., 337:313, 2001.
[49] A. Jaroń-Becker, A. Becker, and F. H. M. Faisal. Signature of molecular orientation andorbital symmetry in strong-field photoelectron angular and energy distribution of diatomic
molecules and small carbon clusters. Laser Phys., 14:179, 2004.
[50] A. Jaroń-Becker, A. Becker, and F. H. M. Faisal. Dependence of strong-field photoelectron
angular distribution on molecular orientation. J. Phys. B: At. Mol. Opt. Phys., 36:L375, 2003.
[51] T. Klamroth. Laser-driven electron transfer through metal-insulator-metal contacts: Timedependent configuration interaction singles calculations for a jelium model. Phys. Rev. B,
68:245421, 2003.
[52] L. Greenman, P. J. Ho, S. Pabst, E. Kamarchik, D. A. Mazziotti, and R. Santra. Implementation of time-dependent configuration-interaction singles method for atomic strong-fieldprocesses. Phys. Rev. A, 74:043420, 2006.
[53] N. Rohringer, A. Gordon, and R. Santra. Configuration-interaction-based time-dependent
orbital approach for ab initio treatment of electronic dynamics in a strong optical laser field.
Phys. Rev. A, 74:043420, 2006.
[54] S. M. Smith, X. Li, A. N. Markevitch, D. A. Romanov, R. J. Levis, and H. B. Schlegel. A
numerical simulation of nonadiabatic electron excitation in the strong field regime: Linear
polyenes. J. Phys. Chem. A, 109:5176, 2005.
[55] J. Zanghellini, M. Kitzler, C. Fabian, T. Brabec, and A. Scrinzi. An MCTDHF approach tomultilectron dynamics in laser fields. Laser Phys., 13:1064, 2003.
[56] M. Nest, T. Klamroth, and P. Saalfrank. The multiconfiguration time-dependent HartreeFock method for quantum chemical calculation. J. Chem. Phys., 122:124102, 2005.
[57] G. Jordan, J. Caillat, C. Ede, and A. Scriniz. Strong field ionization of linear molecules: A
correlated three-dimensional calculation. J. Phys. B: At. Mol. Opt. Phys., 39:S341, 2006.
[58] T. Kato and H. Kono. Time-dependent multiconfiguration theory for electronic dynamics ofmolecules in an intense laser field. Chem. Phys. Lett., 392:533, 2004.
[59] T. Otobe and K. Yabana. Density-functional calculation for the tunnel ionization rate ofhydrocarbon molecules. Phys. Rev. A, 75:062507, 2007.
[60] A. Castro, M. A. L. Marques, J. A. Aloso, G. F. Bertsch, and A. Rubio. Excited states
dynamics in time-dependent density-functional theory: High-field molecular dissociation andharmonic generation. Eur. Phys. J. D, 28:211, 2004.
[61] E. Runge and E. K. U. Gross. Density-functional theory for time-dependent systems. Phys.
Rev. Lett., 52:997, 1984.

140[62] D. A. Telnov and Shih-I Chu. Effect of electron structure and multiphoton dynamical responseon strong-field multiphoton ionization of diatomic molecules with arbitrary orientation: An
all-electron time-dependent density-functional-theory approach. Phys. Rev. A, 79:041401,
2009.
[63] E. P. Fowe and A. D. Bandrauk. Nonlinear time-dependent density-functional-theory studyof ionization and harmonic generation in CO2 by ultrashort intense laser pulses: Orientation
effects. Phys. Rev. A, 81:023411, 2010.
[64] J. Liu, Z. Guo, J. Sun, and W. Liang. Theoretical studies on electronic spectroscopy anddynamics with the real-time time-dependent density-functional theory. Front. Chem. China,
5:11, 2010.
[65] M. A. L. Marques, X. López, D. Varsano, A. Castro, and A. Rubio. Time-dependent densityfunctional approach for biological chromophores: The case of the green fluorescent protein.
Phys. Rev. Lett., 90:258101, 2003.
[66] A. Castro, M. Isla, J. Martı́nez, and J. A. Alonso. Scattering of a proton with the Li4 cluster:Non-adiabatic molecular dynamics description based on time-dependent density-functional
theory. Chem. Phys., 399:130–134, 2012.
[67] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Rubio. Octopus: a first-principles toolfor excited electron–ion dynamics. Comput. Phys. Commun., 151:60–78, 2003.
[68] S. Botti, A. Castro, N. Lathiotakis, X. Andrade, and M. A. L. Marques. Optical and magnetic
properties of boron fullerenes. Phys. Chem. Chem. Phys., 11:4523–4527, 2009.
[69] S. I. Chu. Recent development of self-interaction-free time-dependent density functional theory for nonperturbative treatment of atomic and molecular multiphoton processes in intenselaser fields. J. Chem. Phys., 123:062207, 2005.
[70] X. Chu and S. I. Chu. Role of the electronic structure and multielectron responses in ionization
mechanisms of diatomic molecules in intense short-pulse lasers: An all-electron ab initio study.
Phys. Rev. A, 70:061402(R), 2004.
[71] X. Chu and S. I. Chu. Time-dependent density-functional theory for molecular processes instrong fields: Study of multiphoton processes and dynamical response of individual valence
electrons of N2 in intense laser fields. Phys. Rev. A, 64:063404, 2001.
[72] X. Chu and S. I. Chu. Self-interaction-free time-dependent density functional theory formolecular processes in intense laser fields. In 8th International Conference on Multiphoton
Processes, volume 525 of AIP Conference Proceedings, pages 415–426, Melville, N.Y., 2000.
American Institute of Physics.
[73] M. Oliveira and F. Nogueira. Generating relativistic pseudo-potentials with explicit incorporation of semi-core states using ape, the atomic pseudo-potentials engine. Comput. Phys.
Comm., 178, 2008.
[74] G. B. Bachelet, D. R. Hamann, and M. Schluter. Relativistic norm-conserving pseudopotentials. Phys. Rev. B, 25:2103, 1982.
[75] R. M. Dreizler and E. K. Gross. Density Functional Theory. Springer, 1990.

141[76] C. Fiolhais, F. Nogueira, and M. Margues. A Primer in Density Functional Theory. Springer,
2003.
[77] W. Koch and M. C. Holthausen. A Chemist’s Guide to Density Functional Theory. WileyVCH Verlag GmbH, 2001.
[78] R. G. Parr and W. Yang. Density-Functional Theory of Atoms and Molecules. Oxford
University Press, 1989.
[79] J.M. Combes, P. Duclos, and R. Seiler. The Born-Oppenheimer approximation. In Rigorous
Atomic and Molecular Physics, volume 74 of NATO Advanced Study Institutes Series, pages
185–213. Springer US, 1981.
[80] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:B864, 1964.
[81] M. Levy. Electron densities in search of hamiltonians. Phys. Rev. A, 26:1200, 1982.
[82] E. H. Lieb. Density Functional Methods in Physics. Plenum, New York, 1985.
[83] A. K. Theophilou. The energy density functional formalism for excited states. J. Phys. C,
12:5419, 1979.
[84] O. Gritsenko and E. J. Baerends. The spin-unrestricted molecular Kohn-Sham solution andthe analogue of Koopmanss theorem for open-shell molecules. J. Chem. Phys., 120:8364,
2004.
[85] P. Strange. Relativistic Quantum Mechanics with Application in Condensed Matter andAtomic Physics. Cambridge University Press, 1998.
[86] K. Capelle. A bird’s-eye view of density-functional theory. Braz. J. Phys., 36:1318–1343,
2006.
[87] J. T. Chayes, L. Chayes, and M. B. Ruskai. Density functional approach to quantum latticesystems. J. Stat. Phys., 38:497, 1985.
[88] P. E. Lammert. Coarse-grained v representability. J. Chem. Phys., 125:074114, 2006.
[89] C. A. Ullrich and W. Kohn. Degeneracy in density functional theory: Topology in the ν andn spaces. Phys. Rev. Lett., 89:156401, 2002.
[90] M. A.L. Miguel, O. J.T. Micael, and B. Tobias. Libxc: A library of exchange and correlation
functionals for density functional theory. Comput. Phys. Commun., 183:2272 – 2281, 2012.
[91] E. S. Gustavo and N. S. Viktor.
Chapter 24 - Progress in the development of
exchange-correlation functionals. In Theory and Applications of Computational Chemistry,
pages 669 – 724. Elsevier, 2005.
[92] F. Bloch. Bemerkung zur Elektronentheorie des Ferromagnetismus und der elektrischen
Leitfähigkeit. Z. Phys., 57:545, 1929.
[93] P. A. M. Dirac. Note on exchange phenomena in the thomas atom. Proc. Camb. Phil. Soc.,
26:376, 1930.

142[94] P. Gombás. Die Statistische Theorie des Atoms und Ihre Anwendungen. Springer-Verlag
Wien, 1949.
[95] J. C. Slater. Quantum Theory of Atomic Structure, volume 2. McGraw-Hill, New York, 1960.
[96] P. R. Antoniewicz and L. Kelinman. Kohn-Sham exchange potential exact to first order inρ(K)/ρ0 . Phys. Rev. B, 31:6779, 1985.
[97] L. Kleinman and S. H. Vosko. Gradient expansion of the exchange-energy density functional:
Effect of taking limits in the wrong order. Phys. Rev. B, 37:4634, 1988.
[98] J. A. Chevary and S. H. Vosko. Importance of Pauli-principle restrictions for accurate numerical evaluation of the exchange-only screening function. Phys. Rev. B, 42:5320, 1990.
[99] L. J. Sham. Computational Methods in Band Theory, page 458. Plenum Press, 1971.
[100] D. C. Langreth and J. P. Perdew. Theory of nonuniform electronic systems. I. Analysis of
the gradient approximation and a generalization that works. Phys. Rev. B, 21:5469, 1980.
[101] D. C. Langreth and M. J. Mehl. Beyond the local-density approximation in calculations ofground-state electronic properties. Phys. Rev. B, 28:1809, 1983.
[102] J. P. Perdew. Accurate density functional for the energy: Real-space cutoff of the gradient
expansion for the exchange hole. Phys. Rev. Lett., 55:1665, 1985.
[103] A. D. Becke. Density functional calculations of molecular bond energies. J. Chem. Phys.,
84:4524, 1986.
[104] A. D. Becke. Density-functional exchange-energy approximation with correct asymptotic
behavior. Phys. Rev. A, 38:3098, 1988.
[105] A. E. DePristo and J. D. Kress. Rational function representation for accurate exchange energyfunctionals. J. Chem. Phys., 86:1425, 1987.
[106] J. P. Perdew and Y. Wang. Accurate and simple analytic representation of the electron-gas
correlation energy. Phys. Rev. B, 45:13244, 1992.
[107] A. D. Becke and M. R. Roussel. Exchange holes in inhomogeneous systems: A coordinatespace model. Phys. Rev. A, 39:3761, 1989.
[108] A. D. Becke. Density-functional thermochemistry. V. Systematic optimization of exchangecorrelation functionals. J. Chem. Phys., 107:8554, 1997.
[109] A. D. Becke. Exploring the limits of gradient corrections in density functional theory. J.
Comput. Chem., 20:63, 1999.
[110] A. D. Becke. A new mixing of Hartree-Fock and local density-functional theories. J. Chem.
Phys., 98:1372, 1993.
[111] M. Levy, N. H. March, and N. C. Handy. On the adiabatic connection method, and scalingof electronelectron interactions in the Thomas-Fermi limit. J. Chem. Phys., 104:1989, 1996.
[112] R. van Leeuwen and E. J. Baerends. Exchange-correlation potential with correct asymptotic
behavior. Phys. Rev. A, 49:2421–2431, 1994.

143[113] P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J. Baerends. Molecular
calculations of excitation energies and (hyper)polarizabilities with a statistical average oforbital model exchange-correlation potentials. J. Chem. Phys., 112:1344–1352, 2000.
[114] M. Gell-Mann and K. A. Brueckner. Correlation energy of an electron gas at high density.
Phys. Rev., 106:364, 1957.
[115] Jr. W. J. Carr and A. A. Maradudin. Ground-state energy of a high-density electron gas.
Phys. Rev., 133:A371, 1964.
[116] P. Nozières and D. Pines. Correlation energy of a free electron gas. Phys. Rev., 111:442,
1958.
[117] Jr. W. J. Carr. Energy, specific heat, and magnetic properties of the low-density electron gas.
Phys. Rev., 122:1437, 1961.
[118] D. M. Ceperley and B. J. Alder. Ground state of the electron gas by a stochastic method.
Phys. Rev. Lett., 45:566, 1980.
[119] J. P. Perdew and A. Zunger. Self-interaction correction to density-functional approximations
for many-electron systems. Phys. Rev. B, 23:5048, 1981.
[120] S. H. Vosko, L. Milk, and M. Nusair. Accurate spin-dependent electron liquid correlationenergies for local spin density calculations: a critical analysis. Can. J. Phys., 58:1200, 1980.
[121] N. Troullier and J. L. Martines. Efficient pseudopotentials for plane-wave calculations. Phys.
Rev. B, 43, 1991.
[122] L. Kleinman and D. M. Bylander. Efficacious form for model pseudopotentials. Phys. Rev.
Lett., 48, 1982.
[123] D. Vanderbilt. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism.
Phys. Rev. B, 41:7892, 1990.
[124] G Seifert, D Porezag, and Th Frauenheim. Calculations of molecules, clusters, and solids
with a simplified LCAO-DFT-LDA scheme. Int. J. Quant. Chem., 58:185–192, 1996.
[125] L. Pan, K. T. Taylor, and C. W. Clark. Perturbation theory study of high-harmonic generation. J. Opt. Soc. Am. B, 7:509, 1990.
[126] M. C. Kohler, T. Pfeifer, K. Z. Hatsagortsyan, and C. H. Keitel. Frontiers of atomic highharmonic generation. arXiv, 1201:5094v1, 2012.
[127] R. A. Ganeev. Harmonic generation in laser-produced plasmas containing atoms, ions andclusters: a review. J. Mod. Opt., 59:409, 2012.
[128] M. C. Chen, C. Mancuso, C. Hernandez-Garcia, F. Dollar, B. Galloway, D. Popmintchev,
P. C. Huang, B. Walker, L. Plaja, A. A. Jaroń-Becker, A. Becker, M. M. Murnane, H. C.
Kapteyn, and T. Popmintchev. Generation of bright isolated attosecond soft X-ray pulses
driven by multicycle midinfrared lasers. Proc. Nat. Acad. Sciences, 111:E2361–E2367, 2014.

144[129] D.M. Villeneuve, J.B. Bertrand, P.B. Corkum, N. Dudovich, J. Itatani, J.C. Kieffer, F. Lgar,
J. Levesque, Y. Mairesse, H. Niikura, B.E. Schmidt, A.D. Shiner, and H.J. Wrner. Studying
the electronic structure of molecules with high harmonic spectroscopy. In Attosecond Physics,
volume 177 of Springer Series in Optical Sciences, pages 159–190. 2013.
[130] E. Frumker, N. Kajumba, J. B. Bertrand, H. J. Wörner, C. T. Hebeisen, P. Hockett, M. Spanner, S. Patchkovskii, G. G. Paulus, D. M. Villeneuve, A. Naumov, and P. B. Corkum. Probing
polar molecules with high harmonic spectroscopy. Phys. Rev. Lett., 109:233904, 2012.
[131] J. B. Bertrand, H. J. Wörner, P. Hockett, D. M. Villeneuve, and P. B. Corkum. Revealing thecooper minimum of N2 by molecular frame high-harmonic spectroscopy. Phys. Rev. Lett.,
109:143001, 2012.
[132] S. Haessler, J. Caillat, W. Boutu, C. Giovanetti-Teixeira, T. Ruchon, T. Auguste, Z. Diveki,
P. Breger, A. Maquet, B. Carré, R. Taieb, and P. Salieres. Attosecond imaging of molecularelectronic wavepackets. Nat. Phys., 6:200–206, 2010.
[133] W. Boutu, S. Haessler, H. Merdji, P. Breger, G. Waters, M. Stankiewicz, L. J. Frasinski,
R. Taieb, J. Caillat, A. Maquet, P. Monchicourt, B. Carre, and P. Salieres. Coherent control
of attosecond emission from aligned molecules. Nat. Phys, 4:545, 2008.
[134] W. Li, X. Zhou, R. Lock, S. Patchkovskii, A. Stolow, H. C. Kapteyn, and M. M. Murnane.
Time-resolved dynamics in N2 O4 probed using high harmonic generation. Science, 322:1207–
1211, 2008.
[135] O. Smirnova, S. Patchkovskii, Y. Mairesse, N. Dudovich, D. Villeneuve, P. Corkum, and
M. Y. Ivanov. Attosecond circular dichroism spectroscopy of polyatomic molecules. Phys.
Rev. Lett., 102:063601, 2009.
[136] C. Vozzi, M. Nero, F. Calegari, G. Sansone, M. Nisoli, S. D. Silverstri, and S. Stagira.
Generalized molecular orbital tomography. Nat. phys., 7:822–826, 2011.
[137] H. J. Wörner, J. B. Bertrand, D. V. Kartashov, P. B. Corkum, and D. M. Villeneuve. Following a chemical reaction using high-harmonic interferometry. Nature, 466:604, 2010.
[138] C. Altucci, R. Velotta, E. Heesel, E. Springate, J. P. Marangos, C. Vozzi, E. Benedettti,
F. Calegari, G. Sansone, S. Stagira, M. Nisoli, and V. Tosa. High-order harmonics generationin alkanes. Phys. Rev. A, 73:043411, 2006.
[139] R. Torres, T. Siegel, L. Brugnera, I. Procino, Jonathan G. Underwood, C. Altucci, R. Velotta,
E. Springate, C. Froud, I. C. E. Turcu, S. Patchkovskii, M. Yu. Ivanov, O. Smirnova, andJ. P. Marangos. Revealing molecular structure and dynamics through high-order harmonic
generation driven by mid-IR fields. Phys. Rev. A, 81:051802, 2010.
[140] P. B. Corkum. Plasma perspective on strong field multiphoton ionization. Phys. Rev. Lett.,
71:1994, 1993.
[141] M. Lewenstein, P. Balcou, M. Yu. Ivanov, A. L’Hullier, and P. B. Corkum. Theory of highharmonic generation by low frequency laser fields. Phys. Rev. A, 49:2117, 1994.

145[142] H. Sakai, J. L. Larsen, L. Wendt-Larsen, J. Olesen, P. B. Corkum, and H. Stapelfeldt. Nonsequential double ionization of D2 molecules with intense 20-fs pulses. Phys. Rev. A, 67:063404,
2003.
[143] M. Lein, J. P. Marangos, and P. L. Knight. Electron diffraction in above-threshold ionizationof molecules. Phys. Rev. A, 66:051404, 2002.
[144] T. Popmintchev, M. C. Chen, P. Arpin, M. M. Murnane, and H. C. Kapteyn. The attosecond
nonlinear optics of bright coherent X-ray generation. Nat. Photonics, 4:822–832, 2010.
[145] D. M. Volkov. über eine Klasse von Ls̈ungen der Diracschen Gleichung. Z. Physik, 94:250,
1935.
[146] I. V. Litvinyuk, Kevin F. Lee, P. W. Dooley, D. M. Bylander, D. M. Villeneuve, andP. B. Corkum. Alignment-dependent strong field ionization of molecules. Phys. Rev. Lett.,
90:233003, 2003.
[147] T. Kanai, S. Minemoto, and H. Sakai. Quantum interference during high-order harmonic
generation from aligned molecules. Nature, 435:470–474, 2005.
[148] M. Lein, N. Hay, R. Velotta, J. P. Marangos, and P. L. Knight. Role of the intramolecularphase in high-harmonic generation. Phys. Rev. Lett., 88:183903, 2002.
[149] M. Lein, P. P. Corso, J. P. Marangos, and P. L. Knight. Orientation dependence of high-orderharmonic generation in molecules. Phys. Rev. A, 67:023819, 2003.
[150] S. K. Son, D. A. Telnov, and S. I. Chu. Probing the origin of elliptical high-order harmonic
generation from aligned molecules in linearly polarized laser fields. Phys. Rev. A, 82:043829,
2010.
[151] H. Sakai and K. Miyazaki. High-order harmonic-generation in nitrogen molecules with subpicosecond visible dye-laser pulse. Appl. Phys. B, 61:493–498, 1995.
[152] Y. Liang, S. Augst, S. L. Chin, Y. Beaudoin, and M. Chaker. High harmonic generation in
atomic and diatomic molecular gases using intense picosecond laser pulse-a comparison. J.
Phys. B, 27:5119–5130, 1994.
[153] C. Lyngå, A. L’Hullier, and C.-G Wahlström. High-order harmonic genration in moleculargases. J. Phys. B, 29:3293–3302, 1996.
[154] Y. Mairesse, J. Higuet, N. Dudovich, D. Shafir, B. Fabre, E. Mevel, E. Constant,
S. Patchkovskii, Z. Walters, M. Y. Ivanov, and O. Smirnova. High harmonic spectroscopy ofmultichannel dynamics in strong-field ionization. Phys. Rev. Lett., 104:213601, 2010.
[155] Y. Xia and A. Jaroń-Becker. Multielectron contributions in elliptically polarized high-orderharmonic emission from nitrogen molecules. Opt. Lett., 39:1461–1464, 2014.
[156] Y. Xia and A. Jaroń-Becker. Mollow sidebands in high order harmonic spectra of molecules.
Opt. Exp., 24:4689, 2016.
[157] D. M. Gaudiosi, B. Reagan, T. Popmintchev, M. Grisham, M. Berrill, O. Cohen, B. C. Walker,
M. M. Murnane, H. C. Kapteyn, and J. J. Rocca. High harmonic generation from ions in acapillary discharge plasma waveguide. Opt. Photonics News, 17:44–44, 2006.

146[158] D. Popmintchev, C. Hernández-Garcı́a, F. Dollar, C. Mancuso, J. Pérez-Hernández, M.C.
Chen, A. Hankla, X. Gao, B. Shim, A. Gaeta, M. Tarazkar, D. Romanov, R. Levis, J. Gaffney,
M. Foord, S. Libby, A. Jaroń-Becker, A. Becker, L. Plaja, M. Murnane, H. Kapteyn, andT. Popmintchev. Efficient soft X-ray high harmonic generation in multiply-ionized plasmas:
the ultraviolet surprise. Science, 350:1225, 2015.
[159] P. Wei, X. Yuan, Z. Zeng, Y. Zheng, J. Jiang, X. Ge, and R. Li. Enhanced high-order harmonic
generation from spatially prepared filamentation in argon. Opt. Exp., 29:17229–17236, 2015.
[160] R. A. Ganeev and H. Kuroda. Extremely broadened high-order harmonics generated by the
femtosecond pulses propagating through the filaments in air. Appl. Phys. Lett., 95:201117,
2009.
[161] B. R. Mollow. Power spectrum of light scattered by two-level systems. Phys. Rev., 188, 1969.
[162] D. E. Browne and C. H. Keitel. Resonance fluorescence in intense laser fields. J. Mod. Opt.,
47:1307, 2000.
[163] T. Zuo, S. Chelkowski, and A. D. Bandrauk. Photon-emission spectra of the H2+ molecularion in an intense laser field. Phys. Rev. A, 49:3943, 1994.
[164] C. R. Jacob and M. Reiher. Spin in density-functional theory. Int. J. Quant. Chem., 112:3661,
2012.
[165] X. Zhou, R. Lock, N. Wagner, W. Li, H. C. Kapteyn, and M. M. Murnane. Ellipticallypolarized high-order harmonic emission from molecules in linearly polarized laser fields. Phys.
Rev. Lett., 102:073902, 2009.
[166] X. B. Zhou, R. Lock, W. Li, N. Wagner, M. M. Murnane, and H. C. Kapteyn. Molecular
recollision interferometry in high harmonic generation. Phys. Rev. Lett., 100:073902, 2008.
[167] A. T. Le, R. R. Lucchese, and C. D. Lin. Polarization and ellipticity of high-order harmonics
from aligned molecules generated by linearly polarized intense laser pulses. Phys. Rev. A, 82,
2010.
[168] J. Levesque, Y. Mairesse, N. Dudovich, H. Pépin, J. C. Kieffer, P. B. Corkum, and D. M.
Villeneuve. Polarization state of high-order harmonic emission from aligned molecules. Phys.
Rev. Lett., 99:243001, 2007.
[169] L. Holmegaard, J. L. Hansen, L. Kalhøj, S. L. Kargh, H. Stapelfeldt, F. Filsinger, J. Küpper,
G. Meijer, D. Dimitrovski, M. Abu-samba, C. P. J. Martiny, and L. B. Madsen. Photoelectron
angular distributions from strong-field ionization of oriented molecules. Nat. Phys., 6:428,
2010.
[170] C. Z. Bisgaard, O. J. Clarkin, G. Wu, A. M. D. Lee, O. Gebner, C. C. Hayden, and A. Stolow.
Time-resolved molecular frame dynamics of fixed-in-space CS2 molecules. Science, 323:1464–
1468, 2009.
[171] G. N. Gibson, R. R. Freeman, and T. J. Mcllrath. Dynamics of the high-intensity multiphoton
ionization of N2 . Phys. Rev. Lett., 67:1230, 1991.

147[172] S. L. Chin, Y. Liang, J. E. Decker, F. A. IIkov, and M. V. Ammosov. Tunnel ionization of
diatomic molecules by an intense CO2 laser. J. Phys. B, 25:L249, 1992.
[173] T. D. G. Walsh, J. E. Decker, and S. L. Chin. Tunnel ionization of simple molecules by anintense CO2 laser. J. Phys. B, 26:L85, 1993.
[174] T. Zuo and A. D. Bandrauk. Charge-resonance-enhanced ionization of diatomic molecularions by intense lasers. Phys. Rev. A, 52:R2511, 1995.
[175] T. Seideman, M. Y. Ivanov, and P. B. Corkum. Role of electron localization in intense-field
molecular ionization. Phys. Rev. Lett., 75:2819, 1995.
[176] W. Li, A. Jaroń-Becker, C. W. Hogle, V. Sharma, X. Zhou, A. Becker, H. C. Kapteyn, andM. M. Murnane. Visualizing electron rearrangement in shape and time during the transition
from a molecule to atoms. Proc. Natl. Acad. Sci., 107:20219, 2010.
[177] E. Lötstedt, T. Kato, and K. Yamanouchi. Intramolecular electron dynamics in the ionizationof acetylene by an intense laser pulse. J. Chem. Phys., 138:104304, 2013.
[178] E. P. Fowe and A. D. Bandrauk. Nonperturbative time-dependent density-functional theoryof ionization and harmonic generation in OCS and CS2 molecules with ultrashort intenselaser pulses: Intensity and orientational effects. Phys. Rev. A, 84:035402, 2011.
[179] W. Chu, H. Li, J. Ni, B. Zeng, J. Yao, H. Zhang, G. Li, C. Jing, H. Xie, H. Xu, K. Yamanouchi,
and Y. Cheng. Lasing action induced by femtosecond laser filamentation in ethanol flame for
combustion diagnosis. App. Phys. Lett., 104(9), 2014.
[180] W. Chu, B. Zeng, J. Yao, H. Xu, J. Ni, G. Li, H. Zhang, F. He, C. Jing, Y. Cheng, and Z. Xu.
Multiwavelength amplified harmonic emissions from carbon dioxide pumped by mid-infrared
femtosecond laser pulses. Europhys. Lett, 97(6):64004, 2012.
[181] T. Streibel and R. Zimmermann. Resonance-enhanced ionization mass spectrometry (REMPIMS): Application for process analysis. Annu. Rev. Anal. Chem., 7:361, 2014.
[182] X. Li, J. C. Tully, H. B. Schlegel, and M. J. Frisch. Ab initio Ehrenfest dynamics. J. Chem.
Phys., 123:084106, 2005.
[183] S. I. Chu and D. A. Telnov. Beyond the Floquet theorem: Generalized Floquet formalisms
and quasienergy methods for atomic and molecular multiphoton processes in intense laser
fields. Phys. Rep., 390:1, 2004.
[184] J.-F. Daigle, A. Jaroń-Becker, S. Hosseini, T.-J. Wang, Y.. Kamali, G. Roy, A. Becker, andS. L. Chin. Intensity clamping measurement of laser filamentation in air at 400 and 800 nm.
Phys. Rev. A, 82:023405, 2010.
[185] J. Yao, B. Zeng, H. Xu., G. Li, W. Chu, J. Ni, H. Zhang, S. L. Chin, Y. Cheng, and Z. Xu.
High-brightness switchable multiwavelength remote laser in air. Phys. Rev. A, 84:051802,
2011.
[186] H. Zhang, C. Jing, G. Liand H. Xie, J. Yao, B. Zeng, W. Chu, J. Ni, H. Xu, and Y. Cheng.
Abnormal dependence of strong-field-ionization-induced nitrogen lasing on polarization ellipticity of the driving field. Phys. Rev. A, 88:063417, 2013.

148[187] J. Ni, W. Chu, C. Jing, H. Zhang, B. Zhang, J. Yao, G. Li, H. Xie, C. Zhang, H. Xu, S. L.
Chin, Y. Cheng, and Z. Xu. Identification of the physical mechanism of generation of coherent
N2+ emissions in air by femtosecond laser excitation. Opt. Exp., 21:8746, 2013.
[188] M. R. Miller, Y. Xia, A. Becker, and A. Jaroń-Becker. Laser-driven nonadiabatic electron
dynamics in molecules. Optica, 3:259, 2016.
[189] H. Y. Tan, Y. Sun, W. Lo, S. J. Lin, C. H. Hsiao, Y. F. Chen, C.-M. Huang, W. C. Lin, S. H.
Jee, H. S. Yu, and C. Y. Dong. Multiphoton fluorescence and second harmonic generation
imaging of the structural alterations in keratoconus ex vivo. Invest. Ophthalmol. Vis. Sci.,
47:5251, 2006.
[190] D. Ding, J. Huang, R. N. Compton, C. E. Klots, and R. E. Haufler. CW laser ionization of
C60 and C70 . Phys. Rev. Lett., 73:1084, 1994.
[191] A. Baev, O. Rubio-Pons, F. Gel’mukhanov, and H. Agren. Optical limiting properties of zincand platinum-based organometallic compounds. J. Phys. Chem. A, 108:7406–7416, 2004.
[192] M. Tsubouchi, B. J. Whitaker, and T. Suzuki. Femtosecond photoelectron imaging on
pyrazine: S1 → T1 intersystem crossing and rotational coherence transfer. J. Phys. Chem.
A, 108:6823–6835, 2004.
[193] H. S. Yoo, M. J. DeWitt, and B. H. Pate. Vibrational dynamics of terminal Acetylenes: II.
Pathway for vibrational relaxation in gas and solution. J. Phys. Chem. A, 108:1365–1379,
2004.
[194] J. R. Barker, L. M. Yoder, and K. D. King. Vibrational energy transfer modeling of nonequilibrium polyatomic reaction systems. J. Phys. Chem. A, 105(5):796–809, 2001.
[195] J. Che, J. L. Krause, M. Messina, K. R. Wilson, and Y. Yan. Detection and control ofmolecular quantum dynamics. J. Phys. Chem., 99:14949–14958, 1995.
[196] A. N. Markevitch, S. M. Smith, D. A. Romanov, H. B. Schlegel, M. Yu. Ivanov, and R. J.
Levis. Nonadiabatic dynamics of polyatomic molecules and ions in strong laser fields. Phys.
Rev. A, 63:011402(R), 2003.
[197] A. N. Markevitch, D. A. Romanov, S. M. Smith, H. B. Schlegel, M. Yu. Ivanov, and R. J.
Levis. Sequential nonadiabatic excitation of large molecules and ions driven by strong laser
field. Phys. Rev. A, 69:013401, 2004.
[198] X. Li, S. M. Smith, A. N. Markevitch, D. A. Romanov, R. J. Levis, and H. B. Schlegel. A timedependent Hartree-Fock approach for studying the electronic optical response of molecules inintense fields. Phys. Chem. Chem. Phys., 7:233, 2005.
[199] S. M. Smith, X. Li, A. N. Markevitch, D. A. Romanov, R. J. Levis, and H. B. Schlegel.
Numerical simulation of nonadiabatic electron excitation in the strong field regime. 2. linear
polyene cations. J. Phys. Chem. A, 109:10527, 2005.
[200] S. M. Smith, X. Li, A. N. Markevitch, D. Romanov, R. J. Levis, and H. B. Schlegel. Numericalsimulation of nonadiabatic electron excitation in the strong field regime. 3. polyacene neutrals
and cations. J. Phys. Chem. A, 111:6920, 2007.

149[201] H. B. Smith, D. A. Romanov, X. Li, J. A. Sonk, H. B. Schlegel, and R. J. Levis. Numericalbound state electron dynamics of carbon dioxide in the strong-field regime. J. Phys. Chem.
A, 114:2576, 2010.
[202] S. Chelkowski, A. Conjusteau, T. Zuo, and A. Bandrauk. Dissociative ionization of H2+ in an
intense laser field: Charge-resonance-enhanced ionization, Coulomb explosion, and harmonic
generation at 600 nm. Phys. Rev. A, 54:3235, 1996.
[203] N. Takemoto and A. Becker. Time-resolved view on charge-resonance-enhanced ionization.
Phys. Rev. A, 84:023401, 2011.
[204] R. Bavli and H. Metiu. Laser-induced localization of an electron in double-well quantum
structure. Phys. Rev. Lett., 69:1986, 1992.
[205] G. L. Kamta and A. D. Bandrauk. Effects of molecular symmetry on enhanced ionization byintense laser pulses. Phys. Rev. A, 75:041401, 2007.
[206] S. I. Chu. Generalized Floquet theoretical apapproach to intense field multiphoton andnonlinear optical processes. Adv. Chem. Phys., 73:739, 1989.
[207] A. Jaroń-Becker. Multiphoton processes in intense laser fields. PhD thesis, Warsaw University,
2001.
[208] L. W. Garland, A. Jaroń, J. Z. Kamiński, and R. M. Potvliege. Off-shell effects in laserassisted electron scattering at low frequency. J. Phys. B: At. Mol. Opt. Phys., 35:2861, 2002.
[209] R. M. Potvliege and R. Shakeshaft. In Atoms in intense laser fields, page 373. Academic
Press, New York, 1992.
[210] A. Jaroń, R. M. Potvliege, and E. Mese. Floquet analysis of laser-induced continuum structures. J. Phys. B: At. Mol. Opt. Phys., 33:1487, 2000.
[211] L. B. Madsen, M. Plummer, and J. F. McCann. Multiphoton ionization of H2+ by intenselight: A comparison of Floquet and wave-packet results. Phys. Rev. A, 58:456, 1998.
[212] P. V. Likhatov and D. A. Telnov. Multiphoton ionization of H2+ at critical internuclear
separation: Non-Hermitian Floquet analysis. J. Phys. B: At. Mol. Opt. Phys., 42:105603,
2009.
[213] T. T. Nguyen-Dang, F. Chateauneuf, O. Atabek, and X. He. Time-resolved dynamics oftwo-channel molecular systems in CW laser fields: Wave-packet construction in the Floquet
formalism. Phys. Rev. A, 51:1387, 1995.
[214] T. T. Nguyen-Dang, C. Lefebvre, H. Abou-Rachid, and O. Atabek. Floquet representation ofabsolute phase and pulse-shape effects on laser-driven molecular photon dissociation. Phys.
Rev. A, 71:023403, 2005.
[215] R. M. Potvliege and R. Shakeshaft. Multiphoton processes in an intense laser field: Harmonic
generation and total ionization rates for atomic hydrogen. Phys. Rev. A, 40:3061, 1989.
[216] D. A. Telnov, J. Wang, and S. I. Chu. Two-color phase control of high-order harmonic
generation in intense laser fields. Phys. Rev. A, 52:3988, 1995.

150[217] K. Drese and M. Hlthaus. Floquet theory for short laser pulses. Eur. Phys. J. D, 5:119, 1999.
[218] D. Telnov and S. I. Chu. Floquet formalism of time-dependent density functional theory.
Chem. Phys. Lett., 264:466, 1997.
[219] P. Salek, T. Helgaker, and T. Saue. Linear response at the 4-component relativistic density
functional level: application to the frequency-dependent dipole polarizability of Hg , AuH andP tH2 . Chem. Phys., 311:187, 2005.
[220] G. Floquet. Sur les équations différentielles linéaires à coefficients périodiques. Ann. I’Ecol.
Norm. Sup., 12:47, 1883.
[221] R. M. Potvliege. STRFLO: a program for time-independent calculation of multiphoton
processes in one-electron atomic systems I. Quasienergy spectra and angular distributions.
Comput. Phys. Commun., 114:42–93, 1998.
[222] C. C. Yu, R. M. Pitzer, and C. W. McCurdy. Molecular photoionization cross sections by thecomplex-basis-function method. Phys. Rev. A, 32:2134, 1985.
[223] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin,
V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F.
Smith, S. Zampini, and H. Zhang. PETSc Web page, 2015.
[224] V. Hernandez, J. E. Roman, and V. Vidal. SLEPc: A scalable and flexible toolkit for thesolution of eigenvalue problems. ACM Trans. Math. Software, 31:351–362, 2005.
[225] J. Aguilar and J. M. Combes. A class of analytic perturbations for one-body Schrödinger
Hamiltonians. Commun. Math. Phys., 22:269–279, 1971.
[226] E. Balslev and J. M. Combes. Spectral properties of many-body Schrödinger operators withdilatation-analytic interactions. Commun. Math. Phys., 22:280–294, 1971.
[227] W. P. Reinhart. Complex coordinates in the theory of atomic and molecular structure anddynamics. Ann. Rev. Phys. Chem., 33:223–255, 1982.
[228] T. N. Rescigno, C. W. McCurdy, and A. E. Orel. Extensions of the complex-coordinate
method to the study of resonances in many-electron systems. Phys. Rev. A, 17:1931–1938,1978.
[229] N. Moiseyev. Non-Hermitian quantum mechanics. Cambridge University Press, 2011.
[230] A. Szaba and N. S. Ostlund. Modern Quantum Chemistry: Introduction to advanced
Quantum Chemistry. Dover Publications, 1989.
[231] Peter M.W. Gill. Molecular Integrals over Gaussian basis functions. Adv. Quantum Chem.,
25:141 – 205, 1994.
[232] M. Galassi, J. Davies, J. Theiler, G. Gough, G. Jungman, P. Alken, M. Booth, and F. Rossi.
GNU Scientific Library Reference Manual, volume 83. 2009.
[233] M. Plummer and J. F. McCann. Complex-basis-function Floquet calculation of multiphoton
ionization rates in hydrogen systems. J. Phys. B, 28:4073–4089, 1995.

151[234] M. Morita and S. Yabushita. Photoionization cross sections of H2+ and H2 with complex
Gaussian-type basis functions optimized for the frequency-dependent polarizabilities. J.
Comput. Chem., 29:2471, 2008.
[235] G. D. Doolen. A procedure for calculating resonance eigenvalue. J. Phys. B: At.. Mol. Phys.,
8:525, 1975.
[236] Y. K. Ho. A resonant state and the ground state of positronium hydride. Phys. Rev. A,
17:1675, 1978.
[237] C. W. McCurdy and F. Martı́n. Implementation of exterior complex scaling in B-splines tosolve atomic and molecular collision problems. J. Phys. B: At. Mol. Opt. Phys., 37:917, 2004.
[238] M. Honigmann, H.-P. Liebermann, and R. J. Buenker. Use of complex configuration interaction calculations and the stationary principle for the description of metastable electronic
states of HCl− . J. Chem. Phys., 133:044305, 2010.
[239] T. N. Rescigno. Atomic photoionization by the complex-basis-function expansion method:
Application to ground-state and metastable Mg . Phys. Rev. A, 31:607, 1985.
[240] R. Matsuzaki, S. Asai, C. W. McCurdy, and S. Yabushita. Construction of complex STO-NG
basis sets by the method of least squares and their applications. Theor. Chem. Acc., 133:1521,
2014.
[241] G. Liuzzi, S. Lucidi, and M. Sciandrone. Sequential penalty derivative-free methods fornonlinear constrained optimization. Siam. J. Optimiz., 20(5):2614–2635, 2010.
[242] R. E. Perez, P. W. Jansen, and J. R. R. A. Martins. PyOpt: A Python-based object-oriented
framework for nonlinear constrained optimization. Struct. Multidiscip. O., 45:101–118, 2012.
[243] S. I. Chu. Intense field multiphoton ionization via complex dressed states: application to the
H atom. Phys. Rev. Lett., 39:1195–1198, 1977.
[244] L. Tao, C. W. McCurdy, and T. N. Rescigno. Grid-based methods for diatomic quantum
scattering problems. III. Double photoionization of molecular hydrogen in prolate spheroidal
coordinates. Phys. Rev. A, 82:023423, 2010.
[245] N. Moiseyev, M Chrysos, and R Lefebvre. The solution of the time-dependent Schrödinger
equation by the (t,t’) method: application to intense field molecular photodissociation. J.
Phys. B: At. Mol. Opt. Phys., 28:2599, 1995.
[246] E. J. Baerends, O. V. Gritsenko, and R. Van Meer. The KohnSham gap, the fundamental
gap and the optical gap: the physical meaning of occupied and virtual KohnSham orbitalenergies. Phys. CHem. Chem. Phys., 15:16408–16425, 2013.
[247] S. Baker, J. S. Robinson, C. A. Haworth, H. Teng, R. A. Smith, C. C. Chirilǎ, M. Lein,
J. W. G. Tisch, and J. P. Marangos. Probing proton dynamics in molecules on an attosecond
time scale. Science, 312:424–427, 2006.
[248] P. Balanarayan and N. Moiseyev. Linear Stark effect for sulfur atom in strong high-frequency
laser fields. Phys. Rev. Lettt., 110:253001, 2013.

152[249] A. Baltuška, Th. Udem, M. Uiberacker, M. Hentschel, E. Goulielmakis, Ch. Gohle,
R. Holzwarth, V. S. Yakovlev, A. Scrinzi, T. W. Hänsch, and F. Krausz. Attosecond controlof electronic processes by intense light fields. Nature, 421:611–615, 2003.
[250] A. D. Bandrauk and F. Légaré. Progress in Ultrafast Intense Laser Science VIII, chapter
Chap. 2. Springer, NY, 2012.
[251] A. D. Bandrauk. Molecules in Laser Fields, chapter Chap. 4. Marcel Dekker Inc., NY, 1994.
[252] D. Bauer, D. B. Milosevic, and W. Becker. Strong-field approximation for intense-laseratom
processes: The choice of gauge. Phys. Rev. A, 72:023415, 2005.
[253] A. D. Becke. On the large-gradient behavior of the density functional exchange energy. J.
Chem. Phys., 85:7184, 1986.
[254] R. V. Leeuwen and E. J. Baerends. Exchange-correlation potential with correct asymptotic
behavior. Phys. Rev. A, 49:2421, 1994.
[255] N. H. Burnett, H. A. Baldis, M. C. Richardson, and G. D. Enright. Harmonic generation in
CO2 laser target interaction. Appl. Phys. Lett., 31:172, 1977.
[256] R. L. Carman, D. W. Forlund, and J. M. Kindel. Visible harmonic emission as a way ofmeasuring profile steepening. Phys. Rev. Lett., 46:29, 1981.
[257] R. L. Carman, C. K. Rhodes, and R. F. Benjamin. Observation of harmonics in the visible
and ultraviolet created in CO2 -laser-produced plasmas. Phys. Rev. A, 24:2649, 1981.
[258] Chan L., M. Siemens, M. M. Murnane, H. C. Kapteyn, S. Mathias, M. Aeschlimann, P. Grychtol, R. Adam, C. M. Schneider, J. M. Shaw, H. Nembach, and T. J. Silva. Ultrafast demagnetization dynamics at the M edges of magnetic elements observed using a tabletop
high-harmonic soft X-ray source. Phys. Rev. Lett., 103:257402, 2009.
[259] P. B. Corkum, N. H. Burnett, and F. Brunel. Above-threshold ionization in the longwavelength limit. Phys. Rev. Lett., 62:1259, 1989.
[260] M. Drescher, M. Hentschel, R. Kienberger, M. Uiberacker, V. Yakovlev, A. Scrini, Th. Westerwalbeslon, U. Kleineberg, U. Heinzmann, and F. Krausz. Time-resolved atomic inner-shell
spectroscopy. Nature, 419:803–807, 2002.
[261] D. J. Fraser, MHR Hutchinson, JP Marangos, YL Shao, JWG Tisch, and M. Castillejo. Highharmonic generation in butane and butadiene. J. Phys. B: At. Mol. Opt. Phys., 28:L739–L745,
1995.
[262] E. Gagnon, P. Ranitovic, X. M. Tong, C. L. Cocke, M. M. Murnane, H. C. Kapteyn, and
A. S. Sandhu. Soft X-ray-driven femtosecond molecular dynamics. Science, 317:1374–1378,
2007.
[263] P. Jemmer and P. J. Knowles. Exchange energy in Kohn-Sham density-functional theory.
Phys. Rev. A, 51:3571, 1995.
[264] F. Krausz and M. Ivanov. Attosecond physics. Rev. Mod. Phys., 81:163, 2009.

153[265] W. Li, V. Sharma, C. Hogle, A. Jaroń-Becker, A. Becker, H. Kapteyn, and M. Murnane.
Ultrafast molecular dynamics probed using strong field ionization. In Second International
Conference on Attosecond Physics, 2009.
[266] Z. H. Loh, M. Khalil, R. E. Correa, R. Santra, C. Buth, and S. R. Leone. Quantum stateresolved probing of strong-field-ionized xenon atoms using femtosecond high-order harmonictransient absorption spectroscopy. Phys. Rev. Lett., 98:143601, 2007.
[267] E. A. McLean, J. A. Stamper, B. H. Ripin, H. R. Griem, J. M. McMahon, and S. E. Bodner.
Harmonic generation in Nd:laser-produced plasmas. Appl. Phys. Lett., 31:825, 1977.
[268] L. Miaja-Avila, C. Lei, M. Aeschlimann, J. L. Gland, M. M. Murnane, H. C. Kapteyn., and
G. Saathoff. Laser-assisted photoelectric effect from surfaces. Phys. Rev. Lett., 97:113604,
2006.
[269] R. S. Mulliken. Intensities of electronic transitions in molecular spectra II. Charge-transfer
spectra. J. Chem. Phys., 7:20, 1939.
[270] G. Saathoff, L. Miaja-Avila, M. Aeschlimann, M. M. Murnane, and H. C. Kapteyn. Laserassisted photoemission from surfaces. Phys. Rev. A, 77:022903, 2008.
[271] A. Savin, C. J. Umrigar, and X. Gonze. Relationship of Kohn-Sham eigenvalues to excitationenergies. Chem. Phys. Lett., 288:391, 1998.
[272] J. M. Schins, P. Breger, P. Agostini, R. C. Constantinescu, H. G. Muller, G. Grillon, A. Antonetti, and A. Mysyrowicz. Observation of laser-assisted auger decay in Argon. Phys. Rev.
Lett., 73:2180–2183, 1994.
[273] A. Srinzi, M. Y. Ivanov, R. Kienberger, and D. M. Villeneuve. Attosecond physics. J. Phys.
B, 39:R1, 2006.
[274] M. E. Siemens, Q. Li, R. Yang, K. A. Nelson, E. H. Anderson, M. M. Murnane, and H. C.
Kapteyn. Quasi-ballistic thermal transport from nanoscale interfaces observed using ultrafast
coherent soft X-ray beams. Nature, 9:26–30, 2010.
[275] R. Stowasser and R. Hoffmann. What do the Kohn-Sham orbitals and eigenvalues mean? J.
Am. Chem. Soc., 121:3414, 1999.
[276] A. V. Titov, A. O. Mitrushenkov, and I. I. Tupitsyn. Effective core potential for pseudoorbitals with nodes. Chem. Phys. Lett., 185:330, 1991.
[277] X. M. Tong, Z. X. Zhao, and C. D. Lin. Theory of molecular tunnel ionization. Phys. Rev.
A, 66:033402, 2002.
[278] R. Torres, N. Kajumba, J. G. Underwood, J. S. Robinson, S. Baker, J. W. Gtisch, R. de Nalda,
W. A. Bryan, R. Velotta, C. Altucci, I. C. E. Turcu, and J. P. Marangos. Probing orbital
structure of polyatomic molecules by high-order harmonic generation. Phys. Rev. Lett.,
98:203007, 2007.
[279] D. J. Tozer and N. C. Handy. Improving virtual Kohn-Sham orbitals and eigenvalues: Application to excitation energies and static polarizabilities. J. Chem. Phys., 109:10180, 1998.

154[280] A. J. Verhaef. Optical detection of tunneling ionization. Phys. Rev. Lett., 104:163904, 2010.
[281] J. Yao, G. Li, C. Jing, B. Zeng, W. Chu, J. Ni, H. Zhang, H. Xie, C. Zhang, H. Li, H. Xu,
S. L. Chin, Y. Cheng, and Z. Xu. Remote creation of coherent emissions in air with two-color
ultrafast laser pulses. New J. Phys., 15(2):023046, 2013.
[282] B. Zhang, W. Chu, G. Liand J. Yao, H. Zhang, J. Ni, C. Jing, H. Xie, and Y. Cheng. Real-timeobservation of dynamics in rotational molecular wave packets by use of air-laser spectroscopy.
Phys. Rev. A, 89:042508, 2014.
[283] T. Zuo, S. Chelkowski, and A. D. Bandrauk. Harmonic generation by the H2+ molecular ion
in intense laser fields. Phys. Rev., 48:3837, 1993.