Theoretical Studies of Ultrafast Correlated Electron
Dynamics in Single and Double Photoionizationby
Hongcheng Ni
B.S., Nanjing University, 2008
M.S., University of Colorado 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
2014This thesis entitled:
Theoretical Studies of Ultrafast Correlated Electron Dynamics in Single and Double
Photoionization
written by Hongcheng Ni
has been approved for the Department of Physics

Prof. Andreas Becker

Prof. Margaret Murnane

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
Ni, Hongcheng (Ph.D., Physics)
Theoretical Studies of Ultrafast Correlated Electron Dynamics in Single and Double Photoionization
Thesis directed by Prof. Andreas Becker

Recent advances in laser technology have led to the generation of attosecond laser pulses,
whose duration is in the range of the natural time scale of electron dynamics, and thus make theobservation and even control of electron dynamics in atoms and molecules possible. While thesingle-electron dynamics is more thoroughly studied, the correlated dynamics of two electrons isless understood, especially in the context of resolving the ultrafast temporal information in double
photoionization.
In this thesis, we first study the energy exchange via electron correlation upon photon absorption over large distances in the two-site double photoionization of the helium dimer, which isfound to be a two-step process. In the first step, one electron in one atom absorbs the photon andgets ionized. In the second step, this electron propagates towards the neighboring atom and knocks
out the other electron. We then introduce the Hamiltonian reduction method to further study theeffects of different interactions in the single and double photoionization of the helium dimer.
Next, we analyze the selection rules for the emission of two electrons from the helium atom,
the helium dimer, and general molecules following the absorption of a few photons in an intenselaser field. In particular, the back-to-back emission of the two electrons with equal energy sharing
is either suppressed or not depending on the number of photons absorbed from the field.
Finally, we study the time delay between the single and double photoionization processes.
We first propose a self-consistent-time method to account for the Coulomb-laser coupling effect andobtain the intrinsic photoabsorption time delay measured by the attosecond streak camera. We
then proceed to time resolve the correlated emission of two electrons in the knockout process of thehelium dimer with respect to the first step of single ionization.

Dedication

To my love.

Acknowledgements

Working at JILA and the University of Colorado Boulder has been a great pleasure, especiallybecause of the knowledgeable and friendly people here, without whom my journey through thegraduate school would not have been possible. Here, I would like to thank every one who has lentme kind help, support, and encouragement.
First of all, I am in heavy debt to my advisor, Prof. Andreas Becker. Doing research is likesailing on the vast ocean, which would not be possible if there is not a beacon. Andreas is just likesuch a beacon for me. He has been a borderless source of creative ideas and physical intuitions thathas always been inspiring me. I remember all the times when I come up with some messy (to me)
figures, from which Andreas is able to disentangle the disordered parts and discern the relevant
physics. Watching him drawing figures and doing all the analyses in front of me have be a gracious
treat, from which I learn a lot. What also impresses me a lot is his logic flow in presentations andI can never learn enough from him on these skills. That said, Andreas is a very nice and easygoing
person in life on the other hand. I would especially thank his warm support for me when I wastrying to look for a postdoctoral position. He had been kindly helping me contacting people heknow and giving me suggestions.
Next, I would like to gratefully thank Dr. Agnieszka Jaroń-Becker, my coadvisor. Agnieszka
has a very broad knowledge in both physics and chemistry, from whom I have gained a lot. Whenever we encounter some unfamiliar phenomena, she would always been able to connect that to aphysical process. Her insights have been helping me to digest and build a more comprehensive viewof the physical world. Besides that, she is extremely warm-hearted and is ready to help at any

vitime. When I was looking for the postdoctoral position, she had been looking for every possibility
for me and have talked to me about the research and life in different places for hours and hours. I
could never thank her enough.
I would also give my gratitude to my previous colleagues, Dr. Shaohao Chen, Dr. Antonio
Picón, and Dr. Norio Takemoto, with whom I had enjoyed great benefit when I started my researchin the group. It would not have been possible for me to dig into all the computer codes, numericalmethods, and physical concepts without their help. I would especially thank Shaohao, who had been
disturbed by me for afternoons and afternoons. He had been always patient and led me through the
problems step by step. Antonio had be particularly helpful in sharing cutting papers and explaining
new physical concepts to me. He had also been constantly helping me when I was looking for the
postdoctoral position. In addition, I would gratefully thank all my current colleagues, Jing Su, Dr.
Carlos Hernández-Garcı́a, Yuqing Xia, Daniel Weflen, Andrew Spott, Michelle Miller, and Cory
Goldsmith, for creating such an active and pleasant atmosphere. Working in the Becker group has
helped me shaping my early career that would lighten my future and benefit me for life.
I would further express my appreciation to my collaborators, Dr. Shaohao Chen, Dr. Camilo
Ruiz, and Prof. Reinhard Dörner, without whom most of my work would not have been possible.
It has also been a delight to work with them. I would especially thank Prof. Dörner, who lent mewarm help in writing a very positive recommendation letter, without which it would not have beenpossible for me to get a postdoctoral position.
Moreover, I would like to thank Prof. Margaret Murnane, Prof. Steven Cundiff, Prof. John
Bohn, and Prof. David Jonas for agreeing to be in my dissertation committee. Their insightful
feedback has helped me to shape my work and thesis better. I would especially thank Prof. Murnane
and Prof. Bohn to proof read my thesis.
The JILA computing center and the Janus supercomputing center have also been continuously
supporting me. It would not have been possible to carry out numerical experiments without theirsupport.
Finally and most importantly, I would like to thank my family and all friends that have been

vii
continuously supporting me and encouraging me for years. I have been always on the go and awayfrom home, chasing my dream since high school, which makes me feel guilty. However, the support,
encouragement, and love I get from my family have never ceased or receded. The life in a foreigncountry has also been difficult. Fortunately, I met a lot of friends here that make me feel unlonely.
Especially luckily, I met my lovely wife, Jin Wen, here. She has been the sunshine in my life andhas been backing me up with her endless love to help me go through the rough and bumpy roads
and fight against all the frustrations.
Thank you all, sincerely.

viiiContents

Chapter
1 Introduction

1

1.1

Strong-field single ionization and the Keldysh parameter . . . . . . . . . . . . . . . .

1

1.2

Strong-field double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3

High-order harmonic generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.4

Generation of single isolated attosecond pulse . . . . . . . . . . . . . . . . . . . . . .

7

1.5

Attosecond streak camera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.6

Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 Numerical methods

16

2.1

The split-operator method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.2

The Crank–Nicolson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.3

Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4

The imaginary-time-propagation method . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.5

The spectral method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.6

Absorbing boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.6.1

The mask-function method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.6.2

The exterior-complex-scaling method . . . . . . . . . . . . . . . . . . . . . . . 27

2.7

The space-partitioning method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.8

Calculating single and double-ionization yields . . . . . . . . . . . . . . . . . . . . . 32

2.9

Parallel computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

ix
3 Single and double photoionization of the helium dimer

36

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.2

Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2.1

The planar 4D model of the helium dimer . . . . . . . . . . . . . . . . . . . . 39

3.2.2

Initial states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.3

Photoelectron angular distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.4

The knockout mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.5

The Hamiltonian reduction method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.6

3.5.1

Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.5.2

Single ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.5.3

Double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.5.4

Exchange ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4 Selection rules in the few-photon double ionization process

68

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.2

Selection rules for transitions to two-electron continuum states . . . . . . . . . . . . 69

4.3

4.4

4.2.1

Selection rules based on general studies of the two-electron continuum . . . . 70

4.2.2

Selection rules based on the symmetry of the transition matrix . . . . . . . . 72

4.2.3

Selection rules based on the correlated Keldysh–Faisal–Reiss theory . . . . . 73

Selection rules in the helium atom . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.3.1

The restricted 3D model of the helium atom . . . . . . . . . . . . . . . . . . . 75

4.3.2

One-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.3.3

Two-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.3.4

Three- and four-photon double ionization . . . . . . . . . . . . . . . . . . . . 84

4.3.5

Momentum distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

Selection rules in the helium atom within the planar 4D model . . . . . . . . . . . . 89

x

4.5

4.4.1

The planar 4D model of the helium atom . . . . . . . . . . . . . . . . . . . . 92

4.4.2

One-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

4.4.3

Two-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

4.4.4

Momentum distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Selection rules in the helium dimer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.5.1

One-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

4.5.2

Two-photon double ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

4.5.3

Momentum distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

4.6

General molecules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

4.7

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

5 Photoabsorption and correlation time delay

110

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.2

Photoabsorption time delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.3

5.4

5.2.1

Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.2.2

The self-consistent-time method . . . . . . . . . . . . . . . . . . . . . . . . . 113

Correlation time delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.3.1

Correlation time delay in the helium dimer . . . . . . . . . . . . . . . . . . . 124

5.3.2

Correlation time delay in the planar 4D helium dimer . . . . . . . . . . . . . 126

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

6 Summary and perspective

134

Bibliography

137

xi

TablesTable
3.1

Double ionization probabilities for the different momentum-space wave functions . . 46

5.1

Two-photon ionization time delay for Li2+ . . . . . . . . . . . . . . . . . . . . . . . . 123

5.2

Correlation time delay for different R and λXUV . . . . . . . . . . . . . . . . . . . . . 128

xii

FiguresFigure1.1

Illustration of the Keldysh theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2

The three step model of the HHG process and electron trajectories . . . . . . . . . .

6

1.3

Generation of an attosecond pulse train . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4

The cut-off harmonics method to retrieve a single isolated attosecond pulse . . . . .

8

1.5

The cut-off harmonics method fails for CEP=±π/2 . . . . . . . . . . . . . . . . . . .

9

1.6

Fourier transform of cut-off pulses above cut-off energy . . . . . . . . . . . . . . . . .

91.7(a) Scheme of the conventional streak camera, (b) scheme of the attosecond streak
camera, and (c) principle of the attosecond streak camera. . . . . . . . . . . . . . . . 12

2.1

The spectrum of the Cylindrical 2D hydrogen atom. . . . . . . . . . . . . . . . . . . 24

2.2

The scheme of the absorbing-boundary method. . . . . . . . . . . . . . . . . . . . . . 25

2.3

The typical shape of the mask function. . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4

The space partitioning for an atom.

2.5

The space partitioning for a diatomic system (atomic nature). . . . . . . . . . . . . . 29

2.6

The space partitioning for a diatomic system (molecular nature). . . . . . . . . . . . 30

2.7

The schematic illustration of the way to choose the bounds in the space-partitioning. . . . . . . . . . . . . . . . . . . . . . . . . . . 28

method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1

Scheme of the numerical model for the helium dimer. . . . . . . . . . . . . . . . . . . 40

xiii
3.2

Spatial distributions of the ground state and first excited state of the planar model
helium dimer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.3

Molecular-frame angular distributions of one of the two electrons following double
photoionization of the helium dimer, obtained using different methods. . . . . . . . . 47

3.4

Molecular-frame angular distributions of one of the two electrons following double
photoionization of the helium dimer, with singlet and triplet contributions. . . . . . 48

3.5

Orientation dependence of the double ionization probability of the helium dimer. . . 50

3.6Internuclear-distance dependence of the double ionization probability of the helium
dimer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.7

Temporal evolution of the two-electron probability distribution of the helium dimer.

52

3.8

Scheme of the numerical model for the helium dimer in the case of parallel orientation. 54

3.9

Spatial distributions of the ground state and first excited state of the planar model
helium dimer with different Hamiltonians. . . . . . . . . . . . . . . . . . . . . . . . . 56

3.10 Comparison of photoelectron angular distributions in single ionization of the helium
dimer with different Hamiltonians. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.11 Angular distributions of the primary and knockout electron in double photoionizationof the helium dimer with different Hamiltonians. . . . . . . . . . . . . . . . . . . . . 61
3.12 Angular distributions of the knockout electron in double photoionization of the helium dimer with different Hamiltonians when emission direction of the primary electron is fixed.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.13 Angular distributions of the knockout electron in double photoionization of the helium dimer with different Hamiltonians when emission direction of the primary electron is fixed.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.14 Exchange ionization under electron impact. . . . . . . . . . . . . . . . . . . . . . . . 65
4.1

The restricted 3D model of the helium atom. . . . . . . . . . . . . . . . . . . . . . . 76

4.2

The photoabsorption scheme in the restricted 3D model of the helium atom. . . . . . 77

xiv
4.3

Probability density distribution following single-photon double ionization. . . . . . . 78

4.4

Probability density distribution following single-photon double ionization. . . . . . . 79

4.5

Probability density distribution following single- and two-photon double ionization. . 81

4.6

Probability density distribution following single- and two-photon double ionization. . 82

4.7

Probability density distribution following two-photon double ionization. . . . . . . . 83

4.8

Probability density distribution following three- and four-photon double ionization. . 85

4.9

Correlated momentum distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.10 Power law. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.11 Correlated momentum distribution obtained by Fourier transform. . . . . . . . . . . 90
4.12 Comparison of correlated momentum distributions cut at pρ = 0. . . . . . . . . . . . 91
4.13 Probability density distribution following single-photon double ionization. . . . . . . 94
4.14 Probability density distribution following two-photon double ionization. . . . . . . . 94
4.15 Correlated momentum distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.16 Momentum distributions of electron 2 for single-photon double ionization when theemission direction of electron 1 is fixed in the z direction. . . . . . . . . . . . . . . . 97
4.17 Probability density distribution following single-photon double ionization. . . . . . . 99
4.18 Probability density distribution following two-photon double ionization. . . . . . . . 102
4.19 Correlated momentum distribution for single-photon double ionization. . . . . . . . . 103
4.20 Correlated momentum distribution for two-photon double ionization. . . . . . . . . . 104
5.1

Flowchart of the self-consistent-time (SCT) method. . . . . . . . . . . . . . . . . . . 116

5.2

Iterations in the self-consistent-time method. . . . . . . . . . . . . . . . . . . . . . . 118

5.3

Convergence of the self-consistent-time method. . . . . . . . . . . . . . . . . . . . . . 119

5.4

Robustness of the self-consistent-time method. . . . . . . . . . . . . . . . . . . . . . 119

5.5

Converged photoabsorption time delays when the photoemission takes places at different positions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

5.6

Scheme of the correlation time delay. . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

xv
5.7

Streaking signals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

5.8

The (a) single ionization yield and (b) single ionization rate, as a function of propagation time for different dbig . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

5.9

The single and double ionization time as a function of the partitioning bounds dbig
and dsmall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

Chapter 1

Introduction

The last 30 years have been a prosperous period of time in ultrafast science thanks to theavailability of femtosecond (1 fs = 10−15 s) laser technology, most notably the Ti:sapphire lasers [1],which are capable of producing intense femtosecond laser pulses tunable in the red to near-infrared(NIR) range with intensities exceeding 1015 W/cm2 in the focal spot. The real-time observation andcontrol of atom motions, as well as energy exchange between atoms in molecules and condensed
matter systems have been achieved (refer to Ref. [2] for a review of femtochemistry). Recent
advances in laser technology have led to the generation of attosecond (1 as = 10−18 s) extreme
ultraviolet (XUV) laser pulses. The duration of these pulses is in the range of the natural timescale of electron dynamics, and thus make the observation and even control of electron dynamics
in atoms and molecules possible (refer to Ref. [3] for a recent review of attosecond physics). In thischapter, we give a brief introduction to the background and some basic concepts of attosecond andstrong-field physics and outline the overall structure of this thesis.

1.1

Strong-field single ionization and the Keldysh parameter
The strong-field single ionization of atoms and molecules can generally be divided into three

categories, i.e., multiphoton ionization, tunneling ionization, and above-barrier ionization, withthe help of the Keldysh parameter γ [4]. The strong-field ionization consists of a subtle relation
between photon energy and photon number (or field intensity). For a high photon energy andrelatively low field intensity, the ionization proceeds in a way that the electron absorbs one or a few

2energetic photons and the interaction with the laser field can be considered as perturbative. For
a low photon energy and high field intensity, in contrast, the ionization occurs in a manner thatthe electron absorbs many photons and the interaction is nonperturbative. In the limiting case ofa large number of absorbed photons, the latter process can be considered as tunneling ionization,
which is illustrated in Fig. 1.1. The strong laser field suppresses the Coulomb field of the nucleus
at one side and results in a barrier, through which the bound electron may tunnel out, leading toionization.

x0
Ip

Figure 1.1. Illustration of the Keldysh theory [4].

To clearly understand the relation between the multiphoton and tunneling ionization, we
examine the time scale of the quasi-static tunneling ionization [4]. As shown in Fig. 1.1, the strong
laser field produces a tunneling exit x0 for the electron with bound energy Ip , thus (Hartree atomic
units, e = m = ~ = 1, are used throughout the thesis, unless stated otherwise)− Ip = −Z− Ex0 ≈ −Ex0 ,
x0

(1.1)

where Z is the nucleus charge and E is the field amplitude. Thus
x0 ≈

Ip
.
E

(1.2)

The velocity of the electron before tunneling can be estimated by
ve =

p

2Ee =

p

2Ip ,

(1.3)

where the kinetic energy of the electron Ee = Ip is a result of the virial theorem. Therefore, the3tunneling time can be estimated as
Ip
x0
τ=
≈
ve
E

s

1
1
=
2Ip
2ω

s

Ip
,
2Up

(1.4)

where

Up =

E
2ω

2
(1.5)

is the ponderomotive energy of the electron (cycle-averaged quiver energy of the electron in a laser
field) and ω is the angular frequency of the laser field. On the other hand, one half period of theqIp1
oscillating field is T2 ≈ 2ω
. If τ > T2 , or
2Up > 1, the field oscillation has already reversed itsqI
direction before the electron can tunnel; while if τ < T2 , or 2Upp < 1, the electron may tunnel out
through the barrier. Thus, we can define the Keldysh parameter [4]
s
Ip
γ=
2Up

(1.6)

as a criteria to categorize the strong-field ionization into multiphoton ionization when γ  1 andtunneling ionization when γ  1. In cases when γ ≈ 1, the two mechanisms coexist. The abovebarrier ionization, on the other hand, is just an extreme case of the tunneling ionization, where thelaser field is so strong and suppresses the Coulomb potential so much that the energy level of theinitially bound electron is already above the tunneling barrier.

1.2

Strong-field double ionization
Single ionization in strong fields, as discussed above, can be categorized into multiphoton

ionization (γ  1) and tunneling ionization (γ  1). In this section, we introduce the different
mechanisms in strong-field double ionization, which can generally be divided into three categories [5,6], i.e., sequential double ionization (SDI), nonsequential double ionization (NSDI), and recollision
excitation with subsequent field ionization (RESI). These mechanisms distinguish themselves, inparticular, in the role of the electron correlation.
In SDI, as the name implies, the two electrons are emitted one after the other via interaction
with the field, without much influence on the emission of each other. Thus, it can be simply

4
regarded as two sequential single ionization processes. SDI can occur in both multiphoton andtunneling regimes.
NSDI, on the other hand, is more complex, since the emission of the two electrons is correlated.
In the multiphoton regime, the photon energy can be regarded as being deposited to one of theelectrons initially, which interacts and exchanges energy with the other electron via their Coulomb
correlation and both electrons are emitted in the end. Two different mechanisms are possible here,i.e., knockout [7, 8] and shakeoff [9]. The knockout mechanism applies to the case when the photon
has a relatively low energy and the first electron (primary electron) has a low kinetic energy after
ionization. On its path of emission from the system, however, this electron may knock out theother electron from the system through their mutual interaction. The shakeoff mechanism, on theother hand, applies to the cases when the photon has a high energy and the first electron is stripedfrom the system rapidly before the rest of the system, including the other electron, can adjust to
this change. Now this system needs to relax to the eigenstates of the new Hamiltonian with one
electron removed. Projecting the wave function onto the eigenstates of the new system, with someprobability the electron will be found in the continuum, leading to double ionization. The electron
correlation in the shakeoff process is reflected in the initial two-electron state, which prepares thesystem in a highly correlated fashion. It is also the correlation in the initial state that keeps theenergy conserved for the emission of the second electron. The separation of these two mechanisms
for double photoionization of the helium atom is around an excess electron kinetic energy of 350
eV [10].
NSDI in the tunneling regime is equally interesting. After the first electron tunnels out andgets accelerated first away from the parent ion, it may be driven back by the same laser field when
the field reverses its direction. During the so-called rescattering [11, 12], the electron may collidewith the parent ion and knock another electron out, resulting in double ionization.
RESI, at last, is similar to NSDI in the tunneling regime, but when the rescattered electron
collides with the parent ion, it does not knock another electron out, but promotes it to an excitedstate, from which it can easily tunnel out in the laser field. RESI is considered as belonging to the5third category since the electron correlation is necessary for the excitation of the second electron,
while the subsequent ionization of the second electron may not directly relate to the first electron.
We can summarize all these mechanisms as a diagram below, using the double ionization of
the helium atom as an example:





multiphotonHe+ + e− −−−−−−−→ He2+ + e− + e− (SDI)













multiphoton




−
−
−
−
−
−
−
→

 Knockout



2+

−


He
+
2e
(NSDI)














 Shakeoff

He



tunneling−−−−−−→ He2+ + e− + e− (SDI)













tunneling

+

−


−−−−−−→ He + e He2+ + 2e− (NSDI)
rescattering −−−−−−−→








tunneling He∗+ + e− −−−−−−→ He2+ + e− + e−

1.3

(RESI)

High-order harmonic generation
In NSDI in the tunneling regime, the rescattered electron collides with the parent ion andknocks out another electron. However, instead of knocking out a second electron, the returning
electron may instead recombine to the ground state and emit a photon. This process is called thehigh-order harmonic generation (HHG) [13–15] since one single energetic photon is emitted in thisrecombination process after the electron absorbs many photons and gains a lot of energy from theoscillating laser field.
The HHG process can be described by a three-step model (see Fig. 1.2) [12]: 1. Tunneling:tunneling ionization of electron through the potential barrier suppressed by the intense laser field.
2. Acceleration: electron driven away from and then back to the ion core by the oscillating laser
field. 3. Recombination: returning electron collides with the ion core with some energy and emitsa photon. Or, in the language of wave mechanics, the returning wave packet, together with theparent wave packet, form a rapidly changing dipole moment, and thus emit light radiation.
The tunneling probability of an electron through the suppressed potential barrier is exponen-

6

Figure 1.2. The three step model of the HHG process and electron trajectories (image credit Ref.
[16]).

Figure 1.3. Generation of an attosecond pulse train (image credit Ref. [17]).

7
tially proportional to the laser field magnitude. However, if an electron is released at exactly a local
maximum of the laser field, it returns to the ion core with zero energy, a temporally uniform field
granted. If it is released before the maxima, it will drift away, not being able to return at all. Whenan electron is however released after the maxima, it will then return with some energy, collide with
the ion core, and free the energy in the form of a photon. If it is released at approximately 18◦
after the field extrema, the return energy reaches its maximum, at 3.17Up . Therefore, the cut-offenergy of the photon emitted is Emax = Ip + 3.17Up [12, 18].
For each photon energy lower than the cut-off energy, there must be two contributions to the
HHG spectrum, one for the electron release before 18◦ which results in the so-called long trajectory,
and another for the electron release after 18◦ which leads to the shorter trajectory, as illustrated inFig. 1.2. Since the photon energy spans broadly, from 0 to the cut-off energy, the pulses generated
by the superposition of a broad range of harmonics are very short, indeed of attosecond duration.
However, the three-step process can occur every half-cycle of the laser field, since the fieldreaches a maximum every half-cycle (see Fig. 1.3). Therefore, for each half-cycle of the laser field,
there will be an attosecond pulse generated. That is why in a long driving pulse an attosecond pulse
train is generated much more easily than an isolated attosecond pulse. This half-cycle periodicity
also ensures that there is only odd-order harmonics in the HHG spectrum (if it is generated fromrare-gas atoms). Moreover, since the driving laser field changes in magnitude over time, the cut-offenergy for each half-cycle also changes, i.e., there are different cut-off energies for each half-cycleof the laser field.

1.4

Generation of single isolated attosecond pulse
The attosecond pulse train, while suitable for a number of experiments, is not ideal for time-

resolved experiments since there is an ambiguity in time due to the interaction with multiple pulses
in the pulse train.
Proposed by Christov, Murnane, and Kapteyn in 1997 [19, 20] and experimentally implemented by Krausz and coworkers in 2001 [21–24], the cut-off harmonics method is a straightforward

8
way to retrieve a single isolated attosecond pulse. Although there are different cut-off energies foreach half-cycle of the laser field, the overall cut-off energy occurs only near the envelope peak of
the laser pulse, as is also shown in Fig. 1.3.

Figure 1.4. The cut-off harmonics method to retrieve a single isolated attosecond pulse (image
credit Ref. [25]).

As illustrated in Fig. 1.4, the idea is to reduce the femtosecond driving laser pulse width
to only about 5 cycles which results in the generation of only 3 attosecond bursts at significant
energy. There is however only one burst, which contributes to the overall cut-off energy. Usinga high pass filter (a thin metal film) to filter out the lower harmonics, one is left with the mostenergetic attosecond burst.
The electric field of a pulse can be expressed as E(t) = E0 (t) cos(ωt + ϕ), where E0 (t) is
the pulse envelope and ϕ is the carrier-envelope phase (CEP). To generate stable single attosecond
pulses, it is necessary that the waveform of the driving laser field is stable from pulse to pulse, i.e.,
the CEP should be stable.
In fact, this method not only requires a stable CEP, but works optimally only for ϕ = 0
and fails for ϕ = ±π/2. As shown in Fig. 1.5, two attosecond pulses almost identical in energy

9Figure 1.5. The cut-off harmonics method fails for CEP=±π/2 (image credit Ref. [24]).

Figure 1.6. Fourier transform of cut-off pulses above cut-off energy (image credit Ref. [23]).

10
are generated at ϕ = −π/2, while only one is generated at ϕ = 0. Filtering out components withenergy lower than the overall cut-off energy still retains two pulses at ϕ = ±π/2. This can also
be seen from the Fourier transform of the cut-off pulses in these two cases, as shown in Fig. 1.6.
No matter how the high pass filter is set, the spectrum left is still oscillating, which is still two
attosecond pulses when transformed back.
As outlined, the technique requires the generation of few-cycle laser pulses, which is technologically still challenging. Another disadvantage of this method is that it has a low efficiency. Theenergy remained after the filter is applied is only about 10−5 of the original energy [26].
Over the years, a number of methods have been proposed and implemented to generate
single isolated attosecond pulses, including the polarization gating method [27–31], the two-color
laser method [27, 32], the surface harmonic generation method [33–35], the spatial filtering method[36, 37], and the phase-matching gating method [38–40].

1.5

Attosecond streak camera
With the attosecond XUV laser pulse available, the first question to ask is how to makeuse of it to resolve attosecond electron dynamics. In the traditional femtosecond spectroscopy,
the pump-probe scheme is usually applied, where the pump femtosecond laser pulse initiates theultrafast dynamics and the probe femtosecond laser pulse measures the change in a physical variableinduced by the pump pulse. By varying the relative delay between the pump and probe pulses in acontrolled manner, the ultrafast femtosecond dynamics of processes in atoms, molecules, or solidstate systems can be observed or even manipulated. Although attosecond XUV laser pulses are
currently available, the pump-probe experiments with these pulses are still difficult up to now dueto the low interaction cross section of matter with XUV laser pulses and the low repetition rate of
the attosecond XUV laser pulse at nanojoule pulse energies [41]. However, alternative techniques
like the attosecond streak camera [42] and the attoclock [43] have been used in observing andcontrolling electronic motion with a temporal resolution of less than 100 as, e.g., the real-timeobservation of tunneling ionization in atoms [43, 44], photoelectron emission from solids [45], and11the tracing of electron dynamics in atoms and molecules [46, 47].
To understand how the attosecond streak camera technique works, we first introduce theconventional streak camera [49, 50], as shown in Fig. 1.7(a). The conventional streak camera is
traditionally used to characterize the temporal profile of a short light pulse. Its basic principle isfirst to convert the light signal into an electron signal, which temporal profile is mapped onto ameasurable spatial profile as follows. As shown in the figure, a metal (gray plate) is illuminated by
the short light pulse (cyan curve) with an unknown temporal profile, and electrons (purple dots) get
emitted from the metal surface due to the photoelectric effect. The stronger the light pulse is at a
certain moment the more electrons will be emitted. As a result, the light signal is converted into anelectron signal. Next, the electron signal passes through an external electric field linearly varying intime, which deflects earlier arriving electrons differently from later arriving electrons. In this way,
the temporal profile of the electron signal is converted into the spatial profile when it is collectedon the fluorescent screen, leading to the characterization of the temporal profile of the incident
light pulse. The conventional streak camera can reach a temporal resolution of a few hundred
femtoseconds. For even faster phenomena, like electronic dynamics in atoms and molecules, theconventional streak camera falls short, since increasing the gradient of the linear varying laser field
to improve the temporal resolution becomes increasingly difficult.
The attosecond streak camera [42], as shown in Fig. 1.7(b), shares similar ideas. Its basic
principle is first to convert the light signal into an electron signal, which is then converted into themeasurable momentum profile instead of a spatial profile. As shown in the figure, a gas jet (notshown) of atoms (yellow dots) is illuminated by the ultrashort XUV light pulse (cyan curve), which
ionizes the atoms and produces free electrons. In this way, the temporal profile of the ultrashortlight pulse is converted into an electron signal. In contrast to the conventional streak camera inwhich a linear varying electric field is applied, here an additional few-cycle NIR streaking laser
pulse (red curve), synchronized with the ionizing XUV laser pulse, is applied, which modulates theasymptotic momentum of the ionized electrons according to its instant of transition into continuum.
By measuring the final momentum profile of the electrons collected at the detector, the electron

12(a) Scheme of the conventional streak camera(b) Scheme of the attosecond streak camera(c) Principle of the attosecond streak camera

Figure 1.7. (a) Scheme of the conventional streak camera (image credit Ref. [3]), (b) scheme of
the attosecond streak camera (image credit Ref. [48]), and (c) principle of the attosecond streak
camera (image credit Ref. [3]).

13signal, and in turn the temporal profile of the ultrashort XUV light pulse, can be characterized.
In this case, however, the ”deflection” occurs within half a period of the streaking pulse, which
provides the possibility for measurements in the attosecond regime.
Although the attosecond streak camera was originally introduced for the characterization ofultrashort light pulses, it can also be used to time resolve ultrafast electron dynamics. As shownin Fig. 1.7(c), electrons released at different instants t0 drift differently in the streaking laser field
and have a different shift in the final momentum according to their instant of release t0
pf ≈ p0 − A(t0 )/c,

(1.7)

where p0 is the asymptotic momentum of the electron in the absence of the streaking field, A(t0 )is the vector potential of the streaking field at the instant of release t0 of the electron, and c is thevacuum speed of light. By scanning the relative delay τ between the ionizing XUV laser pulse andthe streaking NIR laser pulse, a streaking trace can be obtained which consists of the momentum
distributions of the electron as a function of the time delay τ . Fitting the streaking trace to thevector potential of the streaking field, we can obtain the time information regarding the electron
release. Eqn. (1.7) is accurate as long as the dipole approximation is valid and there is no otherpotential present. However, in the ionization of an atom, the Coulomb potential of the residual ion
is present, and the electron motion in the coupled fields of the Coulomb potential and the streaking
laser field induces an error to Eqn. (1.7). In Sec. 5.2, we will study this Coulomb-laser coupling
effect and account for the additional time delay induced to obtain the intrinsic time delay in thephotoionization process of atoms and molecules.

1.6

Outline of this thesis
Recent advances in a variety of light sources such as the ultrashort intense lasers [3, 15],

synchrotron radiations [51], and free-electron lasers [52–54], novel measurement techniques like thecold-target recoil-ion momentum spectroscopy [55, 56], and theoretical approaches [57, 58] have
led to a profound understanding of the correlated emission of two electrons in light-induced double

14
ionization of atoms and molecules [51, 59–73]. The role of electron correlation in double ionization,
as compared to the single-electron emission where only electron-nucleus interaction plays a role, ishowever still less understood, especially regarding the recently observed energy exchange between
electrons over distances of several Angstroms [74]. This thesis is thus devoted to the theoreticalstudies and numerical simulations of ultrafast correlated electron dynamics in single and double
photoionization of atoms and molecules, in an effort to better understand the highly correlated
few-body dynamics.
In Chapter 2, we introduce the numerical methods used in this thesis for an ab-initio solutionof the time-dependent Schrödinger equation (TDSE) in the grid representation with the finitedifference method. The numerical methods to obtain the ground state of a system and to propagate
the initial state in a time-dependent potential, the absorbing boundaries, and the space-partitioning
method to obtain single and double ionization contributions are presented. Parallel computation,
which uses multiple concurrent threads for a speedup of the calculations, is also briefly discussed.
In Chapter 3, we first study the energy exchange via electron correlation upon photon absorption over large distances in the two-site double photoionization of the helium dimer, which isfound to be a two-step process in our time-dependent numerical simulations. The double ionization
probability is shown to be closely related to that of the photoemission of an electron from one of
the helium atoms along the internuclear axis. Together with an analysis of the temporal evolutionof the two-electron probability distribution, this provides direct evidence for the two-step knockout
mechanism by which the photon energy is shared between the electrons over distances of several
Angstroms in the dimer. In the first step, as shown, one electron (the primary electron) in oneof the atoms absorbs the photon and gets ionized. In the second step, this electron propagates
towards the neighboring atom and knocks out the other electron (the secondary electron). We thenintroduce the Hamiltonian reduction method, which selectively removes certain interactions fromthe full Hamiltonian, in order to further study the effects of different interactions in the single anddouble photoionization of the helium dimer.
In Chapter 4, we first analyze the selection rules for the emission of two electrons from the15helium atom following the absorption of a few photons in an intense laser field. The rules arise
due to the symmetries of the accessible final states in the two-electron continuum. In particular,
the back-to-back emission of the two electrons with equal energy sharing is either suppressed or
not depending on the number of photons absorbed from the field. We then further generalize theselection rules to the helium dimer and general molecules and conclude that the rules apply as longas the system, excluding the two active electrons, has the central symmetry respect to its center ofmass.
In Chapter 5, we study the time delay between the single and double photoionization processes, namely the correlation time delay. To this end, we first propose a self-consistent-time
method to account for the Coulomb-laser coupling effect and obtain the intrinsic photoabsorption
time delay measured by the attosecond streak camera. We then proceed to time resolve the correlated emission of two electrons in the knockout process of the helium dimer with respect to thefirst step of single ionization using the variable space-partitioning method as well as the attosecond
streaking technique, which is extended to many-electron systems. The streaked variable for the
single ionization is the momentum of the first electron, while the streaked variable for the double
ionization is the center-of-mass momentum of both electrons. The correlation time delays obtainedfor different laser parameters are compared with classical estimates.
At last, a summary is given in Chapter 6.

Chapter 2

Numerical methods

Central to the analysis of time resolving the laser-matter interaction is solving the TDSE of
the system. For complex systems containing many charged particles, solving TDSE directly can be
extremely difficult. In this case, one can use standard chemical methods such as the time-dependent
Hartree–Fock method and the time-dependent density-functional theory. For simpler systems where
one or two electrons are involved in the interaction, we can solve the associated TDSE using abinitio methods, which are more accurate than the chemical methods and are capable of resolving
more details of the physics of a certain phenomenon.
In this chapter, we discuss the numerical methods to solve the TDSE in grid representations,
i.e., using the finite-difference method. The basic parameters to characterize a specific simulation
are the spatial step ∆x (and ∆y, ∆z, etc.), number of grid points Nx (and Ny , Nz , etc.), and thetime step ∆t.

2.1

The split-operator method
In this section, we describe the split-operator method, which reduces the time propagationof the multidimensional electron wave function to different directions. To make the scheme transparent, we illustrate the method with a system containing only one electron, the TDSE of which
reads
i

∂
ψ(r, t) = H(r, t)ψ(r, t),
∂t

(2.1)

17
where the Hamiltonian of the system isH(r, t) =[p − A(t)]2+ V (r),2(2.2)

in which A(t) is the vector potential of the laser pulse and V (r) is the Coulomb potential of theresidual ion. For small time steps ∆t, the evolution of the wave function can be approximated asψ(r, t + ∆t) ≈ e−iH(r,t)∆t ψ(r, t).

(2.3)

In Cartesian coordinates, the Hamiltonian, by convenience, can be decomposed to different
directions as
H(r, t) = Hx (r, t) + Hy (r, t) + Hz (r, t),(2.4)

in which[px − Ax (t)]2 1+ V (r),23[py − Ay (t)]2 1
Hy (r, t) =+ V (r),2
3[pz − Az (t)]2 1
Hz (r, t) =+ V (r),2
3

Hx (r, t) =

(2.5)
(2.6)
(2.7)

where Hx , Hy , and Hz can be regarded as the Hamiltonian in the x, y, and z directions, respectively.
The time-propagation operator can now be approximated ase−iH(r,t)∆t ≈ e−iHx (r,t)∆t/2 e−iHy (r,t)∆t/2 e−iHz (r,t)∆t e−iHy (r,t)∆t/2 e−iHx (r,t)∆t/2 ,(2.8)

which is accurate to the second order in time. This follows straightforward fromex̂+ŷ = eŷ/2 ex̂ eŷ/2 + O3 ,(2.9)

which can be proven using the Baker–Campbell–Hausdorff formula1
1
1log(ex̂ eŷ ) = x̂ + ŷ + [x̂, ŷ] + [x̂, [x̂, ŷ]] − [ŷ, [x̂, ŷ]] + · · · .
21212(2.10)

Using the split-operator method, the propagation in arbitrary dimensions can be reduced
to the propagations in one direction at a time, which greatly improves the speed of numericalcomputations. Thus, for clarity, we restrict our discussions of the numerical methods below to one
dimension.

18

2.2

The Crank–Nicolson method
The Crank–Nicolson method [75] is a unitary, energy conservative, and unconditionally stable

method to solve the TDSE with an accuracy of second order in time and space. Here, the time
evolution operator can be written in the Cayley form as
e−iH(t)∆t ≈1 − iH(t) ∆t
2
1 + iH(t) ∆t
2,(2.11)

which is an unitary transformation accurate to the second order in time. From this it follows




∆t
∆t
1 + i H ψ(t + ∆t) = 1 − i H ψ(t).
22(2.12)

With the split-operator method, the Hamiltonian in the x direction has the form
Hx = A(x, y)∂2∂+ B(x, y)+ V (x, y),2∂x∂x(2.13)

where y stands for all other coordinates.
The Crank–Nicolson method for the coordinate x readsψn+1 (t + ∆t) − 2ψn (t + ∆t) + ψn−1 (t + ∆t)∆t
An (y)
ψn (t + ∆t) + i
2
(∆x)2

ψn+1 (t + ∆t) − ψn−1 (t + ∆t)+Bn (y)+ Vn (y)ψn (t + ∆t)2∆xψn+1 (t) − 2ψn (t) + ψn−1 (t)ψn+1 (t) − ψn−1 (t)∆t
An (y)+ Bn (y)+ Vn (y)ψn (t) ,=ψn (t) − i
2
(∆x)2
2∆x
(2.14)
where n represents the index of the respective grid point in the x direction. This leads to a system

19
of linear equations with tridiagonal matrices

2 X1


1 X2



 0

 .
 ..



 0


0

2 Y1


1 Y2



 0
=
 .
 ..



 0


0



3 X1

0

0

···

0

0

2 X2

3 X2

0

···

0

0

1 X3

2 X3

3 X3

···

..
.

..
.

..
.

..

.

0
..
.

0
..
.

0

0

0

···

1 XN −1

2 XN −1

0

0

0

···

0

1 XN

3 Y1

0

0

···

0

0

2 Y2

3 Y2

0

···

0

0

1 Y3

2 Y3

3 Y3

···

..
.

..
.

..
.

..

.

0
..
.

0
..
.

0

0

0

···

1 YN −1

2 YN −1

0

0

0

···

0

1 YN



0

  ψ1 (t + ∆t) 






0   ψ2 (t + ∆t) 





0   ψ3 (t + ∆t) 




.. 
.
..


. 








3 XN −1  ψN −1 (t + ∆t)ψN (t + ∆t)
2 XN



(2.15)  ψ1 (t) 






0   ψ2 (t) 





0   ψ3 (t) 

,


.. 
.
..


. 








3 YN −1  ψN −1 (t)ψN (t)
2 YN
0where N is the total number of grid points in the x direction andi∆ti∆t
An (y) −Bn (y),
2
2(∆x)
4∆x
i∆t
i∆t
An (y) +Vn (y),
2 Xn = 1 −
2
(∆x)
2
i∆t
i∆t
An (y) +Bn (y),
3 Xn =
2(∆x)24∆xi∆ti∆t
An (y) +Bn (y),
1 Yn = −
2
2(∆x)
4∆x
i∆t
i∆t
An (y) −Vn (y),
2 Yn = 1 +
2
(∆x)
2
i∆t
i∆t
An (y) −Bn (y).
3 Yn = −
2
2(∆x)
4∆x
1 Xn =(2.16)(2.17)(2.18)(2.19)(2.20)(2.21)

There are a number of algorithms to solve this system of equations, such as the forward andbackward substitution [76]. The matrix multiplication of the right-hand side can be carried out
straightforward. For the left-hand side, one gets that ψ1 (t + ∆t) depends on ψ2 (t + ∆t) from thefirst line. After substitution to the second line, one gets that ψ2 (t + ∆t) depends on ψ3 (t + ∆t) and20
so on. Finally, in the last line, one obtains the value of ψN (t + ∆t). Now one just inserts this value
into the (N − 1)th line and gets ψN −1 (t + ∆t). One can continue this procedure until the whole
wave function is determined.

2.3

Boundary conditions
We must pay special attentions to the borders of the grid, especially when using symmetricor antisymmetric wave functions (half of the grid in use), cylindrical coordinates, and sphericalcoordinates.
If we use a symmetric wave function in x, which has the grid points defined in the positive x
1
direction from ∆x
2 through (N − 2 )∆x, the other side in the negative x direction is just a symmetricmirror of the positive side. Further considering ψ − ∆x
= ψ + ∆x(or ψ0 = ψ1 ), we have22∂ψ2 − ψ0ψ2 − ψ1
ψ1 =
=
,
∂x
2∆x
2∆x
∂2
ψ0 − 2ψ1 + ψ2
ψ2 − ψ1
ψ1 =
=
.
2
2
∂x
(∆x)
(∆x)2

(2.22)(2.23)

This leads to coefficients in Eqn. (2.21) as
i∆t
i∆t
i∆t
A1 (y) −B1 (y) +V1 (y),
2
2(∆x)
4∆x
2
i∆t
i∆t
i∆t
A1 (y) +B1 (y) −V1 (y).
2 Y1 = 1 +
2
2(∆x)
4∆x
2

2 X1 = 1 −(2.24)(2.25)For an antisymmetric wave function (ψ0 = −ψ1 ), the respective equations are∂ψ2 − ψ0ψ2 + ψ1
ψ1 =
=
,
∂x
2∆x
2∆x
∂2
ψ0 − 2ψ1 + ψ2
ψ2 − 3ψ1
ψ1 =
=
,
∂x2
(∆x)2
(∆x)2
3i∆t
i∆t
i∆t
A1 (y) +B1 (y) +V1 (y),
2 X1 = 1 −
2
2(∆x)
4∆x
2
3i∆t
i∆t
i∆t
A1 (y) −B1 (y) −V1 (y).
2 Y1 = 1 +
2
2(∆x)
4∆x
2

(2.26)(2.27)(2.28)(2.29)

For the radial coordinate in cylindrical and spherical coordinates one assumes ψ − ∆x=2
ψ + ∆xand the above equations for the symmetric wave functions can be used.
2

21
The angular coordinate in cylindrical coordinates leads to the boundary condition ψN = ψ1 ,which can be achieved by padding the grid points with this periodic boundary condition in mind.
Similar treatment can be used for spherical coordinates.

2.4

The imaginary-time-propagation method
The imaginary-time-propagation method [77] is a way to obtain the eigenstates of the system

studied. Within the imaginary-time-propagation method, the eigenstates are calculated in the orderof their energies. That is to say, the ground state needs to be obtained first, then the first excitedstate, and next the second excited state and so on. First, we start from an arbitrarily chosen initial
guess for the ground-state wave function of the system ψ, which can be expanded in terms of thestill-unknown eigenstates of the system
ψ(0) =

Xcn ψn (0),(2.30)

X(2.31)nthe time-evolved version of which isψ(∆t) =

Xcn ψn (∆t) =n

cn ψn (0)e−iEn ∆t ,nwhere cn is the expansion coefficient and En is the eigenenergy of the n-th eigenstate.
Instead of using the real time ∆t for the time propagation, we use an imaginary ∆t. Thus,
let us substitute ∆t = −i∆τ , where ∆τ is now a real time, we then haveψ(∆t) =

Xcn ψn (0)e−En ∆τ .

(2.32)

nSince the ground state has the lowest energy E0 , the weight of the ground state in the full wave
function grows with time as compared to excited states.
By carrying out the imaginary time propagation multiple times and normalizing the fullwave function after each time step, one will finally obtain the ground state of the system whenconvergence is reached. The criteria for convergence can be the energy of the state, i.e., when thecalculated energy of the previous time step is approximately equal to the energy of the current
step. In this thesis, we use an absolute error of 10−16 a.u. in all results presented.

22
The first excited state can be obtained in the same way, but one needs to project out the
ground state from the full wave function (using Gram–Schmidt orthogonalization) at each time
step. Similarly, higher excited states can be obtained by projecting out all lower states.

2.5

The spectral method
The eigenenergies can be obtained very accurately using the imaginary-time-propagation

method since the convergence criteria based on the eigenenergy can be chosen very strict, and theeigenstates are also reasonably accurate. However, if we propagate the eigenstates without external
field, we can often still see a very small portion of the probability coming out from the boundeigenstates. This indicates that the eigenstates obtained using the imaginary-time-propagation
method are not extremely accurate. This may not be a problem when we study strong ionization,in which the ionization signal overwhelms the error, it may however pose a risk to the cases of weakionization.
To solve this problem, we can propagate the eigenstates in the free field for a long time
until the fake ionization signals have gone out of the grid and normalize the eigenstates again.
Alternatively, we can use the spectral method [78] to improve the eigenstates, whose accuracy
depends on that of the eigenenergies here.
The spectral method is based on the spectral properties of the Hamiltonian. Using thismethod, both eigenenergies and eigenstates can be obtained. Following is a brief review of thespectral method [78] that we only present for the nondegenerate case here. To this end, we write
the solution to the TDSE as a superposition of eigenstates
ψ(r, t) =X

An un (r) exp(−iEn t),(2.33)nwhere An is the coefficient, un (r) is the eigenstate, and En is the corresponding eigenenergy. Let
us define the correlation function J(t) as
J(t) = hψ(r, 0)|ψ(r, t)i =Xn|An |2 exp(−iEn t),(2.34)

23
and its Fourier transform is
J(E) =

X|An |2 δ(E − En ).

(2.35)n

As it is obvious from the equation above, we can get the spectrum of the correlation function with infinitely sharp peaks corresponding to the eigenenergies of the system if we record thecorrelation function J(t) for an infinitely long time. However, in reality this is of course impossible. Instead, we can only record J(t) up to a certain time T , and the resonances have the formsin[(E − En )T /2]/(E − En ), whose sidelobes might confuse the overall spectrum. To overcome thisproblem, we can add a normalized Hanning window function w(t)/T to the correlation function
J(t) =

X|An |2n

wherew(t) =w(t)exp(−iEn t),
T





1 − cos(2πt/T ),(0 6 t 6 T )




0,

(t > T )(2.36)(2.37)and thus
J(E) =

X|An |2 L(E − En ),(2.38)nwhere L(E − En ) is a lineshape function corresponding to the Fourier transform of the windowfunction. Therefore, the eigenenergies of the system can be obtained all at once, from the spectrum
properties of the Hamiltonian corresponding to the system. Shown in Fig. 2.1 is the spectrum of
the hydrogen atom obtained with a cylindrical 2D model.
After the eigenenergies are known, the eigenstates corresponding to the eigenenergies can beobtained. To this end, we multiply both sides of Eqn. (2.33) by w(t) exp(iEt)/T and integrate from0 to T , which yields
Z T
ψ(r, E) =ψ(r, t)0Xw(t)exp(iEt)dt =An un (r)L(E − En ).
Tn(2.39)The eigenstate um (r) corresponding to the eigenenergy Em can be obtained by substituting Em tothe equation above, i.e.,ψ(r, Em ) = Am um (r)L(0) +X
n6=m

An un (r)L(Em − En ) ≈ Am um (r)L(0).

(2.40)

24

5

Spectrum strength (arb. units)

10

0

10

−5

10

−10

10

−15

10

−0.5

−0.4

−0.3
−0.2
Energy (a.u.)

−0.1

0

Figure 2.1. The spectrum of the Cylindrical 2D hydrogen atom. The peaks correspond to theeigenenergies of the atom. Parameters are ∆ρ = ∆z = 0.1, Nρ = Nz = 1000, ∆t = 0.01, recording
time T = 10000.

25
The advantage of the spectral method, as mentioned above, is that all the eigenenergies ofthe system can be obtained at once, as long as the correlation function is recorded long enough sothat its Fourier transform shows sharp peaks corresponding to the eigenenergies. However, to get
an eigenenergy as accurate as the results from the imaginary-time-propagation method is very time
consuming. Thus, if we only need to obtain the ground state and a few excited states, we can use
the imaginary-time-propagation method, which yields very accurate eigenenergies and reasonably
accurate eigenstates, which can then be improved using the spectral method based on the accurate
eigenenergies corresponding to the eigenstates obtained. Otherwise, if we need to obtain highly
excited states, we can use the spectral method to obtain both eigenenergies and eigenstates.

2.6

Absorbing boundaries
The essence of the grid-based numerical method is to initially keep the entire wave function

on the grid. However, when ionization happens, the wave packet may reach the border of the grid.
In order to hold the outgoing parts of the wave function on the grid, the grid space needs to belarge, which is time and memory consuming for a numerical simulation.

Figure 2.2. The scheme of the absorbing-boundary method.

To overcome this difficulty, grids with absorbing boundaries are used. Without absorbing
boundaries, using a border would be equivalent to having an infinitely high wall, which would
induce strong reflections of the wave packet that easily obscure the relevant physics. On the other
hand, with absorbing boundaries, which introduces a smooth transition between the region ofinterest and the border, minimal reflections occur. The scheme of the absorbing-boundary method

26
is illustrated in Fig. 2.2.
There are two major types of absorbing boundaries used in the grid-based numerical simulations, i.e., the mask-function method [79] and the exterior-complex-scaling method [80, 81].

2.6.1

The mask-function method

1.0
0.8
0.6
0.4
0.2
-4

-2

2

4

Figure 2.3. The typical shape of the mask function. The index γ is chosen to be 1/6 and theabsorber (shaded area) spans 10% of the grid size in each direction.

In the mask function method [79], the probability is multiplied by a function in each time
step, which introduces a smooth and gradual change of the probability to 0 at the border of thegrid. Usually, the form of the mask function is






0|
cosγ π2 |x−x,Lf (x) =
1,(|x| > |x0 |)(2.41)(|x| 6 |x0 |)where x0 denotes the border of the absorber region and L is its width. The absorber is usuallychosen to span 10% of the grid size in each direction and common choices of the index γ are 1/6,
1/8, etc. Shown in Fig. 2.3 is the typical shape of the mask function.
The mask-function method is usually good enough to suppress the reflection of the probability
from the border. However, in case of a strong ionization (with an ionization yield larger than 10−2 ),the reflections may become remarkable. In this case, we need to use alternative methods such as
the exterior-complex-scaling method.

27
2.6.2

The exterior-complex-scaling method

The idea of the exterior-complex-scaling method [80, 81] is to rotate the coordinate of theparticle by an angle η (0 < η < π/2) into the complex plane when it is larger than a certain value
x0 (border of the absorber), x → xeiη , or




x,





x → |x0 | + (x − |x0 |)eiη ,







−|x0 | + (x + |x0 |)eiη .

(|x| 6 |x0 |)(x > |x0 |)(2.42)(x < −|x0 |)By means of this transformation, the wave function is kept unchanged when |x| 6 |x0 | andfalls off exponentially when |x| > |x0 | due to the decaying factor in the propagation operator present
in this transformation. A similar method with the same goal is the complex-absorbing-potential
method [82–84], in which an additional complex potential is added to the Coulomb potential for|x| > |x0 |V (x) →V (x),(|x| 6 |x0 |)V (x) + iηW (x).

(|x| > |x0 |)(2.43)

The advantage of the exterior-complex-scaling method, in most cases, is that it leads to abetter absorption than the mask-function method. The advantage of the mask-function method, on
the other hand, is that different ionization yields, such as single and double-ionization yields, canbe obtained separately, when used together with the space-partitioning method presented below.

2.7

The space-partitioning method
The space-partitioning method [85] is a numerical method to separate single and double-

ionization contributions from the whole wave function. The essence of this method lies in the factthat the bound state, the single-ionization wave packet, and the double-ionization wave packet are
separated in the two-particle space. Considering an atomic system for example, both electrons arelocated near the nucleus in the bound states, one of the electrons remains located near the nucleus

28
while the other is far from the nucleus in the event of single ionization, and both electrons are far
from the nucleus in the double ionization process. This observation enables us to spatially separate
the different contributions in the full wave function.

r2

25
20
SI
BD

15

DI

10
5
0

5

10

15

20

25

r1

Figure 2.4. The space partitioning for an atom with dsmall = 6 and dbig = 12.

Specifically, for an atomic system, the different contributions, i.e., the bound state (BD), thesingle-ionization part (SI), and the double-ionization part (DI), can be identified as(1) BD: r1 6 dbig ∩ r2 6 dbig ;(2) SI: (r1 6 dsmall ∩ r2 > dbig ) ∪ (r2 6 dsmall ∩ r1 > dbig );(3) DI: (r1 > dsmall ∩ r2 > dbig ) ∪ (r2 > dsmall ∩ r1 > dbig ),where r1 = |x1 | is the distance of the first electron to the nucleus, r2 = |x2 | is that of the second
electron to the nucleus, and dsmall and dbig are the two bounds to distinguish different contributionsin the full wave function. Such an atomic partitioning for dsmall = 6 and dbig = 12 is illustrated inFig. 2.4.

29
For a diatomic system (two-center system), the partitioning becomes complex, although thebasic principle remains simple, i.e., the two electrons are located near either nuclei but not near
the same nuclei for the bound state, one of the electrons is located near one of the nuclei while theother is far from the other nucleus for the single ionization event, and both electrons are far fromthe nuclei for double ionization.

x2

20
10

-30

-20

SI

10

-10

20

x1

BD
DI

-10

-20

-30
Figure 2.5. The space partitioning for a diatomic system (atomic nature) in the direction of themolecular axis with R = 10, dsmall = 4 and dbig = 7.

Specifically, for such a diatomic system (of atomic nature), the different contributions, i.e.,
BD, SI, and DI, after assuming the middle of the two nuclei as the origin and defining r11 =|x1 − R/2|, r12 = |x1 + R/2|, r21 = |x2 − R/2|, and r22 = |x2 + R/2| with R standing for theinternuclear distance, can be identified as(1) BD: (r11 6 dbig ∩ r22 6 dbig ) ∪ (r12 6 dbig ∩ r21 6 dbig );(2) SI: ((r11 6 dsmall ∩r22 > dbig )∪(r22 6 dsmall ∩r11 > dbig )∪(r12 6 dsmall ∩r21 > dbig )∪(r21 6
dsmall ∩ r12 > dbig )) ∩ (r11 > dbig ∪ r22 > dbig ) ∩ (r12 > dbig ∪ r21 > dbig );30(3) DI: ((r11 > dsmall ∩ r22 > dbig ) ∪ (r22 > dsmall ∩ r11 > dbig )) ∩ ((r12 > dsmall ∩ r21 >dbig ) ∪ (r21 > dsmall ∩ r12 > dbig )),and shown in Fig. 2.5 with R = 10, dsmall = 4 and dbig = 7.

x2

20
10

-30

-20

SI

10

-10

20

x1

BD
DI

-10

-20

-30
Figure 2.6. The space partitioning for a diatomic system (molecular nature) in the direction of
the molecular axis with R = 10, dsmall = 4 and dbig = 7.

If, however, we further identify the wave packet between the two nuclei as bound and oneof the electrons between the nuclei and the other far from the nuclei as single ionization, we caninstead identify the different contributions for such a diatomic system (of molecular nature) as(1) BD: (r11 6 dbig ∩ r22 6 dbig ) ∪ (r12 6 dbig ∩ r21 6 dbig ) ∪ (r1 6 R/2 ∩ r2 6 R/2);(2) SI: ((((r11 6 dsmall ∩r22 > dbig )∪(r22 6 dsmall ∩r11 > dbig )∪(r12 6 dsmall ∩r21 > dbig )∪(r21 6
dsmall ∩r12 > dbig ))∩(r11 > dbig ∪r22 > dbig )∩(r12 > dbig ∪r21 > dbig ))∪((r11 > dbig ∪r22 >dbig ) ∩ (r12 > dbig ∪ r21 > dbig ) ∩ (r1 6 R/2 ∪ r2 6 R/2))) ∩ (r1 > R/2 ∪ r2 > R/2);(3) DI: ((r11 > dsmall ∩ r22 > dbig ) ∪ (r22 > dsmall ∩ r11 > dbig )) ∩ ((r12 > dsmall ∩ r21 >dbig ) ∪ (r21 > dsmall ∩ r12 > dbig )) ∩ (r1 > R/2 ∩ r2 > R/2),

31
which is exemplified in Fig. 2.6 with R = 10, dsmall = 4 and dbig = 7.

Figure 2.7. The schematic illustration of the way to choose the bounds in the space-partitioning
method.

The bounds dsmall and dbig are chosen such that the final single and double ionization yields
do not change significantly upon variation about their values. To this end, we first increase dbig
gradually starting from a small value, and for each dbig we propagate the initial state in the laser
field long after the end of the pulse until the single ionization yield converges. From the change inthe final single ionization yield as a function of dbig , we choose the smallest value of dbig such thatthe yield does not change significantly for larger values. After choosing a value for dbig , we thenincrease dsmall gradually starting from a small value, and for each dsmall we propagate the initialstate in the laser field until the double ionization yield converges. From the double ionization
yield as a function of dsmall , we choose a value for dsmall beyond which the yields does not change
significantly. After these procedures, both values of dsmall and dbig are fixed. Shown in Fig. 2.7 isa schematic illustration of this method, in which b1 < b2 < b3 < b4 < b5 , and the final ionization
yield is essentially converged at b3 .

32

2.8

Calculating single and double-ionization yields
The space-partitioning method, when used together with the mask-function method, provides

a way to calculate single and double-ionization yields by accumulating the probability absorbed by
the boundary. The single and double-ionization rate at a certain time t can be written as
Z
Z
1
1
RSI (t) =|ψ(x, t)|2 dx −M (x)|ψ(x, t)|2 dx,
∆t SI
∆t SI
Z
Z
1
1
2
RDI (t) =|ψ(x, t)| dx −M (x)|ψ(x, t)|2 dx,
∆t DI
∆t DI

(2.44)(2.45)

where M (x) is the mask function, and the integrations are carried over the respective ionization
regions identified by the space-partitioning method. The corresponding ionization yields are
Z T
PSI (T ) =RSI (t)dt,(2.46)0Z T
RDI (t)dt,(2.47)PDI (T ) =0

where the integration time T needs to be large enough so that all the ionized wave packet reaches
the boundary.
When ionization time matters, i.e., when we need to find out when the ionized wave packet
comes out from the system, the wave packet can be considered as single ionization immediately
after it enters the single-ionization region, and the wave packet can be regarded as double ionization
immediately after it enters the double-ionization region, as opposed to the previous case where it
is considered as ionization only when the probability gets absorbed by the boundary. In this case,the ionization yields can be written as
Z T −∆t
Z
0
PSI (T ) =RSI (t)dt +|ψ(x, T )|2 dx,
0SI
Z T −∆t
Z
0
PDI (T ) =RDI (t)dt +|ψ(x, T )|2 dx.
0(2.48)(2.49)DI

The first method, i.e., collecting the electrons at the boundary, is used to calculate the single
and double-ionization yields in this thesis, unless stated otherwise. The second method is onlyused when we need to obtain the ionization time directly from the spatial distribution of the wave
function (refer to Sec. 5.3.2.2).

33

2.9

Parallel computation
Studies of correlated electron dynamics, where two or more electrons are involved, require

higher dimensional calculations and parallelization of the code. There are two methods for code
parallelization, namely shared-memory parallelization and distributed-memory parallelization.
Shared-memory parallelization means that all working threads share the same common memory that resides in the same (physical or virtual) compute node. OpenMP is typically used to implement shared-memory parallelization. The advantage of this method is that the memory is shared
between different working threads and thus minimum communication between different threads isneeded, which makes this method fast and effective. The disadvantage, on the other hand, is thatthe maximum number of threads in a single compute node where the common memory resides
is limited and thus the parallelization is restricted to a few threads, typically 12 or 16 currently,
depending on how many cores or threads the compute node has.
Distributed-memory parallelization means that all working threads carry their own private
memories that do not necessarily locate in the same compute node. MPI is typically used to
implement distributed-memory parallelization. The advantage of this method is that the number ofavailable working threads is not limited to the number of cores or threads in a single compute node,and its number can reach the maximum available threads within a whole cluster. The disadvantage,on the other hand, is that the memory is distributed between different working threads, and if thedifferent threads work together to calculate the same physical quantity, heavy communications
between threads can be induced. While the communication within the same compute node is fast,that between different nodes can be very slow depending on the speed of ethernet and the amountof messages to be passed across compute nodes.
In our case, where we study correlated electron dynamics, different parts of the full wave
function are highly correlated. If we use distributed-memory parallelization, heavy communications are needed which would become a bottleneck in scaling up the performance. So we choose
the shared-memory parallelization for our code. We have also estimated that using distributed-

34
memory parallelization, a minimum performance boost, if any, would be achieved as compared tousing shared-memory parallelization, even if we use a lot more compute nodes. Cost-effectively,
we thus stick to the shared-memory parallelization, i.e., using OpenMP, for the purpose of code
parallelization.
Using OpenMP is very easy, straightforward, and induces minimum changes to the singlethreaded code. The minimum working example below shows the comparison between the singlethreaded code (left column) and OpenMP parallelized multiple-threaded code (right column), where
the number of threads to use can be determined by the environmental variable OMP_NUM_THREADS.
1
2
3
4
5
6

int i ;for ( i =0; i < N ; i ++){// code here
}

1
2
3
4
5
6

int i ;
# pragma omp parallel for private ( i )for ( i =0; i < N ; i ++){// code here}

As we can see, only one line starting with the #pragma compiler directive is added whichtells the compiler to compile with OpenMP (if, say, the argument -fopenmp is passed to the g++compiler). The variables inside the parenthesis of the private keyword are the private variablesof each working thread that are not to be shared between different threads. If operators likeaccumulation are used, the reduction method is necessary to avoid the case when different threads
add to the same variable simultaneously and the variable loses track of some additions. A minimum
example is given as follows. In this case, each thread keeps a private copy of the variable sum andperforms the addition of sum from each thread at the end.
1
2
3
4
5
6
7
8

int i ;
int sum =0;# pragma omp parallel for private ( i ) reduction (+: sum )for ( i =0; i < N ; i ++){// code heresum += i ;}If, by any means, we have to disable parallelization, we only need to set the environmental

35
variable OMP_NUM_THREADS to 1, or, at compile time, omit the argument -fopenmp of the g++compiler. Hence, adding such compiler directive as #pragma does not add much to the working
load and does not induce any harm to the computation since we can turn it off at any time.
Therefore, we strongly suggest to always parallelize the code using OpenMP even if we do not use
it later on.
To reach the maximum performance, we need to use parallelization whenever possible. If
we parallelize the code everywhere else but leave one segment single-threaded, then this segment
may become the bottleneck. Furthermore, to reduce the parallelization initialization and collection
overhead, we usually parallelize the outermost loop if there are more than one loops. At the sametime, the sequence of the loops needs to be arranged in such a way that the rule of strided memory
access is obeyed, i.e., the memory needs to be accessed in a continuous flow, in order to maximize
the performance.

Chapter 3
Single and double photoionization of the helium dimer1

In this chapter, we study the energy exchange via electron correlation upon photon absorption
over large distances in double photoionization of the helium dimer. Results of numerical simulationsof the interaction of a planar 4D helium dimer model with a short light pulse are found to be ingood agreement with recent experimental data for the angular distribution of the emitted electron.
The double ionization probability is closely related to that of the photoemission of an electron fromone of the helium atoms along the internuclear axis. Together with an analysis of the temporal
evolution of the two-electron probability distribution, this provides direct evidence for the knockout
mechanism by which the photon energy is shared between the electrons over distances of several
Angstroms in the dimer.
Based on the results of numerical simulations, we further study the scattering effects in single
and double photoionization of the helium dimer. To this end, we use abridged Hamiltonians in which
different interactions between the electrons and the laser field and/or between the charged particles
are removed from the full Hamiltonian. By comparing photoelectron angular distributions obtainedin respective numerical simulations, we are able to identify the role of the Coulomb interactions
between the electrons and between the electrons and the nuclei for the primary as well as theknockout electron in the processes.
1

The results of this chapter are in part based on “H. Ni, C. Ruiz, R. Dörner, and A. Becker, Phys. Rev. A 88,
013407 (2013)” [86] and “H. Ni and A. Becker, Phys. Rev. A 89, 033402 (2014)” [87].

37

3.1

Introduction
Despite recent studies on correlated electron dynamics [51, 59–74], the role of electron corre-

lation in double photoionization of a molecule, in particular, regarding the recently observed energy
exchange between electrons over distances of several Angstroms [74], is still less understood thanin the atomic case. In this respect rare gas dimers, which are formed via the attractive van der
Waals interaction by two atoms at larger equilibrium distances than those in a typical diatomic
molecule, are interesting targets. Among them, the helium dimer is by far the most extended.
Here, interatomic Coulombic decay (ICD) [88], as well as double photoionization [89], which areboth mediated by electron-correlation effects, have been observed recently.
In double photoionization, the energy of a single photon absorbed from the light field is shared
between two electrons, leading to the correlated emission of both electrons from the target. In arare gas dimer, the two emitted electrons can either originate from the same atom or from different
atoms. Due to the strong localization of electrons in the dimer, the minimum energy required
for the emission of two electrons from different atoms in the dimer is about twice the energy forsingle ionization of the rare gas atom. For He2 , this minimum energy is about 49.2 eV, which isconsiderably smaller than the energy of about 79 eV needed for emission of both electrons from thesame helium atom in the dimer. The threshold for ionization plus excitation (to the n = 2 state)in He is 65.4 eV. Below this threshold ICD as a double ionization channel is closed. By selecting aphoton energy of 49.2 eV < ~ω < 65.4 eV, one can therefore study electron correlation effects indouble photoionization of the helium dimers over a bond length of more than 5 a.u., which are not
ICD related.
Recent experimental data suggest that the emission of the two electrons in this case proceeds
as follows [89]: A primary electron localized at one atom in the helium dimer absorbs the photon
energy from the field, then propagates along the internuclear axis and transfers its energy to asecond electron in the neighboring helium atom. This process is called knockout in view of theclose analogy to a similar mechanism known in double photoionization of atoms [7]. Since the atoms

38
in the helium dimer are well separated it was further argued [89] that in this case the knockout
process can be considered as the photoionization of a helium atom followed by a subsequent electronimpact ionization at the other helium atom in the dimer. This interpretation was supported by acomparison of the experimental data with theoretical results for electron-impact ionization of thehelium atom. In this chapter, we analyze the experimental data and the underlying mechanism for
double photoionization of the helium dimer via time-dependent ab-initio numerical simulations.
However, the role of the interaction between each pair of charged particles during the single
and double photoionization processes is still difficult to analyze. In a numerical simulation, in
principle, it is possible to remove interactions between certain pairs of particles in the Hamiltonianand compare the results to those obtained with the full Hamiltonian. However, using this techniquein an atomic or typical molecular system, the initial state inevitably changes as compared to thefull simulation. Consequently, it is difficult to distinguish between the influence of the initial stateand dynamic effects when analyzing the role of the removed interactions. In contrast, this method
is applicable to rare gas dimers, since they have much larger equilibrium distances than those of
typical diatomic molecules. Due to the large internuclear separation and almost negligible overlap
between the orbitals of the two atoms, changes in the initial state are small when certain interactions
are removed.
As we will further show in this chapter, in the helium dimer it is possible to study therole of individual Coulomb interactions during the single and double photoionization in ab-initio
numerical simulations. To this end, we consider abridged Hamiltonians, in which certain Coulombinteractions are removed without noticeable change in the initial state. We then compare theangular distributions obtained using the abridged Hamiltonian to those obtained with the full
Hamiltonian. In this way we are able to identify the role of each Coulomb interaction during theprocess. Moreover, we are also able to distinguish and compare the angular distributions of the
primary electron and the knockout electron by removing the coupling between one electron and thelaser field from the full Hamiltonian.
The remainder of the chapter is organized as follows. We first briefly outline the planar two-

39
active-electron model of the helium dimer used in the simulations. Then, we analyze the projection
to different final states in the calculations to obtain the photoelectron angular distributions andpresent comparisons with recent experimental data [89]. Based on the good agreement between experimental and theoretical data we proceed to provide evidences towards the knockout mechanism,based on predictions for the dependence of the ionization signals on the orientation of the dimer
axis and the separation of the atoms in the dimer, and a real-time visualization of the knockout
mechanism behind the long-range electron correlation. Afterwards, we introduce and discuss different abridged Hamiltonians in which certain interactions are removed from the full Hamiltonian. A
comparison between the ground states of the full and the abridged Hamiltonians will be given anddiscussed. We then apply the different Hamiltonians in numerical simulations to study the effects
of different interactions in single and double photoionization.

3.2

Numerical simulations
In this section we present the numerical model for the helium dimer used in the simulations,in which one electron at each atom in the dimer is considered as active. We further outline themethods used in the numerical simulations and then present the predictions for the lowest energy
eigenstates based on this numerical model.

3.2.1

The planar 4D model of the helium dimer

A solution of the TDSE including the dynamics of all six charged particles and the interactionof all four electrons with the external light field is currently not conceivable. We therefore propose aplanar two-active electron model of the helium dimer, in which the electrons are located at different
atoms in the dimer and their dynamics is constrained to the same plane, as shown schematicallyin Fig. 3.1. Such a model excludes ICD but it is well suited to treat the knockout mechanism. The
orientation of the dimer axis is chosen to be at an angle θ to the axis of linear polarization of thelight field (here, the z-axis). It is further assumed that the correlated emission of the two electrons
proceeds fast enough (a few tens of attoseconds, 1 as = 10−18 s) to hold the internuclear distance40

Figure 3.1. Scheme of the numerical model for the helium dimer. One electron at each atom is
considered active and restricted to the same plane, with coordinates (x1 , z1 ) and (x2 , z2 ), respectively. The internuclear distance R and dimer orientation θ with respect to the laser polarization
direction ε (z-axis) are fixed during the simulations.

41
between the atoms as well as the orientation of the dimer axis fixed in space during the simulation.
The corresponding four-dimensional (two spatial dimensions for each electron) model Hamiltonianis then given by
p2 + p2z1
p2 + p2z2
H = x1+ x2+ E(t) (z1 + z2 )22
+ VSAE (r11 ) + VSAE (r12 ) + VSAE (r21 ) + VSAE (r22 )

(3.1)

1
1
+ ,
+p
2
2
2
R
(x1 − x2 ) + (z1 − z2 ) + b
where pi = (pxi , pzi ) and ri = (xi , zi ) (with i = 1, 2) are the momentum operators and spatial
coordinates of the two active electrons and E(t) is the laser field, which is assumed to be linearlypolarized in z-direction. R is the internuclear distance that can have different orientations θ withrespect to the polarization direction and
s
2 
2
R
R
rij =
xi + (−1)j sin θ + zi + (−1)j cos θ + a2
2
2

(3.2)

corresponds to the distance between ith electron and jth nucleus (i, j = 1, 2) with a2 = 0.201
and b2 = 0.01 denoting soft-core parameters to smooth the Coulomb singularity in the practicalcomputations. The linearly polarized laser field is given in the form of
2

E(t) = cos



ωt
2N


cos (ωt + φ) ,(3.3)

where ω is the central frequency, N is the number of cycles, and φ is the carrier-envelope phase of
the laser pulse, which is set to −π/2 in the present simulations.
VSAE (r) = −Zc + a1 e−a2 r + a3 re−a4 r + a5 e−a6 r
r

(3.4)

is a single active electron (SAE) potential for the helium atom with Zc = 1.0, a1 = 1.231, a2 = 0.662,
a3 = −1.325, a4 = 1.236, a5 = −0.231, and a6 = 0.480 [90], which we adopt here for the planar case.
The TDSE of this 4D model is solved numerically using the Crank–Nicolson method on a grid with
Nx1 = Nx2 = 200 and Nz1 = Nz2 = 300 grid points, a grid step of ∆x1 = ∆x2 = ∆z1 = ∆z2 = 0.3,
and a time step of ∆t = 0.05.

42
In our simulations we considered double photoionization due to the interaction of the model
helium dimer with a 4-cycle XUV light pulse at a central wavelength of 20 nm (corresponding toa photon energy of 62 eV and a bandwidth of ±15.5 eV) and a peak intensity of 1 × 1014 W/cm2 .
At this bandwidth in a real four-electron system, the ICD channel would be open for the highenergy tail of the pulse. In our two-electron calculations, ICD is however excluded and hence does0|
not obscure our findings. We employed absorbing boundaries of the form cos1/6 π2 |x−x
withL|x| > |x0 |, where x0 denotes the border of the boundary region and L its width. The boundary waschosen to span 10% of the grid size in each direction. In order to calculate the total probabilities
and electron angular distributions, we stopped the simulation before the respective part of the wave
function reached the boundaries.
In order to analyze the correlated electron emission from the helium dimer and to comparewith the experimental data on the photoelectron angular distributions, the contributions to single
and double photoionization of the helium dimer were obtained by partitioning the spatial fourdimensional grid as illustrated in Sec. 2.7 for a diatomic system of molecular nature with dsmall = 4
a.u. and dbig = 7 a.u.

3.2.2

Initial states

The lowest energy eigenstates of the model, which we used as initial states for our numericalsimulations, are obtained via imaginary time propagation. The initial guess state for the imaginary
time propagation was chosen to be a 4D Gaussian wave function. For an internuclear separationof 5.6 a.u., which corresponds to the minimum of the helium dimer potential [91], and orientationof the dimer along the z-axis, spatial distributions of the two lowest energy eigenstates are shownin Fig. 3.2 as functions of z1 and z2 (upper row) and z = z1 − z2 (lower row), respectively. The
distributions are integrated over the other coordinates. Due to the large internuclear distance thestates are, as expected, very close in energy (−1.80740 a.u. and −1.80726 a.u., respectively), withvalues close to other theoretical calculations [92]. Note that the emergence of these two separate
states is due to our two-active electron model that selects two electrons out of a total of four. The

43

0

0

(c) 10(d)10−1Probability

Probability

−1

10

−2

10

−3

10

−4

10

−2

10

−3

10

−4

10

−5

10
−10

10

−5

−5

0
z (a.u.)

5

10

10
−10

−5

0
z (a.u.)5

10

Figure 3.2. Spatial distributions (on a logarithmic scale) of the ground state [panels (a) and (c)]and first excited state [panels (b) and (d)] of the planar model helium dimer. The distributions in
the upper row are integrated over x1 and x2 and shown as a function of z1 and z2 . The distributions
2and shown as a function of z = z1 − z2 .
in the lower row are integrated over x1 , x2 , and Z = z1 +z
2Note the node (minimum) at z1 = z2 or z = 0 in the distributions of the first excited state.

44
ground state is a spin-singlet state which has a multiplicity of 1 while the first excited state is aspin-triplet state with a multiplicity of 3. We performed separate simulations for both states as theinitial state and present below, if not mentioned differently, the sum of the contributions by taking
account of the multiplicity of each state.

3.3

Photoelectron angular distributions
The final state of double photoionization of the dimer represents a state with two electrons in

the continuum of a doubly charged two-center potential. In order to obtain photoelectron angular
distributions, we intend to project the respective spatial part of the wave function at the endof the simulation onto an approximate final-state wave function. Since an analytical solution ofthis four-body problem is not known in Cartesian coordinates, we have tested several two-electron
wave functions, in which different parts of the full Coulomb interactions between the charged
particles have been taken into account. To this end, we express the Coulomb potential between theelectrons and the residual ions as well as between the two electrons in the center-of-mass (COM,X = (x1 + x2 )/2 and Z = (z1 + z2 )/2) and relative coordinates (x = x1 − x2 and z = z1 − z2 ) of
the two electrons as
1
V (x, z; X, Z) = √2
x + z2

s

2 
2
2
X
x
R
z
R
+
VSAE 
X − (−1)i + (−1)j sin θ + Z − (−1)i + (−1)j cos θ  .
2
2
2
2

(3.5)

i,j=1

For |X|, |Z|, R  |x|, |z|, V reduces toV (x, z; X, Z) ≈ Vee (x, z) = √1x2 + z 2,(3.6)

which suggests to approximate the final state as a product of a Coulomb wave in the relative anda plane wave in the COM coordinates, which we denote as ψee . Thus, in this approximation the
Coulomb interaction between the two electrons is retained while the Coulomb interactions between
the electrons and the nuclei are neglected. In order to approximately test the relative importanceof the Coulomb interactions in the final state, we have further considered the case in which the45interaction between the electrons is neglected. As a further approximation we then assumed thatthe total charge of the nuclei is concentrated at the center of the coordinate system. In thisapproximation the full potential reduces toV (x, z; X, Z) ≈ VeN (X, Z) = − √4
X2 + Z2,(3.7)

and the final state ψeN is approximated by a product of a plane wave in the relative and a Coulomb
wave in the COM coordinates of the electrons for a monomer with twice the residual charge. As athird option, we also considered a simple product of two plane wave functions to describe the final
state wave function of the two electrons, which we denote as ψplane .
Since our model is planar, we used a two-dimensional scattering wave function [93, 94] for
the Coulomb wave,
ψ(k, r) =∞Xe−πη/21(−2iρ)|l| iρ iσ|l|√Γ |l| + + iηe e
2
(2|l|)!
πl=−∞ il(φr −φk )1e× 1 F1 |l| + + iη, 2|l| + 1, −2iρ
,
2
2π

(3.8)

where 1 F1 is the confluent hypergeometric function. Here, k and r represent either the relative or
COM momentum and displacement of two charged particles in two dimensions with k, r, φk , and φr
being their magnitudes and phases, respectively. η = µZ1 Z2 /k, where µ is the reduced mass of thetwo particles (1/2 in the case of Vee and 2 in the case of VeN ) and Z1 and Z2 are their charges (−1and −1 in the case of Vee and −2 and 2 in the case of VeN ), ρ = kr, and σ|l| = arg Γ(|l|+ 21 +iη) is thetwo-dimensional Coulomb phase. Alternatively, we also expressed the two-dimensional scatteringwave function in terms of the three-dimensional Coulomb wave function Fl ,ψ(k, r) =

∞
X
l=−∞

r

2
eil(φr −φk )(−i)|l| eiσ|l| F|l|− 1 (η, ρ)
,
2
πρ
2π

(3.9)

which enables us to use available packages like GNU Scientific Library (GSL) for the calculation.
In the practical computations we truncated the sum over l to lmax = 12 for ψee and lmax = 16
for ψeN . While we found that the results for the angular momentum distributions presented below
are converged for projection to ψee , the momentum-space wave function for projection onto ψeN

46
Table 3.1. Double ionization probabilities for the different momentum-space wave functionsRψplaneψeeψeN

5.6
3.5520 × 10−5
3.4636 × 10−5
1.0908 × 10−5101.9182 × 10−5
1.7428 × 10−5
6.8121 × 10−6

was not fully converged. However, we checked that our qualitative conclusions below are correct.
We attribute the difficulty in the convergence to the fact that we used a monomer Coulomb wave
function for a two-center problem in ψeN and that the residual potential is not purely 1/r. The
degree of convergence can indeed be estimated by comparing the double ionization probabilities
obtained for the corresponding momentum-space wave function in each case, given in Tab. 3.1.
Evidently, the probabilities of the momentum-space wave functions obtained by projection ontoψplane and onto ψee agree rather well with each other, while the result for the projection onto ψeNindicates the lack of convergence with respect to lmax .
In Fig. 3.3 we compare results for the molecular-frame angular distributions of one of thetwo electrons in the double photoionization after projections onto ψplane [panels (a) and (b)], ontoψee [panels (c) and (d)], and onto ψeN [panels (e) and (f)]. The distributions are integrated overthe emission angle of the other electron and the energies of the electrons. The theoretical data areshown in each set as a function of the emission angle of the electron. In order to compare withrecently obtained experimental data [89], the theoretical results are averaged over the orientationof the dimer axis and weighted with respect to singlet and triplet state contributions according to
the multiplicities of the states. Furthermore, due to the short duration of the pulse used in thesimulations, the calculated distributions show a slight left-right asymmetry along the polarization
direction. We have removed this asymmetry by taking the mean of the original result and its mirror
image with respect the plane perpendicular to the polarization axis through the center of mass ofthe two nuclei. Finally, the theoretical results are matched to coincide with the experimental
data at the maximum of the distributions. Note that the experimental results are integrated overR = 5.1 ∼ 6.8 a.u. in panels (a), (c), and (e), and R = 9.4 ∼ 10.9 a.u. in panels (b), (d),47

Figure 3.3. Molecular-frame angular distributions of one of the two electrons following double
photoionization of the helium dimer, obtained by projection of the respective part of the wave
function onto different approximate final two-electron states at R = 5.6 a.u. (left column) andR = 10 a.u. (right column). Shown is a comparison between the theoretical results (solid lines)
and the experimental data (solid circles). For the theoretical results in panels (a) and (b), ψplane
is used as the approximate final two-electron state, while for the results in panels (c) and (d), ψee ,and for those in panels (e) and (f), ψeN is used.

48

Figure 3.4. Molecular-frame angular distributions of one of the two electrons following double
photoionization of the helium dimer, obtained by projection of the respective part of the wave
function onto ψee at R = 5.6 a.u. (left column) and R = 10 a.u. (right column). Shown is acomparison between the theoretical results (solid lines) and the experimental data (solid circles).
For the theoretical results in panels (a) and (b), full results are shown, while for the results inpanels (c) and (d), singlet contributions, and for those in panels (e) and (f), triplet contributions
are shown.

49
and (f) (subset of data shown in Fig. 3(b) of Ref. [89]). The comparison shows that within thepresent model the approximation of the final-state wave function as product of a plane wave in
the COM coordinates and a Coulomb wave in the relative coordinates qualitatively appears to be
the most appropriate one of the three considered options. We therefore use ψee to obtain angular
distributions for the remainder of our studies. In order to show the individual contributions fromsinglet and triplet initial states, we further show in Fig. 3.4 the same angular distributions obtainedby projection onto ψee but from singlet and triplet states separately.

3.4

The knockout mechanism
On the basis of the comparison presented in Figs. 3.3 and 3.4, we proceed to obtain in-sights into the mechanism of the correlated electron emission from the helium dimer. As mentioned
above, it was argued based on the experimental data and a comparison with theoretical results fromelectron-impact ionization that the double photoionization proceeds via a knockout-type mechanism. According to this mechanism the double ionization probability should be closely related
to the probability that the primary photoelectron, launched at one helium atom, “hits” the second electron which is localized at the neighboring atom. The results in Figs. 3.5 and 3.6 confirm
this close relation based on our theoretical data, which provides strong evidences for the knockout
mechanism.
First, we compare in Fig. 3.5(b) the probability of double photoionization as a function of
the orientation of the dimer axis with respect to the polarization direction of the field (solid line)with the experimental data (solid circles [89]) and a cos2 δ-distribution (dash-dotted line). The
comparison shows not only an excellent agreement between theoretical predictions and experimental
data but moreover indicates a close relation between the double ionization probability and the
direction of the primary photoelectron momentum. The latter is close to the expected p-wave (orcos2 -) distribution as can be seen from the theoretical results in panel (a) for orientation of dimer
axis parallel (red solid line) as well as perpendicular (blue dashed line) to the polarization directionof the field. The small deviation from the cos2 -distribution in the case of parallel orientation

50

Figure 3.5. Comparison of the lab-frame dipole-transition angular distribution cos2 δ (dasheddotted line) with (a) the theoretical photoelectron angular distributions (AD) for single ionization(SI) of the dimer oriented either parallel (red solid line) or perpendicular (blue dashed line) tothe polarization axis and (b) the probability of double photoionization (DI) as a function of theorientation of the dimer axis (theoretical prediction: red solid line, experimental data: solid circles,
integrated over R = 4.5 ∼ 6.8 a.u. (subset of data shown in Fig. 3(a) of Ref. [89])). The polarization
direction of the field is along the horizontal axis in both panels. The internuclear distance was chosen
to be R = 5.6 a.u. and the parameters of the XUV field are the same as in Fig. 3.4.

−6

DI Probability

5

x 10

4

Data
Fit with a/R+b

3
2
1
5

10
R (a.u.)15

Figure 3.6. Double photoionization probability (open circles) as a function of internuclear distance
R for orientation of the dimer axis parallel to the polarization direction of the field. It scales as 1/R
with the internuclear distance. The fitting parameters are a = 2.88 × 10−5 and b = −4.56 × 10−7 .
The parameters of the XUV field are the same as in Fig. 3.4.

51
indicates the effect of elastic scattering of the photoelectron at the second atom in the single
ionization data, in agreement with similar observations for the neon dimer [95].
Further evidence for the knockout mechanism is given by the dependence of the double photoionization probability on the internuclear distance R. For the knockout mechanism, the probability of the primary photoelectron “hitting” the second electron that produces double ionization
is expected to decrease with an increase of R, since the wavepacket representing the primary photoelectron expands in space. Our theoretical results for the double photoionization probability,
shown in Fig. 3.6, indeed follow closely the expected 1/R trend for the case of parallel orientationof the dimer axis in our planar model. We may emphasize that in the real 3D helium dimer we
expect an even stronger 1/R2 decrease. We may further note that the double-to-single ionization
ratio follows the same trend as the double photoionization probability.
Moreover, our time-dependent simulations provide a “movie” of the scattering process as afunction of time. During the knockout process, it is expected that one of the electrons approaches
the other after the initial absorption of the photon. This can be visualized via the temporal evolutionof the two-electron probability distribution as a function of the relative coordinate z = z1 −z2 along
the polarization direction. The respective results of our numerical simulations for an internuclear
distance of 15 a.u. are shown in Fig. 3.7. The initial state probability distribution, localized at
the two nuclei, has been subtracted in order to enhance the visibility of the small contributions
leading to single and double ionization. The results for the temporal evolution starting from theground state (singlet, panel (a)), and the first excited state (triplet, panel (b)) confirm our previous
conclusions regarding the knockout mechanism. The exchange-asymmetry in the wave function of
the triplet state illustrates itself as a minimum at z = 0, which is particularly apparent when thetwo electrons are close to each other. Due to the integration over x1 and x2 , we observe a minimumin the distributions instead of a node, which would show up at fixed transverse distances.

52

Figure 3.7. Temporal evolution of the two-electron probability distribution of the helium dimerin z = z1 − z2 starting from the ground state [singlet, panel (a)] and the first excited state [triplet,
panel (b)], which exhibit the differences between singlet and triplet scattering. In the simulation,
the dimer axis is oriented parallel to the polarization axis.

53

3.5

The Hamiltonian reduction method
In this section, we consider abridged Hamiltonians and study the Coulomb interactions be-

tween the electrons and the nuclei in the helium dimer by comparing the angular distributions
obtained using the abridged Hamiltonian to those obtained with the full Hamiltonian. To this end,
we first briefly discuss the numerical model used, which just restricts the orientation of the dimerin the planar 4D model in Sec. 3.2 to be parallel to the laser field, since the scattering effects insingle and double photoionization are strongest in the parallel orientation case. We then introduce
and discuss different abridged Hamiltonians in which certain interactions are removed from the full
Hamiltonian. A comparison between the ground states of the full and the abridged Hamiltonians
will be given and discussed. Next, we apply the different Hamiltonians in numerical simulations
to study scattering effects of the photoelectron in single photoionization. Afterwards, we compare
angular distributions of the primary and knockout electron as well as to study the role of Coulombinteractions in the double photoionization process.

3.5.1

TheoryIn this section we outline the four-dimensional model Hamiltonian used for the helium dimerin the simulations, which is a restriction of the planar 4D model of the helium dimer in Sec. 3.2
to the z direction. We then present options to restrict the full model Hamiltonian in numericalsimulations in order to label the electrons and analyze the effect of certain interactions on the single
and double photoionization processes.
As shown in Fig. 3.8, we use the planar 4D model of the helium dimer in Sec. 3.2 but restrict
the orientation of the dimer axis to be parallel to the laser polarization direction, or the z direction,
unless explicitly stated otherwise. In this scheme, nucleus 1 (N1 ) is located at the right side andnucleus 2 (N2 ) is located at the left side. Here we restrict the orientation of the dimer axis tobe parallel to the laser field since the scattering effects in single and double photoionization are

54strongest in this case. The full Hamiltonian of the system is given by
p2 + p2z1
p2 + p2z2
H = x1+ x2+ E(t) (z1 + z2 )22
+ VSAE (r11 ) + VSAE (r12 ) + VSAE (r21 ) + VSAE (r22 )(3.10)

1
1
+p
+ ,
2
2
2
R
(x1 − x2 ) + (z1 − z2 ) + b
where
s
rij =



R 2
x2i + zi + (−1)j+ a2 .
2(3.11)

Figure 3.8. Scheme of the numerical model for the helium dimer in the case of parallel orientation. One electron at each atom is considered to be active and restricted to the same plane, withcoordinates (x1 , z1 ) and (x2 , z2 ), respectively. The internuclear distance R and the dimer axis are
fixed in space. The dimer axis is chosen parallel to the laser polarization direction ε during thesimulations, unless explicitly stated otherwise.

For further analysis we used abridged Hamiltonians in test simulations, in which certain
interaction terms are removed from the full model Hamiltonian, Eqn. (3.10). One goal to use these
reductions was to distinguish the primary electron, which first interacts with the light field, fromthe secondary electron by labeling the two electrons. This can be done by removing the term E(t)z2
from the full Hamiltonian, which gives
p2 + p2z1
p2 + p2z2
HE1 = x1
+ x2
+ E(t)z1
2
2+ VSAE (r11 ) + VSAE (r12 ) + VSAE (r21 ) + VSAE (r22 )+p1(x1 − x2 )2 + (z1 − z2 )2 + b2

+

1
.
R

(3.12)

55
In the corresponding test simulations we can then identify electron 1 as the primary and electron
2 as the secondary or knockout electron in double photoionization and, thus, compare the angular
distributions for these electrons.
The present model of the helium dimer further facilitates the analysis of dynamical effects
due to the interactions between electrons and nuclei. As mentioned above, the two electrons are
initially predominantly located at different atoms in the dimer. This can be also seen from thespatial distributions of the lowest energy eigenstates of the planar helium dimer model, which are
shown in Fig. 3.9 for an internuclear separation of R = 5.6 a.u. corresponding to the minimumof the helium dimer potential [91]. These energy eigenstates were obtained via imaginary time
propagation. The spatial distributions are shown as functions of z1 and z2 and are integrated over
the other coordinates. The ground state [panel (a)] is a singlet state, while the first excited state[panel (b)] is a triplet state.
Due to the large internuclear distance, the interaction of each electron with the distant nucleus
in the initial dimer state is rather small. Hence, in two further abridged Hamiltonians, we neglected
first the interaction of electron 2 with nucleus 1, VSAE (r21 ), to get
p2 + p2z1
p2 + p2z2
HE1,local1 = x1
+ x2
+ E(t)z1
2
2+ VSAE (r11 ) + VSAE (r12 ) + VSAE (r22 )(3.13)

1
1
+p
+ ,
2
2
2
R
(x1 − x2 ) + (z1 − z2 ) + b
and then VSAE (r21 ) as well as VSAE (r12 ) to get
p2 + p2z1
p2 + p2z2
HE1,local2 = x1
+ x2
+ E(t)z1
2
2+ VSAE (r11 ) + VSAE (r22 )(3.14)

1
1
+p
+ .
2
2
2
R
(x1 − x2 ) + (z1 − z2 ) + b
In the first of the two Hamiltonians, HE1,local1 , the symmetry of the electron-nucleus interactions
is removed and, therefore, the electrons are localized at one of the two nuclei, i.e., electron 1 at
nucleus 1 and electron 2 at nucleus 2, as can be seen from the ground state distribution shown

56

Figure 3.9. Spatial distributions (on a logarithmic scale) of the ground state [panel (a)] andfirst excited state [panel (b)] of the planar model helium dimer with full model Hamiltonian andthat of the ground states of the planar model helium dimer with abridged Hamiltonians, HE1,local1[Eqn. (3.13), panel (c)] and HE1,local2 [Eqn. (3.14), panel (d)]. The distributions are integrated overx1 and x2 and shown as functions of z1 and z2 .

57
in Fig. 3.2(c). We note that in this Hamiltonian all the interactions of the primary electron, i.e.,electron 1, which interacts with the light field, are retained. On the other hand, in HE1,local2 theinteraction of the primary electron with the second nucleus (nucleus 2) is removed as well and thecorresponding initial state is represented in Fig. 3.2(d). As we will show below, comparison of theresults from the simulations with the different Hamiltonians enables us to study, e.g., the scatteringof the primary electron at the second nucleus (nucleus 2) in the photoelectron angular distributions.
We note that due to the localization of the electrons at either one of the nuclei, the corresponding ground states are no longer symmetric regarding to the exchange of electrons and theterminology of singlet and triplet states does not apply to the states of the two abridged Hamiltonians. The difference between the ground states of the full and the two abridged Hamiltonians canbe quantified via their overlaps as:
singlet 2
i| = 0.499562;|hψE1,local1 |ψfull(3.15)

triplet 2|hψE1,local1 |ψfulli| = 0.499251;(3.16)

singlet 2
i| = 0.498898;|hψE1,local2 |ψfull(3.17)

triplet 2
i| = 0.498805.
|hψE1,local2 |ψfull(3.18)All overlaps are close to 0.5, indicating the small effect of the interaction between an electron,
localized at one atom in the dimer, and the distant nucleus in the initial state. The value of 0.5
instead of 1.0 arises due to the localization of each electron at one of the two nuclei in the abridged
Hamiltonians.

3.5.2

Single ionization

In this section, we use the model Hamiltonians to study scattering effects in single photoionization of the dimer, or more precisely, the dynamical effects of the photoelectron that is emitted
at one of the atoms and then scattered at the other atom in the dimer. To this end, we obtained
the photoelectron angular distributions by taking a Fourier transform of the contributions to thewave function that were identified as single ionization in the numerical simulations. From the re-58

Figure 3.10. Comparison of photoelectron angular distributions in single ionization of the helium
dimer using the full [panel (a)] and abridged Hamiltonians, HE1,local1 [Eqn. (3.13), panel (b)] andHE1,local2 [Eqn. (3.14), panel (c)]. Results represented by solid lines are obtained for parallelorientation of the dimer axis with respect to the polarization direction of the field, which is chosen
along the horizontal axis. For the full simulations [panel (a)] also shown is the distribution obtainedfor perpendicular orientation (θ = 90◦ , blue dashed line) of the dimer axis and for all simulations
the distributions are compared to the expected cos2 -distribution (dash-dotted lines). For the sakeof comparison, the distributions in panel (a) are scaled such that their areas match, while in panels(b) and (c) the distributions are matched at one point. In all simulations a 4-cycle 20 nm XUVlaser pulse with a peak intensity of 1 × 1014 W/cm2 has been used.

59
sults for the full model Hamiltonian, shown in Fig. 3.10(a), we see, as pointed out earlier, that thephotoelectron angular distributions depend on the orientation of the dimer axis with respect to thepolarization direction of the laser. For an orientation perpendicular to the polarization axis, thedistribution (dashed line) is close to the expected p-wave distribution (cos2 , dash-dotted line) for
single-photon absorption from a s-state. On the other hand, for parallel orientation we observe a
deviation from the cos2 -distribution, and the angular distribution appears to be narrower aroundthe polarization direction with two small wings in the region between 30◦ and 60◦ with respect
to the polarization axis. These deviations in the distributions were previously interpreted as anindication for the scattering of the photoelectron at the other atom in the dimer.
To shed further light on this interpretation, we present in the other panels of Fig. 3.10
results obtained with the abridged Hamiltonians, HE1,local1 [Eqn. (3.13), panel (b)] and HE1,local2[Eqn. (3.14), panel (c)], for orientation of the dimer parallel to the polarization direction. In both
set of simulations only one of the two electrons, namely electron 1, interacted with the field and,moreover, this active electron 1 was located at nucleus 1, which is the nucleus located on the rightin the present figures. In this configuration we can compare photoemission with (emission to theleft) and without (emission to the right) scattering at the second atom, which is located to theleft of the initial position of the active electron. The result in panel (b), in which all Coulombinteractions of the active electron with the other electron and the nuclei are taken into account,confirm that the nonscattered angular distribution is close to the cos2 -distribution (dash-dottedline), while the scattered angular distribution shows the deviation observed in the full simulations.
This demonstrates the effect of scattering on the photoelectron angular distribution. We note thatthe small deviations between the numerical results and the cos2 -distributions in the right part ofthe lobe, most notably in the kink along the polarization direction in Fig. 3(c), reflect the numericalerrors in the present calculations.
We can even further distinguish between the effect of electron-electron and electron-ion scattering. This can be done by comparing the results in panels (b) and (c) of Fig. 3.10. In the latter
simulations we deliberately neglected the interaction between the active electron and the second

60
nucleus, i.e., the nucleus on the left. The comparison shows that the angular distribution for the
scattered photoelectron changes. Emission of the electron along the polarization direction is now
strongly suppressed, while there are two maxima at an angle of about 40◦ with respect to thepolarization axis. Thus, we can conclude that it is the electron-electron scattering which gives riseto the two wings in the full distribution. On the other hand, the electron-nucleus (or electron-ion)scattering causes, due to the attractive nature of the interaction, a narrowing (or focusing) to the
angular distribution along the polarization direction.

3.5.3

Double ionization

In this section we analyze the results of simulations using the full and abridged Hamiltonians
to study the correlated electron emission from the model dimer. To this end, we first investigate
the projection onto different approximate two-electron wave functions to obtain the photoelectron
angular distributions. We then proceed to compare the angular distributions of the primary andthe knockout electrons.
In previous studies it has been found [89] that the double photoionization of the helium
dimer proceeds via the so-called knockout mechanism for low photon energies. According to thismechanism, first the primary electron gets ionized after interaction with the field and then shares
its energy via collisional impact ionization with the second knockout electron. As mentioned above
we can distinguish between the two electrons by restricting the electron-field interaction to one of
the two electrons and locating the electrons at one of the atoms in the dimer.
In Fig. 3.11 we present results of simulations with the different abridged model Hamiltonians,
Eqns. (3.13) and (3.14). In the panels in the left column the angular distributions of the primaryelectron, i.e., the electron interacting with the field, are presented, while in the right column thoseof the knockout electron are shown. In all cases the dimer was oriented along the polarization axis,
which is the horizontal axis in each of the panels. We have integrated the distributions over theemission angle of the other electron and the energies of the two electrons.
Due to the location of the electrons in our models we can study the emission direction of the61

Figure 3.11. Angular distributions of the primary (left column) and knockout electron (rightcolumn) in double photoionization of the helium dimer. Results from simulation with different
abridged Hamiltonians are compared: panels (a) and (b) correspond to HE1,local1 [Eqn. (3.13)] andpanels (c) and (d) correspond to HE1,local2 [Eqn. (3.14)]. In each of the simulations, a 4-cycle 20
nm XUV laser pulse with a peak intensity of 1 × 1014 W/cm2 and parallel orientation of the dimer
axis with respect to the polarization direction of the field (along the horizontal axis of the figure)was considered.

62
electrons with respect to their initial location in the dimer as well as the role of electron-nucleus
and electron-electron interaction during the impact ionization. The result in Fig. 3.11(a), obtainedusing HE1,local1 [Eqn. (3.13)], in which all interactions of the primary electron are retained, showsthat the primary electron is scattered both in forward and backward direction with respect tothe propagation direction (from right to left in the figure) after absorption of the photon. Further
restriction by suppression of the electron-nucleus scattering in the energy transfer process, HE1,local2[Eqn. (3.14)], reveals that the electron-electron scattering favors backward scattering of the primaryelectron, as expected in a head-on collision from a repulsive center [Fig. 3.11(c)].
On the other hand the knockout electron is predominantly ionized along the momentum direction of the incoming primary electron in the impact ionization process, c.f., Fig. 3.11(b) and
Fig. 3.11(d). The contribution from the electron-electron interaction leads to a rather broad angular distribution with respect to the molecular axis [Fig. 3.11(d)], in agreement with a similar
distribution of the primary electron in the opposite direction [Fig. 3.11(c)]. The inclusion of theelectron-nuclear contribution leads to a distribution which is more aligned along the polarization
direction [Fig. 3.11(b)].
Finally, we obtained the angular distributions of the knockout electron when the primaryelectron is fixed along certain emission directions. For the results presented in Fig. 3.12, we considered emission of the primary electron either in backward (with respect to the impact direction,
panels on the left) or forward direction (panels on the right), using abridged Hamiltonians HE1,local1[Eqn. (3.13), panels (a) and (b)] and HE1,local2 [Eqn. (3.14), panels (c) and (d)]. The distributions
show that the two electrons are predominantly emitted back-to-back in particular for backward
scattering of the primary electron. As can be seen in the distributions, we observe for backward
scattering that inclusion of electron-nucleus interaction leads to a stronger alignment of the distribution along the direction of electron impact [c.f., Fig. 3.12(a) and Fig. 3.12(c)]. On the other hand,for forward direction we see pronounced wings in the angular distributions of the knockout electron
when HE1,local1 is used where all interactions of the primary electrons are included [Fig. 3.12(b)].
Similar conclusions hold for other emission angles of the primary electron with respect to63

Figure 3.12. Angular distributions of the knockout electron in double photoionization obtainedusing HE1,local1 [Eqn. (3.13), panels (a) and (b)] and HE1,local2 [Eqn. (3.14), panels (c) and (d)].
The emission of the primary electron is fixed along the direction indicated by the arrow. The fieldparameters were the same as in Fig. 3.11.

64

Figure 3.13. Same as Fig. 3.12, but for different emission direction of the primary electron.

65
the molecular axis and the polarization direction as well, as can be seen from the results presented
in Fig. 3.13. In general, we observe that the electrons are predominantly emitted in opposite halfplanes to the polarization direction. For backward scattering of the primary electron, inclusion of
the electron-nucleus interaction leads to a narrower angular distribution of the knockout electron.

3.5.4

Exchange ionization

Using abridged Hamiltonians, we can further study interesting phenomena like the exchange
ionization under electron impact. To this end, we use the abridged Hamiltonian HE1 [Eqn. (3.12)],in which only electron 1 interacts with the XUV laser pulse. In this scenario, electron 1 is theprimary electron during the knockout double photoionization process. In the second step of thisprocess, electron 1 knocks out electron 2 from the other atom. However, it is also possible thatelectron 1 gets captured by the residual ion from which electron 2 is emitted, a process we denoteas exchange ionization.

Figure 3.14. Exchange ionization under electron impact for (a) singlet and (b) triplet scattering
using the abridged Hamiltonian HE1 [Eqn. (3.12)].

The exchange ionization can be clearly seen in Fig. 3.14, which shows a snapshot of theprobability density distribution as a function of z1 and z2 (integrated over x1 and x2 ) of the helium
dimer at t = 12 a.u. after interacting with a 4-cycle XUV laser pulse with central wavelength of 20

66
nm and peak intensity of 1 × 1014 W/cm2 starting from the (a) singlet and (b) triplet state using
the abridged Hamiltonian HE1 [Eqn. (3.12)]. As can be seen in the figure, the two intense stripes
near z2 = 0 represent the single ionization of electron 1 (primary electron), and the two weakerstripes near z1 = 0 stand for the single ionization of electron 2 (secondary electron). Since only
electron 1 interacts with the laser pulse, only it can be ionized in the first step. Thus electron 2
must be knocked out by electron 1 under its impact in the second step, and electron 1 gets captured
at the same time since the vertical stripes near z1 = 0 represent single ionization instead of double
ionization.

3.6

Summary
In conclusion, we have proposed a planar two-active electron model to analyze the recently

observed correlation effects in double photoionization of the helium dimer. Results of numericalsimulations for the angular distribution of the electron are found to be in good agreement withthe experimental data. Furthermore, theoretical predictions for the double ionization probability
as functions of the orientation of the dimer axis as well as the distance of the atoms in the dimer
reveal direct evidence for the knockout mechanism, in which the photon energy is transferred overdistances of several Angstroms in the dimer to be shared between the electrons. Our conclusions
are further supported by an analysis of the temporal evolution of the two-electron probability
distribution. An analogue situation to the one discussed here for the helium dimer is encounteredin the creation of two-site double core holes of molecules [96, 97] discovered recently. Here thephoton is absorbed at one K-shell, e.g., in the C2 H2 molecule [97]. The two K-shells have almost
negligible overlap and their spatial separation is much bigger than the extension of the individual
K-shell itself. Hence we expect that the conclusions drawn from our calculations are also valid forthe case of single-photon-induced two-site double core hole production.
We have also shown how scattering effects in single and double photoionization of the helium
dimer can be studied in ab-initio numerical simulations. Due to the large internuclear separation in
the dimer, the interaction between the electrons and between each electron and the distant nucleus

67
in the ground state of the dimer is small. This has enabled us to consider abridged Hamiltonians for the dimer model, in which these interactions are neglected without changing the initialstate significantly. Based on the comparison of angular distributions of the primary photoelectron
and the secondary knockout electron obtained with different abridged Hamiltonians in numericalsimulations, we were able to identify the role of the different Coulomb interactions. For single
photoionization our results show that the electron correlation interaction leads to scattering of thephotoelectron off the dimer axis, when the latter is oriented along the polarization direction of the
field. On the hand the electron-nucleus interaction tends to focus the electron along the polarization direction. In double photoionization the primary photoelectron is scattered both in forward
as well as backward direction with respect to its initial momentum direction obtained upon photon
absorption. Keeping the direction of the primary electron fixed we then further found that the twoelectrons are primarily emitted back-to-back in the double photoionization process.

Chapter 4
Selection rules in the few-photon double ionization process1

In this chapter, we analyze the selection rules for the emission of two electrons from thehelium atom and the helium dimer following the absorption of a few photons in an intense laser
field. The rules arise, as generalization of the well-studied one-photon case, due to the symmetries of
the accessible final states in the two-electron continuum. We show, in particular, that an increase of
the number of absorbed photons leads to alternating suppression and non-suppression of the backto-back emission of the two electrons with equal energy sharing. Results of numerical simulations
using models of the helium atom and the helium dimer are in agreement with the theoreticalpredictions. We further generalize the selection rules and find that they apply whenever the system
has central symmetry with respect to its center of mass.

4.1

Introduction
The evolution of three or more Coulomb-interacting particles in the continuum after the

break-up of the atomic or molecular system has been a topic of recurring interest over many years.
Despite the complex dynamics of an interacting few-body system, the final-state distributions of theparticles often show some general structures. In particular, prominent nodes in the distributions
represent restrictions for the correlated momenta of the outgoing particles. These restrictions
arise either due to the Coulomb repulsion between the charged particles or due to symmetries(spin, angular momentum, parity, etc.) of the final state of the few-body process. In this respect,
1

The results of this chapter are in part based on “H. Ni, S. Chen, C. Ruiz, and A. Becker, J. Phys. B 44, 175601(2011)” [98].

69
the break-up of an atomic or molecular system following absorption of a single photon has beenstudied in detail. Selection rules for the correlated electron momenta arising from the symmetryof the final state have been established for single-photon double ionization of the helium atom [99]and the hydrogen molecule [100], as well as the general case of a N -electron break-up following
photoionization of the helium atom [101]. In contrast, the selection rules after absorbing multiple
photons have not been studied yet. However, we may expect that the final two-electron statesin few-photon double ionization possess a few definite symmetries and that the well-establishedselection rules in the single-photon case can be applied to the few-photon interactions as well.
In this chapter we study the restrictions to the position and momentum distributions of theemitted electrons in the continuum following the few-photon double ionization of the helium atomand the helium dimer from the ground state. To this end, we first briefly discuss the application of
the final-state selection rules, derived by Briggs and co-workers [99–101], to the few-photon double
ionization. The selection rules based on the symmetry consideration and the correlated Keldysh–
Faisal–Reiss (KFR) theory are discussed as well. We then make use of a model for the interactionof the helium atom [63] and the helium dimer [86] with a laser field to test our expectations innumerical simulations. In our simulations we will, in particular, analyze double ionization pathways
in which up to four photons are absorbed simultaneously by the two electrons to overcome the double
ionization threshold. Based on these studies, we further generalize the selection rules to arbitrary
molecules.

4.2

Selection rules for transitions to two-electron continuum states
The symmetries of a state in the two-electron continuum can lead to nodes in the correspond-

ing configuration space wave function and the correlated momentum distribution [99]. Preparationof a system in an initial state of a given symmetry (e.g., the ground state) as well as the symmetryof the transition operator (e.g., photon absorption) limit the symmetries of the final states that canbe accessed in a specific process [100]. If in a transition all accessible final states share one or more
symmetries or final states of different symmetries lead to the same nodes, there are zeros in the

70differential cross sections of the corresponding process. These selection rules have been accurately
studied in the theory of single-photon double ionization of atoms and molecules and reproduced inexperiments (for a review, see [51]).

4.2.1

Selection rules based on general studies of the two-electron continuum

In order to identify and apply selection rules for double ionization of the helium atom andthe helium dimer by absorption of a few photons, we first briefly review the relevant symmetriesof the final states in the two-electron continuum, as discussed by Briggs and co-workers (see Ref.
[99, 100] and references therein). For the specification of the nodes in the configuration space wave
function and the correlated electron momentum distributions we employ the position vectors r1
and r2 and the momentum vectors k1 and k2 of the two electrons as well as the Jacobi coordinatesr− = r1 − r2 (relative coordinate) and r+ = (r1 + r2 )/2 (center-of-mass coordinate) and thecorresponding momentum vectors k− = (k1 − k2 )/2 and k+ = k1 + k2 . We will consider double
ionization from the ground state of the helium atom and the helium dimer in a linearly polarizedlaser field with polarization axis aligned along the quantization axis of the system, which is chosen
to be the z-axis. Furthermore, energies of the photons will be small enough such that the photon
absorption can be described by an electric dipole transition. Thus, the quantum numbers M = 0(projection of the total angular momentum to the quantization axis) and S = 0 (total spin) do notchange in the cases considered here.
For single-photon double ionization from the ground state of the helium atom and the helium
dimer, which is of 1 Seven symmetry, the transition into the two-electron continuum leads to exclusive
population of the 1 Podd state. In this case and for linear polarization three selection rules have been
identified leading to nodes in the correlated wave function and momentum distribution [99, 100].
We analyze each of these selection rules in view of their application to the absorption of more than
one photon by the atom:(1) There is a node in the two-electron wave function and a corresponding vanishing contribution to the correlated momentum distributions for both r1 and r2 (k1 and k2 ) perpendicular

71
to the quantization axis of the system (here, the polarization axis of the field) if the final
state is of odd parity. Since in dipole transitions the parity of the state changes, this selection rules applies for absorption of either an even or an odd number of photons from thefield depending on the parity of the initial state of the system. Thus, in the case of double
ionization from the ground state of helium atom we expect the corresponding node to show
up for odd-photon processes.
(2) For final states with even S and odd parity (or odd S and even parity) there appears anode in the configuration space wave function (or momentum distribution) for r1 = −r2(or k1 = −k2 ), which is equivalent to r+ = 0 (or k+ = 0). Since the total spin quantum
number S does not change, while the parity of the state does change in an electric dipole
transition, the selection rules applies again either for an even-photon or an odd-photon
process depending here on the total spin and the parity of the initial state. For few-photon
double ionization from the ground state, the node appears, as in the case of selection rule(1), if the number of absorbed photons is odd.
(3) Another node exists for r1 = r2 (k1 = k2 ) and θ1 = π − θ2 for final states with even S, odd
L, and odd parity (or odd S, even L, and even parity). Concerning the total spin quantum
number and the parity the same considerations as in the case of selection rule (2) apply.
Thus, for double ionization from the ground state (even S, even parity and even L) the
node is expected to be present for absorption of an odd number of photons.
In summary, these selection rules should lead to alternating suppression and non-suppressionof electron emission for certain configurations in the two-electron continuum wave function andmomentum distribution as the number of absorbed photons increases beyond one. Our analysis isbased on the assumption that the electron-field coupling is perturbative, which is justified for theparameter regime considered in the numerical simulations below [102]. We may however note thatsimilar conclusions have been drawn based on a S-matrix analysis of double ionization of helium
atom by simultaneous multiphoton absorption in the nonperturbative intensity and wavelength

72
regime [103, 104], which will be discussed below. Please note that we have restricted our discussion
above to those selection rules, which are relevant for the most significant case of double ionization ofhelium atom from the ground state in a linearly polarized laser pulse. This case will be considered in
the numerical simulations below as well. For an initial state with different symmetry or a different
choice of the field polarization, other or further symmetry rules may have to be considered.

4.2.2

Selection rules based on the symmetry of the transition matrix

In the previous section, we outline the selection rules based on the general studies of thetwo-electron continuum by Briggs and co-workers [99, 100]. Indeed, the selection rules can be well
reproduced by studying the symmetry in the transition matrix, which determines the symmetry inthe two-electron continuum based on the symmetry in the initial state and the transition matrix
itself.
The transition matrix for a one-photon process here can be written as
T1 = hψf |E · (r1 + r2 )|ψi i,(4.1)

where ψi is the initial state, ψf is the final state, r1 and r2 are the position vector of the twoelectrons respectively, and E is the electric field of the laser pulse. If the laser pulse is linearlypolarized in the z direction, then the transition matrix for the one-photon process reduces to
T1 = hψf |E(z1 + z2 )|ψi i = 2Ehψf |Z|ψi i,(4.2)

where Z = (z1 + z2 )/2 is the center of mass of the two electrons in the z direction. From the expression for the transition matrix, it is obvious that the final-state wave function needs to be opposite
in Z-symmetry to the initial-state wave function to make the transition matrix T1 nonzero. If theinitial state is the ground state that is symmetric in Z, then the final state must be antisymmetricin Z, i.e., there will be a node in the configuration space wave function at Z = 0, or a node in thecorrelated momentum distribution at pZ = 0.

73
In contrast, the transition matrix for a two-photon process would beT2 =Xhψf |E · (r1 + r2 )|ψ̃ihψ̃|E · (r1 + r2 )|ψi i = 4E 2 hψf |Z 2 |ψi i,

(4.3)

ψ̃

where ψ̃ is the intermediate state. Here, it is evident that the final-state wave function needs tohave the same Z-symmetry as the initial-state wave function to make the transition matrix T2
nonzero. If the initial state is symmetric in Z, then the final state must be symmetric in Z as well,i.e., there is not necessarily a node in the configuration space wave function at Z = 0 (or pZ = 0 in
the correlated momentum distribution).
For higher-photon processes, the pattern goes on in an alternating manner. An increase of
the number of absorbed photons leads to alternating nodes and antinodes of the two-electron wave
function at Z = 0. If the initial state is symmetric in Z, then odd-photon double ionization leads
to a node at Z = 0, and even-photon double ionization results in an antinode at Z = 0.
Note that the initial-state wave function needs to have a definite Z-symmetry for the selection
rules to apply. In a system whose Hamiltonian is symmetric in Z, it can be approved that the
eigenstates of this system should have a definite symmetry in Z. Therefore, the selection rules
apply to the systems whose Hamiltonian is symmetric in Z. For the Hamiltonian of a two-electron
system, this Z-symmetry indeed only exists when X = (x1 + x2 )/2 = 0 and Y = (y1 + y2 )/2 = 0.
As a result, the restriction in the correlated emission of the two electrons in the configuration
space wave function is indeed at X = Y = Z = 0, which is equivalent to r1 = −r2 in selection
rule (2) above. As a result, an increase of the number of absorbed photons leads to an alternating
suppression and non-suppression of the back-to-back emission of the two electrons with equal energy
sharing.

4.2.3

Selection rules based on the correlated Keldysh–Faisal–Reiss theory

The KFR theory [4, 105, 106] is a S-matrix [107] treatment of the above-threshold-ionization
process [108–110], and the same selection rules can be concluded from the correlated version of the
KFR theory [103, 104]. The fully differential cross section for a N -photon process in this theory

74
can be written as
d3 σ (N )2∝ JN(α0 , k1 , k2 ) ,dΩ1 dΩ2 dE2(4.4)

where α0 is the quiver radius, andJN (α0 , k1 , k2 ) =∞X


JN +2m (α0 · (k1 + k2 )) Jm

m=−∞

Up
ω


,

(4.5)

where Up is the ponderomotive potential, ω is the frequency of the laser field, and J is the Bessel
function. In a linearly polarized laser pulse, the quiver radius is along the polarization direction,or the z direction in the present case. Here, JN reduces to∞X

JN (α0 , k1 , k2 ) =
JN +2m (α0 pZ )Jm

m=−∞

Up
ω


,

(4.6)

which further reduces to
JN (α0 , k1 , k2 ) =∞X


JN +2m (0)Jm

m=−∞

Up
ω


,

(4.7)

at pZ = 0. Note that the Bessel function Jn (0) is non-zero only when n = 0, and that N and m
are integers, we therefore have (at pZ = 0)

JN (α0 , k1 , k2 ) =


 



J−N/2 Uωp , (even N )



0.

(4.8)

(odd N )As a result, the same selection rules as discussed above can be found from the correlated KFRtheory, namely there is a suppression of back-to-back emission of the two electrons with equal
energy sharing when odd number of photons are absorbed, while there is no such a suppression foreven number of photons absorbed.

4.3

Selection rules in the helium atom
The general selection rules above apply for the helium atom, the helium dimer, and homoge-

neous diatomic molecules. To study them more thoroughly, we apply numerical model simulations
to the helium atom and the helium dimer for more details. In this section, we focus on the selection
rules in the helium atom. To this end, we first introduce a restricted 3D model of the helium atomand then numerically study the selection rules in this model.

75
4.3.1

The restricted 3D model of the helium atom

Following earlier achievements towards the solution of the full TDSE for the helium atominteracting with an intense laser pulse in view of an analysis of nonsequential double ionization at
optical and near-infrared laser wavelengths (see [111] and references therein), in recent years muchtheoretical interest has been devoted to the double ionization of atoms by absorption of two photons(see [112–116] for most recent studies and references therein). In the case of the two-photon process,results for the total cross sections as well as for momentum distributions have been reported. In
particular, the passage from a simultaneous absorption of the photons (and simultaneous ejectionof the electrons) to sequential double ionization, in which absorption of the first photon generates
a singly charged ion before absorption of the second photon leads to the removal of the second
electron, has been studied.
We are interested in analyzing the occurrence of nodes for particular configurations of theposition and momentum vectors in the probability density and the momentum distributions innumerical simulations. To this end, we may make use of any (lower-dimensional) model, which
exhibits all or some of the above selection rules. Here, we employ of a model for the helium atom
in which the center-of-mass motion of the two electrons is restricted to the polarization axis of
the linearly polarized external laser field. The Hamiltonian of this two-electron model (in Jacobi
coordinates) is given by [63]H(ρ, z, Z; t) =p2ρ + p2z +p2Z
pZ AZ (t)
1
−
+p
2
4
c
ρ + z2
2

2

(4.9)

−p
−p
,
ρ2 /4 + (Z + z/2)2 + a2
ρ2 /4 + (Z − z/2)2 + a2
where a2 is a soft-core parameter to smooth the Coulomb singularity.
Our model, as shown in Fig. 4.1, restricts the center-of-mass motion of the two electrons along
the polarization direction, but does not restrict the symmetry of the dipole transition operator.
Since a linearly polarized laser field interacts with the center-of-mass of the two electrons along thepolarization axis, the symmetry of the dipole transition operator is present in the model. Thus,
the considerations regarding the symmetries of the accessible final states discussed above apply in

76

Figure 4.1. The restricted 3D model of the helium atom. The blue dot represents the helium
nucleus, and the two green dots stand for the two electrons in the helium atom, the center of massof which is fixed along the polarization direction ε (the z direction) of the laser field.

this model as well. In general, in our model the momentum components of the individual electron
along the polarization direction differ and therefore both equal and unequal energy sharings occur.
However, along Z = 0 and pZ = 0 only equal energy sharing is present in the model. This coincides
with the kinematics for the nodes according to selection rules (2) and (3), which are predicted tooccur for equal energy sharing. Thus, for these two selection rules the model captures the relevant
dynamics. Selection rule (1) can occur for special cases of unequal momentum components in thetransverse direction as well, which is not captured by the present model and is a restriction in viewof the present analysis. Moreover, in our model the node arising due to selection rule (3) coincides
with that due to selection rule (2). Thus, in the remainder of the article, we will mainly discuss
our results in view of selection rule (2).
In the actual computations the TDSE of the two-electron model is solved using the Crank–
Nicolson method. The initial ground state wave function ψ(ρ, z, Z; t = 0) is computed throughimaginary time propagation using a grid size of ∆ρ = ∆z = ∆Z = 0.3, and a2 = 0.135. The
field-free ground state energies of helium atom and helium ion are found to be −2.937 and −1.985
respectively, which are both close to the real values of −2.904 and −2. The propagation of the wave
function in the field is conducted using a box size with Nρ = 300, Nz = 600, and NZ = 300 points
0|in the respective directions. An absorbing boundary of the form cos1/6 π2 |x−xwith |x| > |x0 |,Lspanning 10% of the box size in a certain direction, is used, where x0 denotes the border of the

77boundary region and L is its width.
In this following sections, we present the results of a series of numerical simulations basedon the model presented above. In our calculations we change the wavelength and peak intensityof the laser pulse such that the double ionization process can proceed via simultaneous absorptionof one up to four photons. The wavelengths are chosen such that the probability for sequential
double ionization is small (or negligible) in each of the cases considered, as shown schematically in
Fig. 4.2. In each of the computations we use a sin2 -pulse with a total length of six cycles.
Continuum

0

hν3
hν2

−0.448hν1−1.985He 2+

He* + (2s)hν3He + (1s)hν2hν3−2.937He(1s 2)

Figure 4.2. The photoabsorption scheme in the restricted 3D model of the helium atom. In thescheme, hν1 , hν2 , and hν3 correspond to one-, two-, three-photon double ionization, respectively.
The wavelengths used in these three cases are λ1 = 14.32 nm, λ2 = 28.64 nm, and λ3 = 42.96 nm,
respectively.

4.3.2

One-photon double ionization

First, we consider a central laser wavelength of 14.32 nm and a peak intensity of I0 = 1 × 1014
W/cm2 . At this wavelength the photon energy of 3.182 a.u. exceeds the double ionization thresholdof our model He atom and the two electrons can be emitted after absorption of a single photon.
Fig. 4.3 shows snapshots of the probability density distribution |ψ(ρ, z, Z; t)|2 integrated over ρ (left78

Figure 4.3. Probability density distribution (on a logarithmic scale) as a function of z and Z(integrated over ρ, left hand column) and as a function of ρ and Z (integrated over z, right handcolumn) following single-photon double ionization. In the results, presented in the panels in thelower row, the ground state wave function was removed before calculating the probability density
distribution. The snapshots are taken at t = 10 a.u. after the end of the laser pulses. Laserparameters: wavelength λ = 14.32 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total pulse length
6 cycles.

79

Figure 4.4. Same as Fig. 4.3 but at times t = 20 a.u. (upper row) and t = 30 a.u. (lower row)
after the end of the laser pulses. Shown are the probability density distributions after removal of
the He ground state from the full wave function.

80
hand panels) or z (right hand panels) at time t = 10 a.u. after the end of the pulse. Plots in thecolumn on the left show distributions in the Z −z plane (integrated over the ρ-coordinate). In theseplots the contributions at the center (z ≈ Z ≈ 0) mainly correspond to the remaining population in
the neutral helium atom after the interaction. This can be clearly seen from the comparison of theplots in the upper and lower panel in the left hand column, since the He ground state contribution
is removed in the results in the lower panel. The single ionized population is displayed along the
diagonals, while the contributions to double ionization can be found in the regions in between
the diagonals. We clearly see a node along Z = 0 in the distribution, as predicted by selection
rule (2) and shown by Briggs and co-workers before [99]. In previous analysis of single-photon
double ionization, usually asymptotic final states (in time-independent S-matrix calculations) wereconsidered, the present results offer a time-dependent view. It is seen from the results in Figs. 4.3
and 4.4 (the latter shows the distributions at two later time instants) that the node is present assoon as the pulse is over.
The plots in the right hand columns of Figs. 4.3 and 4.4 offer the complementary view ofthe Z − ρ distributions, integrated over z. While the distributions at small ρ are mainly due tocontributions from the neutral He atom and the He+ ion, the population at large ρ indicates ahighly correlated process, i.e., the double ionization. Again, in all of the plots we see a node at
Z = 0 for the contribution belonging to double ionization, as expected from the selection rule (2)and the earlier work on single-photon double ionization [99].

4.3.3

Two-photon double ionization

Next, we consider two-photon double ionization, for which we expect that the node at Z = 0
is not present according to the selection rules. Simultaneous absorption of two photons does occur
either as a second above-threshold double ionization peak (ATDI) [103, 104] at higher intensities
in the wavelength regime considered above, or at longer wavelengths at which the photon energy
is smaller than the double ionization threshold. In our simulations we have analyzed both cases.
First, we present in Fig. 4.5 results at the same wavelength as above but at a higher peak

81

Figure 4.5. Probability density distribution (on a logarithmic scale) as a function of z and Z(integrated over ρ, left hand column) and as a function of ρ and Z (integrated over z, right handcolumn). The ground state wave function was removed before calculating the probability density
distribution. The snapshots are taken at the end of the pulse (upper row) and t = 10 a.u. after
the pulse (lower row). Laser parameters: wavelength λ = 14.32 nm, peak intensity I0 = 1 × 1016
W/cm2 , and total pulse length 6 cycles.

82
intensity of I0 = 1 × 1016 W/cm2 . While in the z − Z distribution taken at the end of the pulse(panel in upper left corner) the evolving population in the double ionization region does not show
a clear structure, we can distinguish the contributions from single-photon and two-photon double
ionization at some later time after the end of the pulse (t = 10 a.u., panel in lower left corner). Please
note that due to the difference in total kinetic energy the single-photon contribution propagates ata smaller velocity than the two-photon contribution. Consequently, we find the latter contribution
at larger distances from the center of the z − Z distribution (i.e., the nucleus) than the one-photon
distribution. Our expectations for the contribution along Z = 0 are obviously realized in the
numerical results for the z − Z distribution (lower panel in left hand column): The absorption of
two photons results in a maximum in the corresponding part of the distribution, while the nodein the single-photon contribution is still present. In the Z − ρ distributions (panels in right handcolumn) there is no clear separation of the two processes possible and, hence, we observe a small
population along Z = 0 due to the two-photon process.

Figure 4.6. Same as Fig. 4.5, but at I0 = 1 × 1017 W/cm2 and t = 7.5 a.u. after the end of thepulse.

Similar conclusions hold for the results at an even higher intensity of I0 = 1 × 1017 W/cm2 ,shown in Fig. 4.6. Here the contribution from the two-photon absorption (2nd ATDI peak) with

83

Figure 4.7. Probability density distribution (on a logarithmic scale) as a function of z and Z(integrated over ρ, left hand column) and as a function of ρ and Z (integrated over z, right handcolumn). The snapshots are taken t = 30 a.u. after the pulse. Laser parameters: wavelength
λ = 28.64 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total pulse length 6 cycles.

84
a maximum along Z = 0 in the Z − z distribution is much stronger, as expected due to thenonlinearity of the process. A clear separation of the single-photon process from the two-photon
process is however not longer possible on the grid size used in the present simulations.
Next, we double the central wavelength of the laser field to λ = 28.64 nm (photon energy of
1.591 a.u.). At this wavelength the (simultaneous) absorption of two photons is needed to overcome
the double ionization threshold. Please note that the photon energy has been chosen such that acompeting sequential double ionization process requires the absorption of at least three photons
and, hence, the corresponding contribution to the two-electron configuration space wave function isstrongly suppressed. The results in Fig. 4.7 do not show a node along Z = 0 either in the Z − z orin the Z − ρ distribution, in agreement with the selection rule (2) for a two-photon process from theHe ground state. Instead, we observe as in the previous results in Figs. 4.5 and 4.6 a maximum inthe contributions from double ionization along this axis in the Z − z distribution. Please note thatthe maximum indicates a preferential back-to-back emission of the two electrons at the present laserparameters, in agreement with recent observations in experiments on two-photon double ionizationof He [117] and Ne [118], as well as recent theoretical predictions [119].

4.3.4

Three- and four-photon double ionization

Finally, we consider double ionization processes via the simultaneous absorption of more than
two photons in order to verify (or, disprove) the alternating appearance of suppression and nonsuppression of electron emission along the Z = 0 axis in the configuration space wave function for
double ionization as the number of absorbed photons increases. We found that at longer wavelengths
the contributions of competing processes such as sequential double ionization are usually of similar(or even larger) strengths than the probability for the ejection of two electrons via simultaneousabsorption of photons. However, at specific wavelengths (and peak intensities), here 42.96 nm and57.28 nm, we observe double ionization which arises from simultaneous three- and four-photon
absorption. The results of the numerical simulations (see Fig. 4.8) clearly exhibit the suppression(for the odd-number-photon process) and the non-suppression (for the even-number-photon process)85

Figure 4.8. Probability density distribution (on a logarithmic scale) as a function of z and Z(integrated over ρ). The snapshots are taken t = 12.25 a.u. after the pulse. Laser parameters:
wavelength λ = 42.96 nm (left hand column) and 57.28 nm (right hand column), peak intensity
I0 = 5 × 1014 W/cm2 , and total pulse length 6 cycles.

86of electron emission along the Z = 0 axis in the Z − z distributions, as expected from the selection
rule (2).

4.3.5

Momentum distributions

As of now, we have demonstrated the selection rules in the configuration space of the twoelectrons. At the same time, the selection rules also apply to the momentum space distributions. As
discussed above: When there is a node in the configuration space, there should be a corresponding
node in the momentum space; and when there is an antinode in the configuration space, there should
also be a corresponding antinode in the momentum space. Due to the lower dimensionality of ournumerical model we cannot provide quantitative predictions for the full momentum distributions
but test certain aspects, here the existence of nodes.
Fig. 4.9 shows correlated momentum distributions for the two electrons in the pZ − pz plane(upper row, integrated over pρ direction) and along the pZ axis (lower row, integrated over both pz
and pρ ). The momentum distributions are obtained by first separating the doubly ionized part ofthe configuration space wave function from the rest. To this end, we partition the coordinate space(refer to Sec. 2.7) as [63]
r1 6 12 ∩ r2 6 12 :He atom(r1 6 6 ∩ r2 > 12) ∪ (r1 > 12 ∩ r2 6 6) : He+ ion







complementary space :He2+ ion,where r1,2 =(4.10)pρ2 /4 + (Z ± z/2)2 . Then, we project the doubly ionized part onto the final-state

wave function, which we approximate by a product state of a Coulomb wave in the relative coordinates ρ and z and a plane wave in the center-of-mass coordinate Z.
The momentum distributions displayed in the three columns correspond to the probability
density distributions(a) in the lower row of Fig. 4.4, in which the single-photon process dominates,

87

Figure 4.9. Correlated momentum distribution of pz and pZ (upper row) and momentum distribution in pZ (lower row, integrated from upper row over pz ) by projection onto the approximate
final-state wave function of a Coulomb wave in the relative coordinates ρ and z and a plane wave inthe center-of-mass coordinate Z. Laser parameters: wavelength λ = 14.32 nm [column (a) and (b)]and 28.64 nm [column (c)], peak intensity I0 = 1 × 1014 W/cm2 [column (a) and (c)] and 1 × 1016
W/cm2 [column (b)], and total pulse length 6 cycles.

88(b) in the lower row of in Fig. 4.5, in which both one- and two-photon processes are present,and(c) in Fig. 4.7, in which the two-photon process is prominent.
It is obvious from the comparison of the results in Fig. 4.9 with the corresponding spatial distributions that both show the same nodes and antinodes. In particular, the nodes for the odd-number
photon processes along pZ = 0 in the pZ − pz distributions (upper row of Fig. 4.9) can be clearly
seen. We have found that the appearance of the nodes is insensitive to the method by which we
extract the ionized population from the full wave function. It is also independent of how long thewave function is propagated after the end of the pulse and of the choice of the basis set for thetwo-electron continuum wave functions.
−1

Contribution of 2−photon process

10

−2

10

−3

10

−4

10

15

10

16

10
Intensity (W/cm2)

17

10

Figure 4.10. Ratio of the probability at pZ = 0 to the maximum probability in the pZ distributions
as a function of intensity at a wavelength of 14.32 nm.

Comparing the momentum distributions in columns (a) and (b) we find that the contribution of the two-photon process increases with intensity at the wavelength of 14.32 nm. This canbe quantified by the ratio of the probability at pZ = 0 to the maximum probability in the pZ
distributions (lower row in Fig. 4.9), as shown in Fig. 4.10. One can see that the ratio increases

89
almost linearly as a function of intensity. This is expected since the ratio basically represents theratio between the probabilities for one- to two-photon absorption, which should follow the powerlaw I 2 /I 1 = I 1 .
Instead of using the approximate final-state wave function of a Coulomb wave in the relative
coordinates ρ and z and a plane wave in the center-of-mass coordinate Z to get the momentum
distributions, we can also use the approximate final-state wave function of a plane wave, which isequivalent to carrying out a Fourier transform to the respective spatial wave function to obtainthe momentum distributions. Shown in Fig. 4.11 is the same plot as in Fig. 4.9 but is obtained byprojection onto a plane wave (Fourier transform).
The difference between these two methods, i.e., projection onto the approximate final-state
wave function of a Coulomb wave in the relative coordinates ρ and z and a plane wave in the centerof-mass coordinate Z (Fig. 4.9) and projection onto the approximate final-state wave function ofa plane wave (Fig. 4.11), is that electron-electron correlation is taken into consideration in theformer case, and it is not accounted for in the latter case. The electron-electron correlation, or
electron-electron repulsion, will cause the two electrons to stay away from each other, which would
result in a node at ρ = z = 0 (and pρ = pz = 0). This can be clearly seen in Fig. 4.12, in which
the first row is the cut of Fig. 4.9 and the second row is the cut of Fig. 4.11 at pρ = 0. A node atpz = 0 is evident in the first row, in which electron-electron correlation is taken into account, and
it is missing in the second row, where electron-electron repulsion is absent.

4.4

Selection rules in the helium atom within the planar 4D model
In the last section, we studied the selection rules in the helium atom using the restricted 3Dmodel. To illustrate that the selection rules do not depend on the specific model, we introduce aplanar 4D model of the helium atom and then study part of the selection rules within this model.

90

Figure 4.11. Correlated momentum distribution of pz and pZ (upper row) and momentum distribution in pZ (lower row, integrated from upper row over pz ) by projection onto the approximate
final-state wave function of a plane wave (Fourier transform). Same laser parameters are used asin Fig. 4.9.

91

Figure 4.12. Correlated momentum distribution of pz and pZ cut at pρ = 0 by projection onto theapproximate final-state wave function of a Coulomb wave in the relative coordinates ρ and z and aplane wave in the center-of-mass coordinate Z (upper row) and by projection onto the approximate
final-state wave function of a plane wave (Fourier transform, lower row). Same laser parameters
are used as in Figs. 4.9 and 4.11.

92
4.4.1

The planar 4D model of the helium atom

The planar 4D model of the helium atom removes the restriction of the center of mass of thetwo electrons to the polarization direction of the laser field in the restricted 3D model of the atom,
but restrains the two electrons and the helium nucleus to the same plane. After this change, eachelectron has two spatial dimensions and the corresponding four-dimensional model Hamiltonian isthen given by
p2 + p2z1
p2 + p2z2
H(x1 , x2 , z1 , z2 ; t) = x1
+ x2
+ E(t)(z1 + z2 )22
(4.11)
2
2
1
p
p
p
,
−
−
+
(x1 − x2 )2 + (z1 − z2 )2 + b2
x21 + z12 + a2
x22 + z22 + a2
where pi = (pxi , pzi ) and ri = (xi , zi ) (with i = 1, 2) are the momentum operators and spatial
coordinates of the two electrons and a2 and b2 are the soft-core parameters to smooth the Coulomb
singularity. The linearly polarized laser field is given in the form ofE(t) = cos2



ωt
2N


cos(ωt + φ),(4.12)

where ω is the central frequency, N is the number of cycles, and φ is the carrier-envelope phase of
the laser pulse, which is set to −π/2 is the present simulations.
In the actual computations the TDSE of the two-electron model is solved using the Crank–
Nicolson method. The initial ground state wave function ψ(x1 , x2 , z1 , z2 ; t = 0) is computed throughimaginary time propagation using a grid size of ∆x1 = ∆x2 = ∆z1 = ∆z2 = 0.3, a2 = 0.164, andb2 = 0.068, which is found to have an energy of −2.9034. The propagation of the wave functionin the field is conducted using a box with Nx1 = Nx2 = 200 and Nz1 = Nz2 = 300 points in
0|
the respective directions. An absorbing boundary of the form cos1/6 π2 |x−xwith |x| > |x0 |,Lspanning 10% of the box size in a certain direction, is used, where x0 denotes the border of theboundary region and L is its width. The double ionization contributions are obtained with thespace-partitioning method in Sec. 2.7 for an atomic system with dsmall = 7 and dbig = 11.

93
4.4.2

One-photon double ionization

First, we consider a central laser wavelength of 13 nm and a peak intensity of I0 = 1 ×1014 W/cm2 . At this wavelength the photon energy of 3.507 a.u. exceeds the double ionizationthreshold of our model helium atom and the two electrons can be emitted after absorption of a single
photon. Fig. 4.13(a) shows the snapshot of the probability density distribution |ψ(x1 , x2 , z1 , z2 ; t)|2
integrated over x1 and x2 at time t = 4 a.u. after the end of the pulse. However, no node is visible
along z1 + z2 = 0, or Z = 0. This can indeed be easily understood as due to the symmetry of the
Hamiltonian, which determines the symmetry of the initial state. To see this clearly, we express
the Hamiltonian in the center-of-mass (X and Z) and relative coordinates (x and z) of the twoelectrons
p2X + p2Z
1+√2
4x + z 2 + b2
2
2
−p
−p
,
x 2
x 2
z 2
2
(X + 2 ) + (Z + 2 ) + a(X − 2 ) + (Z − z2 )2 + a2

H(x, X, z, Z; 0) =p2x + p2z +(4.13)

where X = (x1 + x2 )/2, Z = (z1 + z2 )/2, x = x1 − x2 , and z = z1 − z2 . Obviously, the Hamiltonian
is not symmetric in Z, since it changes if we let Z → −Z. This means that there is no symmetry
with respect to Z for the initial state either, which thus, according to the arguments in Sec. 4.2.2,
tells us that there should not be a node at Z = 0 for the final state either. However, the node canbe retrieved if we consider the probability distribution at X = 0, x = 0, or z = 0. In the first case,this leads to the following Hamiltonianp2X + p2Z
1+√2
4x + z 2 + b2
22
−p x 2−p x 2,z 2
2( 2 ) + (Z + 2 ) + a( 2 ) + (Z − z2 )2 + a2

H(x, X, z, Z; 0) =p2x + p2z +(4.14)

which is symmetric in Z. The respective probability density distribution is shown in Fig. 4.13(b),which clearly has a node at Z = 0. In the restricted 3D model this restriction is always imposed,
as discussed before.

94

Figure 4.13. Probability density distribution (on a logarithmic scale) as a function of z1 and z2(integrated over x1 and x2 ) following single-photon double ionization. The distribution in panel(b) represents a cut at X = 0. The snapshots are taken at t = 4 a.u. after the end of the laser
pulse. Laser parameters: wavelength λ = 13 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total
pulse length 12 cycles.

Figure 4.14. Probability density distribution (on a logarithmic scale) as a function of z1 and z2(integrated over x1 and x2 ) following two-photon double ionization. The distribution in panel (b)represents a cut at X = 0. The snapshots are taken at t = 4 a.u. after the end of the laser pulse.
Laser parameters: wavelength λ = 26 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total pulse
length 6 cycles.

95
4.4.3

Two-photon double ionization

Next, we consider two-photon double ionization, for which we expect that the node at Z = 0
is not present according to the selection rules. Here, we only consider the case when the photon
energy is smaller than the double ionization threshold and a simultaneous absorption of two photons
is necessary to produce double ionization. To this end, we double the central wavelength of thelaser field to λ = 26 nm (photon energy of 1.753 a.u.). Indeed, the results in Fig. 4.14 do not show
a node along Z = 0 either in the probability distribution integrated over x1 and x2 or in that which
is further restricted to X = 0, in agreement with the selection rule (2) for a two-photon processfrom the helium ground state.

4.4.4

Momentum distributions

The selection rules, just as in the restricted 3D model, also apply to the momentum-space
distributions. When there is a node in the configuration space, there should be a corresponding
node in the momentum space, and when there is an antinode in the configuration space, there
should also be a corresponding antinode in the momentum space.
Fig. 4.15 shows correlated momentum distributions for the two electrons in the pz1 − pz2
plane for the single-photon process (upper row) and for the two-photon process (lower row). The
momentum distributions are obtained by projecting the doubly ionized part onto the final-state
wave function, which we approximate by a product state of a Coulomb wave in the relative coordinates x and z and a plane wave in the center-of-mass coordinates X and Z. Then an integration
over px1 and px2 results in the momentum-space distribution in the pz1 − pz2 plane. Note that the
distributions in the right column further restrict the distributions to pX = 0, just like we restricted
to X = 0 in the spatial distributions in order to retrieve the node in odd-photon double ionization.
Obviously, there is a node present at pZ = 0 for the single-photon process when the pX = 0 condition is enforced, while the node is missing for the two-photon process even when the cut at pX = 0
is done.

96

Figure 4.15. Correlated momentum distribution of pz1 and pz2 for the single-photon double ionization (upper row) and two-photon double ionization (lower row) by projection onto the approximate
final-state wave function of a Coulomb wave in the relative coordinates x and z and a plane wave
in the center-of-mass coordinates X and Z. The distributions in the right column represent a cut
at pX = 0. Laser parameters: wavelength λ = 13 nm (upper row) and 26 nm (lower row), peakintensity I0 = 1 × 1014 W/cm2 , and total pulse length 12 cycles (upper row) and 6 cycles (lower
row).

97

Figure 4.16. Momentum distributions of electron 2 for the single-photon double ionization process
when electron 1 is fixed to be emitted along the +z direction. In panels (b), (c), and (d), themagnitude of the momentum of electron 1 is further restricted to certain values, i.e., 0.4, 0.8, and1.2, respectively. Laser parameters: wavelength λ = 13 nm, peak intensity I0 = 1 × 1014 W/cm2 ,and total pulse length 12 cycles.

98
The momentum-space wave function further provides the opportunity to reveal the selection
rules in a more straightforward way, i.e., by fixing the momentum of one electron and observing the
momentum distribution of the other electron. Shown in Fig. 4.16(a) is the probability distributionof electron 2 when the emission of electron 1 is fixed along the +z direction, which is marked as anarrow, for the single-photon double ionization process. There is no node in the opposite directionof the momentum of electron 1. If we however further select the magnitude of the momentum ofelectron 1, as it has been done for the results presented in the other panels of Fig. 4.16, we findthat, electron 2 cannot be emitted back-to-back to electron 1 for equal energy sharing.
Thus, the selection rules in the few-photon double ionization of the helium atom do notdepend on the specific numerical model we are using, and the selection rules apply in bothconfiguration-space and momentum-space wave functions. Generally, the two electrons cannot
be emitted back-to-back with equal energy sharing in an odd-photon process, but that is allowed
in an even-photon process.

4.5

Selection rules in the helium dimer
In the last sections, we studied the selection rules in the helium atom. In this section, we

focus on the selection rules in the helium dimer using the planar 4D model introduced in Sec. 3.2.

4.5.1

One-photon double ionization

First, we consider a central laser wavelength of 20 nm and a peak intensity of I0 = 1 ×1014 W/cm2 . At this wavelength the photon energy of 2.279 a.u. exceeds the double ionizationthreshold of our model helium dimer and the two electrons can be emitted after absorption of asingle photon. Fig. 4.17(a) and Fig. 4.17(c) shows snapshots of the probability density distribution|ψ(x1 , x2 , z1 , z2 ; t)|2 integrated over x1 and x2 at time t = 4 a.u. after the end of the pulse. Again,
no node is visible along z1 +z2 = 0, or Z = 0. Similarly, this can indeed be easily understood as dueto the symmetry of the Hamiltonian, which determines the symmetry of the initial state. To see
this clearly, we express the Hamiltonian in the center-of-mass (X and Z) and relative coordinates

99

Figure 4.17. Probability density distribution (on a logarithmic scale) as a function of z1 andz2 (integrated over x1 and x2 ) following single-photon double ionization starting from the singlet(upper row) and triplet (lower row) states. The distributions in the right column represent a cut
at X = 0. The snapshots are taken at t = 4 a.u. after the end of the laser pulse. Laser parameters:
wavelength λ = 20 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total pulse length 12 cycles.

100(x and z) of the two electrons
p2 + p2Z
1
H(x, X, z, Z; 0) =p2x + p2z + X+√
2
4
x + z 2 + b2
s


2 
2
x R
z R
X + − sin θ + Z + − cos θ + a2 
+ VSAE 
2
2
2
2

s

2 
2
x R
z R
X + + sin θ + Z + + cos θ + a2 
+ VSAE 
2
2
2
2

s

2 
2
x R
z R
X − − sin θ + Z − − cos θ + a2 
+ VSAE 
2
2
2
2
s


2 
2
R
R
x
z
X − + sin θ + Z − + cos θ + a2  ,+ VSAE 22
2
2(4.15)

where X = (x1 + x2 )/2, Z = (z1 + z2 )/2, x = x1 − x2 , and z = z1 − z2 . As in the planar He case, theHamiltonian is not symmetric in Z, since it changes if we let Z → −Z. This means that there is nosymmetry with respect to Z for the initial state and therefore no node at Z = 0 for the final state
either. However, as before the node can be retrieved if we consider the probability distribution at
X = 0, x = 0, or z = 0. For example, for X = 0, the Hamiltonian becomes
p2 + p2Z
1
H(x, X, z, Z; 0) =p2x + p2z + X+√
4x2 + z 2 + b2

s

2 
2
x R
z R
+ VSAE − sin θ + Z + − cos θ + a2 
2
2
2
2
s


2 
2
x
R
z
R
+ VSAE + sin θ + Z + + cos θ + a2 
2
2
2
2
s


2 
2
x
R
z
R
+ VSAE + sin θ + Z − − cos θ + a2 
2
2
2
2
s


2 
2
z R
x R
+ VSAE − sin θ + Z − + cos θ + a2  ,2
2
2
2(4.16)

which is symmetric in Z. The respective probability density distribution is shown in Fig. 4.17(b)and Fig. 4.17(d), the former of which clearly has a node at Z = 0. Note that the node (minimum)at z = 0 is due to the exchange asymmetry in the two-electron wave function.

101
4.5.2

Two-photon double ionization

Next, we consider two-photon double ionization, for which we expect that the node at Z = 0
is not present according to the selection rules. As before, we consider only the case when the energyof the photon is smaller than the double ionization threshold and a simultaneous absorption of two
photons is necessary to produce double ionization. To this end, we double the central wavelengthof the laser field to λ = 40 nm (photon energy of 1.140 a.u.). The results in Fig. 4.18 do not show anode along Z = 0 either in the probability distribution integrated over x1 and x2 or in that whichis restricted to X = 0, in agreement with the selection rule (2) for a two-photon process from thedimer ground state.

4.5.3

Momentum distributions

The selection rules, just as in the case of the helium atom, also apply to the momentumspace distributions. Fig. 4.19 shows correlated momentum distributions for the two electrons in thepz1 − pz2 plane for single-photon double ionization. The momentum distributions are obtained byprojecting the doubly ionized part onto the final-state wave function, which we approximate by a
product state of a Coulomb wave in the relative coordinates x and z and a plane wave in the centerof-mass coordinates X and Z. Then an integration over px1 and px2 results in the momentum-space
distribution in the pz1 − pz2 plane. Note that for the distributions in the right column, we havefurther restricted the distributions to pX = 0. Again, there is a node present at pZ = 0 for the
single-photon double ionization process when the pX = 0 condition is enforced.
Fig. 4.20, in contrast, shows correlated momentum distributions for the two electrons in thepz1 − pz2 plane for the two-photon double ionization and there is, as expected, no node present at
pZ = 0 for the two-photon double ionization process even when the pX = 0 condition is enforced.

102

Figure 4.18. Probability density distribution (on a logarithmic scale) as a function of z1 and z2(integrated over x1 and x2 ) following two-photon double ionization starting from the singlet (upperrow) and triplet (lower row) states. The distributions in the right column represent a cut at X = 0.
The snapshots are taken at t = 4 a.u. after the end of the laser pulse. Laser parameters: wavelength
λ = 40 nm, peak intensity I0 = 1 × 1014 W/cm2 , and total pulse length 6 cycles.

103

Figure 4.19. Correlated momentum distribution of pz1 and pz2 for the single-photon double
ionization double ionization by projection onto the approximate final-state wave function of a
Coulomb wave in the relative coordinates x and z and a plane wave in the center-of-mass coordinates
X and Z starting from the singlet (upper row) and triplet (lower row) states. The distributionsin the right column represent a cut at pX = 0. Laser parameters: wavelength λ = 20 nm, peakintensity I0 = 1 × 1014 W/cm2 , and total pulse length 12 cycles.

104

Figure 4.20. Correlated momentum distribution of pz1 and pz2 for the two-photon double ionization double ionization by projection onto the approximate final-state wave function of a Coulomb
wave in the relative coordinates x and z and a plane wave in the center-of-mass coordinates X andZ starting from the singlet (upper row) and triplet (lower row) states. The distributions in theright column represent a cut at pX = 0. Laser parameters: wavelength λ = 40 nm, peak intensity
I0 = 1 × 1014 W/cm2 , and total pulse length 6 cycles.

105

4.6

General molecules
In previous sections, we focused on one- and two-center systems. In this section, we discuss thesymmetry and selection rules in general multicenter systems during few-photon double ionization.
To this end, we illustrate the procedure for a homogeneous N -center system. The corresponding
Hamiltonian, using center of mass and relative coordinates of the two electrons, can be written as
p2X + p2Z1+p4(x1 − x2 )2 + (z1 − z2 )2 + b2

s
2 
2

N
X
Z
X
+
VSAE 
+ x − rj sin θj ++ z − rj cos θj + a2 
2
2
j=1
s


2 
2
N
X
X
Z
+
VSAE 
− x − rj sin θj +− z − rj cos θj + a2  .
22

H(x, X, z, Z; 0) =p2x + p2z +(4.17)

j=1

When we change Z to −Z, the Hamiltonian becomesp2X + p2Z1+p24(x1 − x2 ) + (z1 − z2 )2 + b2
s


2 
2
N
X
X
Z
+
+ x − rj sin θj +− z + rj cos θj + a2 
VSAE 
2
2
j=1

s

2 
2
N
X
X
Z
− x − rj sin θj ++ z + rj cos θj + a2  .
+VSAE 
2
2

H 0 (x, X, z, Z; 0) =p2x + p2z +(4.18)

j=1

Clearly, there is no Z symmetry in this Hamiltonian since it changes after Z → −Z. However,
if we set X = 0, x = 0, or z = 0, and at the same time enforce an central symmetry of the systemwith respect to the origin (i.e., the center of mass), then there is Z symmetry in the Hamiltonian.
We illustrate for a four-center system why this central symmetry is necessary. The corresponding

106
Hamiltonian for X = 0 is
p2X + p2Z1+p4(x1 − x2 )2 + (z1 − z2 )2 + b2
s


2
Z
+ VSAE  [x − r1 sin θ1 ]2 ++ z − r1 cos θ1 + a2 
2

s
2

Z
− z − r1 cos θ1 + a2 
+ VSAE  [x + r1 sin θ1 ]2 +
2

s
2

Z
+ VSAE  [x − r2 sin θ2 ]2 ++ z − r2 cos θ2 + a2 
2
s


2
Z
+ VSAE  [x + r2 sin θ2 ]2 +− z − r2 cos θ2 + a2 
2

s

2
Z
+ VSAE  [x − r3 sin θ3 ]2 ++ z − r3 cos θ3 + a2 
2
s


2
Z
+ VSAE  [x + r3 sin θ3 ]2 +− z − r3 cos θ3 + a2 
2

s

2
Z
+ VSAE  [x − r4 sin θ4 ]2 ++ z − r4 cos θ4 + a2 
2

s

2
Z
+ VSAE  [x + r4 sin θ4 ]2 +− z − r4 cos θ4 + a2  ,2H(x, 0, z, Z; 0) =p2x + p2z +(4.19)

107
which, after we change Z to −Z, becomesp2X + p2Z1+p24(x1 − x2 ) + (z1 − z2 )2 + b2
s

2

Z
+ VSAE  [x − r1 sin θ1 ]2 +− z + r1 cos θ1 + a2 
2

s

2
Z
+ VSAE  [x + r1 sin θ1 ]2 ++ z + r1 cos θ1 + a2 
2
s


2
Z
+ VSAE  [x − r2 sin θ2 ]2 +− z + r2 cos θ2 + a2 
2
s


2
Z
+ z + r2 cos θ2 + a2 
+ VSAE  [x + r2 sin θ2 ]2 +
2

s

2
Z
− z + r3 cos θ3 + a2 
+ VSAE  [x − r3 sin θ3 ]2 +
2

s

2
Z
+ VSAE  [x + r3 sin θ3 ]2 ++ z + r3 cos θ3 + a2 
2
s


2
Z
+ VSAE  [x − r4 sin θ4 ]2 +− z + r4 cos θ4 + a2 
2
s


2
Z
+ VSAE  [x + r4 sin θ4 ]2 ++ z + r4 cos θ4 + a2  .
2H 0 (x, 0, z, Z; 0) =p2x + p2z +(4.20)The Hamiltonian remains unchanged if any of the following conditions is satisfied:
r1 = −r2 ,r3 = −r4 ;(4.21)

r1 = −r3 ,r2 = −r4 ;(4.22)

r1 = −r4 ,r2 = −r3 .

(4.23)

That is to say, the system has to have central symmetry relative to the origin (center of mass).
For systems consisting of an odd number of centers, such pairs as shown above do not exist.
In this case, one atom has to be located at the origin. In a three-center system, for example, any

108of the following conditions has to be satisfied:
r1 = 0,r2 = −r3 ;(4.24)

r2 = 0,r1 = −r3 ;(4.25)

r3 = 0,r1 = −r2 .

(4.26)

Then, the central symmetry is still necessary.
Indeed, the system does not have to be homogeneous as discussed above. For example, inthe four-center system above, when the first condition, i.e.,r1 = −r2 ,r3 = −r4(4.27)

satisfies, we only need to require nucleus 1 and 2 to be the same and nucleus 3 and 4 to be thesame separately. Here, again, the central symmetry is still valid.
Furthermore, our analysis above relies on a two-dimensional model for each of the electrons.
However, it can be easily shown that our conclusions do not depend on this 2D analysis and are
valid for real 3D cases as well. Molecules, such as C60 , have central symmetry. Therefore, thefew-photon double ionization of C60 should also follow the same selection rules as in the case of
simple atoms.
Note that the central symmetry is needed when the two active electrons are included as well
as excluded, since both definite symmetries in the initial-state and the final-state Hamiltonians arenecessary so that both initial-state and final-state wave functions possess definite symmetries. For
molecules like C6 H6 or C60 , the p electrons are delocalized, thus the selection rules apply when
two p electrons are ionized if the remaining p electrons redistribute faster than the duration of
the ionizing laser pulse. If the p electrons redistribute slower than the duration of the ionizing
laser pulse, then the selection rules do not apply. In this sense, we may be able to map out theredistribution time of the p electrons by controlling the duration of the ionizing laser pulse. In
contrast, if two s electrons are ionized from C6 H6 or C60 , then the selection rules do not apply
since the s electrons are localized.

109
Other interesting applications of the selection rules include determining the ionizing orbitalsof the electrons in molecules like CO2 . If the two electrons are ionized from the different π bonds inthe molecule, then the selection rules apply; while the selection rules do not apply in other cases,
e.g., when they are ionized from the same π bond or when any of them is ionized from the σ bonds
or any inner s orbital of the individual atoms in the molecule.

4.7

Summary
In summary, we showed that the well-known selection rules for single-photon double ionization

can be generalized to the emission of two electrons following the simultaneous absorption of a few
photons from an intense laser pulse. These selection rules, which arise due to the symmetries of theaccessible final states in the two-electron continuum, correspond to nodes in the configuration space
wave function and the correlated momentum distributions. In particular, the back-to-back emissionof the two electrons with equal energy sharing is alternating either suppressed (forbidden) or nonsuppressed as the number of absorbed photons increases. The theoretical predictions are tested
well in numerical simulations. The numerical results for one- to four-photon double ionization of
the helium atom clearly show the presence of the node in the back-to-back emission in the case of
odd-number photon processes. Moreover, we showed that the presence of the selection rules does
not depend on the numerical model. Furthermore, we have also shown that the same selection rules
exist in the helium dimer, which is a diatomic molecule. Based on the studies of the homogeneous
two-center molecules like the helium dimer, we further studied multicenter molecules in general,and found that the selection rules apply whenever the molecule has the central symmetry.

Chapter 5

Photoabsorption and correlation time delay

In this chapter, we propose a self-consistent-time method to account for the Coulomb-laser
coupling effect and obtain the intrinsic photoabsorption time delay measured by the attosecond
streak camera. We illustrate this method with a one dimensional numerical model of the hydrogenatom. The concept of the method is based on an iteration to obtain the photoabsorption timedelay. Starting from an initial guess for this time delay, we iteratively obtain an improved timedelay by using the classical trajectory calculations and comparison with the vector potential of
the streaking pulse. We demonstrate the convergence and robustness of the self-consistent-time
method, investigate the photoemission at different positions, and study the resonant two-photon
ionization time delay.
Furthermore, we resolve the correlation time delay, defined as the delay between the absorption of the photon by the first electron (single ionization) and the transfer of the second electron
into the continuum (double ionization) in a single-photon double-ionization process of the helium
dimer via the variable space-partitioning method as well as the attosecond streaking technique,which is extended to many-electron systems. The streaked variable for the single ionization is
the momentum of the first electron, while the streaked variable for the double ionization is thecenter-of-mass momentum of both electrons. The correlation time delays obtained for different
laser parameters are compared with classical estimates.

111

5.1

Introduction
A recent report [120] on the photoemission time delay of 21 as from the 2s and 2p shells of

the neon atom using the attosecond streak camera induced a vigorous debate among theoreticians
about the origin of the time delay [121–140]. A large part of them was focused on the discussion
concerning the role of the Wigner-Smith (WS) time delay [141, 142] in the photoionization, whichitself diverges for long-range potentials like the Coulomb potential. Among others, our previous
works [138–140] addressed the concern regarding the divergence of the observed time delays byshowing that it can be written as a piecewise field-free time delays weighted by the ratio of theinstantaneous streaking field strength relative to its strength at the instant of transition of theelectron into continuum. The theoretical analysis was based on the agreement between resultsof ab-initio quantum-mechanical and classical-mechanical simulations of the propagation of theelectron in the combined field of the Coulomb potential of the residual ion and the streaking laser
field for the streaking time delay. This led to the interpretation of the measured time delay as dueto the Coulomb-laser coupling effect, i.e., the coupled dynamics of the electron in the fields of the
Coulomb potential and the streaking laser pulse.
In this chapter, we propose an alternative self-consistent-time (SCT) method to account for
the Coulomb-laser coupling effect and obtain the intrinsic photoabsorption time delay measured
by the attosecond streak camera. We illustrate this method using a one dimensional numericalmodel of the hydrogen atom. The concept of the method is based on an iteration to obtain thephotoabsorption time delay. The convergence and robustness of the method is also demonstrated.
Further, the photoemission at different positions and the resonant two-photon ionization time delay
is studied as well.
Beyond the temporal resolution of single-electron processes, one next challenge in attosecond
science is the real-time observation of energy exchange between electrons via their mutual correlation. The basic process to study energy exchange between electrons is the single-photon double
ionization. At low photon energy, one electron (primary electron) absorbs the photon before its

112
energy is shared with the other electron via the Coulomb interaction VC (refer to Sec. 3.4), andfinally both electrons emit from the target together [see Fig. 5.6(a)]. The correlation time, overwhich the energy is exchanged, can be defined as the delay between the absorption of photon energy by the primary electron and the emission of both electrons, or equivalently, the time delay
between single and double ionization. In this chapter, we further study the correlation time delayin the helium dimer via the variable space-partitioning method as well as the attosecond streaking
technique. The correlation time delays obtained for different laser parameters are compared withclassical estimates.

5.2

Photoabsorption time delay
In this section, we describe the one dimensional numerical model of the hydrogen atom usedin the simulations and develop the SCT method.

5.2.1

Numerical model

To make the SCT scheme transparent, we use a one dimensional numerical model of thehydrogen atom with HamiltonianH=

p2+ [EXUV (t) + Es (t)]x + V (x),
2

(5.1)

where
V (x) = − √Zx2 + a2(5.2)

is the Coulomb potential for the hydrogen atom with effective charge Z = 1 and soft-core parameter
a2 = 2, and EXUV and Es are the ionizing XUV laser field and streaking IR laser field, respectively,
both of which are given in the form ofE(t) = cos2



ωt
2N


cos(ωt + φ),(5.3)

where ω is the central frequency, N is the number of cycles, and φ is the carrier-envelope phase of
the laser pulse, which is set to −π/2 in the present simulations.

113
The TDSE is numerically solved using the Crank–Nicolson method on a grid with N = 200000
points, a grid step of ∆x = 0.05, and a time step of ∆t = 0.01. The ground state energy of themodel hydrogen atom is found to be ε0 = −0.500009 under current conditions. In the simulations,the interaction with a 16-cycle XUV laser pulse at a central wavelength of 30 nm and a peakintensity of 5 × 1014 W/cm2 and a 3-cycle IR laser pulse at a central wavelength of 2400 nm and apeak intensity of 1 × 1011 W/cm2 is used. A streaking pulse of 2400 nm instead of 800 nm is usedbecause we want to keep the XUV pulse shorter than one half period of the streaking pulse. We0|employ absorbing boundaries of the form cos1/6 π2 |x−xwith |x| > |x0 |, where x0 denotes theLborder of the boundary region and L is its width. The boundary was chosen to span 10% of thegrid size in each direction.
To obtain the streaking trace, we propagate the ionizing wave packet over long times andlarge distances until both laser fields cease and before the respective ionizing wave packet reaches
the boundary. The momentum spectrum is obtained by applying the Fourier transform to thecorresponding ionizing wave packet, which is well separated from the remaining bound state. We
verified that the error due to the finite grid step and projection onto the plane wave is negligible.

5.2.2

The self-consistent-time method1

In this section, we present the framework of the SCT method. First, we detail the SCT scheme
and illustrate it with a flowchart. Second, we show the convergence properties and robustness of
the SCT method. Third, we study the cases where photoemission takes places at different positions
other than the center of the hydrogen atom. Last, we study the resonant photoionization process,where an intrinsic photoabsorption time delay exists.

5.2.2.1

Scheme of the self-consistent-time method

To extract the final intrinsic photoabsorption time delay, we first obtain the streaking tracefrom the quantum-mechanical simulations. To this end, we scan the delay τ between the ionizing
1

The results presented in this section are currently prepared for publication (H. Ni, J. Su, A. Becker, and A.
Jaroń-Becker).

114
XUV laser pulse and the streaking IR laser pulse, and for each τ we obtain the momentum distribution of the ionizing wave packet after it has been propagated far away from the remaining boundstate. Afterwards, the expectation value of the momentum distribution is obtained for each τ , whichconstitutes the quantum streaking trace p0 (τ ). To minimize the error caused by finite propagation
distances, we correct the quantum streaking trace p0 (τ ) by calculating the asymptotic momentum
via energy conservation based on the expectation value of the position of the corresponding ionizing
wave packet.
Next we obtain the quantum time delay ∆t0 based on the quantum streaking trace p0 (τ )by fitting it to the vector potential As of the streaking laser field with fitting parameter α andtemporal offset ∆ts ,piter (τ ) = p0 − αAs (τ + ∆titer )/c,

(5.4)

where p0 is the asymptotic momentum of the electron in the absence of the streaking laser fieldand c is the vacuum speed of light. This apparent measured quantum time delay ∆t0 can wellbe used as the initial guess for the intrinsic photoabsorption time delay, although it may include
other contributions. At this moment, the Coulomb-laser coupling effect has not been taken into
account from the apparent measured photoabsorption time delay and no SCT iteration has been
implemented yet, and thus the superscript iter = 0, which tracks the number of iterations done, isused in notations p0 (τ ) and ∆t0 .
To obtain the momentum shift in the quantum streaking trace induced by the Coulomb-laser
coupling effect, we classically propagate the electron in the combined Coulomb-laser fields using
the Runge–Kutta method. To this end, we assume the photoemission takes place at the center ofthe hydrogen atom (xi = 0) with initial momentum
pi = sgn(p0 )qq
p20 − 2V (xi ) = sgn(p0 ) p20 − 2V (0),(5.5)

where p0 is the asymptotic momentum of the electron in the absence of the streaking laser field, as
mentioned before.
We propagate the electron classically in the combined field of the Coulomb potential of the

115
residual ion and the laser pulses, starting at instant
iter
titer
i (τ ) = τ + ∆t

(5.6)

and end the calculation at the same time instant tf as the quantum-mechanical calculation. Afteriterthis we note the final position xiterf (τ ) and final momentum pf (τ ) of the electron. The parametricdependence of ti on the XUV-IR delay τ implies that we need to repeat this calculation for each τ .
Now we repeat the classical propagation of the electron first in the case where there is only thestreaking laser field during the same time span from titeri (τ ) to tf , and then in the case where there
is only the Coulomb potential within the same spatial extension from xi to xiterf (τ ), separately, after
which the momentum shift in the quantum streaking trace induced by the Coulomb-laser coupling
effect can be obtained by iter  iter  iter∆pitercoupling (τ ) = pf,IR+Coulomb (τ ) − pi − pf,IR (τ ) − pi − pf,Coulomb (τ ) − pi ,(5.7)

or equivalently,
q
iter
iter
(τ
)
+
A[t
(τ
)]
−
sgn(p
)
p20 − 2V [xiter(τ)=
p∆piter0if,IR+Coulomb
couplingf (τ )],(5.8)

by taking advantage of the conservation of energy and canonical momentum and the fact thatA(tf ) = 0 since we end our propagations well after the end of the streaking pulse.
The improved streaking trace, where the Coulomb-laser coupling effect has been accounted
for, can then be obtained as
piter+1 (τ ) = p0 (τ ) − ∆piter
coupling (τ ).

(5.9)

Fitting this improved streaking trace to the vector potential of the streaking field again, an improved
photoabsorption time delay ∆titer+1 results, which can be used as the next guessed photoabsorption
time delay.
We repeat the iteration until ∆titer+1 and ∆titer are approximately equal, the condition ofwhich in the present simulations is∆titer+1 − ∆titer 6 1 × 10−6 .

(5.10)

116

Start

Quantum streaking trace p0 (τ ) (iter=0)

Quantum apparent time delay ∆t0 (iter=0)

Begin SCT calculation from titeri (τ ) [Eqn. (5.6)]
Coupling-induced shift ∆piter
coupling (τ ) [Eqn. (5.8)]
Improved streaking trace piter+1 (τ ) [Eqn. (5.9)]

Improved photoabsorption time delay ∆titer+1∆titer+1 ≈ ∆titer ?

Noiter + +Yes
Stop
Figure 5.1. Flowchart of the self-consistent-time (SCT) method.

117
As soon as this condition is fulfilled, convergence within the chosen error bound is reached, and thefinal intrinsic photoabsorption time delay is achieved. We summarize the iteration procedure as a
flowchart in Fig. 5.1.

5.2.2.2

Properties of the self-consistent-time method

We illustrate the iteration process by showing the streaking traces for the case of the one
dimensional model hydrogen atom during the iteration. Shown in Fig. 5.2 are the streaking tracesof the model hydrogen atom when exposed to a 16-cycle 30 nm XUV pulse with a peak intensity of
5 × 1014 W/cm2 and a 3-cycle 2400 nm IR streaking pulse with a peak intensity of 1 × 1011 W/cm2 .
The (blue) open circles are the expectation values of the momentum at different XUV-IR delay τand the solid (red) lines are their fits to the vector potential of the streaking field. Small insets
in the figures show the enlargements of the central part of the corresponding streaking trace fromτ = −20 to 20. Shown in panel (a) is the streaking trace before the SCT process, shown in panel(b) is the streaking trace after the first SCT iteration, and shown in panel (c) is the streaking trace
after the last SCT iteration.
As can be seen in the insets of the figures, there is an apparent measured photoabsorption
time delay for iter = 0, when no Coulomb-laser coupling effect has been taken into account, since
the central part of the trace tilts towards one side. After one self-consistent iteration, this apparent
time delay is already much reduced and hard to see, since the central trace is almost balanced on
both sides. For the next three iterations, the apparent time delay reduces slightly to finally reachconvergence within the present error bound for iter = 4.
Now, we illustrate the convergence process of the self-consistent iterations by extracting thetime delay through fitting the streaking traces in Fig. 5.2 to the vector potential of the streaking
laser field according to Eqn. (5.4). Shown in Fig. 5.3 are the fitted photoabsorption time delays foreach iteration. The solid (red) line with open circles stands for the fitted time delay for emission
to the positive direction, and the dashed (blue) line with open squares denotes that for emission
to the negative direction. As it is obvious from the figure, the convergence is rapidly reached after

118
iter = 0
(a)

2
1.5

Momentum (a.u.)

1
0.5
0
−0.5
−1
−1.5
−2

(b)

−400

−200

0
Delay (a.u.)iter = 1

200

400

−400

−200

0
Delay (a.u.)iter = 4

200

400

−400

−200

0
Delay (a.u.)

200

400

2
1.5

Momentum (a.u.)

1
0.5
0
−0.5
−1
−1.5
−2

(c)

2
1.5

Momentum (a.u.)

1
0.5
0
−0.5
−1
−1.5
−2

Figure 5.2. Iterations of the streaking trace of the one dimensional model hydrogen atom when
exposed to a 16-cycle 30 nm XUV pulse with a peak intensity of 5 × 1014 W/cm2 and a 3-cycle 2400
nm IR streaking pulse with a peak intensity of 1 × 1011 W/cm2 . The (blue) open circles are theexpectation values of the momentum at different XUV-IR delay τ and the solid (red) lines are theirfits to the vector potential of the streaking field. Small boxes in the figures are the enlargementsof the central part of the corresponding streaking trace from τ = −20 to 20. Shown in panel (a)is the streaking trace before any SCT iteration, shown in panel (b) is the streaking trace after the
first SCT iteration, and shown in panel (c) is the streaking trace after the last SCT iteration.

119
0.5

∆ t (a.u.)

0

0

−0.5

−1
Positive direction
Negative direction
−1.5
0

1

2
Iteration

3

4

Figure 5.3. Convergence of the self-consistent iterations of the photoabsorption time delay of the
one dimensional model hydrogen atom when exposed to a 16-cycle 30 nm XUV pulse with peakintensity of 5 × 1014 W/cm2 and a 3-cycle 2400 nm IR streaking pulse with a peak intensity of
1 × 1011 W/cm2 . The solid (red) line with open circles stands for the fitted time delay for emission
to the positive direction, and the dashed (blue) line with open squares denotes that for emission
to the negative direction. The inset shows the extracted time delays for iterations 1 through 4 ona different scale.

100
Positive direction
Negative direction80∆ t (a.u.)

60
40
20
0
−20
0

1

2

3

4

5

Iteration

Figure 5.4. Robustness of the self-consistent-time (SCT) method of the photoabsorption timedelay of the one dimensional model hydrogen atom when exposed to a 16-cycle 30 nm XUV pulse
with 5 × 1014 W/cm2 peak intensity and a 3-cycle 2400 nm IR streaking pulse with 1 × 1011 W/cm2
peak intensity. The solid (red) line with open circles stands for the fitted time delay for emission
to the positive direction, and the dashed (blue) line with open squares denotes that for emission to
the negative direction.

120
only 4 iterations.
Next, we illustrate the robustness of the SCT method by deliberately using an obviously
wrong initial time delay. Thus, instead of using the fitted photoabsorption time delay from thequantum streaking trace p0 (τ ) as the initial guess for the time delay, we use a largely exaggerated
wrong time delay as initial guess, and observe whether and how the convergence of the SCT method
takes place. Shown in Fig. 5.4 are the fitted photoabsorption time delay for each iteration when an
initial guess of 100 a.u. is taken. The solid (red) line with open circles stands for the fitted timedelay for emission to the positive direction, and the dashed (blue) line with open squares denotesthat for emission to the negative direction. Obviously, the iteration still rapidly converges to thecorrect intrinsic photoabsorption time delay despite of the large deviation of the initial guess fromthe correct value.

5.2.2.3

Different photoemission positions

Up to now, we have only considered the case when the photoemission takes place at thecenter of the atom, i.e., xi = 0. In this section, we remove this restriction by assuming different
photoemission positions of the electron in the atom and repeat the SCT iterations for these different
initial positions.
Shown in Fig. 5.5 are the converged photoabsorption time delays of the one dimensionalmodel hydrogen atom when exposed to a 16-cycle 30 nm XUV pulse with 5 × 1014 W/cm2 peakintensity and a 3-cycle 2400 nm IR streaking pulse with 1 × 1011 W/cm2 peak intensity when thephotoemission takes places at different positions xi . The solid (red) line stands for the fitted timedelay for emission to the positive direction, and the dashed (blue) line denotes that for emission
to the negative direction. As obvious from the figure, there is no unique solution for the photoemission position and photoabsorption time. For each assumed photoemission position xi , there isa corresponding photoabsorption time ti = τ + ∆t that fulfills all the requirements.
However, if we demand that photoionization for emission to positive and negative directions
has to occur at the same position due to the symmetry of the initial state, i.e., xpos
= xnegii , we find

121

1
Positive direction
Negative direction∆ t (a.u.)

0.5

0

−0.5

−1
−10

−5

0
xi (a.u.)

5

10

Figure 5.5. Converged photoabsorption time delays of the one dimensional model hydrogen atomwhen exposed to a 16-cycle 30 nm XUV pulse with 5 × 1014 W/cm2 peak intensity and a 3-cycle
2400 nm IR streaking pulse with 1 × 1011 W/cm2 peak intensity when the photoemission takes
places at different positions xi . The solid (red) line stands for the fitted time delay for emission to
the positive direction, and the dashed (blue) line denotes that for emission to the negative direction.

122
that there is an unique solution, namely ∆t = −0.00291755 ≈ 0. On the other hand, if we demandthat the photoionization for emission to positive and negative directions has the same time delay,i.e., ∆tpos = ∆tneg , we find that there is also an unique solution, namely xi = 0.000706299 ≈ 0.
While it is impossible to distinguish between all these non-unique solutions, it is straightforwardfrom the physical sense that, xi ≈ 0 and ∆t ≈ 0 is the most probable solution.
That said, we can still seek to do an ensemble average over the initial spatial distributionof the ground state of the hydrogen atom to take into account the possibility of photoemission at
different positions. To this end, we associate a weight for the converged photoabsorption time delay
at each xi according to the probability of the initial state at that point. For the ground state of our
1D model, the averaged photoabsorption time delay can be shown to be ∆tpos
ave = −0.0029164 foremission to the positive direction and ∆tneg
ave = −0.0031470 for emission to the negative direction.
If we further average over the emission to both directions, we have ∆tave = −0.0030317.

5.2.2.4

Resonant two-photon ionization time delay

In the previous sections, we have considered direct photoionization of the hydrogen atom
where no intermediate state plays a role. Here, photoabsorption occurs instantaneously, and no
intrinsic photoabsorption time delay is present. However, one may consider the case when anintermediate state plays a role, e.g., the resonant two-photon ionization, in which the absorption of
the first photon promotes the electron to an excited state, and the absorption of the other photon
ionizes the electron. If the electron absorbs the two photons in sequence, there might be a timedelay between these steps [143].
In this section, we consider the resonant two-photon ionization of the Li2+ ion [143]. The
reason to choose this ion is that it has larger energy gaps than the hydrogen atom and we canensure that only one excited state participates in this process. For comparison, we also consider
the case of nonresonant two-photon ionization of the same ion.
To this end, we use the same Coulomb potential as in Eqn. (5.2), but with an effective charge
Z = 3 and soft-core parameter a2 = 0.15. The ground-state energy, the first-excited-state energy,

123
Table 5.1. Two-photon ionization time delay for Li2+
Case
Nonresonant
Resonant∆tnon-SCT
0.00509.203688∆tSCT−0.000059
9.135658

and the second-excited-state energy of the model ion are found to be ε0 = −5.31689, ε1 = −2.31053,
and ε2 = −1.29547, respectively, using the same grid parameters. In the simulations, the interaction
with a 40-cycle XUV laser pulse at a central frequency of 3.51403 for the nonresonant case anda central frequency of 3.00636 for the resonant case (resonant with the first excited state startingfrom the ground state) and a peak intensity of 1 × 1013 W/cm2 and a 3-cycle IR laser pulse at a
central wavelength of 4800 nm and a peak intensity of 1 × 1011 W/cm2 is used.
Listed in Tab. 5.1 are the two-photon ionization time delays in both cases of nonresonant andresonant two-photon ionization. The values are averaged over emissions to two directions. Shownare the time delays obtained from both the previous non-SCT method [138–140] and the SCT
method. As can be seen from the table and it has been concluded by us recently [138–140], there
is no intrinsic photoabsorption time delay if it is a nonresonant two-photon ionization process, inwhich no intermediate states are involved. On the other hand, if the XUV pulse is tuned to be inresonance with the first excited state of the model ion, there is an intrinsic time delay of about 9
a.u. This suggests that there might be a time delay between the absorption of the first and thesecond photon, i.e., the electron may stay at the first excited state for a while before absorbing
another photon.
As for the accuracy of these two methods, the SCT method can be deemed as more accurate
as it can be seen from the obtained time delay in nonresonant two-photon ionization. If we assume
the accurate value of this time delay is 0, then the SCT method has an error of almost two orders
lower than the previous non-SCT method. In the case of resonant two-photon ionization, the nonSCT method then has a relative error of 0.74%. Although it is not as accurate as the SCT method,
the non-SCT method is still quite accurate considering this small error. For the detailed studies of

124
the resonant two-photon ionization time delay, interested readers are referred to Ref. [143].

5.3

Correlation time delay
In this section, we study the correlation time delay, which is the time delay for the twoelectrons to interact and exchange energy and be transferred to the continuum relative to thesingle ionization process, which is the first step of the knockout double photoionization process.
To this end, we first discuss the method that employs the attosecond streak camera to streak the
momentum of the primary electron in the first step, and streak the center-of-mass momentum of
the two electrons in the second step, then we show the results of the correlation time delay.

5.3.1

Correlation time delay in the helium dimer

To time-resolve correlated electronic dynamics, it has been recently proposed [144] to extend
the attosecond streak camera technique [42]. In the case of single ionization, the asymptotic
momentum p0 of the electron freed at time instance ts is shifted by a value proportional to thevector potential A of the infrared pulse at the release time ts
pf = p0 − A(ts )/c.

(5.11)

By varying the time delay between the triggering attosecond XUV pulse and the streaking infrared
pulse, the time of ionization of the electron is mapped onto its final momentum distribution. We
apply the attosecond streak camera technique to many-electron dynamics by noting that it is thecenter-of-mass electron momentum that couples to the external field [63, 145]. Thus, Eqn. (5.11)
can be readily extended to many-electron dynamics by replacing the single-electron momentum by
the center-of-mass momentum of the many-electron system:
Pf = P0 − ZA(tc )/c,

(5.12)

where tc denotes the time of correlated emission of two or more electrons from the atom or molecule,and Z stands for the number of electrons in the many-electron ionization process. We may further

125
note that the center-of-mass electron momentum is often observed in experiments using the so-called
COLTRIMS technique via the momentum of the recoil ion [55, 56].

Figure 5.6. Scheme of the correlation time delay τ . Vc is the Coulomb correlation between thetwo electrons.

Theoretical time-delay analysis for correlated electron emission using the streaking techniquehave been previously done for the helium atom. Schemes based on measurements of the angle
between the two escaping electrons [144] or the absolute quantum phase [146] have been analyzed.
However, according to our analysis [147] the time delay between single and double ionization for
two electrons in the ground state of helium atom is very small and hence difficult to resolve by
the attosecond streak camera technique. Alternative scenarios, e.g., the correlated dynamics in the
excited helium atom, which has been considered as well [144], are experimentally more difficult
to realize since those are easily being ionized by the streaking pulse itself. Here, we propose analternative to study the time delay of correlated electron emission in the helium dimer instead.
Due to the large separation between the atoms, dimers are ideal systems to study electron
correlation over long distances [86, 89]. At low photon energies, double photoionization of the dimer
proceeds by energy sharing over distances of approximately the internuclear distance of the dimer
via the knockout mechanism [86, 89] [cf. Fig. 5.6(b)]. We therefore expect that in the dimer thetime delay between photon absorption by the primary electron and transfer of part of the energy
to the other electron is much larger than in the helium atom.

126
5.3.2

Correlation time delay in the planar 4D helium dimer

In order to study this problem, we employ the 4D planar model of the helium dimer with
two active electrons located at separate atoms in the dimer (see Sec. 3.2).

5.3.2.1

Correlation time delay using attosecond streak camera

We employ the attosecond streaking camera technique to resolve the time delay within the
double photoionization process. As mentioned above, we identify this delay as the temporal difference between the time instants of single and double ionization of the dimer. The time instants
can be identified by streaking the single electron and the center-of-mass momenta of the electrons
along the polarization direction of the streaking pulse, respectively [see Eqns. (5.11) and (5.12)].
The corresponding streaking signals are obtained by varying the time delay between the ionizing
XUV pulse with respect to a 800 nm infrared streaking pulse (peak intensity of 5 × 1012 W/cm2 , asin2 -shaped envelope, and a total of 3 cycles) and are shown in Fig. 5.7. To save computation time,
only a scan of the time delay between the two pulses at the central streaking cycle is done. The
white solid lines passing through the middle of the streaking signals correspond to the expectation
values of the respective streaked momentum, which are fitted to the shape of the vector potentialof the infrared streaking pulse as p(t) = a + bAIR (t + τstreak ), where a, b, and τstreak are fitting
parameters. One clearly sees a small, but observable shift between these two traces, which is theapparent measured time delayDISIτstreak = τstreak− τstreak.

(5.13)

Note that the multiple peaks in the pZ momentum distributions is due to the two-center interferences [64, 148–150] in the helium dimer, which may result in an inaccuracy in the obtained
expectation value of the momentum.
As discussed above, this apparent time delay consists of two parts. The first part is theinherent (or intrinsic) time delay to the single and double ionization process, and the other part isthe propagation time delay due to the Coulomb-laser-coupling effect because of the presence of the127

Figure 5.7. Central-cycle streaking signals for single and double ionization contributions of theparallel oriented helium dimer starting from the singlet ground state at R = 5.6 with a XUV laser
pulse at a central wavelength of λXUV = 20 nm and a peak intensity of 1 × 1014 W/cm2 .

128
Table 5.2. Correlation time delay for different R and λXUV
R (a.u.)
5.6
5.6
10.0

λXUV (nm)20
1820

τclassical (as)
81.5
75.0
146

τevolution (as)
96.8
90.3
172

τstreak (as)
78.9
60.0
82.2

τextrapolate (as)
99.9
86.5
147

streaking laser field itself. This Coulomb-laser coupling effect is difficult to quantify in a two-centersystem and we may therefore, for now, assume that the coupling leads to similar propagation timedelays for single and double ionization. Under this assumption, the effect would play a small role
only in the difference
SI
SI
DI
DI
SI
DI
)
− τcoupling) − (τstreak− τcoupling= (τstreak− τcorrelation
τcorrelation = τcorrelation(5.14)SIDI),− τcoupling= τstreak − (τcoupling

and we getτcorrelation ≈ τstreak .

(5.15)

Shown in Tab. 5.2 is the correlation time delay as a function of the wavelength of the XUVlaser pulse and the internuclear distance of the dimer as compared to the classical time delay
estimated byRτclassical ≈ q,
2(ωXUV − Ip+ )(5.16)

where the denominator is just the classical velocity of the primary electron after its ionization, andanother estimation τevolution of the correlation time delay by observing the temporal evolution of thetwo-electron wave function as shown in Fig. 3.7. In this table, the streaking results are obtainedby collecting the electrons emitted in every direction in space since the primary electron may
get scattered to different directions from the neighboring atom and the center-of-mass momentaof the two electrons in the double ionization contribution can also be in different directions. The
streaking results are averaged over the emission in +z and −z directions and we have also taken into
account contributions from the singlet and triplet states by averaging according to their respective
multiplicities.

129
As it is obvious from the table, the result for the configuration with R = 10 has a quite large
discrepancy with the estimates, in part due to the difficulty in obtaining an accurate expectation
values of pZ and in part due to the difficulty in accounting for the Coulomb-laser-coupling induced
propagation time delay in the continuum for this two-center system.

5.3.2.2

Correlation time delay using the variable space-partitioning method

Since it has been difficult to obtain the intrinsic correlation time delay from the attosecond
streaking technique, here we instead use a more direct method to obtain the correlation time delaythat we denote as the variable space-partitioning method, which is similar to the space-partitioning
method that is used to obtain the single and double ionization contributions to the full wave
function.
In this method the configuration space is partitioned into different regions like for the spacepartitioning method in Sec. 2.7, and the probability entering the single ionization region is immediately regarded as single ionization, and the probability entering the double ionization region isimmediately regarded as double ionization, as discussed in the second method in Sec. 2.8. Witha pair of bounds dsmall and dbig to partition the different regions, the single and double ionization
yields are obtained as a function of propagation time according to the second method in Sec. 2.8.
The ionization yields are differentiated with respect to the propagation time which yields ionization rates. The single and double ionization times are then defined as the maxima of the respective
ionization rates.
In order to account for the variation due to the choice of the boundaries, we vary the pair ofbounds dsmall and dbig and extrapolate the respective ionization times to dsmall = 0 and dbig = 0.
To this end, we start from the base pair of bounds dsmall = 4 and dbig = 7, and either increasedbig when dsmall is fixed or increase dsmall when dbig is fixed. In the first case, we increase dbig
from 7 to 12 with step 1 when dsmall is fixed at 4. Then we plot the single ionization yield andthe corresponding ionization rate for each pair of bounds. Shown in Fig. 5.8 are the (a) single
ionization yield and (b) single ionization rate, as a function of propagation time for different dbig ,

130

−4

Ionization yield

(a) 4

x 10

dbig=7
dbig=83

dbig=9dbig=102

dbig=11
dbig=12

1

0
0

10

20

30

Time (a.u.)−5

Ionization rate

(b) 8

x 10

dbig=7
dbig=86

dbig=9dbig=104

dbig=11
dbig=12

2

0
0

10

20

30

Time (a.u.)
Figure 5.8. The (a) single ionization yield and (b) single ionization rate, as a function of propagation time for different dbig , of the parallel oriented singlet helium dimer at an internuclear distanceof R = 5.6 when it is exposed to a linearly-polarized (in the direction of the dimer axis) 4-cycle
XUV laser pulse with a central wavelength of 20 nm and a peak intensity of 1 × 1014 W/cm2 . Notethat time 0 is defined as the beginning of the laser pulse here.

131of the parallel oriented singlet helium dimer at an internuclear distance of R = 5.6 when it is
exposed to a linearly-polarized (in the direction of the dimer axis) 4-cycle XUV laser pulse with acentral wavelength of 20 nm and a peak intensity of 1×1014 W/cm2 . It can be clearly seen from thefigure that, as dbig gets bigger, the ionization yield curve shifts to later times, and the corresponding
peak of the ionization rate also shifts to later times. We can define the peak of the ionization rate
as the respective ionization time and extrapolate to dbig = 0. Similarly, in the second case, we
increase dsmall from 4.5 to 7 with step 0.5 when dbig is fixed at 7, and the corresponding double
ionization yield and rate can be obtained. Again, we identify the peak of the double ionization rate
as the double ionization time and extrapolate to dsmall = 0.
Shown in Fig. 5.9 is the single and double ionization time as a function of the partitioning
bounds dbig and dsmall . The values of tSI and tDI are obtained by using a linear fit to the data points
and extrapolation to dbig = 0 and dsmall = 0 respectively. In each panel, both singlet and triplet
states are considered by obtaining the ionization time via averaging the ionization yield according
to the multiplicities of the singlet and triplet states. The single ionization time as a function of
dbig is shown in the left column and the double ionization time as a function of dsmall is shown in
the right column. Shown in the first row are the results for application of a 20 nm XUV laser pulse
to a helium dimer at an internuclear distance of R = 5.6, shown in the second row are those for a
18 nm XUV laser pulse and a helium dimer at R = 5.6, and shown in the third row are those for a
20 nm XUV laser pulse and a helium dimer at R = 10. Other laser parameters are the same as inFig. 5.8.
Note that no streaking laser field is used, and thus there is no Coulomb-laser-couplinginduced time delay here. Therefore, the values of tSI and tDI shown in the figure correspond to theintrinsic photoabsorption time. In Tab. 5.2 we list the correlation time delay, defined as the timedelay between single and double ionization, obtained from extrapolation using the variable spacepartitioning method τextrapolate and the time delay obtained from other methods. The results agree
quite well, except those from the streaking method, in strong support of the knockout mechanism,for different internuclear distances R and XUV wavelengths λXUV . In particular, the time delay

132

12
11
10

Single ionization time(c) 13
12

8

9

dbig

10

11

11
10
9
8
7

(e) 14
13

8

9

dbig

10

11

11
10

9

dbig

10

11

12

5

5.5
6
dsmall

6.5

7

6.5

7

6.5

7

R=5.6, λ=18nm, DI
tDI=7.8262

16

15

14

(f) 21

R=10, λ=20nm, SI
tSI=4.9233

8

15

13
4.5

12

12

9
7

16

(d) 17R=5.6, λ=18nm, SI
tSI=4.2505

R=5.6, λ=20nm, DI
tDI=8.981

17

14
4.5

12

Dingle ionization time

9
7Single ionization time

Dingle ionization time13(b) 18R=5.6, λ=20nm, SI
tSI=4.8505

Dingle ionization timeSingle ionization time(a) 145

5.5
6
dsmall

R=10, λ=20nm, DI
tDI=11.0129

20

19

18

17
4.5

5

5.5
6
dsmallFigure 5.9. The single and double ionization time as a function of the partitioning bounds dbig
and dsmall . The values of tSI and tDI are obtained by using a linear fit to the data points andextrapolate to dbig = 0 and dsmall = 0 respectively. In each panel, both singlet and triplet states
are considered by obtaining the ionization time via averaging the ionization yield according to themultiplicities of the singlet and triplet states. The single ionization time as a function of dbig isshown in the left column and the double ionization time as a function of dsmall is shown in the rightcolumn. Shown in the first row are the results for application of a 20 nm XUV laser pulse to ahelium dimer at an internuclear distance of R = 5.6, shown in the second row are those for a 18
nm XUV laser pulse and a helium dimer at R = 5.6, and shown in the third row are those for a
20 nm XUV laser pulse and a helium dimer at R = 10. Other laser parameters are the same as inFig. 5.8. Note that time 0 is defined as the beginning of the laser pulse here.

133
increases with either increase of the XUV wavelength or the internuclear distance R, consistent
with the classical knockout mechanism.

5.4

Summary
As a summary, we propose a SCT method to account for the Coulomb-laser coupling effectand obtain the intrinsic photoionization time delay measured by the attosecond streak camera.
This method is based on an iteration method. Convergence is reached when the time delay remains
unchanged (within a given error bound) from one iteration to the next. We illustrated this method
using a one dimensional numerical model for the hydrogen atom and the Li2+ ion. This SCT
method features a rapid convergence, a robustness regarding initial guesses, and a better accuracy
than the previous non-SCT method.
Furthermore, we studied the correlation time delay, which is the time delay for the twoelectrons to interact and exchange energy and is defined as the delay between the absorption of
the photon (single ionization) and the emission of both electrons (double ionization), during singlephoton double-ionization of the helium dimer using the attosecond streaking technique, which isextended to many-electron systems. The streaked variable for the single ionization is the momentumof the first electron, while the streaked variable for the double ionization is the center-of-mass
momentum of both electrons. The correlation time delays obtained for different laser parameters
are compared with classical estimates and the estimates from the temporal evolution of the twoelectrons. However, rather poor agreement between the results using the attosecond streak camera
and the estimates was reached, due to the difficulty in accounting for the propagation time delay
induced by the Coulomb-laser coupling effect in the two-center system. To circumvent this problem,
we then used a more direct variable space-partitioning method to obtain the correlation time delay,
from which a good agreement with the estimates was reached.

Chapter 6

Summary and perspective

This thesis is devoted to the theoretical studies and numerical simulations of ultrafast correlated electron dynamics in single and double photoionization of atoms and molecules. To thisend, in Chapter 1, we gave a brief introduction to the background and some basic concepts ofattosecond and strong-field physics, including the strong-field single and double ionization, HHG,
and the attosecond streak camera technique. Next, in Chapter 2, we introduced the numericalmethods used in this thesis to solve the TDSE in the grid representation. The numerical methods
to obtain the ground state of a system and to propagate the initial state in a time-dependent potential, the absorbing boundaries, and the space-partitioning method that can be used to obtainsingle and double ionization contributions were presented. Strategies for parallel computation werealso briefly discussed.
In Chapter 3, we studied the double photoionization of the helium dimer using the planar
4D model, in which the electron interacts and exchanges energy via their Coulomb correlation
upon photon absorption over large distances. We first studied the different methods to obtain
the photoelectron angular distribution and found that the theoretical results agree well with theexperimental data if the momentum distribution is obtained by projection of the respective spatial
wave function onto an approximate two-electron wave function that is a product of a Coulomb
wave in their relative coordinates and a plane wave in their center-of-mass coordinates. We thenstudied the orientation and internuclear-distance dependence of the double-ionization probability
and revealed that the mechanism of the low-photon-energy two-site double photoionization of the135helium dimer is the two-step knockout process. The knockout process was also evident from ananalysis of the temporal evolution of the two-electron probability distribution.
We then introduced the Hamiltonian-reduction method to further study the effects of different
interactions in the single and double photoionization of the helium dimer. First, we removed theinteraction of the laser field with one electron. In this way, we were able to distinguish the primaryand secondary electron in the knockout process. Next we studied the scattering of the primaryelectron by the neighboring atom by removing the interaction between the secondary electron anda nucleus and locating the two electrons at certain nuclei. We further found that the electronnucleus interaction tends to focus the angular distribution to the direction of the dimer axis while
the electron-electron interaction tends to diverge the distribution. Exchange ionization, in which
the primary electron knocks out the other electron but gets captured itself, was also observed byremoving the interaction of the laser field with one electron.
In Chapter 4, we analyzed the selection rules for the the emission of two electrons from thehelium atom, the helium dimer, and general molecules following the absorption of a few photons
in an intense laser field. We first examined the selection rules starting from the general studiesof the two-electron continuum, the symmetry of the transition matrix, and the correlated KFRtheory. Then, we introduced the restricted 3D model of the helium atom, and numerically studied
the apparent nodes in the configuration and momentum space probability distributions and foundthat, in particular, the back-to-back emission of the two electrons for equal energy sharing is either
suppressed or not depending on the number of photons absorbed from the laser field. In addition,
we switched to a different planar 4D model of the helium atom and numerically studied part of theselection rules again and concluded that the presence of the selection rules do not depend on thespecific numerical model used. Furthermore, we studied the selection rules in the helium dimer andrevealed that the same selection rules applied to the helium dimer as well. As a step forward, we
generalized the selection rules to general multicenter molecules and showed that the rules apply as
long as the molecule has an central symmetry with respect to its center of mass.
In Chapter 5, we studied the photoabsorption time delay in the photoionization within a

136model hydrogen atom and the correlation time delay in the double photoionization of the helium
dimer. A SCT method was proposed to account for the Coulomb-laser coupling effect and obtain
the intrinsic photoabsorption time delay measured by the attosecond streak camera. The SCT
method started from an initial guess of the time delay, from which a streaking trace resulted, anda new iterated time delay can in turn be obtained. This process was continued until convergence
was reached and the time delay intrinsic to the photoabsorption process was retrieved. The SCT
method was demonstrated to be fast and robust in convergence. This method was also shown to
work in the case of two-photon resonant photoionization process.
Finally, we attempted to resolve the correlation time delay in the double photoionization
process of the helium dimer using a planar 4D model. The attosecond streaking technique wasused for this purpose, which was extended to many-electron systems. However, the agreement of
the time delay to classical estimates and the estimates from the temporal evolution of the twoelectron wave function was insufficient due to the difficulty in accounting for the Coulomb-laser
coupling effect and obtaining the intrinsic correlation time delay in this two-center system. As
an alternative, we used a variable space-partitioning method, which directly monitors the ionizing
wave packet, and were able to obtain a correlation time delay that was close to the estimates.

Bibliography[1] P. F. Moulton, Journal of the Optical Society of America B 3, 125 (1986).
[2] A. H. Zewail, Angewandte Chemie International Edition 39, 2586 (2000).
[3] F. Krausz and M. Ivanov, Reviews of Modern Physics 81, 163 (2009).
[4] L. V. Keldysh, Soviet Physics JETP 20, 1307 (1965).
[5] R. Dörner, H. Schmidt-Böcking, T. Weber, T. Jahnke, M. Schöffler, A. Knapp, M. Hattass,
A. Czasch, L. Ph. H. Schmidt, and O. Jagutzki, Radiation Physics and Chemistry 70, 191(2004).
[6] W. Becker, X. Liu, P. J. Ho, and J. H. Eberly, Reviews of Modern Physics 84, 1011 (2012).
[7] J. A. R. Samson, Physical Review Letters 65, 2861 (1990).
[8] A. Kheifets, Journal of Physics B: Atomic, Molecular and Optical Physics 34, L247 (2001).
[9] D. N. Fittinghoff, P. R. Bolton, B. Chang, and K. C. Kulander, Physical Review Letters 69,
2642 (1992).
[10] T. Schneider, P. L. Chocian, and J.-M. Rost, Physical Review Letters 89, 073002 (2002).
[11] K. J. Schafer, B. Yang, L. F. DiMauro, and K. C. Kulander, Physical Review Letters 70,
1599 (1993).
[12] P. B. Corkum, Physical Review Letters 71, 1994 (1993).
[13] J. L. Krause, K. J. Schafer, and K. C. Kulander, Physical Review Letters 68, 3535 (1992).
[14] J. J. Macklin, J. D. Kmetec, and C. L. Gordon, Physical Review Letters 70, 766 (1993).
[15] T. Popmintchev, M.-C. Chen, P. Arpin, M. M. Murnane, and H. C. Kapteyn, Nature Photonics 4, 822 (2010).
[16] M. F. Kling and M. J. Vrakking, Annual Review of Physical Chemistry 59, 463 (2008).
[17] S. T. Cundiff, Nature Physics 3, 16 (2007).
[18] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Physical Review
A 49, 2117 (1994).

138[19] I. P. Christov, M. M. Murnane, and H. C. Kapteyn, Physical Review Letters 78, 1251 (1997).
[20] K. J. Schafer and K. C. Kulander, Physical Review Letters 78, 638 (1997).
[21] M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec,
P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz, Nature 414, 509 (2001).
[22] M. Drescher, M. Hentschel, R. Kienberger, G. Tempea, C. Spielmann, G. A. Reider, P. B.
Corkum, and F. Krausz, Science 291, 1923 (2001).
[23] A. Baltuska, T. Udem, M. Uiberacker, M. Hentschel, E. Goulielmakis, C. Gohle,
R. Holzwarth, V. S. Yakovlev, A. Scrinzi, T. W. Hansch, and F. Krausz, Nature 421,
611 (2003).
[24] R. Kienberger, E. Goulielmakis, M. Uiberacker, A. Baltuska, V. Yakovlev, F. Bammer,
A. Scrinzi, T. Westerwalbesloh, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz,
Nature 427, 817 (2004).
[25] H. C. Kapteyn, M. M. Murnane, and I. P. Christov, Physics Today 58, 39 (2005).
[26] J.-F. Hergott, M. Kovacev, H. Merdji, C. Hubert, Y. Mairesse, E. Jean, P. Breger, P. Agostini,
B. Carré, and P. Salières, Physical Review A 66, 021801 (2002).
[27] P. B. Corkum, N. H. Burnett, and M. Y. Ivanov, Optics Letters 19, 1870 (1994).
[28] M. Kovacev, Y. Mairesse, E. Priori, H. Merdji, O. Tcherbakoff, P. Monchicourt, P. Breger,
E. Mével, E. Constant, P. Saliéres, B. Carré, and P. Agostini, The European Physical Journal
D - Atomic, Molecular, Optical and Plasma Physics 26, 79 (2003).
[29] I. J. Sola, E. Mevel, L. Elouga, E. Constant, V. Strelkov, L. Poletto, P. Villoresi, E. Benedetti,
J.-P. Caumes, S. Stagira, C. Vozzi, G. Sansone, and M. Nisoli, Nature Physics 2, 319 (2006).
[30] G. Sansone, E. Benedetti, F. Calegari, C. Vozzi, L. Avaldi, R. Flammini, L. Poletto, P. Villoresi, C. Altucci, R. Velotta, S. Stagira, S. De Silvestri, and M. Nisoli, Science 314, 443(2006).
[31] P. Tzallas, E. Skantzakis, C. Kalpouzos, E. P. Benis, G. D. Tsakiris, and D. Charalambidis,
Nature Physics 3, 846 (2007).
[32] C. Siedschlag, H. G. Muller, and M. Vrakking, Laser Physics 15, 916 (2005).
[33] G. D. Tsakiris, K. Eidmann, J. Meyer-ter Vehn, and F. Krausz, New Journal of Physics 8,
19 (2006).
[34] F. Quéré, C. Thaury, P. Monot, S. Dobosz, P. Martin, J.-P. Geindre, and P. Audebert,
Physical Review Letters 96, 125004 (2006).
[35] B. Dromey, M. Zepf, A. Gopal, K. Lancaster, M. S. Wei, K. Krushelnick, M. Tatarakis,
N. Vakakis, S. Moustaizis, R. Kodama, M. Tampo, C. Stoeckl, R. Clarke, H. Habara, D. Neely,
S. Karsch, and P. Norreys, Nature Physics 2, 456 (2006).
[36] C. A. Haworth, L. E. Chipperfield, J. S. Robinson, P. L. Knight, J. P. Marangos, and J. W. G.
Tisch, Nature Physics 3, 52 (2007).

139[37] V. T. Platonenko, V. V. Strelkov, and F. V. Ignatovich, Quantum Electronics 29, 601 (1999).
[38] I. Christov, H. Kapteyn, and M. Murnane, Optics Express 7, 362 (2000).
[39] I. Thomann, A. Bahabad, X. Liu, R. Trebino, M. M. Murnane, and H. C. Kapteyn, Optics
Express 17, 4611 (2009).
[40] M.-C. Chen, C. Hernández-Garcı́a, C. Mancuso, F. Dollar, B. Galloway, D. Popmintchev, P.C. Huang, B. Walker, L. Plaja, A. Jaroń-Becker, A. Becker, T. Popmintchev, M. M. Murnane,and H. C. Kapteyn, arXiv:1401.0240 [physics] (2013).
[41] L. Gallmann, C. Cirelli, and U. Keller, Annual Review of Physical Chemistry 63, 447 (2012).
[42] J. Itatani, F. Quéré, G. L. Yudin, M. Y. Ivanov, F. Krausz, and P. B. Corkum, PhysicalReview Letters 88, 173903 (2002).
[43] P. Eckle, A. N. Pfeiffer, C. Cirelli, A. Staudte, R. Dörner, H. G. Muller, M. Büttiker, and
U. Keller, Science 322, 1525 (2008).
[44] M. Uiberacker, T. Uphues, M. Schultze, A. J. Verhoef, V. Yakovlev, M. F. Kling, J. Rauschenberger, N. M. Kabachnik, H. Schröder, M. Lezius, K. L. Kompa, H.-G. Muller, M. J. J.
Vrakking, S. Hendel, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz, Nature
446, 627 (2007).
[45] A. L. Cavalieri, N. Müller, T. Uphues, V. S. Yakovlev, A. Baltuška, B. Horvath, B. Schmidt,
L. Blümel, R. Holzwarth, S. Hendel, M. Drescher, U. Kleineberg, P. M. Echenique, R. Kienberger, F. Krausz, and U. Heinzmann, Nature 449, 1029 (2007).
[46] E. Goulielmakis, Z.-H. Loh, A. Wirth, R. Santra, N. Rohringer, V. S. Yakovlev, S. Zherebtsov,
T. Pfeifer, A. M. Azzeer, M. F. Kling, S. R. Leone, and F. Krausz, Nature 466, 739 (2010).
[47] M. Odenweller, N. Takemoto, A. Vredenborg, K. Cole, K. Pahl, J. Titze, L. P. H. Schmidt,
T. Jahnke, R. Dörner, and A. Becker, Physical Review Letters 107, 143004 (2011).
[48] E. Goulielmakis, M. Uiberacker, R. Kienberger, A. Baltuska, V. Yakovlev, A. Scrinzi, T. Westerwalbesloh, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz, Science 305, 1267
(2004).
[49] D. Bradley, B. Liddy, and W. Sleat, Optics Communications 2, 391 (1971).
[50] M. Y. Schelev, M. C. Richardson, and A. J. Alcock, Applied Physics Letters 18, 354 (1971).
[51] J. S. Briggs and V. Schmidt, Journal of Physics B: Atomic, Molecular and Optical Physics
33, R1 (2000).
[52] W. Ackermann et al., Nature Photonics 1, 336 (2007).
[53] T. Shintake, H. Tanaka, T. Hara, T. Tanaka, K. Togawa, M. Yabashi, Y. Otake, Y. Asano,
T. Bizen, T. Fukui, S. Goto, A. Higashiya, T. Hirono, N. Hosoda, T. Inagaki, S. Inoue,
M. Ishii, Y. Kim, H. Kimura, M. Kitamura, T. Kobayashi, H. Maesaka, T. Masuda, S. Matsui,T. Matsushita, X. Maréchal, M. Nagasono, H. Ohashi, T. Ohata, T. Ohshima, K. Onoe,
K. Shirasawa, T. Takagi, S. Takahashi, M. Takeuchi, K. Tamasaku, R. Tanaka, Y. Tanaka,
T. Tanikawa, T. Togashi, S. Wu, A. Yamashita, K. Yanagida, C. Zhang, H. Kitamura, and
T. Ishikawa, Nature Photonics 2, 555 (2008).

140[54] P. Emma, R. Akre, J. Arthur, R. Bionta, C. Bostedt, J. Bozek, A. Brachmann, P. Bucksbaum,
R. Coffee, F.-J. Decker, Y. Ding, D. Dowell, S. Edstrom, A. Fisher, J. Frisch, S. Gilevich,
J. Hastings, G. Hays, P. Hering, Z. Huang, R. Iverson, H. Loos, M. Messerschmidt, A. Miahnahri, S. Moeller, H.-D. Nuhn, G. Pile, D. Ratner, J. Rzepiela, D. Schultz, T. Smith,P. Stefan, H. Tompkins, J. Turner, J. Welch, W. White, J. Wu, G. Yocky, and J. Galayda,
Nature Photonics 4, 641 (2010).
[55] R. Dörner, V. Mergel, O. Jagutzki, L. Spielberger, J. Ullrich, R. Moshammer,
H. Schmidt-Böcking, Physics Reports 330, 95 (2000).

and[56] J. Ullrich, R. Moshammer, A. Dorn, R. Dörner, L. P. H. Schmidt, and H. Schmidt-Böcking,
Reports on Progress in Physics 66, 1463 (2003).
[57] T. N. Rescigno, M. Baertschy, W. A. Isaacs, and C. W. McCurdy, Science 286, 2474 (1999).
[58] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martı́n, Reports on Progress inPhysics 64, 1815 (2001).
[59] R. Wehlitz, F. Heiser, O. Hemmers, B. Langer, A. Menzel, and U. Becker, Physical Review
Letters 67, 3764 (1991).
[60] O. Schwarzkopf, B. Krässig, J. Elmiger, and V. Schmidt, Physical Review Letters 70, 3008
(1993).
[61] F. Maulbetsch and J. S. Briggs, Journal of Physics B: Atomic, Molecular and Optical Physics
26, 1679 (1993).
[62] A. S. Kheifets and I. Bray, Journal of Physics B: Atomic, Molecular and Optical Physics 31,
L447 (1998).
[63] C. Ruiz, L. Plaja, L. Roso, and A. Becker, Physical Review Letters 96, 053001 (2006).
[64] D. Akoury, K. Kreidi, T. Jahnke, T. Weber, A. Staudte, M. Schoffler, N. Neumann, J. Titze,
L. P. H. Schmidt, A. Czasch, O. Jagutzki, R. A. C. Fraga, R. E. Grisenti, R. D. Muino,
N. A. Cherepkov, S. K. Semenov, P. Ranitovic, C. L. Cocke, T. Osipov, H. Adaniya, J. C.
Thompson, M. H. Prior, A. Belkacem, A. L. Landers, H. Schmidt-Bocking, and R. Dorner,
Science 318, 949 (2007).
[65] R. Dörner, H. Bräuning, O. Jagutzki, V. Mergel, M. Achler, R. Moshammer, J. M. Feagin,
T. Osipov, A. Bräuning-Demian, L. Spielberger, J. H. McGuire, M. H. Prior, N. Berrah, J. D.
Bozek, C. L. Cocke, and H. Schmidt-Böcking, Physical Review Letters 81, 5776 (1998).
[66] C. Siedschlag and T. Pattard, Journal of Physics B: Atomic, Molecular and Optical Physics
38, 2297 (2005).
[67] W. Vanroose, D. A. Horner, F. Martı́n, T. N. Rescigno, and C. W. McCurdy, PhysicalReview A 74, 052702 (2006).
[68] S. Baier, C. Ruiz, L. Plaja, and A. Becker, Physical Review A 74, 033405 (2006).
[69] J. Colgan, M. S. Pindzola, and F. Robicheaux, Physical Review Letters 98, 153001 (2007).

141[70] Y. Nabekawa, H. Hasegawa, E. J. Takahashi, and K. Midorikawa, Physical Review Letters
94, 043001 (2005).
[71] M. Kurka, J. Feist, D. A. Horner, A. Rudenko, Y. H. Jiang, K. U. Kühnel, L. Foucar, T. N.
Rescigno, C. W. McCurdy, R. Pazourek, S. Nagele, M. Schulz, O. Herrwerth, M. Lezius,
M. F. Kling, M. Schöffler, A. Belkacem, S. Düsterer, R. Treusch, B. I. Schneider, L. A.
Collins, J. Burgdörfer, C. D. Schröter, R. Moshammer, and J. Ullrich, New Journal of
Physics 12, 073035 (2010).
[72] A. Rudenko, Y. H. Jiang, M. Kurka, K. U. Kühnel, L. Foucar, O. Herrwerth, M. Lezius,
M. F. Kling, C. D. Schröter, R. Moshammer, and J. Ullrich, Journal of Physics B: Atomic,
Molecular and Optical Physics 43, 194004 (2010).
[73] A. Becker, R. Dörner, and R. Moshammer, Journal of Physics B: Atomic, Molecular andOptical Physics 38, S753 (2005).
[74] T. Hartman, P. N. Juranić, K. Collins, B. Reilly, N. Appathurai, and R. Wehlitz, PhysicalReview Letters 108, 023001 (2012).
[75] J. Crank and P. Nicolson, Advances in Computational Mathematics 6, 207 (1996).
[76] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes, 3rded. (Cambridge University Press, 2007).
[77] R. Kosloff and H. Tal-Ezer, Chemical Physics Letters 127, 223 (1986).
[78] M. D. Feit, J. A. Fleck Jr., and A. Steiger, Journal of Computational Physics 47, 412 (1982).
[79] J. L. Krause, K. J. Schafer, and K. C. Kulander, Physical Review A 45, 4998 (1992).
[80] C. W. McCurdy, M. Baertschy, and T. N. Rescigno, Journal of Physics B: Atomic, Molecularand Optical Physics 37, R137 (2004).
[81] F. He, C. Ruiz, and A. Becker, Physical Review A 75, 053407 (2007).
[82] U. V. Riss and H.-D. Meyer, Journal of Physics B: Atomic, Molecular and Optical Physics
26, 4503 (1993).
[83] R. Santra and L. S. Cederbaum, Physics Reports 368, 1 (2002).
[84] J. G. Muga, J. P. Palao, B. Navarro, and I. L. Egusquiza, Physics Reports 395, 357 (2004).
[85] J. S. Parker, B. J. S. Doherty, K. J. Meharg, and K. T. Taylor, Journal of Physics B: Atomic,
Molecular and Optical Physics 36, L393 (2003).
[86] H. Ni, C. Ruiz, R. Dörner, and A. Becker, Physical Review A 88, 013407 (2013).
[87] H. Ni and A. Becker, Physical Review A 89, 033402 (2014).
[88] T. Havermeier, T. Jahnke, K. Kreidi, R. Wallauer, S. Voss, M. Schöffler, S. Schössler, L. Foucar, N. Neumann, J. Titze, H. Sann, M. Kühnel, J. Voigtsberger, J. H. Morilla, W. Schöllkopf,
H. Schmidt-Böcking, R. E. Grisenti, and R. Dörner, Physical Review Letters 104, 133401(2010).

142[89] T. Havermeier, T. Jahnke, K. Kreidi, R. Wallauer, S. Voss, M. Schöffler, S. Schössler, L. Foucar, N. Neumann, J. Titze, H. Sann, M. Kühnel, J. Voigtsberger, A. Malakzadeh, N. Sisourat,
W. Schöllkopf, H. Schmidt-Böcking, R. E. Grisenti, and R. Dörner, Physical Review Letters
104, 153401 (2010).
[90] X. M. Tong and C. D. Lin, Journal of Physics B: Atomic, Molecular and Optical Physics 38,
2593 (2005).
[91] R. A. Aziz and M. J. Slaman, The Journal of Chemical Physics 94, 8047 (1991).
[92] W. Cencek and K. Szalewicz, International Journal of Quantum Chemistry 108, 2191 (2008).
[93] L. S. Davtyan, G. S. Pogosyan, A. N. Sisakyan, and V. M. Ter-Antonyan, Theoretical andMathematical Physics 74, 157 (1988).
[94] X. L. Yang, S. H. Guo, F. T. Chan, K. W. Wong, and W. Y. Ching, Physical Review A 43,
1186 (1991).
[95] T. Jahnke, A. Czasch, M. Schöffler, S. Schössler, M. Käsz, J. Titze, K. Kreidi, R. E. Grisenti,
A. Staudte, O. Jagutzki, L. P. H. Schmidt, S. K. Semenov, N. A. Cherepkov, H. SchmidtBöcking, and R. Dörner, Journal of Physics B: Atomic, Molecular and Optical Physics 40,
2597 (2007).
[96] L. S. Cederbaum, F. Tarantelli, A. Sgamellotti, and J. Schirmer, The Journal of Chemical
Physics 85, 6513 (1986).
[97] P. Lablanquie, T. P. Grozdanov, M. Žitnik, S. Carniato, P. Selles, L. Andric, J. Palaudoux,
F. Penent, H. Iwayama, E. Shigemasa, Y. Hikosaka, K. Soejima, M. Nakano, I. H. Suzuki,
and K. Ito, Physical Review Letters 107, 193004 (2011).
[98] H. Ni, S. Chen, C. Ruiz, and A. Becker, Journal of Physics B: Atomic, Molecular and Optical
Physics 44, 175601 (2011).
[99] F. Maulbetsch and J. S. Briggs, Journal of Physics B: Atomic, Molecular and Optical Physics
28, 551 (1995).
[100] M. Walter and J. S. Briggs, Physical Review Letters 85, 1630 (2000).
[101] A. W. Malcherek and J. S. Briggs, Journal of Physics B: Atomic, Molecular and Optical
Physics 30, 4419 (1997).
[102] M. G. Makris, P. Lambropoulos, and A. Mihelič, Physical Review Letters 102, 033002 (2009).
[103] A. Becker and F. H. M. Faisal, Physical Review A 50, 3256 (1994).
[104] A. Becker and F. H. M. Faisal, Physical Review A 51, 3390 (1995).
[105] F. H. M. Faisal, Journal of Physics B: Atomic and Molecular Physics 6, L89 (1973).
[106] H. R. Reiss, Physical Review A 22, 1786 (1980).
[107] A. Becker and F. H. M. Faisal, Journal of Physics B: Atomic, Molecular and Optical Physics
38, R1 (2005).

143[108] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman, Physical Review Letters42, 1127 (1979).
[109] J. H. Eberly, J. Javanainen, and K. Rzażewski, Physics Reports 204, 331 (1991).
[110] D. B. Milošević, G. G. Paulus, D. Bauer, and W. Becker, Journal of Physics B: Atomic,
Molecular and Optical Physics 39, R203 (2006).
[111] J. S. Parker, B. J. S. Doherty, K. T. Taylor, K. D. Schultz, C. I. Blaga, and L. F. DiMauro,
Physical Review Letters 96, 133001 (2006).
[112] H. Bachau, Physical Review A 83, 033403 (2011).
[113] A. Palacios, D. A. Horner, T. N. Rescigno, and C. W. McCurdy, Journal of Physics B:
Atomic, Molecular and Optical Physics 43, 194003 (2010).
[114] L. Tao, C. W. McCurdy, and T. N. Rescigno, Physical Review A 82, 023423 (2010).
[115] R. Nepstad, T. Birkeland, and M. Førre, Physical Review A 81, 063402 (2010).
[116] D. A. Horner, T. N. Rescigno, and C. W. McCurdy, Physical Review A 81, 023410 (2010).
[117] A. Rudenko, L. Foucar, M. Kurka, T. Ergler, K. U. Kühnel, Y. H. Jiang, A. Voitkiv,
B. Najjari, A. Kheifets, S. Lüdemann, T. Havermeier, M. Smolarski, S. Schössler, K. Cole,M. Schöffler, R. Dörner, S. Düsterer, W. Li, B. Keitel, R. Treusch, M. Gensch, C. D. Schröter,
R. Moshammer, and J. Ullrich, Physical Review Letters 101, 073003 (2008).
[118] R. Moshammer, Y. H. Jiang, L. Foucar, A. Rudenko, T. Ergler, C. D. Schröter, S. Lüdemann,
K. Zrost, D. Fischer, J. Titze, T. Jahnke, M. Schöffler, T. Weber, R. Dörner, T. J. M. Zouros,
A. Dorn, T. Ferger, K. U. Kühnel, S. Düsterer, R. Treusch, P. Radcliffe, E. Plönjes, and
J. Ullrich, Physical Review Letters 98, 203001 (2007).
[119] S. X. Hu, J. Colgan, and L. A. Collins, Journal of Physics B: Atomic, Molecular and Optical
Physics 38, L35 (2005).
[120] M. Schultze, M. Fiess, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A. L.
Cavalieri, Y. Komninos, T. Mercouris, C. A. Nicolaides, R. Pazourek, S. Nagele, J. Feist,
J. Burgdorfer, A. M. Azzeer, R. Ernstorfer, R. Kienberger, U. Kleineberg, E. Goulielmakis,
F. Krausz, and V. S. Yakovlev, Science 328, 1658 (2010).
[121] A. S. Kheifets and I. A. Ivanov, Physical Review Letters 105, 233002 (2010).
[122] I. A. Ivanov, Physical Review A 83, 023421 (2011).
[123] I. A. Ivanov, Physical Review A 86, 023419 (2012).
[124] C.-H. Zhang and U. Thumm, Physical Review A 82, 043405 (2010).
[125] C.-H. Zhang and U. Thumm, Physical Review A 84, 033401 (2011).
[126] J. C. Baggesen and L. B. Madsen, Physical Review Letters 104, 043602 (2010).
[127] M. D. Śpiewanowski and L. B. Madsen, Physical Review A 86, 045401 (2012).

144[128] S. Nagele, R. Pazourek, J. Feist, K. Doblhoff-Dier, C. Lemell, K. Tőkési, and J. Burgdörfer,
Journal of Physics B: Atomic, Molecular and Optical Physics 44, 081001 (2011).
[129] S. Nagele, R. Pazourek, J. Feist, and J. Burgdörfer, Physical Review A 85, 033401 (2012).
[130] R. Pazourek, J. Feist, S. Nagele, and J. Burgdörfer, Physical Review Letters 108, 163001(2012).
[131] R. Pazourek, S. Nagele, and J. Burgdörfer, Faraday Discussions 163, 353 (2013).
[132] M. Ivanov and O. Smirnova, Physical Review Letters 107, 213605 (2011).
[133] L. R. Moore, M. A. Lysaght, J. S. Parker, H. W. van der Hart, and K. T. Taylor, PhysicalReview A 84, 061404 (2011).
[134] J. M. Dahlström, A. L’Huillier, and A. Maquet, Journal of Physics B: Atomic, Molecularand Optical Physics 45, 183001 (2012).
[135] J. M. Dahlström, D. Guénot, K. Klünder, M. Gisselbrecht, J. Mauritsson, A. L’Huillier,A. Maquet, and R. Taı̈eb, Chemical Physics 414, 53 (2013).
[136] J. Su, H. Ni, A. Becker, and A. Jaroń-Becker, Physical Review A 87, 033420 (2013).
[137] J. Su, H. Ni, A. Becker, and A. Jaroń-Becker, Journal of Modern Optics 60, 1484 (2013).
[138] J. Su, H. Ni, A. Becker, and A. Jaroń-Becker, Physical Review A 88, 023413 (2013).
[139] J. Su, H. Ni, A. Becker, and A. Jaroń-Becker, Chinese Journal of Physics 52, 404 (2014).
[140] J. Su, H. Ni, A. Becker, and A. Jaroń-Becker, Physical Review A 89, 013404 (2014).
[141] E. P. Wigner, Physical Review 98, 145 (1955).
[142] F. T. Smith, Physical Review 118, 349 (1960).
[143] J. Su, H. Ni, A. Jaroń-Becker, and A. Becker, Time delays in two-photon ionization, submitted(2014).
[144] A. Emmanouilidou, A. Staudte, and P. B. Corkum, New Journal of Physics 12, 103024(2010).
[145] C. Ruiz and A. Becker, New Journal of Physics 10, 025020 (2008).
[146] A. S. Kheifets, I. A. Ivanov, and I. Bray, Journal of Physics B: Atomic, Molecular and Optical
Physics 44, 101003 (2011).
[147] C. Ruiz and A. Becker, unpublished .
[148] K. Kreidi, D. Akoury, T. Jahnke, T. Weber, A. Staudte, M. Schöffler, N. Neumann, J. Titze,
L. P. H. Schmidt, A. Czasch, O. Jagutzki, R. A. Costa Fraga, R. E. Grisenti, M. Smolarski,
P. Ranitovic, C. L. Cocke, T. Osipov, H. Adaniya, J. C. Thompson, M. H. Prior, A. Belkacem,
A. L. Landers, H. Schmidt-Böcking, and R. Dörner, Physical Review Letters 100, 133005(2008).

145[149] M. S. Schöffler, K. Kreidi, D. Akoury, T. Jahnke, A. Staudte, N. Neumann, J. Titze, L. P. H.
Schmidt, A. Czasch, O. Jagutzki, R. A. Costa Fraga, R. E. Grisenti, M. Smolarski, P. Ranitovic, C. L. Cocke, T. Osipov, H. Adaniya, S. Lee, J. C. Thompson, M. H. Prior, A. Belkacem,
T. Weber, A. Landers, H. Schmidt-Böcking, and R. Dörner, Physical Review A 78, 013414(2008).
[150] K. Kreidi, D. Akoury, T. Jahnke, T. Weber, A. Staudte, M. Schöffler, N. Neumann, J. Titze,
L. P. H. Schmidt, A. Czasch, O. Jagutzki, R. A. C. Fraga, R. E. Grisenti, R. D. Muiño,
N. A. Cherepkov, S. K. Semenov, P. Ranitovic, C. L. Cocke, T. Osipov, H. Adaniya, J. C.
Thompson, M. H. Prior, A. Belkacem, A. Landers, H. Schmidt-Böcking, and R. Dörner, The
European Physical Journal Special Topics 169, 109 (2009)