Strong-field interaction with bicircular and ellipticallypolarized laser pulses
byY. Gebre
M.S., University of Colorado Boulder, 2022
B.S., Illinois Institute of Technology, 2018
B.S., Illinois Institute of Technology, 2018

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
2023
Committee Members:
Andreas Becker, ChairJohn Bohn
Jose P. D’Incao
Margaret Murnane
Ralph Jimenez

ii
Gebre, Y. (Ph.D., Physics)
Strong-field interaction with bicircular and elliptically polarized laser pulses
Thesis directed by Prof. Andreas Becker

Ultrafast physics, which encompasses phenomena occurring at the attosecond scale, provides
the foundation for understanding electron dynamics in atoms, molecules, and materials. To explore
this realm, ultrashort laser pulses are utilized as powerful tools for discerning dynamics and unraveling the underlying physics. By varying the laser polarization and employing pulse combinations,a wide range of intriguing ultrafast phenomena can be probed. In this thesis numerical simulations,specifically of the interaction of atoms with circularly and elliptically polarized pulses, are used toshed light on the fascinating dynamics exhibited by atoms when subjected to intense ultrafast laser
irradiation.
We first provide an overview of the dynamics that form the basis of the research conducted
and the numerical methods employed to model the atom-laser interaction through the solution of
the time-dependent Schrödinger equation. Building upon this foundation we then delve into theinvestigation of atoms interacting with bichromatic circularly polarized laser pulses. The distribution of population among various excited states is examined, and a new mechanism is proposed
to explain the results. Next, we shift our focus to the analysis of photoelectron spectra generated
by the ionization of atoms interacting with circularly and elliptically polarized pulses. The study
specifically considers the effects of different initial magnetic states on the resulting spectra, shedding
light on the underlying electron emission process. Lastly, we present two new numerical methods:
A Monte Carlo simulation technique is introduced as a means to study electron dynamics andionization in laser fields of long wavelengths. Second, an ab-initio solution for diatomic systems,
based on the single-active-electron and Born-Oppenheimer approximations, is developed, applied,
and tested.

Dedication

To my family, friends and Tesla.

iv

Acknowledgements

Throughout my journey as a graduate student, I have come to deeply appreciate the invaluable role played by the mentors and guides who have shaped my academic career. I am immensely
grateful to the supportive community at JILA, whose unwavering assistance has enabled me toundertake the research presented in this work.
First and foremost, I express my heartfelt gratitude to my advisor, Dr. Andreas Becker, whose
instrumental guidance and unwavering support have been pivotal to my success. His mentorshiphas shaped my research and made me a much better scientist. I would also like to acknowledge theinvaluable mentorship provided by Dr. Joel Venzke during the early years of my graduate studies.
Additionally, I extend my sincere thanks to Cory Keasling, Jim McKown, and Eric Alvarado for their consistent and reliable support in addressing my computational questions andneeds throughout this journey.
I am also indebted to my professors and friends from my undergraduate years, whose influence
has been instrumental in my achievements. I would like to give special thanks to my undergraduate
research advisor, Bryce Littlejohn, for fostering an exceptional research environment that nurtured
the growth of his students. I am also grateful to Pranava Teja Surukuchi, whose mentorship served
as a shining example of the kind of graduate student I aspired to become.
Lastly, my deepest gratitude goes to my family for their unwavering love and support. Their
presence has provided me with a sense of belonging and has made this journey possible.

v

Contents

Chapter
1 Introduction and background1

1.1

Attosecond time scale and dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Atomic units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.1

Above threshold ionization and AC stark shift . . . . . . . . . . . . . . . . .

5

1.3.2

Multi-photon and tunneling ionization regimes . . . . . . . . . . . . . . . . .

8

1.3.3

Strong field ionization rates . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.4

Rydberg state excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5

Raman transitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.6

Outline of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2 Time-dependent Schrödinger equation

17

2.1

Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.2

Ab-initio solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2.1

Basis and coordinate system. . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.2.2

Kinetic energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2.3

Single active electron potential . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2.4

Laser coupling in length gauge . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.2.5

Laser coupling in velocity gauge . . . . . . . . . . . . . . . . . . . . . . . . . 25

vi

2.3

2.4

2.2.6

Bound and continuum states . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.2.7

Time propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.2.8

Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Quantum trajectory Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . 32
2.3.1

Adiabatic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.3.2

Non-adiabatic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.3.3

Time integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Laser pulses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3 Pathways to excitation of atoms with bicircular laser pulses

46

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.2

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3

3.2.1

Excitation with corotating pulses . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.2.2

Excitation with counterrotating pulses . . . . . . . . . . . . . . . . . . . . . . 54

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4 Transitions between Rydberg states in two-color corotating circularly polarized laser pulses 59
4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.2

Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.3

4.2.1

Population distribution in time-delayed bicircular corotating pulses . . . . . . 63

4.2.2

Transitions between Rydberg states . . . . . . . . . . . . . . . . . . . . . . . 65

4.2.3

Higher-order transitions in (ω, pω) corotating circularly polarized pulses . . . 68

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5 Photoelectron emission spectra of atoms in intense short circularly and elliptically polarizedlaser pulses

72

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.2

Absorption pathways and photoelectron spectra . . . . . . . . . . . . . . . . . . . . . 73

vii
5.3

Emission angles of photoelectrons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6 Analysis of photoelectron energy spectra using Quantum Trajectory Monte Carlo simulations

85

6.1

QTMC with PPT ionization rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

6.2

Photoelectron spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6.3

Angular distribution of photoelectrons . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6.4

Advantages and limitations of QTMC simulations . . . . . . . . . . . . . . . . . . . . 93

7 Ab-initio solution for diatomic systems

95

7.1

Additional approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

7.2

Matrix elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.2.1

Spherical part of the potential . . . . . . . . . . . . . . . . . . . . . . . . . . 98

7.2.2

Cylindrical part of the potential . . . . . . . . . . . . . . . . . . . . . . . . . 101

7.3

Bound states for diatomic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

7.4Time propagation and observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7.5

Current performance and future work . . . . . . . . . . . . . . . . . . . . . . . . . . 111

8 Summary and outlook

113

Bibliography

116

Appendixviii

TablesTable
1.1

Hartree atomic units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7.1

Parameters for molecular single active electron potential obtained using analytical

2fits of potentials for hydrogen and oxygen molecules [200]. . . . . . . . . . . . . . . . 106

ix

FiguresFigure1.1Time and length scales of ultrafast processes studied in atomic, molecular and solidstate physics (taken from [1]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2(a) Single-photon ionization, where a single energetic photon leads to ionization.(b)Ionization through the absorption of multiple photons.(taken from [11]). . . . . . . .

1.3

2

3

Strong field ionization: (a) tunneling ionization. (b) above barrier ionization (takenfrom [12]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.4

Diagram of above threshold ionization (taken from [12]). . . . . . . . . . . . . . . . .

5

1.5

Photoelectron energy spectra illustrating above threshold ionisation (ATI) of hydrogen at a laser intensity I = 1014 W/cm2 and frequency ℏω = 2eV (taken from [14]). .

1.6

6

Illustration of energy level shift demonstrating a linear increase in ionization potential
with intensity. The diagram illustrates that at intensity I1 , ionization can occur
through the absorption of four photons, whereas at higher intensities I2 and I3 , theelectron does not have enough energy to ionize (taken from [16]). . . . . . . . . . . .

1.7

7

Various photoionization regimes are shown for a given laser pulse. The photon
energies corresponding to the laser wavelengths of 800 nm and 400 nm are marked
by vertical dotted red and blue lines, respectively. The dashed gray line represents
where γ = 1 applies to the ionization of an atom or molecule with an ionizationenergy of 7.54eV creating the boundary between the multiphoton regime (γ ≫ 1)
and the tunneling regime (γ ≪ 1)(taken from [12]). . . . . . . . . . . . . . . . . . . .

8

x
1.8

Diagrams for two-photon Raman transitions. Λ−type Raman transition on the left,
V −type on the right.

1.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Diagrams for three photon Raman transitions. Λ−type Raman transition on the left
and V −type in the middle and S-type on the right.

2.1. . . . . . . . . . . . . . . . . . 15

Regions of laser parameters where the dipole and non-relativistic approximations arevalid (taken from [52]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.2

Illustration of the polarization direction of a right-handed circularly polarized light(RCPL) pulse, as it changes as a function of time (taken from [66]). . . . . . . . . . 37

2.3

Comparison of the exit radius for the TIPIS model and the non-adiabatic model [63].
The Figure shows the exit radius as a function of the electric field strength, both
expressed in a.u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.4

Illustration of the reject-sampling method for sampling a distribution. The Figure
shows a scatter plot of events in the X-Y space, where red dots represent rejected
samples (Yi > P (Xi )) and green dots represent accepted samples (Yi ≤ P (Xi )).

2.5. . 42

Temporal profile of a laser pulse at 200 nm wavelength, 1×1014 W/cm² intensity, and15 cycles. The red line depicts the envelope function, while the blue line represents
the electric field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.1

Excited state distribution as function of n (vertical axis) and ℓ (horizontal axis) for(a) m = −1, (b) m = −2, (c) m = −3 and (d) m = −4 at the end of 20 (at 800 nm)
cycle pulses (40 cycle at 400 nm) with sin squared envelope and total peak intensityof 1 × 1014 W/cm2 for corotating laser pulses of equal intensity (taken from [174]).

3.2

49

Absorption pathways in corotating laser pulses at frequencies ω and 2ω starting froma m = 0-state. Without lack of generalization it is assumed that both pulses haveleft-handed helicity. Absorption of a photon at frequency ω and at frequency 2ω isrepresented by a red and blue arrow, respectively. The numbers in the boxes denotethe minimum number of photons to reach a certain level (taken from [174]).

. . . . 50

xi
3.3

Excited state distribution as function of orbital angular quantum number ℓ summed
over n ≥ 4 and m at (a) I400 = 5 × 1013 W/cm2 , I800 = 5 × 1012 W/cm2 , (b)I400 = 5 × 1013 W/cm2 , I800 = 1 × 1013 W/cm2 , (c) I400 = 5 × 1013 W/cm2 ,I800 = 5 × 1013 W/cm2 , and (d) I400 = 1 × 1013 W/cm2 , I800 = 5 × 1013 W/cm2 .
Pulse durations: 20 cycles at 400 nm, 10 cycles at 800 nm (taken from [174]).

3.4. . . 51

Excited state distribution as function of n (vertical axis) and m (horizontal axis)
summed over ℓ. Laser parameters: 20 (800 nm) cycle pulses with sin squared envelope
and total peak intensity of 1 × 1014 W/cm2 for corotating laser pulses of equal
intensity (taken from [174]).

3.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Distribution in magnetic quantum states for (a) ℓ = 2 and (b) ℓ = 3 and different
peak intensities of the 800 nm pulse. I400 = 5 × 1013 W/cm2 and other parameters
are as in Fig. 3.4 (taken from [174]).

. . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.6

Same as Fig.3.1 but for counterrotating laser pulses (taken from [174]).

. . . . . . . 54

3.7

Absorption pathways in counterrotating laser pulses at frequencies ω and 2ω startingfrom a m = 0-state. Without lack of generalization it is assumed that the pulse at
frequency ω has left-handed helicity, while the second harmonic pulse has righthanded velocity. Other symbols as in Fig. 3.2 (taken from [174]).

3.8. . . . . . . . . . 55

Excited state distributions as a function of m, summed over n ≥ 4 and ℓ (top),
and as a function of n and m, summed over ℓ (bottom), for the interaction with aleft-handed circularly polarized laser pulse at 800 nm (20 cycles) and a right-handed
circularly polarized laser pulse at 400 nm (40 cycles). Both pulses have the same
peak intensity of 5 × 1013 W/cm2 (taken from [174]).

3.9. . . . . . . . . . . . . . . . . 56

Orbital angular momentum distributions in excited states induced by counterrotating
laser pulses at 400 nm (20 cycles) and 800 nm (10 cycles) at peak intensities of (a)I400 = 5 × 1013 W/cm2 , I800 = 5 × 1012 W/cm2 , (b) I400 = 5 × 1012 W/cm2 ,I800 = 5 × 1013 W/cm2 , and (c) I400 = 5 × 1013 W/cm2 , I800 = 5 × 1013 W/cm2(taken from [174]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

xii
4.1

Schemes for Λ− (upper row), V- (middle row) and S-transitions (lower row) between
Rydberg atomic states at central frequencies ω (red arrows) and 2ω (blue arrows)changing the magnetic quantum number by ∆m = −1 (left) or ∆m = +1 (right).
The two applied corotating circularly polarized fields are assumed to be left-handed
polarized (taken from [179]).

4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Comparison of excited state distributions as a function of magnetic quantum numberm (summed over n and ℓ with n ≥ 4) for excitation with (a) circularly polarized laser
pulse at 267 nm (20 cycle, 5 × 1013 W/cm2 ), (b) circularly polarized laser pulse at
534 nm (10 cycle, 1 × 1014 W/cm2 ), and (c) corotating bicircular laser pulse (takenfrom [179]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.3

Excitation probability in states with magnetic quantum number m (summed over n
and ℓ with n ≥ 4) for two corotating circularly polarized pulses at 267 nm (20 cycle,
5 × 1013 W/cm2 ) and 534 nm (10 cycle, 1 × 1014 W/cm2 ) as a function of time delay
between the pulses. At negative (positive) time delays the 267 nm (534 nm) pulse
precedes (taken from [179]).

4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Excitation probability as function of magnetic quantum number (summed over n andℓ with n ≥ 4, upper row) and as function of principal and orbital angular momentum
quantum numbers (summed over m, lower row) for two corotating circularly polarized
pulses at 267 nm (20 cycle, 5 × 1013 W/cm2 ) and 534 nm (10 cycle, 1 × 1014 W/cm2 ).
Results are obtained for initial states prepared in n0 = 7 and ℓ0 = 3, m0 = −3 (a,e), ℓ0 = 4, m0 = −4 (b, f), ℓ0 = 5, m0 = −5 (c, g), and ℓ0 = 6, m0 = −6 (d, h)(taken from [179]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.5

Comparison of normalized excitation probabilities in states with magnetic quantum
number m (summed over n and ℓ with n ≥ 4) for two corotating circularly polarized
pulses at 267 nm (30 cycle, 5 × 1013 W/cm2 ) and 534 nm (15 cycle, 1 × 1014 W/cm2 ),800 nm (10 cycle, 1 × 1014 W/cm2 ), and 1068 nm (7.5 cycle, 1 × 1014 W/cm2 ) (takenfrom [179]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

xiii
4.6

Excitation probability as function of magnetic quantum number (summed over n andℓ with n ≥ 4) and as function of principal and orbital angular momentum quantum
numbers (summed over m, lower row) for two corotating circularly polarized pulses
at 267 nm (30 cycle, 5 × 1013 W/cm2 ) and 800 nm (10 cycle, 1 × 1014 W/cm2 ).
Results are obtained for initial states prepared in n0 = 10 and ℓ0 = 3, m0 = −3(a), ℓ0 = 5, m0 = −5 (b), ℓ0 = 7, m0 = −7 (c), and ℓ0 = 9, m0 = −9 (d) (takenfrom [179]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.1

Photon absorption pathways for a right-handed circularly polarized pulse with initialstates of m = −1 (red) and m = 1 (blue). Direct path for corotating electrons
indicated by the red line, direct path for counterrotating electrons indicated by theblue line, and off-set path for counterrotating electrons marked by the light-blue line. 73

5.2

Comparison of photoelectron energy spectra for neon atom, initially in l = 1, m = 1(corotating) state, indicated by the red curve and atoms, initially in l = 1, m = −1
(counterrotating) state, indicated by the blue curve. The top panel represents thetotal spectrum, while the bottom panel shows the spectrum originating from thedirect paths only. The results were obtained using a 10-cycle right-handed circularly
polarized laser pulse with an intensity of 3 × 1014 W/cm2 . . . . . . . . . . . . . . . . 75

5.3

Same as Fig. 5.2 but for an 800 nm pulse. . . . . . . . . . . . . . . . . . . . . . . . . 76

5.4

Distribution of population in the orbital angular and magnetic states of neon atom
after interaction with a 400 nm circularly polarized laser pulse. Laser parameters
are the same as in Fig. 5.2. Left for initial counterrotating state with l = 1, m = −1
and right for corotating state with l = 1, m = 1. . . . . . . . . . . . . . . . . . . . . . 77

xiv
5.5

Photoelectron spectrum for a neon atom interacting with an elliptically polarizedlaser pulse. The top row corresponds to results for the initial (counterrotating)state with l = 1, m = −1, and the bottom row corresponds to those for the initial(corotating) state with l = 1, m = 1. The left column displays the angle-resolved
emission, the middle column shows the emission as a function of angle for each ATI
peak, and the right column represents the spectrum generated with the three highest
contributing components. The pulse parameters are as follows: wavelength of 400
nm, intensity of 3 × 1014 W/cm², and ellipticity of 0.7. . . . . . . . . . . . . . . . . . 79

5.6

The angle corresponding to the maximum emission of electrons for neon atom interacting with an elliptically polarized pulse. The laser parameters are the same as inFig. 5.5 for ellipticities of 0.8 (blue) and 0.9 (red). . . . . . . . . . . . . . . . . . . . 80

5.7

Illustration of the dependence of angle of electron emission based on the model given
by Eq. (5.3). Top panel shows the phase of the three contributing components while
the bottom panel shows the resulting spectra. . . . . . . . . . . . . . . . . . . . . . . 82

5.8

Phase of the three biggest contributing components of the photoelectron spectrum ofneon interacting with an elliptically polarized pulse. The plots on the left correspondto the corotating case, while the plots on the right correspond to the counterrotatingcase. The top row represents an ellipticity of 0.7, and the bottom row represents anellipticity of 0.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.1

The photoelectron spectrum of a neon atom interacting with an elliptically polarized
pulse at 400 nm is shown in the Figure. The plot on the left corresponds to an initial(corotating) state with quantum numbers l = 1 and m = 1, while the plot on theright corresponds to an initial (counterrotating) state with quantum numbers l = 1
and m = −1. The laser intensity is 9 × 1014 W/cm2 , the pulse has an ellipticity of
0.7 and a duration of 10 cycles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6.2

Same as Fig. 6.1 but for a laser pulse with a wavelength of 600 nm. . . . . . . . . . 90

xv
6.3

Angle-resolved spectrum of a neon atom subject to interaction with an ellipticallypolarized pulse. The plots on the left-hand side display results from the QTMC
simulation, while the plots on the right-hand side depict the results derived fromthe ab-initio solution. The top plots correspond to the l = 1, m = −1 initial state,
whereas the bottom plots represent the l = 1, m = −1 initial state. The parameters
are: Laser intensity of 3 × 1014 W/cm2 , ellipticity of 0.7, and pulse duration of 10
cycles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6.4

Same as Fig. 6.3 but for 600 nm laser pulse. . . . . . . . . . . . . . . . . . . . . . . . 92

7.1

Illustration depicting the orientation of a diatomic molecule, where the axes o andρ represent the coordinate along the molecule axis and the radial distance perpendicular to the molecular axis, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 96

7.2

Lowest-energy orbitals for diatomic (hydrogen) molecule, where 1σg (lower part) and
1σu (upper part) orbitals are formed by combining two 1s hydrogen atom orbitals(taken from [202]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

7.3

Hydrogen molecule orbitals from TISE calculations, where left is the 1σg - and rightis the 1σu -orbital, respectively. The color corresponds to |Ψ|2 eiϕ where ϕ is the phase.107

7.4

Same as Fig. 7.2 but for 1πg (upper part) and 1πu (lower part) orbitals, which areformed by combining two 2p hydrogen atom orbitals (taken from [202]). . . . . . . . 108

7.5

Same as Fig. 7.3 but for 1πg and 1πu orbital. The color corresponds to |Ψ|2 eiϕ where
ϕ is the phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.6

photoelectron spectra for hydrogen molecule interacting with a 10 cycle laser pulseof 200 nm (left) and 400 nm (right) and a peak intensity of 1.5 × 1014 W/cm2 . . . . 110

7.7

Hydrogen molecule bound state calculation and time propagation scaling as a function of utilized cores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Chapter 1

Introduction and background

1.1

Attosecond time scale and dynamics
An often cited illustration highlighting the significance of ultrafast time scales is the use of ahigh-speed camera in 1878 to capture an image of a horse running on a track. The purpose of thephotography was to settle a debate regarding whether all four of a horse’s feet left the ground while
galloping. The outcome revealed that indeed all four legs are lifted. This story serves to emphasize
the need for tools that operate at the time scale of the phenomena we wish to investigate in order tocomprehend physical processes. Natural phenomena occur across a broad range of time scales. Fig.
1.1 presents a graph illustrating the time scales of various processes studied in atomic, molecular,
and solid-state physics [1]. To examine the behavior of electrons within atoms, molecules, and othersystems, it is imperative to possess temporal resolution that captures their dynamics. Electron
dynamics occur on the scale of attoseconds (1 as = 10−18 seconds) and ultrafast atomic physics
investigates dynamics at this temporal scale [1–4]. As the motions of nuclei, atoms, and molecules
are negligible at this time scale, our focus can solely be on the dynamics of electrons. Ultrafast laser
pulses are one of the tools employed to explore this time scale [5, 6]. Attosecond-scale laser pulses
can be employed to illuminate electron dynamics and provide an understanding of the underlying
physics governing processes such as tunneling, multiphoton ionization, high harmonic generation,and related phenomena.

2

Figure 1.1: Time and length scales of ultrafast processes studied in atomic, molecular and solidstate physics (taken from [1]).

1.2

Atomic units
The units used in atomic calculations are atomic units, abbreviated as a.u. We will be using

the Hartree atomic units for all the work presented in this thesis. In atomic units the reduced
Planks constant (ℏ), electron mass (m), and the charge of the electron (e) are set to 1. The othernecessary units are defined so that there is consistency for dimensionless physical constants toretain their values. The Hartree atomic units for relevant quantities are presented in Table 1.1. In
Quantitylengthmass
energyvelocitytime
angular momentum

Unit
Bohr radius, ao
electron mass, me
twice of ionization potential of H
electron velocity in 1st Bohr orbit
period of Bohr orbitalℏValue in SI
0.05292 × 10−9 m
9.1093837 × 10−31 kg
27.2110 eV
2.18 × 108 cm/s
24.18884 as
1.05457 × 10−34 J · s

Table 1.1: Hartree atomic units

atomic units, the speed of light is given by c = 1/α, where α is the fine structure constant. For thelaser intensity we express it in W/cm2 (watts per centimeter2 ), where 1 a.u = 3.51 × 1016 W/cm2 .
Atomic units enable us to represent quantities in their natural scales, which proves particularly
advantageous for the numerical methods discussed in subsequent Chapters. By utilizing atomic

3units, we avoid the need to work with extremely small or large numbers when expressing various
quantities of interest.

1.3

Ionization
The interaction between light and matter and the ionization process has been a topic of study

for over a century, with significant breakthroughs in the early 1900s. One of the most significant
contributions was made by Albert Einstein in 1905 when he proposed the photoelectric effect [7].
According to this effect, electrons can absorb photons of energy larger than a threshold frequency,
and this energy helps them overcome the Coulomb potential and escape the atoms, molecules, or
solids. This phenomenon is now known as single-photon ionization, and the energy of the ionizedelectron is given by the equation E = ℏω − Ip , where ℏ is the reduced Planck constant, ω isthe frequency of the absorbed photon, and Ip is the ionization potential. Later, in 1931, Maria
Goeppert Mayer predicted that ionization could also occur through the absorption of two or more
photons [8]. This phenomenon is known as multi-photon ionization, and it was confirmed with
the invention of the laser, which provided enough intensity to allow for the absorption of multiple
photons [9, 10] (see Fig. 1.2).

Figure 1.2: (a) Single-photon ionization, where a single energetic photon leads to ionization.(b)Ionization through the absorption of multiple photons.(taken from [11]).

4
In addition to multi-photon ionization, there is another process in strong-field physics bywhich electrons can escape their parent nucleus, known as tunneling. When subjected to a strong
field, the atomic potential becomes distorted, creating a barrier through which the electron cantunnel. This barrier is created by the superposition of the atomic potential and the electric field, as
demonstrated in Fig. 1.3(a). The height and width of the barrier are determined by the strength ofthe electric field applied. As the strength of the field increases, the potential barrier becomes smaller
and narrower, allowing for an increase in the rate of ionization. For extremely large intensities, thepotential barrier becomes lower than the energy level of the bound electron, as shown in Fig. 1.3(b).
This phenomenon is known as over-the-barrier ionization, in which the electron can escape without
having to tunnel. Multi-photon absorption and tunneling are not mutually exclusive explanationsof ionization. The particular process that takes place relies on the properties of the target atomand the field applied to it. Further details regarding this subject will be explored in section 1.3.2.

Figure 1.3: Strong field ionization: (a) tunneling ionization. (b) above barrier ionization (takenfrom [12]).

5
1.3.1

Above threshold ionization and AC stark shift

When an electron is exposed to an intense field, it can absorb more photons than necessary
to escape the binding potential of the parent nucleus, resulting in a phenomenon known as abovethreshold ionization (ATI) [13, 14]. The additional energy gained by the electron translates to its
kinetic energy. The phenomenon is illustrated in Fig. 1.4, where the electron absorbs two additional
photons beyond what is required for ionization. This effect can be observed in the photoelectron
spectrum of electrons ionized by intense laser fields. Since the energy of the photon is quantized,
each additional absorption process results in a gain of ℏω in the energy of the ionized electrons. As
a consequence of this process, peaks appear in the photoelectron spectra that are separated by ℏω,where each peak corresponds to an additional photon absorption during ionization.

Figure 1.4: Diagram of above threshold ionization (taken from [12]).

An example of a (computed) photoelectron spectrum is shown in Fig. 1.5. The Figure shows theenergy spectrum of electrons ionized from hydrogen atom. We see well defined ATI peaks separatedby ℏω for each additional photon process.

6

Figure 1.5: Photoelectron energy spectra illustrating above threshold ionisation (ATI) of hydrogenat a laser intensity I = 1014 W/cm2 and frequency ℏω = 2eV (taken from [14]).

As laser technology progressed, it was noticed that the lower peaks in the photoelectron
spectrum became less prominent and eventually vanished as the laser intensity and pulse duration
increased. This phenomenon is due to the Stark-shift and occurs when the intense electric fieldof the laser alters the energy levels of atoms [15–19]. Fig. 1.6 demonstrates the shift in energy
levels as a function of laser intensity. This dynamic shift in energy is more pronounced for energy
states near the continuum, while deeply bound states experience minimal changes in their energy.
The magnitude of the Stark-shift for the loosely bound states is approximately determined by the
ponderomotive energy, expressed as
Up =

Eo2
,
4ω 2

(1.1)

where Eo is the electric field amplitude. The ponderomotive energy is the cycle-averaged kinetic
energy of the quiver motion of the electron in the presence of an oscillating electric field. Whenconsidering both the Stark-shift and ATI effects, the energy of the emitted electrons is described

7by
En = nℏω − Ip − Up .

(1.2)

The first term in this equation represents the energy of the absorbed photons, while Ip is theionization potential of the atom. The last term represents the Stark-shift of the states near thecontinuum, which shifts the ionization energy to an effective ionization energy of Ief f = Ip + Up .
The effect of the Stark-shift on the spectrum was only recognized with the advent of short laser
pulses. Earlier experiments used longer pulses in which the electron had the opportunity to exit
the focal point before the intensity diminished considerably. Once outside the focus, the electron
experiences a gradient force that accelerates it and imparts energy equal to the ponderomotiveenergy. Thus, the energy loss due to the Stark-shift was gained by the quiver motion in the field
and the shift in the spectrum was not present.

Figure 1.6: Illustration of energy level shift demonstrating a linear increase in ionization potential
with intensity. The diagram illustrates that at intensity I1 , ionization can occur through theabsorption of four photons, whereas at higher intensities I2 and I3 , the electron does not haveenough energy to ionize (taken from [16]).

8
1.3.2

Multi-photon and tunneling ionization regimes

The phenomenon of ionization is a complex process influenced by various factors, such as theproperties of the atom being targeted and the strength of the electric field applied. Fig. 1.7 shows
the distinct ranges in which each ionization process occurs, emphasizing the significance of the
Keldysh parameter in determining whether multi-photon ionization or tunneling prevails. Keldysh
introduced the dimensionless Keldysh parameter in 1965 [20], which incorporates the ionization
potential of the target atom as well as the intensity and frequency of the electric field. The Keldysh
parameter is given by the equation,

γ=

ω

s
p
2Ip
Ip
=
.
Eo
2Up

(1.3)

Figure 1.7: Various photoionization regimes are shown for a given laser pulse. The photon energies
corresponding to the laser wavelengths of 800 nm and 400 nm are marked by vertical dotted red
and blue lines, respectively. The dashed gray line represents where γ = 1 applies to the ionizationof an atom or molecule with an ionization energy of 7.54eV creating the boundary between themultiphoton regime (γ ≫ 1) and the tunneling regime (γ ≪ 1)(taken from [12]).

For a given target atom, a field with a long wavelength and high intensity corresponds to a Keldysh
parameter of γ ≪ 1, which is referred to as the adiabatic limit. In this regime, the width of the9barrier is small, and the field is varying slowly, allowing the electron to tunnel through the barrier
before the field changes direction. In this regime, the field is considered to be static during thetunneling process, hence the name adiabatic. On the other hand, for a weaker field and short wavelength, a Keldysh parameter γ ≫ 1 is observed. In this regime, the electron needs only to absorb
a few photons to overcome the binding potential and ionize vertically. The fields will be too weak
to allow effective tunneling, and the vertical channel becomes the dominant process for ionization.
However, electrons can also be ionized with a mixture of the two processes in the parameter spacereferred to as the non-adiabatic tunneling regime. This occurs when γ ∼ 1, where the electron
tunnels through the barrier, but the field dynamics affect the ionization process. This regime is ofparticular interest as it presents an opportunity to study the interplay between multi-photon ionization and tunneling. The difference between adiabatic tunneling and non-adiabatic tunneling is
explained in more detail in the following sections. It is important to note that the ionization process
can also be influenced by other factors such as the details of the atomic potential, the polarizationof the laser field, and the temporal structure of the pulse. Therefore, a complete understanding ofionization requires the consideration of these additional factors as well.

1.3.3

Strong field ionization rates

The phenomenon of strong field ionization has been a topic of extensive research and development since Keldysh’s influential work in 1965. Over the years, many models have been proposed
to explain the underlying physics of this process, including the Ammosov-Delone-Krainov (ADK)model [21, 22], the Perelomov, Popov, and Terent’ev (PPT) model [23, 24], as well as more recent models like those proposed by Yudin and Ivanov [25], or the quantitative tunneling formula(QTF) [26]. In this discussion, we will focus on two models that have been particularly influential
in the field of strong field ionization. The first model we will discuss is the ADK formula, which
provides a means of predicting the ionization rate in the adiabatic limit discussed above. The
ADK model takes into account the properties of the target atom and the electric field to provide a

10formula for the ionization rate in this regime. The second model we will discuss is the PPT model,
which predicts strong field ionization in the non-adiabatic regime. In this regime, the electric field
has a weaker intensity and shorter wavelength, causing the electron to absorb a few photons during
the tunneling process. The PPT model takes into account the dynamic effects of the electric field
during the tunneling process, leading to a more accurate prediction of the ionization rate in thisregime. Both the ADK and PPT models assume the strong-field approximation (SFA) [20, 27–29],
which neglects the Coulomb potential of the nuclei once the electron has transitioned into a continuum state. The SFA has been shown to be a reasonable approximation for many experimental
conditions and has greatly simplified the theoretical description of strong field ionization.
The ADK formula for strong field ionization is given by
WADK (to , p⃗i ) ∝ exp −κ (p⊥ip − p⊥ )
2κ3
exp −
,
3|E(to )||E(to )|

(1.4)

where to is the time of ionization, p⃗i is the initial momentum vector, p⊥ip is the initial momentum
perpendicular to the field in the plane of polarization, p⊥ is the initial momentum perpendicularp
to the field in the plane of propagation and, κ = 2Ip . In accordance with the adiabatic limit,
the momentum parallel to the direction of the field is assumed to be zero. The ADK formula is alimiting case of the PPT formula which can be applied in the adiabatic and non-adiabatic regimes.
Here we give the ionization formula as stated in Ref. [30] that applies the PPT formula for a shortrange potential and different initial states. The ionization rate in the PPT model is given as
Z π
WP P T (to , p⃗i ) ∝

ω

π
−ω

2
iS(v(t))ϕlm (v(t)) e

dt ,

(1.5)

where v(t) = p + A(t), p is the conserved canonical momentum, A(t) is the vector potential and11
ϕlm (v(t)) is a function of the initial state (φ(r)) given by
1
v(t)2 + 2Ip φ̃ (v(t))
2Z
1
φ̃ (v(t)) =e−ip·r φ(r) dr.
(2π)3/2ϕlm (v(t)) =(1.6)

S is the action of the system and is given by
Z t
S=−

1
[p + A(t)]2 + Ip dt.
2
0

(1.7)

Using the saddle point approximation, we can simplify the expression for the ionization rate. In Eq.
(1.5), the integral contains an exponential term with a negative argument, causing it to rapidly
approach zero as the magnitude of the action integral increases. Hence, we anticipate that the
integrand for the ionization rate reaches its highest value at the minimum of the action. The saddle
point approximation assumes that the primary contribution to the integral in Eq. (1.5) arises fromthe critical points of the exponential’s argument. Consequently, the integral can be approximated bysumming the contributions from these critical points, which are commonly known as saddle points.
This approximation holds for laser parameters where ω ≪ Ip and ω ≪ Up [30]. Applying the saddle
point approximation to the ionization rate of the PPT model yields the following expression [30]:s
WP P T (to , p⃗i ) ∝X

ϕlm (v(ti )) exp [iS(v(ti ), ti )]

i

2

2π
S ′′ (v(ti ), ti )

,

(1.8)

where ti are the saddle points given by the saddle point equation1[p + A(ti )]2 + Ip = 0.
2

(1.9)

The set of critical points obtained from the saddle point equation is used to compute the ionization
rate. In section 2.3.2, we demonstrate the procedure for solving the saddle point equation for alaser pulse with arbitrary polarization. Additionally, we provide a physical interpretation of these12
saddle points and present an application of the PPT formula for investigating ionization dynamics.

1.4

Rydberg state excitation
Rydberg states are highly excited states of an atom, in which the outermost electron is on

average at a very large distance from the nucleus. Rydberg states have a relatively large size andexhibit unique physical properties, such as long lifetimes and strong dipole moments. Electrons
can be excited into Rydberg states through the absorption of one or more photons. The excitedstates play a role in many fundamental processes in strong field atomic physics like non-sequential
double ionization, above-threshold ionization (ATI), and high harmonic generation (HHG). There
has been a surge of interest in understanding the mechanism for exciting Rydberg states throughintense laser pulses, with researchers exploring the distribution of excited state populations acrossdifferent quantum numbers [31–43]. Studies, such as those conducted in [38, 39, 44], revealed acorrelation between the closing of a channel caused by the AC-Stark shift and the distribution of
excited state population over the principal quantum numbers. This effect suggests that as the excitation probability adjusts, the channel closing occurs, requiring an additional photon for ionization,
particularly at threshold intensities. Recent research has also focused on understanding the distribution of the Rydberg state population with respect to the angular momentum quantum number.
For linearly polarized pulses, numerical calculations align well with semiclassical estimations in
predicting the angular momentum state with the largest population [38, 39, 45]. Furthermore, Floquet theory has demonstrated that the parity of populated Rydberg states is dictated by Np − 1,
where Np is the number of photons necessary for ionization [40, 44]. Another work [46] on the
Rydberg state population has shown that ultrashort laser pulses can excite Rydberg states via theabsorption of different numbers of photons, leading to a mix of even and odd number of photon
processes. When the intensity of the laser pulse is low, the distribution of population across theangular momentum states can be attributed to the parity effect resulting from the selection rulesof multiphoton absorption. Additionally, the suppression of pathways due to population trappingin low excited states also contributes to this pattern, as long as the laser pulse is not too short.

13
At high intensity the separation in population of odd vs. even angular momentum quantum statescannot be achieved even for long pulses. The loss of the parity effect in the distribution at longpulses with high peak intensity is due to interference effects in excitation or frustrated ionizationin the tunneling regime.
In Chapter 3, we extend the investigation of the population distribution of Rydberg states bystudying the interaction between a hydrogen atom and bichromatic circularly polarized laser pulses.
It has been shown that the relative helicity of the two bicircular fields influences the level of ionization observed [47–49]. The suggested explanation for this phenomenon is that the counterrotating
pulses generate an ionization pathway that leads to resonance-enhanced ionization. However, theauthors of this study did not examine the distribution of excited states across quantum numbers.
By conducting such an analysis, we can gain an additional understanding of the importance of
excited states in the ionization pathways. This is because the excitation that occurs in a resonant
multiphoton process is expected to be influenced by the selection rules governing the interaction
between spin-angular momentum and a bicircular pulse. We investigate the distribution of excitedstates across quantum numbers and provide an explanation for the observed distribution using
selection rules and classical cutoff predictions.

1.5

Raman transitions
One of the ultrafast dynamics explored in this thesis pertains to the transfer of population

between different atomic states. When an atomic or molecular system interacts with a laser field,
electrons can undergo transitions between states by either absorbing or emitting photons. By
considering the polarization of the photon and the initial state of the electron, we can exert some
control over the states that become populated. A mechanism commonly employed for this purpose
is a Raman transition [50]. In a Raman transition, a pair of continuous wave lasers are adjusted tocouple two distinct states to an intermediate state. The first laser pulse facilitates the transfer of the
population from the initial state to the intermediate state, while the second laser pulse transfers the
population from the intermediate state to the desired final state. Notably, it has been demonstrated

14
that this process allows for efficient population control, ensuring that the intermediate state does
not sustain any steady-state population [50, 51].
Raman transitions can occur through V − or Λ− type configurations, depending on thechosen intermediate state. Fig. 1.8 provides a simplified illustration of both the V −transitionand Λ−type Raman transitions. In a Λ−type Raman transition, the first photon facilitates the
population transfer from the ground state to the excited intermediate state, and the second photon
subsequently transfers the excited population to a lower-lying final state. The V −transition follows
a similar process, except that the intermediate state lies below the initial state.

E2⟩

E3⟩

E1⟩
E1⟩

E3⟩

E2⟩

Figure 1.8: Diagrams for two-photon Raman transitions. Λ−type Raman transition on the left,
V −type on the right.

The transfer of population between intermediate and final states is subject to restrictions on
the quantum numbers, as explained in Chapter 2, based on selection rules. For circularly polarized
pulses, the change in the magnetic quantum number for each step is determined by the polarizationof the pulse, following the selection rules. This restriction on quantum numbers presents a means
to transfer populations among magnetic states. In Chapter 4, we propose that Raman transitions
involving three photons serve as the underlying process responsible for the observed distributionof quantum numbers in the Rydberg states of hydrogen when interacting with bicircular pulses.
Fig. 1.9 provides a simplified visual representation of these Raman transitions, specifically involving
three photons. An important aspect of the three-photon Raman transition is that, upon completing

15
a full cycle of the three-photon process, there is a net change in the magnetic number between
the initial and final states. This mechanism enables the transfer of population among Rydberg
states that cannot be accessed through absorption alone. Further elaboration on this process
and its implications for the distribution of populations among different quantum numbers will be
extensively discussed in Chapter 4.

E2⟩
Δm = − 1

E1⟩

E3⟩

E3⟩
E1⟩

Δm = + 2

E4⟩
E2⟩

E3⟩

Λ type Raman transition

ΛE1⟩

type Raman transition

E2⟩
S type Raman transition

Figure 1.9: Diagrams for three photon Raman transitions. Λ−type Raman transition on the left
and V −type in the middle and S-type on the right.

1.6

Outline of thesis
The thesis is organized as follows. In Chapter 2, we provide an overview of two numericalapproaches for solving the Time-Dependent Schrödinger Equation (TDSE) to investigate the interaction between atomic systems and laser pulses. The first approach is the ab-initio method, which
relies on first principles for solving the TDSE. We delve into the intricacies and implementation
details of this method. The second approach we introduce is the quantum trajectory Monte Carlo(QTMC) method, a Monte Carlo simulation technique that combines analytical ionization rates
with quantum trajectories. Chapter 3 presents the findings of a study involving bi-circular laser
pulses interacting with hydrogen atoms. We specifically examine the population of Rydberg states

16
and the distribution of the excited state population across different quantum numbers. Moving
on to Chapter 4, we demonstrate how Raman transitions, involving three-photon processes, canexplain the observed population distribution among various magnetic states for atoms interactingwith corotating bi-circular laser pulses. In Chapter 5, we investigate the photoelectron spectrum
resulting from the interaction of atoms with short circularly and elliptically polarized pulses. We
analyze the differences observed in the spectrum for atoms initially in different magnetic states.
In Chapter 6, we employ the QTMC method to study the photoelectron spectrum of Neon interacting with an elliptically polarized pulse. We compare these outcomes with the ab-initio solution
to validate the accuracy of the QTMC method. Finally, Chapter 7 extends the ab-initio method
to explore the dynamics of diatomic molecules within the Born-Oppenheimer approximation. We
present preliminary results for the ab-initio solution of the TDSE for the hydrogen molecule in the
SAE approximation.

Chapter 2

Time-dependent Schrödinger equation

The interaction of atomic systems with electromagnetic fields can be described using thetime-dependent Schrödinger equation (TDSE). The TDSE is given byi∂Ψ = ĤΨ∂t(2.1)

with Ĥ = Ĥo + Vint ,where Ψ is the wave function, Ĥ is the Hamiltonian of the system and i is the imaginary number.
The Hamiltonian consists of the atomic system Hamiltonian (Ĥo ) and the interaction potential(Vint ). The atomic system Hamiltonian for a system of Ne electrons and Nn fixed nuclei can be
written as
Ĥo =

Ne
X
i=1

N

n
X
p̂i 2
Zn
−
2
ri,n

n=1

!
+

Ne
X

1
,
r
e
,e
1
2
e <e
1

(2.2)

2

where p̂i is the momentum operator of the ith electron and Zn is the charge of the nth nucleus. The
rei ,ej term corresponds to the distance between the ith and j th electron, and ri,n is the distance
between the ith electron and the nth nuclei. The interaction potential term of the Hamiltonian
represents the interaction of the atomic system with the laser pulse. The interaction term will be
discussed in the next section where we make some approximations to obtain the form that we use in
our work. There are currently no known analytical solutions for the TDSE at the laser parameters
used in our work. Instead, we use numerical solutions to solve the TDSE and to study electron
dynamics driven by the laser fields. Due to the large-scale computation necessary for numerical18
solutions, the TDSE can only be solved numerically for systems with at most two active electrons.
In most cases, we focus on processes that can be modeled by simulating the dynamics of a single
active electron and assuming a mean-field approach for the rest of the system. In the next section,
we present the methodology behind our approach and the approximations used to implement a
numerical solver of the TDSE.

2.1

Approximations

Figure 2.1: Regions of laser parameters where the dipole and non-relativistic approximations arevalid (taken from [52]).

To find a numerical solution of the TDSE we make some approximations. The first one is totreat the electromagnetic field classically using Maxwell’s equations. This approximation holds forlaser pulses with high intensity or large numbers of photons. For the intensities used in our work
(1013 − 1014 W/cm2 ), this approximation holds. Second, we use the dipole and non-relativistic
approximation [52]. In this approximation, we neglect the spatial variation of the laser pulse andassume that it is uniform across the size of the atom. This assumption is valid when the wavelengthof the laser pulse is much larger than the radius of the atom, i.e. λ ≫ ratom . We also neglect the19magnetic component of the field as it has a negligible contribution to the interaction for the laserparameters used in our study. The parameter space for the laser field where these approximations
are valid is referred to as the dipole oasis and is shown in Fig. 2.1.
With the dipole and non-relativistic approximation and a semi-classical treatment of the laser
pulse, the laser interaction term of the Hamiltonian can be written as (either in length or in velocity
gauge):

Vint =

X−E(t) · riiVint =

X(2.3)−p̂i · A(t),iwhere E and A are the electric field and vector potential of the laser field, respectively. To model
systems other than the hydrogen atom, an effective potential must be used. An effective potential
is necessary to include the (static) effects of other electrons in our model without solving the
TDSE directly for more than one electron. We include the effects of other electrons by using thesingle active electron approximation (SAE). In the SAE approximation, all electrons except one areconsidered to be static or frozen in place. The laser is assumed to only interact with the ’active’
electron and the dynamics of that electron are studied. With the SAE approximation, the TDSE
can be written as


∂
1 2
i ψ(r, t) = − ∇ + VSAE (r) + Vint (r, t) ψ(r, t)
∂t
2

(2.4)

where ψ is the wave function of the active electron and Vint is the interaction term described in
Eq. (2.3) for a single electron. The SAE potential is given by

VSAE (r) =−Z−Zc e−cr X−−
an e−bn r .
r
r
n

(2.5)

The constants an , bn , and c are fit parameters. The SAE potentials in Eq. (2.5) are determined byfitting the functional form to a potential obtained from density functional theory (DFT) [53]. The
TDSE Eq. (2.4) is used to model all the laser-atom interactions that we consider in this study.

20

2.2

Ab-initio solution
Ab-initio methods are ones that solve the TDSE directly from first principles. As mentioned

above the TDSE considered for the problems in this thesis does not have an analytical solution andsome approximations are made to solve it numerically. In the ab-initio approach, only a simple
system with one, at most two active electrons can be solved. In this work, only systems with
one active electron will be studied. Here the ab-initio method used to solve the TDSE for all the
calculations in our study is outlined.

2.2.1

Basis and coordinate system

Our numerical solution uses a finite difference method to represent the wave function on agrid in a three-dimensional coordinate system. To exploit the spherical symmetry of the potential in
Eq. (2.5), spherical coordinates (r, ϕ, θ) are used to represent the wave function. In this coordinate
system, we look for a separable solution for the wave function and choose the following basis

ψ(r, t) =X χl (r)lmrYlm (θ, ϕ) =X|lm⟩ ,(2.6)

lm

where χl (r) denotes the radial portion of the wave function, Ylm (θ, ϕ) corresponds to the sphericalharmonics, and l, and m represent the orbital, and magnetic quantum numbers, respectively. We
express the radial component on a grid, while the angular components are expressed using a spherical harmonic basis. The operators in the Hamiltonian are represented by matrices in our chosen
basis, and the matrix elements are computed by evaluating the inner product as

Hlml′ m′ = ⟨lm| H l′ m .

(2.7)

Since the radial portion is represented on a grid, the integral above is done over the solid angle.
The spatial derivatives in the Hamiltonian can be represented using the finite difference method.

21
The first and second-order derivatives in the Hamiltonian can be discretized on the grid asdf (x − 2h) − 8f (x − h) + 8f (x + h) − f (x + 2h)
f (x) ≈
dx
12h

(2.8)

d2
−f (x − 2h) + 16f (x − h) − 30f (x) + 16f (x + h) − f (x + 2h)
f (x) ≈
.
dx2
12h2

(2.9)

The accuracy of the method can be calculated using Taylor expansion of the finite element termswith respect to f (x). For the central difference formula that are shown in Eqs. (2.8) and (2.9),
the accuracy is the grid spacing to the fourth power or O(h4 ). In the upcoming sections, we willcompute the matrix elements for each term in the Hamiltonian. The procedure involves obtaining
an expression for the operators in spherical coordinates and then evaluating the inner product withour chosen basis.

2.2.2

Kinetic energy

The kinetic energy term in the Hamiltonian is the Laplace operator scaled by a constant,which can be written in spherical coordinates as
 2



∇2
∂2
1
∂
1 ∂
∂
1
T =−
= − 2 r 2r +sin(θ)+.
22r∂r
sin(θ) ∂θ∂θ
sin2 (θ) ∂ϕ2(2.10)

The angular part of the kinetic energy term is the angular momentum operator L and we can write
Eq. (2.10) as
T =

L2
1 ∂2
r
−
.
2r ∂r22r2(2.11)

Spherical harmonics are used to expand our wave function which are eigenvectors of the angularmomentum operator. When L acts on the basis we get the following well-known eigenvalues

L2 l′ m′ = l′ (l′ + 1) l′ m′ .

(2.12)

22
Using Eqs. (2.11) and (2.12) the matrix element for the kinetic energy operator is given byTlml′ m′ =∇2 ′ ′
lm
lm
2

Z
T

lml′ m′=χl (r) ∗
Ylm
r




=



1 ∂2
l′ (l′ + 1) ′ ′
lm
r
−
lm
2r ∂r22r2

1 ∂2l′ (l′ + 1)r−
2r ∂r22r2

χl′ (r)Yl′ m′ dΩ.
r(2.13)(2.14)

Using the orthogonality of the spherical harmonics, Eq. (2.14) simplifies toχl (r)
Tlml′ m′ = δll′ mm′
r



1 d2
l′ (l′ + 1)
−
2 dr2
2r2



χl′ (r),r(2.15)

where the derivative will be represented on the grid using the finite difference discretization scheme,
as discussed in section 2.2.1. The matrix has a tri-diagonal and five-diagonal structure respectively
for second and fourth-order finite difference representation of the derivative.

2.2.3

Single active electron potential

The matrix element for the SAE potential given by Eq. (2.5) is much simpler since the singleactive electron potential used in the present work has spherical symmetry. The matrix elements
for the SAE potential is given by*
Vlml′ m′ =Z
Zc e−cr X
lm − −−an e−bn r l′ m′
r
r
n

+
.

(2.16)As before we use the orthogonality of the spherical harmonics and Eq. (2.16) simplifies toχl (r)
Vlml′ m′ = δll′ mm′rz Zc e−cr X− −−an e−bn rrrn!

χl′ (r).
r(2.17)

23
2.2.4

Laser coupling in length gauge

The laser coupling term in length gauge in the dipole and non-relativistic approximation is
written as
Vint (r, t) = −E(t) · r.

(2.18)The matrix elements are more easily calculated if we use a Cartesian coordinate system to expand
the dot product and then convert it back to spherical coordinates when taking the inner productwith our basis. When expanding the dot product we haveVint (r, t) = −Ez z − Ex x − Ey y,(2.19)

where Ei is the electric field in the i-th coordinate. Eq. (2.19) can be expressed in sphericalcoordinates utilizing the expansion of the Cartesian coordinates in spherical harmonics given byr

4π
rY00 (θ, ϕ)
3
r
2π
y=i
r (Y1−1 (θ, ϕ) + Y11 (θ, ϕ))
3
r
2π
x=
r (Y1−1 (θ, ϕ) − Y11 (θ, ϕ)) .
3z=(2.20)

Using Eq. (2.20) we can calculate the matrix elements of the laser coupling term for each axis ofpolarization. For polarization in the z-axis, the matrix elements are*r

Dlml′ m′ = lm −EzrDlml′ m′ = − Ez

4π
3

4π
rY00 l′ m′
3+

Ylm Y00 Yl′ m′ dΩS(2.21)
Z

χl (r) χl′ (r)
r
r
r


.

24
This can be simplified further by noting that the integral of three spherical harmonics can be
written in terms of the Clebsch-Gordan coefficient asZ
4πsYl∗1 m1 Yl∗2 m2 Yl3 m3 = (−1)m1(2l1 + 1)(2l2 + 1)
⟨l1 , 0, l2 , 0|l3 , 0⟩ ⟨l1 , m1 , l2 , m2 |l3 , m3 ⟩ .
4π(2l3 + 1)

(2.22)

When applying Eq. (2.22) to our expression, the matrix element for the laser coupling with a laser
field polarized along the z-axis can be written asrm

Dlml′ m′ = −Ez (−1)(2l + 1)(2l′ + 1)
l, 0, 1, 0|l′ , 0
12πl′ , −m, 1, 0|l′ , m′χl (r) χl′ (r)
r
r
r


.
(2.23)The same procedure can be used to calculate the laser coupling term for fields polarized in thex-axis and y-axis. Using Eqs. (2.20) and (2.22) the matrix elements for x-axis polarization are
given byr

(2l + 1)(2l′ + 1)
l, 0, 1, 0|l′ , 0 ×
12π χl (r) χl′ (r)
′
′
′
′
′
′
l , −m, 1, −1|l , m − l , −m, 1, 1|l , mr.
rr
Dlml′ m′ = −Ex (−1)m(2.24)

Similarly for y-axis polarization we getr(2l + 1)(2l′ + 1)
l, 0, 1, 0|l′ , 0 ×
12π χl (r) χl′ (r)
′
′
′
′
′
′
l , −m, 1, −1|l , m + l , −m, 1, 1|l , mrrr
Dlml′ m′ = −iEy (−1)m(2.25).
The matrix elements show how different angular momentum states are coupled to each other.
The coupling determines the transitions between states that are allowed in our system. This iscommonly referred to as the selection rules. We can get the selection rules for the system by using
the properties of the Clebsch-Gordan coefficients. For fields polarized in the z-axis, the conditionm = m′ and |l′ − l| = 1 must be satisfied to have a non-zero matrix element. For light polarized in

25
the x- and y-axis, the selection rules are |m′ − m| = 1 and |l′ − l| = 1. The selection rules allow us
to know the elements of the laser coupling matrix that are zero, eliminating the need to compute
those. This greatly reduces the amount of computation needed when using ab-initio calculations.

2.2.5

Laser coupling in velocity gauge

The laser coupling term in velocity gauge within the dipole and non-relativistic approximation
is given by
D(r, t) = −p̂ · A(t).

(2.26)In Cartesian coordinates, this can be expressed as

D(r, t) = iAx

∂
∂
∂
+ iAy
+ iAz .
∂x∂y∂z

(2.27)

The calculation for the matrix elements can be simplified by expressing the laser coupling operator
in terms of other operators that have simple matrix elements. We do this by using the commutation
relations for the momentum operator given bypˆi = [∇2 , ri ].

(2.28)

The commutation relation allows us to express the laser coupling in terms of the Laplace andposition operator. We have calculated the matrix elements for the Laplace and position operator
in section 2.2.2 and 2.2.4 respectively and the results are used here. For a field linearly polarized
in the z-axis the matrix elements can be calculated as follows:D(r, t) = iAz∂= iAz [∇2 , z]∂z
Dlml′ m′ =lm iAz

∂ ′ ′
lm
∂z

(2.29)

Dlml′ m′ = lm iAz [∇2 , z] l′ m′ .

(2.30)(2.31)

26
We can now take the inner product to calculate the matrix elements in our basis. Using the resultsfrom sections 2.2.2 and 2.2.4 we haverDlml′ m′ = iAz 2

4π
3



−L2 Y10 (θ, ϕ) + Y10 (θ, ϕ)L2 ′ ′∂lm Y10 (θ, ϕ)
+
lm
∂r
2r





∂
−l′ (l′ + 1) + l(l + 1) ′ ′Dlml′ m′ = iAz 2 lm
lm .
+
∂r
2r
r
(2l + 1)(2l′ + 1)
Dlml′ m′ = iAz 2(−1)ml, 0, 1, 0|l′ , 0 l′ , −m, 1, 0|l′ , m′ ×
12π


∂
−l′ (l′ + 1) + l(l + 1)
χl (r)
χl′ (r),+∂r
2r(2.32)(2.33)(2.34)The matrix elements for the other two polarization axes are calculated using the same procedure.
Here we simply state the results. The matrix elements for x-axis polarization are given byr(2l + 1)(2l′ + 1)
l, 0, 1, 0|l′ , 0l′ , −m, −1, 0|l′ , m′ −
12π



∂
−l′ (l′ + 1) + l(l + 1)′′′
l , −m, 1, 0|l , m χl (r)+χl′ (r)∂r
2r
m

Dlml′ m′ = iAx 2(−1)(2.35)

and for y-axis polarization the matrix elements arer(2l + 1)(2l′ + 1)
l, 0, 1, 0|l′ , 0l′ , −m, −1, 0|l′ , m′ +
12π



∂
−l′ (l′ + 1) + l(l + 1)′′′+
χl′ (r)l , −m, 1, 0|l , m χl (r)∂r2rm

Dlml′ m′ = iAy 2(−1)(2.36)

The first order derivative in the velocity gauge laser coupling term is represented using finite
difference as explained in section 2.2.1. As anticipated due to the gauge invariance of the TDSE,
the Clebsch-Gordan coefficients reveal that the selection rules for both the length gauge and thevelocity gauge are identical. It is important to note that while the ab-initio method explained
here maintains gauge invariance, the numerical properties of the method differ depending on theinteraction Hamiltonian employed to represent the atom-laser coupling. For a converged result
all calculations of the observables are the same irrespective of the gauge. However, the velocity

27
gauge exhibits faster convergence when considering the expansion required in spherical harmonics.
This discrepancy is attributed to the canonical momentum of the Hamiltonian, which is expressed
differently for the length and velocity gauges. The canonical momentum is given byΠ = mṙ for length gauge and(2.37)
Π = mṙ − qA for velocity gauge
In the velocity gauge, the electron’s canonical momentum incorporates the momentum gained fromthe laser field and subtracts it from the kinetic energy momentum component. As a result, themomentum in the velocity gauge is typically much smaller compared to the length gauge. Thisdifference leads to the requirement of fewer terms in the expansion in the velocity gauge, as thedegree of the spherical harmonics is linked to the electron’s momentum.

2.2.6

Bound and continuum states

The states of any atomic (or other) system can be distinguished into bound and continuum
states. Bound states have energies that are discrete and negative while continuum states haveenergies that are positive and form a continuous spectrum. The states correspond to the eigenstatesof the field-free Hamiltonian (Ĥo ). The eigenstates of the Hamiltonian are calculated by solving
the eigenvalue problemĤo ψi = Ei ψi(2.38)

1Ĥo = − ∇2 + VSAE (r),2(2.39)

where ψi is an eigenstate and Ei is the corresponding energy. An important detail in determining
the bound states is the boundary condition that we impose on the states. We defined the states ofthe system before asψ(r, t) =X χl (r)lmrYlm (θ, ϕ).

(2.40)

28
To find the states of the system, we only need to calculate the radial portion of the wave function.
The angular portion is specified by the spherical harmonic functions. To obtain the radial part we
write the equation above asrψ(r, t) =Xχl (r)Ylm (θ, ϕ).

(2.41)

lm

This form puts a constraint on the radial function χ(r) = 0 at r = 0. This gives us the boundary
condition necessary to find a unique solution to the time-independent Schrödinger equation. In
our ab-initio method, bound states are calculated using the Krylov-Schur method provided byan algorithm in the SLEPc library [54]. The matrix elements for the kinetic energy and SAE
potential do not couple states with different orbital quantum numbers and are independent of
the magnetic quantum number. This allows the calculations for eigenstates with different orbital
quantum numbers to be carried out independently and the states can be labeled using the principle
and orbital quantum numbers. For the basis used in our solution, Eq. (2.38) can be written asĤo χl (r) = Ei χl (r)(2.42)

where χl (r) are the radial functions given in Eq. (2.6). The calculations for the continuum states
can be carried out in a similar manner. However, the finite representation of our wave function on
the grid will discretize the energies of the continuum states. With this approach, the continuum
states that can be calculated are limited by the size of the grid and the grid spacing used torepresent the wave functions. Another approach is to use the shooting method [55, 56] to calculate
the continuum states. With the shooting method, a continuum state can be determined for any
choice in energy allowing for a preselected range and resolution in the energies of the states thatare needed. We implement the shooting method as outlined in Ref. [56] to calculate the continuum
states used in our solution of the TDSE. The shooting method is a valuable technique employedin our numerical solution, allowing us to convert a boundary value problem into an initial value
problem. We use ϕkl (r) to denote the continuum radial functions where k and l represent the energy
and angular momentum quantum numbers, respectively. The process begins at the boundary point29
r0 = 0 and proceeds toward the edge of the grid as we propagate the solution. At the boundary
point r0 = 0, the value of the wave function is established based on the boundary condition we
impose on the radial function, as explained above. More specifically, we set ϕkl (r0 ) = 0. To initiate
the iterative procedure of the shooting method, we assign an arbitrary non-zero value to the wave
function at r1 = dr, denoted as ϕkl (r1 ). Specifically, we set ϕkl (r1 ) to 1. This assumption is valid
and provides the correct solution up to a constant factor. Subsequently, we normalize the wave
function, which adjusts ϕkl (r1 ) to the accurate value. Subsequent points on the grid are computed
using the following iterative formula:
2 l(l + 1)
ϕkl (ri+1 ) = ϕkl (ri ) dr
+ 2V (ri ) − 2E + 2 − ϕkl (ri−1 ).
r2(2.43)

Once we have obtained the wave function values using the shooting method, we can normalize
the wave function by utilizing the asymptotic solution for the Coulomb potential. The asymptotic
solution is given by the Coulomb wave function:lπ Z
ϕkl (r ≫ 1) = sin kr −+ ln(2kr) + δkl ,2
k(2.44)

where δkl represents the phase shift calculated in Ref. [56]. As long as our grid is sufficiently large
and the solution obtained through the shooting method matches the Coulomb wave with a phaseshift, it is valid to use the asymptotic solution for normalization. The normalization factor for thecontinuum states is given by:








1


Akl =  s2

ϕ′ kl(r) |ϕkl (r)|2 + (k+
Z
)
kr

.

(2.45)

r=rmaxBy incorporating this factor, we obtain the accurate radial functions for the continuum states.
The phase shift is important when utilizing these continuum waves to compute observables and is30
calculated by matching the phase of our radial function with the asymptotic solution. The phaseis given by






δkl = arg iϕ′ (r)iϕkl (r) + k+kl Zkr(2kr)iZ/k− kr +

lπ 

2

.

(2.46)

r=rmax

To find a detailed explanation for the normalization factor and phase equations please refer toRef. [56]. To validate our work, we have compared the states obtained using the shooting methodwith those obtained through Hamiltonian diagonalization and we obtained the same states. The
shooting method is applicable to atomic systems with spherically symmetric potentials, and we have
employed it for the ionization processes we have studied in the work presented in this thesis. Froma computational standpoint, the shooting method is favored over diagonalization. This preference
arises due to the extensive time and memory resources required for diagonalizing large matrices.

2.2.7

Time propagation

The Crank-Nicolson method is employed to propagate the wave function as a function of time.
By utilizing this method, the wave function undergoes a unitary time evolution, guaranteeing thepreservation of its norm throughout the propagation process. The total Hamiltonian is employed
for the propagation, and each subsequent time step is calculated accordingly. The propagation stepis given byψ(r, t + ∆t) ≈ e−iĤ∆t ψ(r, t),(2.47)

wheree−iĤ∆t ≈1 − i∆2 Ĥ
1 + i∆
2 Ĥ

.

(2.48)

Using Eqs. (2.47) and (2.48) the second order Crank-Nicolson scheme for the time propagation canbe written as




∆
∆
1 + i Ĥ ψ(r, t + ∆t) = 1 − i Ĥ ψ(r, t).
22(2.49)

31
The selection rules outlined in section 2.2.4 result in sparse matrices when applying the CrankNicolson method to the system of equations. These sparse matrices enable the utilization of iterativealgorithms, facilitating efficient time propagation. To achieve this, we employ the PETSc library[57, 58] and its implementation of the Generalized Minimal Residual Method (GMRES) for solving
linear equations during time propagation. The PETSc library simplifies the parallelization of
the Crank-Nicolson method, allowing our calculations to be executed on multiple cores with aperformance that scales linearly.
During time propagation, a portion of the wave function might reach the end of the grid
allocated to represent the wave function. This results in the reflection of outgoing wave packets
that reach the end of the grid. Reflections from the edge cause nonphysical interference effects,
leading to numerical errors in our solution. To remove the reflections, we implement the exterior
complex scaling (ECS) method [59] on the edge of the grid. The ECS method rotates a portion of
the grid into the complex plane leading to an exponential decay of the wave function at the edgeof the grid. We implement the ECS using the procedure outlined in Ref.s [59, 60].

2.2.8

Observables

The time-propagated wave function is used to calculate all the relevant observables for ourstudy. The value of an observable is determined by taking the expectation value of the corresponding
operator with the time-propagated wave function. For operator Ô the expectation value is given
byObservable =D

E2
Ψ Ô Ψ
,(2.50)

where Ψ is the time-propagated wave function. The observables studied in this thesis are thepopulation of bound states, ionization probability, and photoelectron spectrum. The population
for a given state is calculated by projecting the time-propagated wave function on the corresponding

32
bound or continuum state. The population in the bound states is given by

Pϕn = |⟨ϕn |Ψ⟩|2 ,(2.51)

where ϕn is the nth bound state. Since the Crank-Nicolson method preserves the norm of the wave
function, the ionization probability can be calculated by subtracting the population in the boundstates from unity as
Pionization = 1 −X|⟨ϕn |Ψ⟩|2 .

(2.52)n

The photoelectron spectrum can be obtained from the time-propagated wave function as
2

1 X
F (k, ϕ, θ) = 2kZe−iδkl (i)l ϕkl (r)∗ ψ(r, t) dx

∗(ϕ, θ)
Yl,m(2.53)l,m∗ (ϕ, θ) are the spherical harmonics, δ is the phase shift, and ϕ (r) is a continuum state
where Yl,mklkl

with momentum k. The methods for calculating the continuum states and phase shift are outlined
in section 2.2.6. Another approach for calculating the photoelectron spectrum is described in thenext section.

2.3

Quantum trajectory Monte Carlo method
Another approach for calculating the photoelectron spectrum is using quantum trajectoryMonte Carlo (QTMC) simulations. QTMC makes use of the two-step model for simulating theionization process. The two-step model breaks the ionization process into two stages. In the firststage, the electron tunnels through the potential barrier of the atom or molecule. The tunneling is
described in a quantum mechanical framework using the strong field approximation [28], outlined
in section 1.3.3. In the second stage, the ionized electron is treated classically as it interacts withthe laser field and the atomic potential.
The method of QTMC uses the ionization rates discussed in section 1.3.3 to provide a probability distribution for the initial conditions of a Monte Carlo simulation. Ionization rates provide a

33
distribution over ionization time and momentum for the emerging electrons. For each set of ionization time and momentum, a corresponding initial position is then calculated to form a set of initialconditions. Each set of initial conditions will be classically propagated in time using Newton’s
equations of motion. The final momentum at the end of the propagation is then used to calculate
the photoelectron spectrum. During propagation, Feynman’s path-integral approach is used toaccount for the phase of the trajectory. Including the phase in the model distinguishes it from theso-called ’classical trajectory Monte Carlo’ models. Models of QTMC differ from one another basedon the probability distribution used to generate the initial conditions. In the following sections, we
describe two QTMC models and the numerical details involved in the simulations.

2.3.1

Adiabatic model

The ionization rate used for the adiabatic model is the Ammosov, Delone, and Krainov (ADK)
tunneling rate [21, 22]. This model works in the adiabatic limit, where the field is considered to bestatic during the ionization process. Details about the ADK model are explained in section 1.3.3.
The ADK ionization rate is given by
κ (v⊥ip − v⊥ )2κ3
WADK (to , ⃗v ) ∝ exp −exp −,3|E(to )||E(to )|


(2.54)

where to is the time of ionization, ⃗vi is the initial velocity vector, v⊥ip is the initial velocity perpendicular to the field in the plane of polarization, v⊥ is the initial velocity in the direction ofpthe electric field’s propagation, and κ = 2Ip . In the ADK model, the tunneled electrons have a
Gaussian distribution for the transversal velocity perpendicular to the instantaneous laser field. In
accordance with the adiabatic limit, the velocity component parallel to the instantaneous laser field
is approximated to be zero. We use the ADK distribution to derive a probability distribution for
the initial momentum and the time of ionization. The initial value for the phase space also requires
the electron’s exit point (initial position after tunneling). The exit point for each trajectory is
related to its ionization time and initial velocity. In the adiabatic limit, the TIPIS model (tunnel

34
ionization in parabolic coordinates with induced dipole and Stark shift) [61–63] can be used to
calculate the exit point. The TIPIS model uses parabolic coordinates to separate the Schrödinger
equation for an atom in the presence of a static field. This results in a one-dimensional problem
along the axis of the field where the electron can ionize. The problem can be solved analytically to
get the position after tunneling. The exit point from the TIPIS model is given byIp +rexit =q
Ip2 − 4βE(to )2E(to)

p
2Ip
with β = 1 −.
2(2.55)

With the exit radius given by the Eq. (2.55), the exit point in Cartesian coordinates can be writtenasxi =−Ei (t)rexit ,
|E(t)|(2.56)

where xi is exit position in the ith coordinate. The tunneling exit point is located at the opposite
side of the direction of the instantaneous field where the barrier is lowest. The initial velocity andionization time (Eq. (2.54)) along with the exit position from the TIPIS model provide the initialcondition for the adiabatic QTMC simulation.

2.3.2

Non-adiabatic model

To include non-adiabatic effects in QTMC, the PPT model [23, 24, 64] of ionization is used.
The PPT and other non-adiabatic models are described in section 1.3.3 in detail. The ionization
rate used in the non-adiabatic model is given bys
W (ts , ⃗v ) ∝ eiS(ts ,⃗v (ts ))

2

2π
′′
S (ts , ⃗v (ts )),(2.57)

where ts is the complex transition time and ⃗v is the initial velocity. S is the action defined in Eq.
(1.7). As explained in section 1.3.3, the ionization rate is evaluated at the complex transition time
which satisfies the saddle point equation. Similar to the adiabatic model discussed in section 2.3.1,
the non-adiabatic ionization rate is used to generate a set of initial conditions for the Monte Carlo35
simulation. Unlike the adiabatic model, the distribution depends on the saddle points (complex
transition times), and the velocity parallel to the instantaneous field is not assumed to be zero. In
this section, we outline how the saddle points can be calculated for a general elliptically polarized
pulse. We then show how the saddle points are used to calculate the velocity parallel to theinstantaneous field and the exit position of the electron.
The complex transition times can be split into real and imaginary parts as ts = t0 + iti .
Where to is the ionization time and ti is thought of as the tunneling time [20]. The saddle points
can be calculated analytically as a function of known quantities, namely the ionization time andinitial perpendicular momentum [65]. The saddle point equation can be written as
1[p + A(ts )]2 + Ip = 02(2.58)

where p is the conserved canonical momentum and A is the vector potential. We choose a laser
pulse polarized in the x − y plane with the vector potential given by

A(t) = Ao (t) sin (ωt)x̂ + ϵAo (t) cos (ωt)ŷ,(2.59)

where Ao (t) is the pulse envelope. The saddle point equation for an elliptically polarized pulse isgiven by1
(px + Ao (to ) sin (ωto ) cosh (ωti ) + iAo (t0 ) cos (ωto ) sinh (ωti ))22
1+ (py + ϵAo (to ) cos (ωto ) cosh (ωti ) − iϵAo (to ) sin (ωto ) sinh (ωti ))221+ p2z + Ip = 0.
2(2.60)

Since we want to solve for the initial velocity, we must relate the velocity to the conserved
canonical momentum. At the time of ionization, the relationship between the conserved momentum

36
and initial velocity is given by p = v − A, which for the laser pulse introduced above becomespx = vx − Ao (t) sin ωtpy = vy − ϵAo (t) cos ωt(2.61)

pz = vz .
Before substituting the above expression in Eq. (2.60), the problem can be simplified by performinga change of basis. We transform the problem from Cartesian coordinates into a coordinate systemthat rotates with the laser field. The basis in the rotating coordinate system has components of[v|| , v⊥ip , v⊥ ], where v|| and v⊥ip are the velocity components parallel and perpendicular to theinstantaneous laser field in the plane of polarization and v⊥ is the momentum component in the
direction of propagation. The equations for the change of basis are given byvx = v|| cos β − v⊥ip sin β(2.62)
vy = v|| sin β + v⊥ip cos β
where β = tan−1 (ϵ tan (ωto )), represents the angle between the major axis of polarization and theinstantaneous laser field. In Fig. 2.2, a circularly polarized pulse is shown to illustrate the conceptof changing the basis to the rotating frame. The arrows indicate the direction of polarization (v|| ),while v⊥ip represents a vector component orthogonal to v|| within the plane of polarization. The
motivation behind this change of basis is to reduce the number of unknowns in the saddle point
equation. By utilizing the new basis, we can express the parallel component of velocity in terms of
the saddle points and other velocity components.

37

Figure 2.2: Illustration of the polarization direction of a right-handed circularly polarized light(RCPL) pulse, as it changes as a function of time (taken from [66]).

In the new basis, the canonical momentum relation given in Eq. (2.61) can be expressed asfollows:"v||v⊥ ϵ tan ωto#−p+ Ao (to ) sin (ωto )1 + ϵ2 tan2 ωto"#v|| ϵ tan (ωto )v⊥+p− ϵAo (to ) cos (ωto )py = sgn(Ex ) p
1 + ϵ2 tan2 ωto1 + ϵ2 tan2 ωtopx = sgn(Ex ) p

1 + ϵ2 tan2 ωto(2.63)

p z = vz
where sgn(Ex ) is the sign of the electric field in the major axis (x−axis in our definition of the laser
pulse). The next step is to substitute the expression for the canonical momentum in Eq. (2.60) andsolve for the saddle points. This can be done numerically with a root-finding algorithm as suggested
in Ref. [67]. But in practice, this method is computationally too expensive for calculations that
require large amounts of trajectories. It is also unnecessary as an analytical solution for the saddle
point can be derived after some manipulation of the saddle point equation. The expression for the38
saddle points for an elliptically polarized laser pulse is given by [65]
#
"

r



2

2

γ

ω
1

,
a ω v⊥ − ϵ + (a4 − ϵ2 ) 1 + eff 4 4 ϵ a v⊥ − ϵ ± a22
a −ϵcosh ωti =

Eo

Eo

a




 h 
i


2

aω
 1 1 − aω v⊥ + a22 1 + γeff
/
2
1
−
v
,
⊥
2
2

ϵEo

ϵ

a

ϵEo

a2 ̸= |ϵ|a2 = |ϵ|(2.64)
pwhere γeff = ω 2(Ip + vz2 /2)/Eo is the effective Keldysh parameter [68] and a = |E(to )|/Eo is thenormalized instantaneous electric field.
Once the saddle points have been obtained, it becomes feasible to calculate both the ionization rate and the remaining initial conditions, including the parallel momentum and the exit
points. To solve for the parallel momentum (v|| ) we utilize expressions for the saddle point andthe ionization time given above. Through algebraic manipulations, the expression for the initial
parallel momentum can be expressed as follows [65]:v|| =(1 − ϵ2 )Eo sin ωto cos ωto(cosh ωti − 1).
aω(2.65)The ionization rate can be calculated by performing the integral in Eq. (1.7) for the action. For
an elliptically polarized pulse the ionization rate is given byp2
Eo
W (to , ti , vo ) = exp −2 ( + Ip + Up )ti + 2px 2 sin (ωto ) sinh (ωti )
2
ω
 2

2
1 − ϵ FoϵEocos (2ωto ) sinh (2ωti ) ,−2py 2 cos (ωto ) sinh (ωti ) +
ω
4ω 3

2

(2.66)

2)Eo
where Up = (1+ϵis the ponderomotive energy and p is the magnitude of the canonical momen4ω 2

tum. For the expression of the ionization rate, we revert back to Cartesian coordinates for a more
simplified expression. In practice when implementing the QTMC calculation we use the rotating
frame and apply Eq. (2.63) to derive the Cartesian components of the momentum for calculating
the ionization rate for a set of initial conditions.
The exit point can be calculated by considering the sub-barrier trajectory [69] where the39electron tunnels through the barrier from the complex transition time to the time of ionization.
The trajectory can be represented by an integral given as
Z tro (t) =
p + A(t′ ) dt′ .

(2.67)

ts

The result of performing the integral in Eq. (2.67) and taking the real part gives the exit point foreach coordinate in terms of known quantities. The exit point in Cartesian coordinates is given byEocos ωto (1 − cosh ωti )ω2ϵEoyo (to ) = 2 sin ωto (1 − cosh ωti )ω
xo (to ) =(2.68)zo (to ) = 0

Figure 2.3: Comparison of the exit radius for the TIPIS model and the non-adiabatic model [63].
The Figure shows the exit radius as a function of the electric field strength, both expressed in a.u.

The exit radius calculated from the adiabatic and non-adiabatic models is shown in Fig. 2.3 as afunction of the electric field strength. We see that the exit radii provided by both models are in good
agreement with each other, within 2.2% [63]. The major difference in the prediction of the model
comes from the parallel component of the initial momentum that leads to different photoelectron

40
spectra [65, 67, 69–74]. With the exit points calculated, all the necessary initial conditions for the
non-adiabatic Monte Carlo simulations are obtained.

2.3.3

Time integration

Once a set of initial conditions is obtained from the ionization models, trajectories are propagated in time to obtain the final momentum. The propagation is done using Newton’s equation of
motion for the combined potential of the ion and the electric field using the following expression:d2 r(t)Zr(t)= −E(t) − 3 .
2
dt
r (t)

(2.69)

We split the time propagation into two intervals. The first interval is for the length of the pulse(to , τf ), and the second is for the propagation after the pulse has ended (τf , ∞). We calculate
the propagation in the first interval numerically using an ODE solver that implements the RungeKutta algorithm with adaptive time steps. In the second interval, the particle moves only under the
influence of the Coulomb potential. The trajectory follows a hyperbolic motion that can be solved
analytically [75, 76] to give the asymptotic momentum. The values of the position and momentum
after the first time interval are used to calculate the asymptotic momentum that will be measured
at a detector. The phase space coordinates at the end of the pulse are written as [r(τf ), p(τf )].
The asymptotic momentum for the particles will then be given byp⃗asy = p

⃗ f ) × ⃗a(τf )) − ⃗a(τf )p(L(τ1 + p2 L2 (τf )(2.70)

where p⃗asy is the asymptotic momentum, p is the magnitude of the momentum at the end of thepulse, L is the angular momentum and ⃗a is the Runge-Lenz vector given by⃗ f ) − ⃗r(τf ) .
⃗a(τf ) = p⃗(τf ) × L(τr(τf )(2.71)

41
To model the interference effects in quantum mechanics, the phase accumulated by the particles
along their trajectories is included [70, 71, 77]. This is done using Feynman’s path integral, where
the phase accumulated by the electron is given by
Z ∞ 2p⃗ (t) 2Z
Φ(to , ro , vo ) = −po , ro + Ip to −
dt
−
2
r
to

(2.72)

where po is the initial momentum, ro is the initial position and to is the time of ionization. The
phase given in Eq. (2.72) is derived in detail in Ref. [71].
The photoelectron spectrum is calculated by binning the trajectories according to their final
momentum. Each cell in momentum space contains trajectories with [ki , ki +∆ki ], where i = x, y, z.
Results in the same bin are added coherently and the photoelectron spectrum is given byF (⃗k) =

Nk q
X

2
j j
w(tjo , ⃗voj )eiΦ(to ,⃗vo ),(2.73)

j=1

where w(tjo , ⃗voj ) is the probability distribution of the initial condition, Φ(tjo , ⃗voj ) is the phase andq
Nk is the total number of trajectories in the k-bin. The w(tjo , ⃗voj ) term in Eq. (2.73) weights
each trajectory based on the probability associated with the initial momentum and ionization
time. Typically, we utilize uniform sampling of initial time and momentum while incorporating
the ionization rate as a weight to account for the likelihood of electron ionization. However, thisapproach is computationally demanding as a significant number of initial values result in negligible
weights and do not contribute to the spectrum. Hence, it is unnecessary to compute trajectories forthese specific initial values. To effectively address this issue, we employ a method outlined in theinitial part of this Chapter. The method involves sampling the initial distribution using ionization
rates. However, a significant challenge arises in generating a probability distribution based on theseionization rates. To tackle this challenge, we utilize the reject-sampling algorithm [78]. Here, we
provide a description of the algorithm for generating a one-dimensional distribution, although themethod can be easily generalized to higher dimensions. The fundamental idea is that instead of

42
directly sampling a random variable in one dimension, we can achieve this by conducting uniform
random sampling in two dimensions and then filtering out only the samples that lie within thearea beneath the graph of the density function. This process is illustrated in Fig. 2.4, where thered dots represent the initial uniform distribution and the black line denotes the density function.
The green points are the accepted points that are distributed according to the density function.
As the number of points used to implement this method becomes exceedingly large, the initialconditions will be generated in accordance with the probability function. We apply this method inhigher dimensions for the set of initial conditions that our ionization rate depends on. If we use theionization rate to generate the set of initial conditions, then the photoelectron spectrum is givenbyF (⃗k) =

Nt
X

2

e

iΦ(tjo ,⃗voj ).

(2.74)

j=1

Figure 2.4: Illustration of the reject-sampling method for sampling a distribution. The Figure showsa scatter plot of events in the X-Y space, where red dots represent rejected samples (Yi > P (Xi ))and green dots represent accepted samples (Yi ≤ P (Xi )).

43
By employing the reject-sampling method, the convergence rate is greatly improved as it
generates trajectory samples based on their probability of occurrence. As a result, we obtain aset of initial conditions for trajectories that have a high probability to occur. This leads to asignificant reduction in the number of time propagations done for trajectories that do not make a
substantial contribution to the spectrum. In general, the reject-sampling method typically requires
approximately 106 to 107 trajectories to compute the photoelectron spectrum, whereas the uniform
sampling method necessitates around 109 trajectories. The exact number of trajectories necessary
depends on the specific laser parameters employed in the atom-laser interaction. Further details
on this matter will be provided in Chapter 6.

2.4

Laser pulses
In both the ab-initio solution and the QTMC simulations, we simulate the interaction between

atoms and laser pulses. In this section, we explain how we define the laser pulses that we use in our
simulations. Laser pulses release light in short bursts lasting from attoseconds to milliseconds. As
the pulses are time-limited, they contain photons in a frequency range around the central frequency.
Defining the laser field can be done either by using its electric field or vector potential. For pulsed
lasers with an envelope shape, defining the laser pulse through its electric field results in a non-zero
vector potential. To ensure that the potential disappears at the end of the pulse, as stated inRef. [79], we first define the vector potential and then derive the electric field from it. The vector
potential is defined as

A(t) = √Ao[cos(ωA t + ϕ)x̂ + ϵ sin(ωA t + ϕ)ŷ] f (t)1 + ϵ2(2.75)

where Ao is the electric potential amplitude, ωA is the central frequency and ϕ is the carrier-envelope
phase. The laser field is polarized in the x − y plane and propagates in the ẑ direction. ϵ is theellipticity of the field that determines the polarization state of the laser light. The intensity of laserqpulses is connected to the amplitude of the electric potential via the equation Ao = IoIω , where44
Io = 3.51 × 1016 W/cm2 . The pulse envelope function f (t) governs the pulse shape, and there
are various envelope functions available. Two commonly used envelope functions are provided as
examples in Eq. (2.76). An example of a sin-squared laser pulse is depicted in Fig. 2.5 where weshow both the field (blue) and the envelope (red) of the pulse.

f (t) =


sin2 ( πtτ ),sin-squared envelope(2.76)

2exp − ln(2) ( 2tτ ) ,Gaussian envelope

Figure 2.5: Temporal profile of a laser pulse at 200 nm wavelength, 1 × 1014 W/cm² intensity, and15 cycles. The red line depicts the envelope function, while the blue line represents the electricfield.

The electric field can be calculated from the vector potential using the relation

E(t) = −

1∂
A(t).
c ∂t

(2.77)

45
For a simple linearly polarized laser pulse such as

A(t) = Ao f (t) cos(ωA t + ϕ)ẑ,(2.78)

the electric field is given byE(t) = −1∂∂f (t)
Eocos(ωA t + ϕ)A(t) = Eo f (t) sin(ωA t + ϕ) −
ẑ.
c ∂t
ωA
∂t

(2.79)

where Eo = ωA Ao . A consequence of defining the electric field in this manner is that the central
frequencies of the electric field and vector potential are not the same [80]. This is due to the second
term in Eq. (2.79) leading to a frequency shift defined by |ωE − ωA |. As shown in Ref. [80], theratio of the frequencies is given by
1+
ωE
≈
ωA

sR ∞
p
2
1 + 4(πN )−2ωA−∞ (t − to ) f (t)dt
R∞.
, where N =2π−∞ f (t)dt(2.80)In the above equation N is the number of cycles that are with one standard deviation of the peakof the pulse and to is given byR∞
to = R−∞∞tf (t)dt−∞ f (t)dt.

(2.81)The ratio of the frequencies given by Eq. (2.80) depends on the pulse shape and number of cyclesof the pulse. It is, however, independent of other laser parameters. The frequency shift becomes
negligible as the number of cycles of the pulse increases. To properly model the laser couplinginteraction as described in section 2.2, one needs to use the corrected frequency for the vector
potential. Further information about the frequency shift and its effects on the solution to the
TDSE can be found in Ref. [81].

Chapter 3

Pathways to excitation of atoms with bicircular laser pulses

3.1

Introduction1
Recently, studies on the interaction of atoms and molecules with intense fields generated

by the superposition of two circularly polarized laser pulses have seen an upsurge in activity instrong-field experiment and theory. For the most part, the renewed interest results from the capability to control the polarization of emitted light in high-order harmonic generation with such
pulses. The physical principle has been proposed and applied first two decades ago [82, 83]. Morerecently, efficient phase matching of circularly polarized high-order harmonic beams in the EUV
and soft X-ray regime using bichromatic beams with counterrotating circular polarization has been
demonstrated [84–89]. Since then, much experimental and theoretical work on high-harmonic generation [90–114], ionization and photoelectron momentum distributions [47–49, 115–134], double
ionization [135–144] and other strong-field processes [145–147] driven by bichromatic circularly
polarized laser pulses has been performed. One interesting aspect in these kind of strong-fieldinteractions is the control of ionization via the helicity of the applied bichromatic pulses [47–49].
Such studies complement related work on the dependence of the ionization rate by a one-color
circularly polarized pulse on the relative helicity between the pulse and the electron in the atomic
orbital [30, 148–165].
For bichromatic circularly polarized laser pulses it has been observed that the probability to
1

Part of the material presented in this Chapter has been published in J. Venzke, Y. Gebre, A. Becker, and A.
Jaroń-Becker, Physical Review A 101, 053425 (2020). Y. Gebre and J. Venzke contributed equally to the work
presented in the publication.

47
ionize an atom is significantly enhanced if the two fields are counterrotating as compared to corotating fields [48]. The experimental observations were interpreted as due to the increased density of
excited states accessible for resonant enhanced multiphoton ionization in the case counterrotating
fields. Results of numerical solutions of the time-dependent Schrödinger equation in Ref. [48] did
confirm a close relation between the ratios of total excitation and ionization probabilities for counterrotating and corotating circularly polarized laser pulses. However, the results for excitation of
the atom were not further resolved by distributions over the quantum numbers (principal, angularmomentum, magnetic). Such analysis potentially can shed further light on the role of excited statesin the pathways to ionization since excitation in a resonant multiphoton process should rely on the
spin-angular momentum selection rules for the absorption of circularly polarized photons (∆l = ±1
and ∆m = ±1).
More generally, analysis of the role of strong-field excitation has recently experienced a renaissance [166–169] following earlier work [170–172]. Concerning the distribution in the excitedstates with respect to the quantum numbers studies for the interaction of atoms with linearlypolarized pulses have been performed only. Theoretical studies have considered the distributionof the population as a function of the principal and/or the angular momentum quantum number [38–40, 44, 81, 173]. In applications of Floquet theory for a monochromatic laser field [44]and numerical calculations for laser pulses with trapezoidal [40] and Gaussian or sin-squared envelopes [81] it has been analyzed how the parity of the populated angular momentum states in such
pulses relates to the selection rules for linearly polarized pulses. In view of the recent experimental
observations discussed above, we extend these studies to interaction of atoms with bichromatic circularly polarized laser pulses. Such study provides the interesting opportunity to not only resolve
the excited state distribution with respect to the principal or angular orbital momentum quantum
number, but in particular to consider the role of the magnetic quantum number as well. For ourstudies we make use of results of numerical solutions of the time-dependent Schrödinger equation(TDSE) for the interaction of the hydrogen atom with intense bichromatic circularly polarized laser
pulses.

48

The rest of this Chapter is organized as follows: The results of the numerical calculations arepresented and analyzed in section 3.2, first for corotating and then for couter-rotating pulses. In
section 3.3 we summarize the insights gained into the excitation pathways in bircircular pulses.

3.2

Results
In this section we present the results for the distributions, first for the corotating and thenfor the counterrotating case, which provide insights into selection rules and excitation pathways inbichromatic multiphoton processes. The interaction with the bicircular laser pulse is implemented
via the total vector potential as:
A(t) = Aω (t) + A2ω (t)

(3.1)

where

2



AΩ (t) = A0,Ω sin

πt
τΩ



× [sin (Ωt) x̂ + ϵΩ cos (Ωt)ŷ]

(3.2)

√

Ω
for Ω = ω and 2ω, respectively. A0,Ω = c ΩIΩ , τΩ = 2πNΩ , and c is the speed of light, where IΩ

is the peak intensity and NΩ denotes the number of cycles. ϵΩ = ±1 denotes the helicity of thefundamental and 2nd harmonic pulse, respectively.

3.2.1

Excitation with corotating pulses

Selection rules for (single-)photon absorption from circularly polarized light are given by∆l = ±1 and ∆m = ±1, where the change in the magnetic quantum number is positive (negative) if
the helicity of the light is right- (left-)handed. Extending the concept to multiphoton absorption, thesimultaneous change in both quantum numbers puts distinct constraints on the parity and helicityof the accessible excited states in the atoms upon absorption of multiple photons. Specifically, it
is expected that states in which ℓ and m are either both even or both odd are being populated

49
during the interaction with the field. This selection holds for the interaction with a single circularly
polarized pulse as well as for the case of a superposition of two (or more) of such fields, independentof the relative helicity of the two pulses.

Figure 3.1: Excited state distribution as function of n (vertical axis) and ℓ (horizontal axis) for(a) m = −1, (b) m = −2, (c) m = −3 and (d) m = −4 at the end of 20 (at 800 nm) cycle pulses(40 cycle at 400 nm) with sin squared envelope and total peak intensity of 1 × 1014 W/cm2 forcorotating laser pulses of equal intensity (taken from [174]).

In Fig. 3.1 we show examples of the population in the excited states of hydrogen atom as afunction of n and ℓ for various m values at the end of the interaction with bichromatic corotating lefthanded circularly polarized pulses. The results clearly confirm the expected population distribution
in states with either odd or even parity for a given value of m according to the selection rules upon
multiphoton absorption. The present results have been obtained for interaction with equal peakintensities I400 = I800 = 5 × 1013 W/cm2 .

50

Figure 3.2: Absorption pathways in corotating laser pulses at frequencies ω and 2ω starting froma m = 0-state. Without lack of generalization it is assumed that both pulses have left-handed
helicity. Absorption of a photon at frequency ω and at frequency 2ω is represented by a red andblue arrow, respectively. The numbers in the boxes denote the minimum number of photons toreach a certain level (taken from [174]).

Due to the correlation in changes of m and ℓ the observed pattern is independent of total
peak intensity, ratio of peak intensities and pulse duration, as long as the dipole approximation
holds. In the present study we have verified this up to intensities of 1 × 1014 W/cm2 . This isdifferent from the case of linear polarization [81], where selective population concerning the parityof the populated excited states is observed for long pulses and low peak intensities only. In thatcase the restriction to a given m-channel and a broad energy spectrum (for short pulse durations)or a significant Stark shift of the excited states (at high peak intensities) leads to a mixing of
population over the states with odd and even orbital angular momentum quantum numbers.
In corotating bicircular laser pulses all the photons have the same spin (either +1 or −1),
consequently the magnetic quantum number always changes either by ∆m = +1 or by ∆m = −1
upon absorption of each photon. For our studies we have chosen left-handed helicity for both pulses
and, hence, only excited states with negative m can be populated upon absorption of photons fromthe ground state with m = 0 (c.f., Fig. 3.2). Therefore, as already mentioned in Ref. [118], only
Rydberg states with high orbital angular quantum number ℓ are accessible. For example, for
excitation of Rydberg states (with n ≥ 4) in the hydrogen atom, the absorption of at least 4

51photons in laser field at 400 nm or at least 8 photons at 800 nm is required. Thus, Rydberg stateswith ℓ < 4 (and m > −4) cannot be populated just via photon absorption alone.

Figure 3.3: Excited state distribution as function of orbital angular quantum number ℓ summed
over n ≥ 4 and m at (a) I400 = 5×1013 W/cm2 , I800 = 5×1012 W/cm2 , (b) I400 = 5×1013 W/cm2 ,I800 = 1 × 1013 W/cm2 , (c) I400 = 5 × 1013 W/cm2 , I800 = 5 × 1013 W/cm2 , and (d) I400 = 1 × 1013
W/cm2 , I800 = 5 × 1013 W/cm2 . Pulse durations: 20 cycles at 400 nm, 10 cycles at 800 nm (takenfrom [174]).

Accordingly, the angular momentum distribution in the Rydberg states is controlled via therelative intensity of the two fields at the fundamental and second harmonic frequency. This isdemonstrated in Fig. 3.3, where the excited state distribution as a function of ℓ, summed for n ≥ 4
and all m, is shown. For large ratio of I400 /I800 = 10 (panel (a)) the distribution is centered, as
expected, about ℓ = 4. As the intensity ratio decreases, high orbital angular momentum states get
increasingly populated due to the impact of the laser pulse at 800 nm.

52

Figure 3.4: Excited state distribution as function of n (vertical axis) and m (horizontal axis)
summed over ℓ. Laser parameters: 20 (800 nm) cycle pulses with sin squared envelope and total
peak intensity of 1 × 1014 W/cm2 for corotating laser pulses of equal intensity (taken from [174]).

Another interesting feature in Fig. 3.3 is that the population in angular momentum states withℓ < 4 increases significantly when the intensities of the two pulses are similar. Further insight canbe gained by the distribution over the magnetic quantum number, which is displayed in Fig. 3.4(c)for the case of equal intensities. It is clearly seen that Rydberg states with magnetic quantum
numbers between m = 0 and m = −3 are populated. In view of the number of photons needed
to reach the Rydberg levels, the population in these states cannot be explained by absorption ofphotons only.
Instead, we propose the following mechanism: Initially, Rydberg states with ℓ ≥ 4 are populated via the absorption pathways shown in Fig. 3.2. Then a redistribution of population occurs
via Raman-type Λ-transitions (c.f., [175, 176]). In the present bichromatic laser field the Λ-processleads to a change in the magnetic quantum number, if photons from both fields are involved. For
the absorption of one 400 nm photon and emission of two 800 nm photons the magnetic quantum
number between initial and final state changes by ∆m = +1.
The order of absorption and emission may vary, i.e., the redistribution process can either
proceed via the continuum (absorption first, Fig. 3.4(a)) or via a lower excited state (emission first,
Fig. 3.4(b)). A larger change in m is achieved either via a sequence of these Λ-processes or by

53
higher order processes (e.g., absorption of two 400 nm photons followed by emission of four 800 nm
photons leading to ∆m = +2). We note that similarly the absorption of two photons at 800 nm
and the emission of a 400 nm photon will lead to a change of ∆m = −1 in the present set-up and,hence, contribute to population of states with higher ℓ and m.

Figure 3.5: Distribution in magnetic quantum states for (a) ℓ = 2 and (b) ℓ = 3 and different peakintensities of the 800 nm pulse. I400 = 5 × 1013 W/cm2 and other parameters are as in Fig. 3.4(taken from [174]).

Our interpretation is further supported by the results in Fig. 3.5, which shows how the
population in states of certain quantum numbers for (a) ℓ = 2 and (b) ℓ = 3 changes as function of
the relative intensity of the two pulses. It is clearly seen that the population in these quantum states,
which are not accessible via direct absorption of photons from the ground state, increases as theintensity of the pulse at 800 nm increases. Thus, these results provide further indications that thepresence of the redistribution process depends on the impact of both pulses and its effectivenessincreases with increase of the total intensity, in agreement with our interpretation of a Λ-type
process. Further analysis of this new phenomenon will be presented in the next Chapter 4.

54
3.2.2

Excitation with counterrotating pulses

Figure 3.6: Same as Fig.3.1 but for counterrotating laser pulses (taken from [174]).

As discussed in the previous subsection, the selection rules by which only states with ℓ andm either both even or both odd hold independent of the relative helicity of the two pulses. Thisis confirmed by the results that we obtained for the interaction with two counterrotating pulses at
equal intensities and equal pulse duration presented in Fig. 3.6. Depending on whether m is evenor odd, the distribution over the orbital angular momentum shows population in states with evenor odd parity. As in the case of corotating pulses, the observed pattern is found independent of
total peak intensity, ratio of peak intensities and pulse duration.

55

Figure 3.7: Absorption pathways in counterrotating laser pulses at frequencies ω and 2ω startingfrom a m = 0-state. Without lack of generalization it is assumed that the pulse at frequency ω hasleft-handed helicity, while the second harmonic pulse has right-handed velocity. Other symbols as
in Fig. 3.2 (taken from [174]).

Since in counterrotating bicircular laser pulses photons of the two fields have opposite spin,
starting from the ground state with m = 0, excited states with both positive and negative magnetic
quantum numbers can be populated. The absorption pathways for the set-up chosen in the present
studies, namely right-handed helicity for the 800 nm pulse and left-handed helicity for the second
harmonic, are shown in Fig. 3.7. As can be seen from the Figure, the magnetic quantum number reflects the difference between the number of 400 nm photons and that at 800 nm absorbed.
Furthermore, it can be seen that for a given total photon energy absorbed states with magnetic
quantum numbers separated by ∆m = ±3 are populated.

56

Figure 3.8: Excited state distributions as a function of m, summed over n ≥ 4 and ℓ (top), and as
a function of n and m, summed over ℓ (bottom), for the interaction with a left-handed circularly
polarized laser pulse at 800 nm (20 cycles) and a right-handed circularly polarized laser pulse at
400 nm (40 cycles). Both pulses have the same peak intensity of 5×1013 W/cm2 (taken from [174]).

These features are clearly present in the population distributions as function of m, summed
over n and ℓ (top), and of n and m, summed over ℓ (bottom), in Fig. 3.8, which show the results
for counterrotating pulses of equal peak intensity. In the Rydberg manifold (n ≥ 4) the highest
populated m-states differ by ∆m = ±3, other states show some but lower population as the manifold
is AC-stark shifted during the interaction with the pulses. In view of the nonlinearity of multiphoton
processes, it is likely that the states showing the largest probability are being populated near the
peak of the pulses at which the highest total intensity is present. Overall, the strongest populationis seen for states with negative magnetic quantum numbers, leading to the conclusion that it most
likely that either five (for excited states with m = −5) or two (for excited states with m = −2)
more 800 nm photons with left-handed helicity than 400 nm photons with right-handed helicity arebeing absorbed.

57

Figure 3.9: Orbital angular momentum distributions in excited states induced by counterrotating
laser pulses at 400 nm (20 cycles) and 800 nm (10 cycles) at peak intensities of (a) I400 = 5 × 1013
W/cm2 , I800 = 5 × 1012 W/cm2 , (b) I400 = 5 × 1012 W/cm2 , I800 = 5 × 1013 W/cm2 , and (c)I400 = 5 × 1013 W/cm2 , I800 = 5 × 1013 W/cm2 (taken from [174]).

The distributions in Fig. 3.8 do not extend much beyond |m| = 5, which is consistent withthe results shown in Fig. 3.9. Thus, there appears to be a highest orbital angular momentum
number ℓmax beyond which the population in the states drops off quickly. This is in agreementwith previous studies for Rydberg state excitation [40, 81] and low energy angular momentum
distributions [177]. In Ref. [177] a random walk analysis of the absorption pathways between theaccessible quantum states is used to obtain the classical orbital angular momentum for a electron√with zero energy in a laser field has been derived as L = 2Zα0 where Z is the charge of theresidual ion and α0 is the quiver radius. Relating classical orbital angular momentum and theorbital angular momentum quantum number by ℓ ≈ L − 1/2 we have estimated the maximum

58ℓ gained in the bicircular counterrotating pulse. The estimates, shown by the solid red lines inFig. 3.9, are in good agreement with the cut-offs seen in the numerical results. We note that the
random walk analysis of Ref. [177] can be applied in the case of counterrotating pulses, since in eachabsorption step ∆m = ±1 and hence, in general, ∆ℓ = ±1 is possible. In contrast, for corotating
pulses the changes in magnetic and angular quantum number are determined in each absorption
step (∆m = −1, ∆ℓ = +1) and hence a random walk analysis is not applicable and a cut-off cannot
be derived.

3.3

Summary
We have studied the distributions over the orbital angular momentum (ℓ) and magnetic (m)quantum numbers in Rydberg states due to the interaction with bichromatic circularly polarizedlaser pulses. Multiphoton selection rules lead to population of states in which ℓ and m are either
both even or both odd, independent of relative helicity, peak intensity and pulse duration. In thecase of corotating pulses the results show that the distribution over the magnetic quantum number
can be controlled via the intensities of the two pulses. Furthermore, we propose that the states
are populated via direct absorption from the ground state and via Λ-type transitions between
Rydberg states of different ℓ and m, involving two photons at the fundamental wavelength andone photon at the second harmonic, which we will further analyze in the next Chapter 4. For
bicircular laser pulses with opposite helicities Rydberg states with magnetic quantum numbersthat differ by ∆m = ±3 are predominantly populated. The pattern allows to gain insights into therelative number of photons absorbed from the two fields. The distribution is however restrictedby the maximum orbital angular momentum quantum number that can be estimated by classicalconsiderations.

Chapter 4

Transitions between Rydberg states in two-color corotating circularly polarizedlaser pulses

4.1

Introduction1
In the previous chapter, we have discussed the role excited states play in the ionization pro-

cess for co- and counterrotating bichromatic circularly polarized laser pulse. The relative helicityof the two applied circularly polarized laser pulses has a significant effect on the manifold of accessible excited states during the interaction. Initially, it has been shown [178] that opposite photon
polarization in counterrotating fields increases the probability of resonant-enhanced ionization ascompared to the case of corotating fields. In Chapter 3 we demonstrated how the quantum selection rules for (multi-)photon absorption results in the selective population of excited states inatoms. Surprisingly, the results of the numerical simulations of the time-dependent Schrödinger
equation in the case of corotating fields however also revealed population in states with orbital
angular momentum and magnetic quantum numbers that are not accessible via sole absorption ofphotons from the two fields. We proposed that these states are populated via Raman-like transitions involving the absorption and emission of photons from both laser fields [174]. Let us assume
without loss of generality that the two corotating fields have central frequencies ω and 2ω and areboth left-handed circularly polarized, giving rise to selection rules of ∆m = −1 for the absorptionand ∆m = +1 for the emission of a photon from either one of these fields. As shown in Fig. 4.1(right column) the absorption of one photon at frequency 2ω (blue arrows) along with the emission
1

Part of the material presented in this Chapter has been published in Y. Gebre, J. Venzke, A. Jaroń-Becker, and
A. Becker, Physical Review A 103, 013101 (2021).

60
of two photons at ω (red arrows) in either a Λ-, V- or S-transition results in a total change ofmagnetic quantum of ∆m = +1. Analogous, photon emission at 2ω along with photon absorptionat ω gives rise to ∆m = −1 (Fig. 4.1, left column). Consequently, Λ-, V- and S-transitions can leadto a redistribution of population in Rydberg states, also beyond the manifold of states accessible
by absorption only.

Figure 4.1: Schemes for Λ− (upper row), V- (middle row) and S-transitions (lower row) between
Rydberg atomic states at central frequencies ω (red arrows) and 2ω (blue arrows) changing themagnetic quantum number by ∆m = −1 (left) or ∆m = +1 (right). The two applied corotating
circularly polarized fields are assumed to be left-handed polarized (taken from [179]).

Changing and controlling the population in atomic energy levels via Λ- and V-transitions are
well-known quantum mechanical processes which involve either one or two optical fields. Exem-61
plary applications based on this phenomenon are coherent population trapping [180,181], inversionfree lasing [182], stimulated Raman adiabatic passage [183, 184], electromagnetic induced transparency [185, 186], interference stabilization [187, 188] and population trapping [189, 190]. In allthese applications the absorption and emission steps each involve one photon, i.e. in total a secondorder process. In contrast the redistribution in the Rydberg states by corotating bicircular (ω, 2ω)laser fields relies on third-order processes requiring three photons, two at the longer and one at the
shorter wavelength, leading to the transitions shown in Fig. 4.1.
The results in the previous Chapter indicated a redistribution during the population of Rydberg states by bicircular laser pulses, which can be controlled via the relative intensity of thetwo fields. In this Chapter we extend the analysis of the redistribution schemes. For our studies
we make again use of results of numerical solutions of the time-dependent Schrödinger equation(TDSE) for the interaction of the hydrogen atom with intense bichromatic circularly polarized laser
pulses. This Chapter is organized as follows: We establish the effects of the simultaneous presenceof both fields on the excitation of atoms in corotating bicircular pulses at frequencies ω and 2ω
in section 4.2.1. By selecting the initial state of the atom in section 4.2.2 we analyze which kindof transitions are effective in the redistribution of population in the Rydberg states. In section
4.2.3 we provide indications that redistribution via even higher-order transitions occur in corotating circularly polarized fields with larger difference of the central frequencies. We end with a briefsummary in section 4.3.

4.2

Results and Discussion
The interaction with the bicircular laser pulse is implemented, as in the previous Chapter,via the total vector potential as:
A(t) = Aω (t) + A2ω (t)

(4.1)

62
where

2



AΩ (t) = A0,Ω sin

πt
τΩ



× [sin (Ωt) x̂ + ϵΩ cos (Ωt)ŷ]

(4.2)

√

Ω
for Ω = ω and 2ω, respectively. A0,Ω = c ΩIΩ , τΩ = 2πNΩ , and c is the speed of light, where IΩ

is the peak intensity and NΩ denotes the number of cycles. ϵΩ = ±1 denotes the helicity of thefundamental and 2nd harmonic pulse, respectively.

Figure 4.2: Comparison of excited state distributions as a function of magnetic quantum numberm (summed over n and ℓ with n ≥ 4) for excitation with (a) circularly polarized laser pulse at 267
nm (20 cycle, 5 × 1013 W/cm2 ), (b) circularly polarized laser pulse at 534 nm (10 cycle, 1 × 1014
W/cm2 ), and (c) corotating bicircular laser pulse (taken from [179]).

63
4.2.1

Population distribution in time-delayed bicircular corotating pulses

In the interaction with circularly polarized laser pulses the selection rules for the orbital
angular momentum and magnetic quantum numbers limit the manifold of excited states accessible
by absorption of photons [174, 178]. In multiphoton absorption in each transition the magnetic
quantum number is changed either by ∆m = +1 (right-handed circular polarization) or ∆m = −1
(left-handed circular polarization). In Fig. 4.2 we compare the distribution in the excited statesof the hydrogen atom as a function of the magnetic quantum number (summed over n and ℓ withn ≥ 4) by single circularly polarized pulses at 267 nm (a) and 534 nm (b) with that due to theinteraction with the superposition of the two pulses (c). The pulse duration of the two pulses werechosen to be the same and in the superposition the maxima of the field were chosen to coincide.
The distributions induced by the individual pulses are limited to narrow ranges in m, inagreement with the excitation via absorption of three photons at 267 nm and six or seven photons
at 534 nm. The total photon energies at central frequency are 13.92 eV (3 photons at 267 nm
and 6 photons at 534 nm) and 16.24 eV (7 photons at 534 nm). Assuming that the Rydberg
states approximatively shift with the ponderomotive energy, which is 0.333 eV (267 nm, 5 × 1013
W/cm2 ) and 2.663 eV (534 nm, 1 × 1014 W/cm2 ) at peak intensity, respectively, the maxima of
the population distributions are in agreement with an excitation near the center of the pulse.
In contrast, the distribution is much broader for the interaction with the bicircular laser pulse.
The states with magnetic quantum numbers between m = −3 and m = −7 are accessible via the
combined absorption of photons from both fields. However, the results also reveal populationin states with m > −3, which cannot be populated just via photon absorption alone. We havepreviously proposed [174] that a redistribution of population between Rydberg states with different
magnetic quantum numbers via Raman-type transitions involving three photons (c.f., Fig. 4.1) are
effective when both pulses are present.

64

Figure 4.3: Excitation probability in states with magnetic quantum number m (summed over n
and ℓ with n ≥ 4) for two corotating circularly polarized pulses at 267 nm (20 cycle, 5 × 1013
W/cm2 ) and 534 nm (10 cycle, 1 × 1014 W/cm2 ) as a function of time delay between the pulses.
At negative (positive) time delays the 267 nm (534 nm) pulse precedes (taken from [179]).

The interpretation that the population in the states outside of m = −3 are only populated
when both pulses overlap in time is confirmed by the results in Fig. 4.3, which show the population
distributions as a function of magnetic quantum number and the time delay between the two pulses.
The pulse parameters were kept the same as for the simulations in Fig. 4.2(c). When the pulses do
not overlap, i.e. for delays |∆t| > 11 fs, the population is concentrated in m = −3 in agreementwith an excitation by the 267 nm laser pulse without impact of the preceding (positive time delays)
or subsequent (negative time delays) 534 nm laser pulse. In contrast, once the pulses overlap intime there appears population in states with m-values larger and smaller than -3. Furthermore,
the range in m over which the population extends as well as the magnitude of population outsideof m = −3 itself is as larger as more the pulses overlap.

65

Figure 4.4: Excitation probability as function of magnetic quantum number (summed over n andℓ with n ≥ 4, upper row) and as function of principal and orbital angular momentum quantum
numbers (summed over m, lower row) for two corotating circularly polarized pulses at 267 nm (20cycle, 5 × 1013 W/cm2 ) and 534 nm (10 cycle, 1 × 1014 W/cm2 ). Results are obtained for initialstates prepared in n0 = 7 and ℓ0 = 3, m0 = −3 (a, e), ℓ0 = 4, m0 = −4 (b, f), ℓ0 = 5, m0 = −5 (c,g), and ℓ0 = 6, m0 = −6 (d, h) (taken from [179]).

4.2.2

Transitions between Rydberg states

The results so far indicate that the simultaneous presence of both pulses is required for the
population of excited states with magnetic quantum numbers other than m = −3 (and m = −6, −7).
Next, we analyze if any of the transitions for the redistribution between Rydberg states, shownin Fig. 4.1, are effective mechanisms in the bicircular pulse. To this end, we consider that thehydrogen atom is initially prepared in a specific Rydberg state. While the processes can proceed
via virtual states the transition probabilities are larger when intermediate real states are involved.
For two of the transitions (V and S) there is first an emission of one or two photons at the smaller
photon energy and, hence, intermediate bound states at lower energies must be accessible to makethe transition effective. Furthermore, we expect that there remains some population in these lower
lying states at the end of the pulse.
For the V-mechanism (Fig. 4.1(b)) the transition to the lowest intermediate state is associatedwith an energy exchange of ∆E = −2ω ≈ −4.64 eV (for the bicircular pulse used in the present
analysis), and ∆m = +1 (for total change of ∆m = −1) or ∆m = +2 (for total change of∆m = +1). Assuming that the initial Rydberg state shifts with the ponderomotive energy during

66
the pulse while the lowest lying intermediate states may shift less in energy, the states at n = 2
level in the hydrogen atom must be accessible for the V-transition to be effective. Similarly,
the S-transition to the lower intermediate state is associated with ∆E = −ω ≈ −2.32 eV and∆ℓ = ∆m = +1 (for total change in m of +1 or -1). Therefore, it is most likely that the transition
occurs via states at the n = 3 level. Due to the restriction in orbital angular momentum (ℓ < n )and magnetic quantum numbers (|m| < n) the accessibility of the intermediate lower lying states
depends on the quantum numbers of the initially prepared Rydberg state, which allows us to test
the presence of the redistribution mechanisms, proposed in Fig. 4.1.
To this end we have performed a series of calculations in which initially the hydrogen atom
is prepared in specific Rydberg states, namely (n0 = 7, ℓ0 = 3, m0 = −3) to (n0 = 7, ℓ0 = 6,
m0 = −6) changing ℓ0 by 1 and m0 by -1. The distributions as function of magnetic quantum
numbers and as function of principal and orbital angular momentum quantum numbers at the endof the interaction with the same bicircular pulse as in Fig. 4.2(c) are presented in Fig. 4.4. The
population in states with magnetic quantum numbers larger than that of the initial state clearly
depends on the choice of the initial state. While for ℓ0 = 3, m0 = −3 there is a redistribution overa range of states up to m = 1, the range gets smaller as ℓ0 increases (smaller m0 ) of the initialstate. The ranges of the distributions indicate that for most of the cases there is a sequence oftransitions occurring. In general, the type of transition (Λ, V, S) may change from step to step inthe sequence.
For the initial state ℓ0 = 3, m0 = −3 (Fig. 4.4 a, e) the total change in m is four, indicating
a sequence of up to four transitions. Any of the final states, except those with m = 1, can beaccessed by either one of the three transitions shown in Fig. 4.1. For the V-transition the process
likely proceeds via n = 2-states with ℓ = −1, 0, 1 and this transition is therefore not resonant for
the final step into the m = 1-states, since the required intermediate state with m = 2 does notexist in the n = 2-manifold of states. The other two pathways (Λ and S) are both allowed foreach transition in the sequence. Another indication that the V- and S-pathway are effective in theredistribution are the remaining populations in the intermediate (n = 2, ℓ = 1) and (n = 3, ℓ = 2)

67
states of these pathways at the end of the pulse.
For the next initial state (ℓ0 = 4, m0 = −4; Fig. 4.4 b, f) the range of populated m-statesindicates again a sequence of up to four transitions. In this case it is however unlikely that thefirst transition in the sequence to m = −3 occurs via the V- or S-pathways since the required
intermediate states with m = −3 (in the n = 3 manifold) and m = −2-states (in the n = 2
manifold) do not exist. This agrees with the reduced overall probability of redistribution over the
Rydberg states with larger magnetic quantum number than the initial state as well as the lower
population in the intermediate states (n = 2, 3) required for the V- and S-transitions in the sequence
beyond the first step.
The trend of reduced probability in the redistribution to states with larger magnetic quantum
numbers continues as the quantum numbers of the initial state are changed to larger absolute values
in ℓ0 and m0 . In the final example considered (ℓ0 = 6, m0 = −6; Fig. 4.4 d, h) the transfer of
population to the m = −5 state can occur via the Λ-pathway only. This pathway however appears tobe effective for just one step in the change of magnetic quantum number with rather low probability
only. This may indicate that the V- and S-transitions via bound states with lower energy are overall
more probable than the Λ-transition via the continuum states in the redistribution of population
between Rydberg states.
The last part of our interpretation about the relative efficiency of the different types oftransitions agrees with the observations for the redistribution into states with smaller magnetic
quantum number than that of the initial state. Except for the (ℓ0 = 3, m0 = −3) initial state,
the V- and S-transitions are not open for a change of ∆m = −1 since the required intermediate
states in the n = 2 and n = 3 manifolds do not exist. Due to the restriction in pathways weobserve population in just one (or two) states with smaller magnetic quantum number. Since theΛ-transition is the only ∆m = −1 pathway for all initial states with ℓ0 ≥ 4, the results confirm
that the probability for this transition is small and independent of the values of the orbital angularmomentum and magnetic quantum numbers.
There are a few more interesting features in the (n, ℓ)-distributions in Fig. 4.4 (lower row),

68
which we will discuss now. These features are not closely related to the main focus of this Chapter,
namely the three-photon Λ-, V- and S-transitions involving both fields. First, we note that the
final distributions in the ℓ0 -channel are spread over all available principal quantum numbers n. We
interpret this result primarily as an indication of the presence of one-pulse Λ-transitions throughthe continuum, consisting of the absorption of a photon followed by the emission of one photon fromone of the two fields. In these transitions the changes in orbital angular and magnetic quantum
number are given by ∆ℓ = 0, ±2 and ∆m = 0. Therefore, the population in the ℓ0 + 2-channels
in some of the final distributions is another indicator for the presence of these kind of transitions.
This process was previously analyzed in the context of interference stabilization [187, 188] andpopulation trapping [189,190] in strong fields. In the final distributions we further note populationin the ground state. Each of the initial Rydberg states considered in this set of calculations is
within the manifold of states that can be reached from the ground state via sole absorption ofphotons from the two fields. Therefore, it is likely that the ground state population is due to theinverse process, namely the deexcitation from the initial Rydberg state via the stimulated emissionof photons .

4.2.3

Higher-order transitions in (ω, pω) corotating circularly polarized pulses

The population of Rydberg states with magnetic quantum numbers in corotating bicircular(ω, 2ω) fields that are not accessible via absorption of photons from the initial states raises thequestion if similar distributions occur in corotating circularly polarized (ω, pω) fields with p > 2
as well. To answer this question we have performed additional calculations for the superpositionof one pulse at 267 nm and a second pulse at 800 nm and 1068 nm, respectively. In each case thepulse durations of the two pulses were chosen to be same and the center of the pulses to coincide.
The pulses at the longer wavelength had the peak intensity of 1 × 1014 W/cm2 , while the peakintensity of the 267 nm was 5 × 1013 W/cm2 , as in the (ω, 2ω) studies above. In test calculations
we have found that sole application of each of the pulses at the long wavelengths does not produce
any significant population in the excited states.

69

Figure 4.5: Comparison of normalized excitation probabilities in states with magnetic quantum
number m (summed over n and ℓ with n ≥ 4) for two corotating circularly polarized pulses at 267
nm (30 cycle, 5 × 1013 W/cm2 ) and 534 nm (15 cycle, 1 × 1014 W/cm2 ), 800 nm (10 cycle, 1 × 1014
W/cm2 ), and 1068 nm (7.5 cycle, 1 × 1014 W/cm2 ) (taken from [179]).

Figure 4.6: Excitation probability as function of magnetic quantum number (summed over n
and ℓ with n ≥ 4) and as function of principal and orbital angular momentum quantum numbers(summed over m, lower row) for two corotating circularly polarized pulses at 267 nm (30 cycle,
5 × 1013 W/cm2 ) and 800 nm (10 cycle, 1 × 1014 W/cm2 ). Results are obtained for initial states
prepared in n0 = 10 and ℓ0 = 3, m0 = −3 (a), ℓ0 = 5, m0 = −5 (b), ℓ0 = 7, m0 = −7 (c), andℓ0 = 9, m0 = −9 (d) (taken from [179]).

70
From the comparison of the results for the normalized excitation probabilities as function of
the magnetic quantum number at the end of the pulse with that of the (ω, 2ω) case in Fig. 4.5 it isseen that the range of states as function of magnetic quantum number and the relative redistribution
probability increases with change from p = 2 to p = 4. Thus, higher-order transitions involving
the emission and absorption of more than three photons are likely effective in the excitation ofhydrogen atom with (ω, pω) fields. The number of possible pathways increases with increase of pand changes in the magnetic quantum number up to ∆m = ±(p − 1) are possible in the individualtransitions.
As in the case of the (ω, 2ω)-field, we have performed additional series of calculations by
preparing the hydrogen atom in specific Rydberg states. The final distributions over the magnetic
quantum numbers for the interaction with the superposition of corotating left-handed circularly
polarized fields at 267 nm and 800 nm are shown in Fig. 4.6. The distributions are obtained for
initial states ranging from (n0 = 10, ℓ0 = 3, m0 = −3) to (n0 = 10, ℓ0 = 9, m0 = −9) changingℓ0 by 2 and m0 by -2. The distributions show the same trends as those in Fig. 4.4. The rangeof populated magnetic quantum states decreases and the overall redistribution probability overthe Rydberg states decreases as ℓ0 increases (smaller m0 ). These features indicate that a similarinterpretation holds as for the (ω, 2ω)-field. As ℓ0 increases, transitions that involve photons fromboth fields and proceed via lower lying states are stepwise excluded, which provides an explanation
for the observed features.

4.3

Summary
Using the solutions of the time-dependent Schrödinger equation we have further analyzed thatthe simultaneous interaction with both fields leads to a distribution in states over a broad rangeof magnetic quantum numbers, in contrast to narrow distributions if only one of the two fields arepresent. Since the range extends beyond that accessible via absorption of photons from the two
fields, we have proposed that higher-order Raman (Λ, V, S) transitions, involving absorption andemission of at least three photons from the two fields, are effective mechanisms in the redistribution

71
of population between Rydberg states.
By selecting specific Rydberg states in the numerical calculations it has been found that the
distributions over the orbital angular momentum and magnetic quantum numbers depend on theinitial state. The changes in the distributions have been interpreted as indications of the presenceof the higher-order Raman transitions. Specifically, the elimination of the V- and S-pathways via
lower-lying states leads to a significant decrease in the overall redistribution probability and therange of populated ℓ- and m-states in the final distribution.
While the presence of the proposed transitions is supported by the present results, we notethat the overall excitation and redistribution probabilities are rather low. This is related to thefact that the competing process of ionization, i.e. transitions to the continuum, has a much higher
probability. Indeed, in all calculations presented in this work the ionization probability is typicallyby a couple of orders of magnitude larger than the total redistributed excitation probability. It is
likely that the ratio between excitation and ionization can be controlled to some degree via thevarious laser parameters such as the intensity ratio (as shown in the previous Chapter), the relative
wavelengths, the relative carrier-to-envelope phase, and the pulse durations of the two fields as wellas the time delay in between them. Although the efficiency of the control may be limited in view of
the number of competing pathways (ionization, direct excitation, Λ−, V- and S-transitions), we notethat the impact of the carrier-to-envelope phase on bound-bound transitions [191,192] and coherent
control protocols to prepare specific Rydberg states [193, 194] have been studied recently. To thisend, it can be useful to further analyze qualitatively and quantitatively each of the transitions
separately using other theoretical approaches, e.g. via lowest order perturbation theory.
Finally, we have shown that redistribution mechanisms are also present in (ω, 3ω)- and(ω, 4ω)-fields. The range of populated m-states increases as larger is the difference in the central
frequencies of the two fields. This is in agreement with the increase in change of magnetic quantum
numbers in higher-order Raman transitions in (ω, pω)-fields with p > 2. Results of calculations
starting from specific initial Rydberg states support the conclusions drawn for the (ω, 2ω)-case.

Chapter 5

Photoelectron emission spectra of atoms in intense short circularly andelliptically polarized laser pulses

5.1

Introduction
The spin polarization of photoelectrons serves as a valuable tool for understanding the struc-ture of atoms and molecules. By utilizing intense laser pulses, photoelectrons can be generated
through the interaction with atoms, molecules, or solids. To achieve spin-polarized electrons,
circularly polarized laser pulses are employed for electron emission. A crucial aspect in generating these electrons lies in the selective emission of electrons with a specific rotation of theirinitial state relative to the helicity of the circularly polarized laser pulse [148, 154]. The rotation of the electron’s initial state is determined by its magnetic quantum number. Numerous
studies have investigated the relationship between the helicity of the initial state and the ionization process, revealing helicity-dependent enhancements of ionization across a wide range of
wavelengths [30, 108, 148, 154, 163, 164, 195, 196]. This study focuses on the realm of multi-photon
ionization, where initial states with the opposite helicity to the applied field are found to be preferentially ionized. A mechanism underlying this phenomenon, that has been recently proposed byour group, is the existence of additional pathways for ionization specific to electrons originatingfrom these states [163,164]. Building upon this previous work, the objective of the present researchis to extend the investigation beyond ionization rates and explore other observables associated withthe relative helicity of the initial state.
It is important to clarify the terminology we will be using in this study. Here, we will refer to73electrons with initial states matching the helicity of the pulse as ”corotating,” while those with the
opposite helicity will be referred to as ”counterrotating.” This distinction differs from the previous
sections, where we utilized the terms ”corotating” and ”counterrotating” to describe the relative
helicity of the two circularly polarized pulses employed in the TDSE calculations for those Chapters.

5.2

Absorption pathways and photoelectron spectra
l=0 l=1

l=2

l=3

l=4

l=m+2
l=m
l=m

m = −1 | +1
Figure 5.1: Photon absorption pathways for a right-handed circularly polarized pulse with initialstates of m = −1 (red) and m = 1 (blue). Direct path for corotating electrons indicated by the
red line, direct path for counterrotating electrons indicated by the blue line, and off-set path forcounterrotating electrons marked by the light-blue line.

As described in the introduction, the enhanced ionization observed for the counterrotating
initial state can be attributed to the availability of additional absorption pathways for that state.
This concept is depicted in Fig. 5.1, which illustrates the absorption pathways for a right-handed
field. This can be explained using quantum selection rules (see also section 2.2.4) which dictatethat for a right-handed circularly polarized pulse, each photon absorption causes a change in theorbital angular quantum number (∆l = ±1) and the magnetic quantum number (∆m = +1).
Consequently, the available states for the electron to occupy as it transitions to an excited state

74
are constrained. Considering the counterrotating initial state (mo = −1, lo = 1), the absorptionof one photon leads to m = 0, which, with ∆l = ±1, can transition to either the l = 0 or
l = 2 states. Conversely, for the corotating state (mo = −1, lo = 1), a transition occurs to them = 2 state, and with the restriction on l, it is only possible to further transition to the l = 2 state.
Thus, the counterrotating state benefits from additional transition pathways during each successive
absorption, resulting in ionization enhancement.
The relative helicity of the field and the initial state also influences another observable:
the photoelectron spectrum. As demonstrated in a study reported in [197], it has been observedthat the shape of the spectrum differs for corotating and counterrotating initial states. In thisinvestigation, we aim to determine if the photoelectron spectrum, in addition to ionization rates,
can serve as a tool to examine the hypothesis of additional pathways for counterrotating initialstates. To facilitate this analysis, we refer back to the diagram in Fig. 5.1 and introduce labels to
distinguish the various pathways. For the corotating initial state, only one pathway is evident, as
illustrated by the red arrows in the right diagram of Fig. 5.1. Conversely, for the counterrotating
case, we observe a pathway denoted by blue arrows that pass through states where l = m + 2, as
well as another pathway marked by light-blue arrows where l = m. We will refer to the pathways
marked in blue and red as ”direct” pathways, while the pathway indicated by light-blue arrows forthe counterrotating case will be labeled as the ”off-set” pathway. The ”off-set” pathway represents
an additional option available solely in the counterrotating case, and we hypothesize that it is
responsible for the observed disparity in the photoelectron spectrum. Furthermore, we proposethat if we examine the spectrum derived from the direct pathway for both the corotating andcounterrotating cases, they should exhibit the same distribution.

75

Figure 5.2: Comparison of photoelectron energy spectra for neon atom, initially in l = 1, m = 1(corotating) state, indicated by the red curve and atoms, initially in l = 1, m = −1 (counterrotating)state, indicated by the blue curve. The top panel represents the total spectrum, while the bottom
panel shows the spectrum originating from the direct paths only. The results were obtained using
a 10-cycle right-handed circularly polarized laser pulse with an intensity of 3 × 1014 W/cm2 .

To verify this hypothesis, we conducted simulations of the interaction between a neon atom
and a circularly polarized laser pulse, considering the initial states of l = 1, m = 1 and l = 1, m =−1. The pulse duration was set to 10 cycles, a wavelength of 400 nm was chosen. As described in
section 2.2.8, the photoelectron spectrum can be obtained from the time-propagated wave function
as
2

1 X
F (k, ϕ, θ) = 2
k

Z

−iδkl

e

l

∗

(i) ϕkl (r) ψ(r, t) dx



∗
Yl,m
(ϕ, θ)(5.1)

l,m∗ (ϕ, θ) are the spherical harmonics, δ is the phase shift, and ϕ (r) is a continuum state
where Yl,mklkl

with momentum k. To calculate the photoelectron spectrum, the summation is performed over the

76various orbital angular and magnetic components of the wave function. However, by restricting
this summation to only include the l and m values that correspond to the direct pathway, we canisolate and analyze the contribution of this specific pathway to the overall spectrum. Applying
this approach to our TDSE calculations, we obtained simulation results presented in Fig. 5.2. In
the top plot of the Figure, we show the spectra for corotating (red) and counterrotating (blue)
helicity of the initial state. As anticipated, we get different shapes for the spectra in the two cases.
In the bottom plot, we show as comparison and for further analysis the portion of the spectra
originating from the direct pathway for both the corotating and counterrotating cases. Indeed, we
observe a high degree of similarity in the spectra for the direct pathway. This finding is in line with
our hypothesis and provides support for the idea that the discrepancy in the spectra between thetwo cases can be attributed to the presence of an additional pathway exclusively available to thecounterrotating case.

Figure 5.3: Same as Fig. 5.2 but for an 800 nm pulse.

77
In order to test that this explanation for the spectra holds for other laser wavelengths as well,
we performed another calculation for a circularly polarized laser pulse with a wavelength of 800
nm. The resulting spectra are shown in Fig. 5.3. We indeed make the same observations from theresults. For the total spectra, the shapes differ for the co and counterrotating initial state (panel
at top), while when we consider the contribution from the direct pathway we get similar shapes forboth spectra (panel at bottom).

Figure 5.4: Distribution of population in the orbital angular and magnetic states of neon atom
after interaction with a 400 nm circularly polarized laser pulse. Laser parameters are the same as
in Fig. 5.2. Left for initial counterrotating state with l = 1, m = −1 and right for corotating statewith l = 1, m = 1.

In addition to analyzing the observable spectra, we can gain more insight by looking at the population in the different l and m components of the time propagated wave function. In Fig. 5.4 we plot
the components for the counter (left) and corotating (right) initial states for the TDSE simulation
for the 400 nm laser pulse considered previously. We plot the components related to the direct
pathway in red and those for the off-set one in blue. We see that the relative amount of populationin the different components of l are similar for the direct pathway, while those for the off-set one
have a different distribution. This reinforces our conclusion that the pathway explanation of theionization process is a valid interpretation and the discrepancy in the observables of the co andcounterrotating cases can be attributed to the additional pathways available in the counterrotating

78
case. We also note that, for the present laser parameters, the major contribution to the spectrum
comes from the l = m components for both cases. This is due to the strong coupling between thesestates.

5.3

Emission angles of photoelectrons
Additionally to the energy distribution, we have also examined the angular distribution of

the ionized electrons. When dealing with a circularly polarized pulse consisting of multiple cycles,
we observe uniform emission across all angles within the plane of polarization. To investigate thedisparity between corotating and counterrotating scenarios, we therefore introduce some ellipticity
to the pulse. The dissimilarity in angular resolved spectra between the two helicities has been
explored in a previous study reported in Ref. [197]. Our focus in this research is a phenomenon
that has gained attention in recent years, namely the variation in the emission angle of the spectra
as a function of energy [67, 69, 71, 72]. To address this, we conducted numerical simulations using
our ab-initio method, which involved the interaction of a neon atom with an elliptically polarizedlaser pulse. Similarly to our previous calculations, we performed two sets of calculations, one for an
initial (counterrotating) state of l = 1, m = −1 and another for a (corotating) state of l = 1, m = 1
for a 400 nm pulse with an ellipticity of 0.7. The results are illustrated in Fig. 5.5, where theresults presented in the top panel correspond to the l = 1, m = −1 case, and the results at thebottom represent the l = 1, m = 1 case. The plots on the left reveal that the angle under which the
maximum of the distribution for each ATI peak occurs changes as the energy varies. To demonstrate
this further, the plots on the right depict the spectra for each ATI peak as a function of angle.
We can observe that the emission angle for the counterrotating case changes more significantlywith each successive peak as compared to the corotating case. Previous studies have indicatedthat the origin of this rotation is associated with non-adiabatic tunneling, where electrons acquire
additional momentum in the parallel component during tunneling due to the interaction with thelaser field [67, 72].

79

Figure 5.5: Photoelectron spectrum for a neon atom interacting with an elliptically polarized laser
pulse. The top row corresponds to results for the initial (counterrotating) state with l = 1, m = −1,
and the bottom row corresponds to those for the initial (corotating) state with l = 1, m = 1. The
left column displays the angle-resolved emission, the middle column shows the emission as a functionof angle for each ATI peak, and the right column represents the spectrum generated with the threehighest contributing components. The pulse parameters are as follows: wavelength of 400 nm,
intensity of 3 × 1014 W/cm², and ellipticity of 0.7.

In addition, we present the results of similar TDSE simulations with ellipticities of 0.8 and
0.9, as shown in Fig. 5.6. This Figure illustrates the variation in the angle of maximum emission
for each ATI peak. Consistent with our previous findings, we observe a change in the emission
angle as a function of energy. Furthermore, it is evident that the counterrotating scenario exhibits
a greater degree of change in the emission angle compared to the corotating case, independent of
the ellipticity of the field. These results emphasize that the observations discussed in this section
are not specific to a particular ellipticity. Within this section, our aim is to provide further insight
into the distinction between the corotating and counterrotating scenarios regarding the extent ofvariation in emission angle for successive ATI peaks. Our particular focus lies in exploring thecorrelation between the emission angle and the different orbital and magnetic state componentsthat contribute to the photoelectron spectrum.

80

Figure 5.6: The angle corresponding to the maximum emission of electrons for neon atom interactingwith an elliptically polarized pulse. The laser parameters are the same as in Fig. 5.5 for ellipticities
of 0.8 (blue) and 0.9 (red).

To focus our analysis, we first identify the orbital angular and magnetic states that provide the mostsignificant contributions to the observed spectra. In both the corotating and counterrotating cases,
we identify the three states that contribute most to each ATI peak. We exclusively plot the spectra
produced by considering only these states in the column on the right of Fig. 5.5. Remarkably,
we observe a highly accurate reproduction of the spectral features that we see in the full spectra(panels in the middle). Upon investigation, we find that the three most contributing components
for the nth ATI peak can be represented as (ln , mn = ln ), (ln − 2, mn − 2), and (ln − 4, mn − 4).
Where (ln , mn = ln ) is the state occupied by absorbing the minimum number of photons to ionize
an electron with energy equal to the nth ATI peak. For the laser parameters used in our TDSE
calculation ln = n + 7. Note that in the case of the counterrotating state, this channel is related tothe additional or off-set pathways.
This identification simplifies our analysis significantly, as the spherical harmonics associatedwith these states in the plane of polarization (θ = π2 ) can be expressed as Yl,m=l ∝ eimϕ . By
considering the interference of these three dominant states, we attempt to calculate the observed

81
spectrum. The spectrum for the nth peak can be described as follows:Fn ∝ eiσ1 eimn ϕ + eiσ2 ei(mn −2)ϕ + eiσ3 ei(mn −4)ϕ

2

(5.2)

where eiσ0 is the phase associated with the projection of the wave function on to the continuum
states as shown in Eq. (5.1). After some simplification the spectrum produced from the interference
will be given by

Fn ∝ 3 + 2 cos(σ12 + 2iθ) + +2 cos(σ23 + 2iθ) + 2 cos(σ13 + 4iθ)

(5.3)

where σij = σi − σj is the phase difference. The expression presented in Eq. (5.3) implies that thephase difference between the three states plays a crucial role. To test this model, we set σ1 = 0
and modify the relative phase between the states in Eq. (5.3) to observe its impact on the change
in the angle of maximum emission. We plot the results of this model in Fig. 5.7. In the top
plot, we depict the phases (σi , i = 1, 2, 3) for the three contributing states, while the bottom plot
showcases the resulting spectrum obtained from Eq. (5.3). It is evident that as the phase difference
between the states increases, the amount of variation in the emission angle as a function of energy
also increases. This finding strongly suggests that the emission angle is determined by the relative
phase of the different components, and in this particular case, it can be explained by considering
only the three most contributing states.

82

Figure 5.7: Illustration of the dependence of angle of electron emission based on the model given
by Eq. (5.3). Top panel shows the phase of the three contributing components while the bottom
panel shows the resulting spectra.

If our assumption is accurate, we would anticipate that the relative phase between the threecomponents obtained from our TDSE calculation exhibit the same behavior as described earlier.
Specifically, for the counterrotating case, we would expect a larger phase difference between thestates compared to the corotating case, in order to account for the observed angular spectrum. To
verify the validity of our assumption, we plot the phases of the three most influential components
in Fig. 5.8. From the result we observe the expected behavior, thereby validating our assumption.
To ensure the applicability of the simplified model across various ellipticities, we plot the phases
for ellipticities of 0.7 and 0.8 in Fig. 5.8. This analysis confirms that the simplified model remains
effective for different values of ellipticity.

83

Figure 5.8: Phase of the three biggest contributing components of the photoelectron spectrumof neon interacting with an elliptically polarized pulse. The plots on the left correspond to the
corotating case, while the plots on the right correspond to the counterrotating case. The top row
represents an ellipticity of 0.7, and the bottom row represents an ellipticity of 0.8.

In summary, this section we focused on analyzing the photoelectron energy spectra of atoms
subjected to intense short circularly and elliptically polarized laser pulses. For the circularly polarized case, our investigation centered on explaining the differences in the photoelectron spectrum
between initial magnetic states with the same and opposite helicity as the laser pulse. Our findings
revealed that the dissimilarity in the spectra originated from the availability of additional absorption pathways specific to the counterrotating initial state. By plotting the spectrum resulting fromthe pathway accessible to both cases, we obtained a similar spectrum, indicating the influence of
these additional pathways. To explore the variations in angular emission between the corotating
and counterrotating cases, we introduced ellipticity to the laser pulse, resulting in angular structure.
Notably, we observed that the emission angle changed in relation to the energy of the photoelectron.
To provide a comprehensive explanation of the observed spectra and emission angles, we developed
a simplified model that incorporated the dominant orbital angular and magnetic states. We foundthat the phase difference between these states was responsible for the observed spectrum, with the84
counterrotating case exhibiting a larger phase difference, leading to a more significant change in theemission angle. Overall, our study sheds light on the intricate dynamics reflected in photoelectron
spectra in intense laser pulses and highlights the role of initial magnetic states and their phasedifferences. These insights contribute to our understanding of laser-atom interactions and havepotential implications in various fields, such as laser physics and spectroscopy.

Chapter 6

Analysis of photoelectron energy spectra using Quantum Trajectory Monte
Carlo simulations

In section 2.3, we explained how the photoelectron spectrum can be calculated using quantum
trajectory Monte Carlo (QTMC) simulations. The method involves a two-step process to describe
the ionization processes. First, the probability of an electron tunneling through the atomic potential
barrier is modeled using ionization rates predicted in the quantum mechanical framework. Second,
the motion of the ionized electron in the combined potential of the parent ion and laser field is
modeled using Newtonian mechanics. If the simulation includes the phase accumulated during the
trajectory, it is called a ”quantum” or ”semi-classical” simulation. In this work, we refer to themas quantum trajectories. To model the tunneling, we utilize the PPT ionization rates described in
section 2.3.2, which results in a non-adiabatic model of ionization [23, 24, 64]. The non-adiabatic
models incorporate the electron’s oscillation while being subjected to the field’s influence during
the tunneling process. For more information on the adiabatic and non-adiabatic models, pleaserefer to sections 2.3.1 and 2.3.2, respectively. In this section, we present some results from QTMC
simulations and compare them to results from our numerical TDSE solutions described in section
2.2.

6.1

QTMC with PPT ionization rate
The SFA approximation is usually the starting point for modeling the tunneling process forstrong-field ionization. In this approximation the transition from an initial bound state into the86continuum is given by
MSF A (⃗k) = −iZ ∞−∞dt ψ⃗k |VL (t)| ψi (t) ,(6.1)

where |ψi (t)⟩ is the initial bound state, VL is the laser-atom interaction in length gauge, and ψ⃗k is
the Volkov state. As describe in section 1.3.3, the integral can be simplified using the saddle point
approximation. Following the work in Refs. [69, 198], the transition probability is given by

MSF A (⃗k) =XV⃗j exp [iS(v(tj ), tj )]
k

j

∂2 ⃗
S(k, tj ) .
∂t2

(6.2)

In most applications of this formula in QTMC simulations, the dependence on the initial state andderivative of the action are ignored [65,67,69–72]. Instead the exponential term is used to calculate
the ionization rate as
W (to , p⃗i ) = MSF A (⃗k)

2

∝ e−2iS(v(ti ),ti ) .

(6.3)

The utilization of the rate in Eq. (6.3) yields a reasonable agreement with other numerical solutionsof the TDSE [67, 71]. However, as discussed in the previous Chapter, certain effects rely on theinitial state of the ionized electrons and this is the topic we want to explore. In the context of
QTMC, our focus lies in discerning variations in the characteristics of the ionized electrons basedon their initial state, more specifically on the dependence on the magnetic quantum number. With
this goal in mind, we incorporate the full transition probability, given in Eq. (6.2), to calculate
the ionization rate. This ionization rate, referred to as the PPT formula [23, 24, 64] provides theionization rate ass
WP P T (to , p⃗i ) ∝X
iϕlm (v(ti )) exp [iS(v(ti ), ti )]

2

2π
′′
S (v(ti ), ti )

,

(6.4)

where p⃗i is the initial momentum, to is the ionization time, and v(t) = p + A(t). The ionization
rate is calculated as a sum over the saddle points (ti ), as described in section 1.3.3. The saddle
points are calculated from the saddle point equation (1.9). The detailed derivation of this initial87
state term can be found in the work conducted by Barth and Smirnova [30]. Here, we present the
final result and explain its implications on the QTMC simulation. The derivation assumes a laser
polarized in the x − y plane and can be readily adapted for other orientations by referring to themethodology outlined in Ref. [30]. The PPT ionization rate’s initial state term can be expressed
as follows:(l − |m|)! |m|ϕlm (v(ti )) ∝P
(l + |m|)! l

ik
p z
2Ip

!2

2

eimϕv (ti ) ,(6.5)

where Plm are associated Legendre polynomials and l and m are the orbital and magnetic quantum
numbers of the initial state. The ϕv (ti ) is the azimuthal angle associated with the vector v described above. It is important to note that the exponential term in Eq. (6.5) does not equal unity
for m ̸= 0. This is due to the complex values associated with the saddle points, which result incomplex arguments for the exponential function. The exponential term is what leads to different
ionization characteristics between states with positive and negative magnetic numbers. To examine
the influence of states with different magnetic numbers on the photoelectron spectrum, we specifically consider the initial states characterized by l = 1, m = 1, and l = 1, m = −1. The exponential
term can be expanded for m = ±1 as

eimϕv (ti )2= |cos ϕv (ti ) + i sgn(m) sin ϕv (ti )|2 .

(6.6)

With the exponential term given above, the PPT ionization rate is given bys
WP P T (to , p⃗i ) ∝X

iS(v(ti ),ti )

e

i

2

2π
S ′′ (v(ti ), ti )|ϕlm (v(ti ))|2
(6.7)

(l − |m|)! |m|P
ϕlm (v(ti )) =(l + |m|)! l

ik
p z
2Ip

!2
|cos ϕv (ti ) + i sgn(m) sin ϕv (ti )|2 .

The ionization rate obtained above will be employed to generate the distribution for the initialconditions utilized in the QTMC simulation presented in the subsequent section. The difference in

88
the ionization characteristics between the two considered states can be attributed to the sgn(m)term in Eq. (6.7). The influence of this term on the ionized electron spectrum has been extensively
studied in Ref. [30]. However, our focus lies on the QTMC approach to replicate the distinct
characteristics of the photoelectron spectrum and angle-resolved spectrum originating from ionizing
electrons from these two states.
Different methods are employed to calculate the photoelectron spectrum for the ab-initio
and QTMC approaches. In the ab-initio method, the wave function is computed at the end of thepulse, enabling the calculation of the photoelectron spectrum. This is achieved by projecting thetime-evolved wave function onto continuum states. In the ab-initio approach, the photoelectron
spectrum is determined as follows (see section 2.2.8):
2

1 X
F (k, ϕ, θ) = 2kZe−iδkl (i)l ϕkl (r)∗ ψ(r, t) dx

∗(ϕ, θ) .
Yl,m

(6.8)

l,m

Here, Y l,m (ϕ, θ) represents the spherical harmonics, δkl corresponds to the phase shift, and ϕkl (r)represents a continuum state with momentum k. The continuum states are determined either by
diagonalizing the field-free Hamiltonian or by implementing the shooting method, as described in
section 2.2.6. The photoelectron spectrum for the QTMC is calculated by binning the trajectories
according to their final momentum. We use the phase accumulated by the electrons along theirpath to incorporate interference effects. As described in section 2.3.3 the spectrum is then givenbyF (⃗k) =

Nk
X

2

e

iΦ(tjo ,⃗voj )

,

(6.9)

j=1

where Φ(tjo , ⃗voj ) is the phase of each trajectory described in section 2.3.3 and Nk is the total numberof trajectories in the k-th bin. By employing different computational methods to compute thespectrum, we can effectively evaluate the Quantum Trajectory Monte Carlo (QTMC) approachand compare the results with those obtained from the ab-initio method. In the following section,
we present the photoelectron spectrum calculations derived from both the ab-initio and the QTMC

89
methods.

6.2

Photoelectron spectrum
In this section, we compare the results for the photoelectron spectrum produced from the

QTMC simulations and the ab-initio TDSE solutions. We calculated the photoelectron spectrumof a neon atom starting in different initial states all subjected to an elliptically polarized pulse. The
reason for choosing an elliptically polarized pulse is explained in the next section. In Fig. 6.1, we
present the outcomes of ab-initio and QTMC calculations for neon atom interacting with a 400 nm
pulse. We see that there is a very good agreement between the QTMC and the ab-initio results forboth the m = 1 and m = −1 states. It is important to note that the laser parameters utilized inthis calculation place us in the non-adiabatic regime, with a Keldysh parameter γ of approximately
2 (for a definition and discussion of the Keldysh parameter, please refer to section 1.3.3). These
parameters were chosen to emphasize the ability of non-adiabatic QTMC simulations to produce
correct results in the regime where the ADK (adiabatic) formulas fail.

Figure 6.1: The photoelectron spectrum of a neon atom interacting with an elliptically polarized
pulse at 400 nm is shown in the Figure. The plot on the left corresponds to an initial (corotating)state with quantum numbers l = 1 and m = 1, while the plot on the right corresponds to an initial(counterrotating) state with quantum numbers l = 1 and m = −1. The laser intensity is 9 × 1014
W/cm2 , the pulse has an ellipticity of 0.7 and a duration of 10 cycles.

90
The performance of QTMC simulations improves when the Keldysh parameter decreases,
which can be achieved by increasing either the intensity or the wavelength of the laser pulse. In
the next example, we choose to increase the wavelength from 400 nm to 600 nm while keeping theintensity constant. Fig. 6.2 displays the photoelectron spectrum for Neon when interacting with
an elliptically polarized pulse of wavelength of 600 nm. We see that indeed the QTMC simulation
reproduces the photoelectron spectrum again well.

Figure 6.2: Same as Fig. 6.1 but for a laser pulse with a wavelength of 600 nm.

It is crucial to highlight that in order to compare the two methods accurately, the observed spectrum
needs to be normalized. This normalization is necessary because certain constant factors are omitted
from the PPT model. Our focus lies not in the total number of ionized electrons, but rather in the
distribution of ionized electrons based on their energy and angle.

6.3

Angular distribution of photoelectrons
In addition to examining the photoelectron spectrum, our analysis encompasses a comparisonof the angular distribution of ionized electrons. Throughout this section and the preceding one, we
employed elliptically polarized pulses to compare the results obtained from QTMC simulations andab-initio calculations. When using circularly polarized pulses, distinguishing the angle-resolved

91
spectrum between the m = 1 and m = −1 states becomes challenging. This difficulty arises
because both cases yield concentric circles within the plane of polarization. Although the relative
height of each circle differs, we prefer elliptically polarized pulses as they generate distinct angularemission characteristics, enabling easy differentiation of ionization from the two states. Employing
the identical parameters as in the previous section, we utilize an elliptically polarized pulse oriented
in the x − y plane. Fig. 6.3 compares the angle-resolved spectrum for 400 nm pulse for a neon
atom with initial states of l = 1, m = 1 and l = 1, m = −1, respectively. The angle ϕ is the polar
angle in the x − y plane. We can see that we have good agreement between the two methods.

Figure 6.3: Angle-resolved spectrum of a neon atom subject to interaction with an ellipticallypolarized pulse. The plots on the left-hand side display results from the QTMC simulation, while
the plots on the right-hand side depict the results derived from the ab-initio solution. The top
plots correspond to the l = 1, m = −1 initial state, whereas the bottom plots represent the l = 1,m = −1 initial state. The parameters are: Laser intensity of 3 × 1014 W/cm2 , ellipticity of 0.7, andpulse duration of 10 cycles.

As discussed in Chapter 5, the emission angle of the photoelectrons changes as a function of

92
energy [72]. The non-adiabatic models of QTMC simulations are shown to reproduce this effectobserved both experimentally and with ab-initio calculations. As discussed in the previous Chapter,
we see that the amount of rotation depends on the magnetic quantum number of the initial state.
The results from our calculations show that the QTMC simulations reproduce this effect well,something that previous models, that ignore the initial term, were not able to do. Thus, we seethat in addition to the energy spectrum, the angular resolved spectrum also accurately accounts
for the effect of the initial state.

Figure 6.4: Same as Fig. 6.3 but for 600 nm laser pulse.

As before we also include a calculation for a 600 nm pulse for comparison of the two methods
in the smaller Keldysh parameter regime. The 600 nm angle resolved spectrum is shown in Fig. 6.4.
A similar good agreement as before for the 400 nm spectrum is found. A crucial aspect in calculating
the angular resolved spectra is the requirement for a larger number of trajectories compared to thephotoelectron spectrum in order to resolve the peaks. The photoelectron spectrum benefits from93
integrating all trajectories within the same momentum range, whereas the angle resolved spectralack this advantage. To conduct a comprehensive three-dimensional comparison of the angularemission, an even larger number of trajectories would be necessary.

6.4

Advantages and limitations of QTMC simulations
QTMC simulations offer several advantages compared to the ab-initio method. First, one

notable advantage is the constant scaling of QTMC with respect to additional dimensions in the
problem. When solving Newton’s equations in three dimensions for the trajectories, the computational cost is less than three times that of solving it for a single dimension. On the other hand, the2ab-initio method scales as rpoints × lmax
, where lmax represents the highest degree of the sphericalharmonics used to describe the system and rpoints is the number of points used to represent theradial grid. In the case of longer wavelengths, lmax can reach values as high as 90 in the present
calculations, resulting in significantly increased memory and time requirements for computing solutions of the three-dimensional TDSE. Another benefit arises from the linear scaling attained throughparallelization in QTMC simulations. Due to the independence of all trajectories, parallelization
becomes straightforward and can be effectively implemented across separate compute nodes. In
contrast, the ab-initio solution can achieve near-linear scaling but is constrained to a single computenode, as communication between different cores is necessary for each iteration of the solution. It
is further important to highlight that the independence of the trajectories in QTMC simulations
enables the potential to improve unconverged results by incorporating additional data. On theother hand, in the case of the unconverged ab-initio solution, the result obtained is indescribableand cannot be progressively improved through iterative methods.
Like other Monte Carlo methods, QTMC simulations have inherent limitations due to theirstatistical nature. One significant limitation that affects the spectrum is the challenge of resolving
fine details and detecting relatively weak signals. In particular, fine details in the spectrum produced
by longer wavelength pulses are affected. This occurs because the ATI peaks in the produced
spectrum are separated by the energy of a single photon. When dealing with longer wavelengths

94(lower frequency), the photon energy becomes relatively small, resulting in less distinct separation
between the peaks and details of the spectrum. Consequently, in the above results for 400 nm and
600 nm, QTMC simulations were calculated with 4 × 108 and 8 × 108 trajectories, respectively,to achieve a similar level of agreement with the ab-initio method. For the results of the energy
spectrum presented above we can see that we have small peaks in the tail of the spectrum. In
the QTMC simulations, these small peaks have fewer trajectories ending with a final momentum
within that particular range. Consequently, to precisely discern the structure of these peaks, asignificant number of trajectories is required in order to sufficiently sample instances within that
range. The aforementioned limitations can be overcome by employing an adequate number oftrajectories. In this context, ”enough trajectories” refers to two aspects. First, it entails achieving
numerical convergence when our objective is to accurately depict the relative scale between peaks
or other features. Second, it involves attaining a sufficient quantity of trajectories to discern the
distinct features of interest.
The QTMC simulation can be a valuable complement to the ab-initio method when investigating ionization phenomena. The ab-initio method is particularly effective when dealing withscenarios where the parameter γ is much greater than 1. This parameter corresponds to laser pulses
with short wavelengths and relatively weaker intensities. In such cases, the ab-initio method offers
quick numerical convergence and requires modest memory resources. However, as we shift towards
longer wavelengths and higher intensity laser pulses, the ab-initio method becomes computationallydemanding. This is due to the necessity of employing extensive expansions in spherical harmonics
and large grids to represent the problem adequately. Conversely, the QTMC simulation performs
exceptionally well in the range where the ab-initio method becomes computationally burdensome,as its range of validity is in the regime of γ ⪅ 1. The scaling of the method is also linear with wavelength and thus has no limitation on modeling pulses with very long wavelengths. By combining
these two methods, we gain the capability to model atom-laser interactions across a wide spectrumof laser parameters, enabling a thorough analysis of ionization dynamics.

Chapter 7

Ab-initio solution for diatomic systems

Ab-initio solutions can be applied not only to atomic systems but also to diatomic systems
in order to solve the TDSE. In a previous Chapter (section 2.2), we described an ab-initio method
for solving the TDSE of atomic systems, and in Chapters 3, 4, 5, and partially 6, we utilized thismethod to investigate certain ultrafast dynamics. Diatomic systems are more complex than atomic
systems, posing greater challenges in obtaining their corresponding TDSE solutions. To address
this issue, we propose an approach in this Chapter to extend the methodology outlined in section
2.2, specifically for diatomic systems. We focus on homo-nuclear diatomic systems and, similar
to before, consider only one active electron. We first provide the numerical details for our abinitio solution to the single-active-electron TDSE for diatomic systems interacting with an intenselaser pulse. We then present some preliminary results of the diatomic ab-initio TDSE solution
by conducting calculations for bound states and generating photoelectron spectra for hydrogenmolecules interacting with ultrafast laser pulses.

7.1

Additional approximations
Similar to the ab-initio solution for atoms discussed in section 2.2, the TDSE solution fora diatomic system requires several approximations. In the case of atoms, the approximations weused were the classical treatment of the field, dipole and non-relativistic approximations, and theuse of SAE potentials. The classical treatment of the field is still valid for the diatomic system
as it relates to the laser pulse and not the target system. Although diatomic systems are larger

96
than their atomic counterparts, the wavelength of a typical laser system is still much larger than
the size of the molecule and so the spatial variation of the field can be ignored. Thus, the dipole
approximation is still valid. The SAE potential serves the same function as in the atomic casewhere we use an effective potential to include the effect of the non-active electrons in our model.
For the diatomic system, it is calculated in a similar manner as described in section 2.2.3. The
diatomic potential is obtained from DFT and fitted into a functional form [199]. The functional
form used is given by

VSAE (r; Ri ) =

2
X
i=1

−

Z
Zc e−c|r−Ri | X
an e−bn |r−Ri | + de−oco −ρcρ ,−−|r − Ri ||r − Ri |n(7.1)

where c, an , and bn are fit parameters, Zc is the shielded charge, r is the position of the electron and
Ri is the position of the ith nucleus. As shown in Fig. 7.1, the variables o and ρ correspond to thecoordinate along the axis of the molecule and the radial distance perpendicular to the molecularaxis respectively. As in the case of the atomic SAE potential, the functional form in Eq. (7.1) has
Coulomb, Yukawa, and exponential terms centered on each nucleus. In addition to this, we also
have a cylindrical term that is oriented along the axis between the two nuclei. More details about
the diatomic SAE potential can be found in Refs. [199, 200].

O

ρ

α
Z

Figure 7.1: Illustration depicting the orientation of a diatomic molecule, where the axes o andρ represent the coordinate along the molecule axis and the radial distance perpendicular to themolecular axis, respectively.

97
In addition to these approximations, we make use of the Born-Oppenheimer approximation[201] for the ab-initio solution for a diatomic system. In the Born-Oppenheimer approximation, themotion of the electron is considered to be independent of the motion of the nuclei. This is because of
the relatively slow motion of the nuclei compared to that of the electrons. The approximation allows
us to solve for the electronic wave function independently of the nuclear dynamics. The TDSE forthe electronic component of a diatomic system within the Born-Oppenheimer approximation canbe written as


∂
1 2
i ψ(r, t; R1 , R2 ) = − ∇ + VSAE (r; R1 , R2 ) + Vint (r, t) ψ(r, t; R1 , R2 ),
∂t
2

(7.2)

where ψ(r, t; R1 , R2 ) is the electronic wave function, Vint is the laser coupling term as described in
Chapter 2 and VSAE is the diatomic SAE potential. With the approximations mentioned above,
the TDSE can be solved numerically. The procedure for formulating an ab-initio solution is similar
to the case for atoms. The main difference lies in the SAE potential structure. In the next sections,
we describe the numerical method and use the results from section 2.2 when appropriate.

7.2

Matrix elements
In order to solve the TDSE using the ab-initio method, it is necessary to expand the wave

function and the Hamiltonian in a selected basis. In section 2.2, we explained this approach indetail, using a basis consisting of radial functions and spherical harmonics. We represent theradial portion on a grid and use spherical harmonics to describe the angular components of oursystem. To reuse the results from the solution of the atomic system, we decided to utilize the same
basis as before, which can be denoted as |lm⟩, where l and m represent the orbital and magnetic
quantum numbers, respectively. In section 2.2, we computed the matrix elements for all terms inthe Hamiltonian using this basis. As the kinetic energy and laser coupling terms are identical forboth atomic and diatomic systems, we can use the same matrix elements for those two operators.
The kinetic energy and laser coupling matrix elements can be found in sections 2.2.2 and 2.2.4,

98respectively. To calculate the matrix elements for the diatomic potential, we need to compute theinner product of the potential with the same basis. The resulting matrix elements can be expressed
as follows:
SAE
Vlml′ m′ =

X

lm |VSAE (r; Ri )| l′ m′ .

(7.3)

i

The potentials are categorized into two groups according to their symmetry. The Coulomb, Yukawa,
and exponential terms, which exhibit spherical symmetry about the position of the nuclei, comprisethe first three terms of the potential. We label this group as the spherical part of the potential. The
remaining term we will refer to as the cylindrical term. In the subsequent sections, we demonstrate
how to compute the matrix elements for both the spherical and cylindrical components of the
diatomic SAE potential.

7.2.1

Spherical part of the potential

The spherical part of the diatomic SAE potential is given by

V sph (|r − Ri |) =

2
X

−

i=1

Z
Zc e−c|r−Ri | X−+
an e−bn |r−Ri | .
|r − Ri ||r − Ri |n(7.4)

The first step in calculating the matrix elements is to express the potential in spherical coordinates.
This can be done by expanding each term in a basis of radial functions and spherical harmonics asfollows:∞l

X X 4π rl
1<=
Y m∗ (θi , ϕi )Ylm (θ, ϕ),l+1 l|r − Ri |2l + 1 r>(7.5)

l=0 m=−l∞l

X X
e−c|r−Ri |= 8c
Il (cr< )Kl (cr> )Ylm∗ (θi , ϕi )Ylm (θ, ϕ)|r − Ri |(7.6)

l=0 m=−l−c|r−Ri |

e

=

∞ X
l 
X


8 (Il (cr< )Kl (cr> )) + 8c Il′ (cr< )Kl (cr> )+l=0 m=−l

(7.7)

′
Il (cr< )Kl (cr> ) Ylm∗ (θi , ϕi )Ylm (θ, ϕ)99
where θi and ϕi are the angular coordinates of the nuclei and Il and Kl are the modified sphericalBessel functions of the first and second kind, respectively. r> and r< correspond to the greater andlesser value between r and Ri . Due to the shared symmetry of the three terms, the expansions all
have the same functional form. The expansions have a radial function labeled by l and sphericalharmonics products for the angles representing the position of the electron and the position of thenuclei. Substituting the expansions, Eqs. (7.5), (7.6) and (7.7), into the spherical potential reduces
it to

V sph =∞ X
L
X

FL (r; R1 )YLm (θ, ϕ)YLm∗ (θR1 , ϕR1 )+L=0 m=−L

(7.8)

∞ X
L
X

FL (r; R2 )YLm (θ, ϕ)YLm∗ (θR2 , ϕR2 ),L=0 m=−L

where FL is a radial function given by

FL (r; R) = −Z

L
X
r<
4π
−
8cZ
I
(cr
)K
(cr
)
−
8a
(IL (bn r< )KL (bn r> )) +
o
<
>
n
L
L
L+1
(2L + 1) r>
n

(7.9)

bn IL′ (bn r< )KL (bn r> )r< + IL (bn r< )KL′ (bn r> )r> .
Combining the sum over the quantum numbers, the spherical potential is given by

V sph =∞ X
L
X

[FL (r; R1 )YLm∗ (θR1 , ϕR1 ) + FL (r; R2 )YLm∗ (θR2 , ϕR2 )] YLm (θ, ϕ).

(7.10)

L=0 m=−L

In order to simplify the calculation, we can orient the molecule along the z-axis without loss ofgenerality. Additionally, we can position the center of the coordinate system at the midpoint
between the molecules, ensuring that both nuclei have the same radial position. With this chosen
position and orientation, the positions of the nuclei in spherical coordinates will be as follows:R, 0, 0
R1 =
2


R
R2 =
, π, 0 ,2(7.11)

100
where R is the internuclear distance. For this configuration, the spherical potential can be writtenasV sph =∞ X
L
X

FL (r;

L=0 m=−L

R m)Y (θ, ϕ) [YLm∗ (0, 0) + YLm∗ (π, 0)] .
2 L(7.12)The expression for the potential can be further simplified by substituting the value of the sphericalharmonics at the coordinates of the nuclei which are given by

Ylm∗ (π, 0) =





0,

m ̸= 0q(−1)l 2l+1
4π , m = 0




0,
m ̸= 0m∗
Yl (0, 0) = q 2l+14π , m = 0 .

(7.13)As indicated by Eq. (7.13), the summation over the magnetic quantum number reduces to theterms where m is equal to zero. By substituting the values for the spherical harmonics, the sphericalcomponent of the SAE potential can be expressed in a simplified form as

V

sph

=2

∞
X

r

L

"
#
r
r
2L + 1R 0
2l + 12l + 1l
FL (r; )YL (θ, ϕ) (−1)
+
.
4π
2
4π
4π

(7.14)

When looking at the term in the bracket in Eq. (7.14), we recognize that only the even terms
contribute to the sum over the angular momentum quantum number l and we can simplify theexpression further. This simplification arises due to the cylindrical symmetry of the diatomic
potential. In the subsequent section we will observe similar behavior for the final term of the SAE
potential. After implementing these simplifications, the spherical potential can be expressed as

V

sph

=2

∞
X
L=0,2,4r

2L + 1R
FL (r; )YL0 (θ, ϕ).
4π2(7.15)

101
With this simplified form, the matrix elements can be computed by taking the inner product of theexpression provided in Eq. (7.15), and we get
∞
X

*

+


2L + 1R0′ ′lm 2
YL (θ, ϕ) l m
FL r;
4π
2
L=0,2,4
r


Z
∞
X
R
2L + 1′m∗
0=2
Yl (θ, ϕ)FL r;
Ylm′ (θ, ϕ)YL (θ, ϕ) dΩ
4π
2
4π

sph
Vlml
′ m′ =

r(7.16)

L=0,2,4

Taking the radial function out of the integral, we are left with an integral of three spherical harmonics over a solid angle.

sph
Vlml
′ m′ = 2∞X
L=0,2,4rZ
2L + 1R′
0
FL r;
Ylm∗ Ylm′ YL dΩ.
4π
24π(7.17)The integral over three spherical harmonics can be expressed in terms of Wigner 3j symbols using
the following relation


Z

r

′

0
m
Ylm∗ Ylm′ YL = (−1)

′L  ll L
(2l + 1)(2l′ + 1)(2L + 1)  l


.
4π
0 0 0
−m m′ 0l′

(7.18)

With this, the matrix elements for the spherical part of the diatomic SAE potential are given by





′
′
l L
R (2L + 1) p
 l l L  l
sph
m
FL r;
Vlml
(2l + 1)(2l′ + 1) .
′ m′ = 2(−1)
24π′L=0,2,40 0 0−m m 0∞X(7.19)
7.2.2

Cylindrical part of the potential

To calculate the matrix elements of the cylindrical term, we need a different approach. Since
we have chosen to align the molecule along the z-axis, the cylindrical part of the SAE potential isgiven by2

2

V cyl = de−(zcz ) −(ρcρ )(7.20)

102
where ρ =p
x2 + y 2 . Similar to the spherical terms discussed in the previous section, we aim torepresent the potential using spherical coordinates. This is achieved by factorizing it into a radialfunction and spherical harmonics. With this goal in mind, we complete the square of the terms in
the exponent to write the potential as

2

2

2

2

2

2

2

2

V cyl = de−(zcz ) −(ρcρ ) −(zcρ ) +(zcρ ) = de−(rcρ ) e−z (cz −cρ ) ,where r =p

(7.21)

ρ2 + z 2 . Expanding the z-axis in terms of spherical harmonics allows us to represent

the expression in Eq. (7.21) using spherical coordinates. The cylindrical potential can then be
written as
V

cyl

q
4π
0 2 2
2
−(rcρ )2 −( 3 rY1 ) (cz −cρ )

= de

e

r
,

where z =

4π 0
rY .
3 1

(7.22)The expression can be simplified by combining constants as

V

cyl

= de

−(rcρ )2 −(c̄rY10 )2er2, where c̄ =4π 2(c − c2ρ ).
3 z(7.23)

To get the functional form we are seeking, we expand the second exponential term in Eq. (7.23)
using a Taylor series. With this expansion the cylindrical potential is given by

V

cyl

= de

−(rcρ )2

∞
X

(−1)n

n=0

(c̄rY10 )2n.
n!

(7.24)

This can be further simplified by noting that the product of two spherical harmonics can be expressed in terms of a sum of spherical harmonics. The relation is given by

r
1
2
Ylm
(θ, ϕ)Ylm(θ, ϕ) =1
2

X
l3 ,m3




l
l
l
l
l
l
2
3  m
(2l1 + 1)(2l2 + 1)(2l3 + 1)  1 2 3   1

 Yl3 3 (θ, ϕ).
4π
0 0 0m1 m2 m3
(7.25)

Each term in the Taylor series can be recursively written as a product of two spherical harmonics.
This allows us to express all terms in the series as a sum over spherical harmonics. Using the

103
relation given in Eq. (7.25), the cylindrical potential can be written as

V cyl = de−(rcρ

)2

∞
X

Pn
(−1)n (c̄r)2n0
L=0,2,... aL YL

n!

n=0

.

(7.26)

To determine the value of aL , we apply Eq. (7.25) recursively and aggregate the coefficients
corresponding to each spherical harmonic degree for every term present in the Taylor series. With
the cylindrical potential written in the form given by Eq. (7.26), the matrix elements can now be
calculated. The matrix elements are given by
*
cyl
Vlml
′ m′ =

lm de−(rcρ

)2

∞
X

(c̄r)2n
(−1)n

Pn

n=0

+

0
L=0,2,... aL YL

l′ m′

n!

(7.27)cyl−(rcρ )2
Vlml
′ m′ = de∞X(−1)n (c̄r)2n
n=0

n!

n
X
L=0,2,...

Z
aL
4π

′

0
Ylm∗ Ylm′ YL dΩAs before we apply Eq. (7.25) to calculate the integral over the spherical harmonics and get thefollowing expression for the matrix elements of the cylindrical potential2cyl−(rcρ )
Vlml
′ m′ = de∞X
n=0(−1)n (c̄r)2nn!

nX
L=0,2,...


r
aL

l′

(2l + 1)(2l′ + 1)(2L + 1)  l

4π
−m m′


L
 . (7.28)
0After obtaining the matrix elements for both the spherical and cylindrical components of the
SAE potential, we can use those to get the complete expression for the matrix elements. The total

104
matrix elements for the SAE diatomic potential in the |l, m⟩ basis can be written as′L  ll L
R (2L + 1) pl
Vlml′ m′ = 2(−1)mFL r;(2l + 1)(2l′ + 1) 

+
2
4π
′
L=0,2,40 0 0−m m 0



r
′
′
n
∞
X (−1)n (c̄r)2n X
l L
(2l + 1)(2l′ + 1)(2L + 1)  l l L  l
2de−(rcρ )
aL



n!
4π
n=0
L=0,2,...
0 0 0−m m′ 0
∞
X







l′

(7.29)
where

FL

R
r;
2


= −Z

L
X
r<
4π
− 8cZo IL (cr< )KL (cr> ) − an 8
(IL (bn r< )KL (bn r> )) +L+1
(2L + 1) r>n(7.30)

bn IL′ (bn r< )KL (bn r> )r< + IL (bn r< )KL′ (bn r> )r> .


Eq. (7.29) expresses the matrix elements through two infinite series: one for the Taylor series andanother for the angular momentum quantum number. In practice, limiting the number of terms
in the Taylor series index n to about ten, allows for the precise representation of the cylindricalpotential within machine precision. The other infinite sum involves the index L for the expansionof the spherical terms of the potential. This sum can be truncated by using the following propertiesof the Wigner 3j symbol.

 


′
0,
 l l L
=





0 0 0
0,




l′
 l

−m m′


 


L 0,
=



0
0,

for L < |l − l′ | or L > l + l′
for l + l′ + L is odd(7.31)
for L < |l − l′ | or L > l + l′
for m ̸= m′

With the above relations, we can write the matrix elements for the potential as

105

′

l+l
X

R (2L + 1) p
(2l + 1)(2l′ + 1)+)2
4π
L=|l−l′ |



!
r
′
′
nX
n
max
n
2n
′
X
l L
(−1) (c̄r)(2l + 1)(2l + 1)(2L + 1)  l l L  l
2de−(rcρ )
aL


.
n!
4π
′
′
n=0
L=|l−l |0 0 0−m m 0m

Vlml′ m′ = δm,m′ δ(l+l′ )even (−1)

2FL (r;(7.32)

7.3

Bound states for diatomic systems
There are some notable differences between the matrix elements of the diatomic and atomic

systems. In the diatomic case, there is a coupling between different angular momentum states.
This results from the break of the spherical symmetry of the system. So, the states of the diatomic
systems cannot be labeled using the l quantum number. The principal and magnetic quantum
numbers are still good quantum numbers and can be used to label the states. We can get an
additional label for the states by considering the parity operator Π̂ which commutes with the
Hamiltonian. The parity operator is related to the symmetry of the system under parity inversion.
For a wave function Ψ(r),Π̂Ψ(r) = (−1)p Ψ(−r)(7.33)

where p = 0 corresponds to even symmetry and p = 1 to odd symmetry. The states exhibiting evenand odd symmetry are alternatively called ”gerade” and ”ungerade” states, respectively. To identify
the symmetry of the states, one can examine the matrix elements of the diatomic potential andobserve the parity of the spherical harmonics. Specifically, the even states contain contributions
from spherical harmonics components with even values of l, such as l = 0, 2, 4..., while the odd
states contain contributions from those with odd values of l, such as l = 1, 3, 5... The states will
be labeled using the principal and magnetic quantum numbers followed by an underscript of u andg for the gerade and ungerade states respectively. We use the spectroscopic notation to designate
the states with m = 0, 1, 2... with the labels λ = σ, π, δ...

106
As stated above the fit parameters in the SAE potential given in Eq. (7.1) are obtained byfitting the functional form to a potential obtained using density functional theory (DFT). For the
diatomic system, we have obtained the fit parameters for hydrogen and oxygen molecules [200].
The fit parameters are given in the table below.

H2
O2

R
Co
1.4314 12.3000 1

Diatomic SAE potential fit parameters
Zc
c
co
cρ
d
a1
b1
1
0.316 13.841 18.637 −0.162 1.205 0.497
15
1.326 3.248 2.778 1.643 10.804 0.466

Table 7.1: Parameters for molecular single active electron potential obtained using analytical fitsof potentials for hydrogen and oxygen molecules [200].

To validate the approach for the bound states, we compare our results with those frommolecular orbital theory (MOT). The states from our method are obtained by solving the timeindependent Schrödinger equation (TISE) via diagonalization of the field-free Hamiltonian and wecompared the calculated states to the MOT orbitals. Given that the wave functions exist in three
dimensions, we take a slice across the x − z plane ensuring that we are correctly depicting thesymmetry of the state. To better compare the results with the MOT orbitals, we incorporate thephase of the wave function into the magnitude, thereby displaying their symmetry.
MOT is a model used to describe the electronic structure of molecules in terms of their constituent atoms and the interactions between them. The molecular orbitals are formed by combining
the atomic orbitals of the constituent atoms. An example of molecular orbitals from MOT for the
first two states of the hydrogen molecule are shown in Fig. 7.2. The 1σg and 1σu states are formedby combining two hydrogen 1s orbitals. It is worth noting that the MOT orbitals in this thesis areused as an illustration and are not depicted to scale.

107

Figure 7.2: Lowest-energy orbitals for diatomic (hydrogen) molecule, where 1σg (lower part) and
1σu (upper part) orbitals are formed by combining two 1s hydrogen atom orbitals (taken from [202]).

Figure 7.3: Hydrogen molecule orbitals from TISE calculations, where left is the 1σg - and right is
the 1σu -orbital, respectively. The color corresponds to |Ψ|2 eiϕ where ϕ is the phase.

The orbitals from the TISE solution presented in Fig. 7.3 match the MOT orbitals, having thesame symmetries of the phase and the shape of the orbitals. We also make the same comparison
for the 1πg and 1πu orbitals. These orbitals from MOT and TISE solutions are shown in Fig. 7.4
and Fig. 7.5, respectively.

108

Figure 7.4: Same as Fig. 7.2 but for 1πg (upper part) and 1πu (lower part) orbitals, which areformed by combining two 2p hydrogen atom orbitals (taken from [202]).

Figure 7.5: Same as Fig. 7.3 but for 1πg and 1πu orbital. The color corresponds to |Ψ|2 eiϕ where
ϕ is the phase.

The comparison of TISE solutions and MOT is primarily qualitative in nature. However,
MOT can be utilized as a means to confirm the expected structure of TISE states. The MOT
orbitals can also help in labeling and understanding the energy ordering of the states calculated by

109
the ab-initio solution.

7.4

Time propagation and observables
Besides analyzing the bound states of diatomic systems, we also investigated their dynamics,

specifically by solving the full TDSE that entails propagating the wave function using the entire
Hamiltonian. Similar to our atomic solution (as explained in section 2.2.7), we employ the CrankNicolson method for the unitary time evolution of the wave function. As stated before, for bothdiatomic and atomic systems, we utilize the same matrix elements for the kinetic energy and lasercoupling terms.
Observables such as population in excited states and ionization are calculated in the samemanner as in the atomic case. These observables are given by
Pϕn = |⟨ϕn |Ψ⟩|2 ,Pionization = 1 −X(7.34)2|⟨ϕn |Ψ⟩| ,nwhere Ψ and ϕn are the time propagated wave function and bound states of the diatomic system,
respectively. Another observable discussed previously in section 2.2.8 is the photoelectron spectrum.
The photoelectron spectrum is calculated by projecting the time propagated wave function on tocontinuum states of the system and is given by

1 X
F (k, ϕ, θ) = 2kZe−iδkl (i)l ϕkl (r)∗ ψ(r, t) dx

2,(7.35)

l,m

where δkl is the phase shift, and ϕkl (r) is a continuum state with momentum k. In atomic systems, there are two approaches to determine the continuum state: using the shooting method
or diagonalizing the field-free Hamiltonian. However, for diatomic systems, the shooting method
is not applicable since the potential lacks spherical symmetry, requiring the Hamiltonian to bediagonalized instead.

110
Examples of photoelectron spectra of the hydrogen molecule resulting from the interaction ofa hydrogen molecule with a laser pulse of 1.5 × 1014 W/cm2 intensity is shown in the Figures below.
We performed two calculations, each consisting of 10 cycles, with laser wavelengths of 200 nm and
400 nm, respectively. Fig. 7.6 shows the peaks in the above-threshold ionization (ATI) processes.
As expected, the energy separation between each peak is equivalent to the energy of a single
photon. Thus, each peak corresponds to ionization by the absorption of an additional photon. The
photoelectron spectrum serves as a good test for the time propagation and the calculated continuum
states from diagonalization. Along with the bound state calculations, the results demonstrate thesuccessful implementation of our method to solve the TDSE for a diatomic system within the SAE

Yield

potential approximation.

Electron energy

Electron energy

Figure 7.6: photoelectron spectra for hydrogen molecule interacting with a 10 cycle laser pulse of
200 nm (left) and 400 nm (right) and a peak intensity of 1.5 × 1014 W/cm2 .

111

7.5

Current performance and future work
The parallel implementation of the TDSE solver for diatomic systems has been successfully

achieved on high-performance computing systems. To achieve parallelization, the code utilizes the
SLEPc and PETSc libraries. Initial tests demonstrate linear scaling with the number of cores used,
as shown in Fig. 7.7, for both bound state calculations and wave function time propagation. It
is important to note that these calculations are computationally more demanding as compared tothose for atomic systems due to the inclusion of additional off-diagonal matrix elements resulting
from the coupling of angular momentum states. The off-diagonal elements lead not only to alarger memory requirement but also to more communication between different cores during matrix
operations. This results in an increase in time for diagonalizing the Hamiltonian and for solving
the linear system of equations for the time propagation.

Figure 7.7: Hydrogen molecule bound state calculation and time propagation scaling as a functionof utilized cores.

In future work, we plan to further develop our ab-initio method by employing B-splines to represent

112
the radial part of the wave function instead of the finite difference method. B-splines offer acomplete basis and enable non-uniform grid sampling, which is crucial for achieving fine samplingnear the nuclei and coarser sampling in the outer regions of the grid. This results in fewer points
needed to represent the radial part of the wave function. This improvement applies not only to
diatomic systems but also to ab-initio methods in general. However, the diatomic system solution
stands to benefit significantly in terms of speed and memory requirements from these changes.
The results presented in this section illustrate the feasibility of using ab-initio solutions tosolve the TDSE for diatomic systems. This approach leads to numerous opportunities for advancing
research in the field of diatomic systems beyond the examples studied here. For instance, by obtaining the SAE potentials of other diatomic systems such as the nitrogen molecule, their dynamics
can be simulated. Moreover, non-homogeneous diatomic systems can be investigated by slightly
modifying the matrix computations outlined in section 7.2. The versatility of this method makes it
a promising tool for investigating the properties and dynamics of various diatomic systems, paving
the way for further exploration and discoveries in the field.

Chapter 8

Summary and outlook

In this thesis, we have investigated strong field interaction of atoms with circularly and elliptically polarized pulses. In Chapter 1 we gave an introduction to ultrafast physics and provided
background for the phenomena that we studied in the following Chapters. In Chapter 2 we discussed the time-dependent Schrödinger equation (TDSE), which is used to model the atom-laserinteraction. We introduced two numerical methods for the solution of the TDSE and discussed
the relevant approximations necessary for each method. The first was an ab-initio solution thatsolves the TDSE from first principles using the single-active-electron (SAE) approximation. The
second method is the quantum trajectory Monte Carlo (QTMC) method, a Monte Carlo simulation technique that combines analytical ionization rates with quantum trajectories. In subsequent
Chapters, we used these methods to study some aspects of ultrafast electron dynamics.
In Chapter 3, we examined the distribution of orbital angular momentum (l) and magnetic(m) quantum numbers in Rydberg states excited by the interaction with bichromatic circularly
polarized laser pulses of frequency ω and 2ω. We showed how selection rules determined the population of states with both l and m being either both even or both odd, irrespective of the laserparameters. For corotating pulses, we showed how adjusting the pulse intensities enabled control
over the distribution of the magnetic quantum number. For counterrotating pulses, we explained
the limitation of the maximum orbital angular momentum by using classical considerations. Furthermore, we analyzed the origin of the predominant population of Rydberg states with ∆m = ±3
based on selection rules for photon absorption.

114
In Chapter 4 we expanded on the work presented in the previous Chapter and proposed amechanism for the observed redistribution of population in the Rydberg states amongst different
magnetic states. We demonstrated that higher-order Raman transitions are effective for the observed redistribution. We also predicted that the mechanism is present for corotating pulses ofhigher frequency differences (ω and nω, n > 2).
In Chapter 5, we analyzed the difference in photoelectron energy spectra of atoms interactingwith intense circularly and elliptically polarized laser pulses for different initial magnetic states. We
demonstrated that an additional absorption pathway specific to the counterrotating initial state is
the source of the spectrum difference. We showed that the observed changes in angular emission
in ellipticity polarized pulses are due to interference caused by phase differences in orbital angularand magnetic states. To this end, we introduced a simple model with results that support ourclaim and explain that the counterrotating case exhibits a larger phase difference, leading to amore significant change in emission angle.
In Chapter 6, we applied the QTMC method to compute the photoelectron energy spectrum
generated by the interaction of atoms with circularly and elliptically polarized pulses. To accountfor the initial state information, we utilized the full PPT ionization rate within the QTMC method.
We compared our findings with the results obtained through the ab-initio method, demonstrating
the efficacy of the QTMC method in accurately generating the spectra. Further improvements to themethod can be made by using ionization models that do not rely on the saddle-point approximation.
With this, the method can be applied to a much wider range of laser parameters and provide evenbetter agreement with the ab-initio method.
Finally, in Chapter 7 we extended the ab-initio solution discussed in Chapter 2 to diatomic
systems. We showed that we can calculate bound states of diatomic systems in the SAE and BornOppenheimer approximation. We demonstrated a full solution to the TDSE by calculating thephotoelectron spectrum for hydrogen molecules. Future advancements can be made by adopting
a more efficient basis for representing the radial part of the wave function, such as the B-spline
method. This improvement would lead to faster computations and reduced memory usage when115
solving the TDSE for diatomic systems. Consequently, it would enable the exploration of a broader
range of captivating dynamics in these systems.

Bibliography[1] F. Krausz and M. Ivanov. Attosecond physics. Rev Mod. Phys. 81, 163 (2009).
[2] Corkum P. and F. Krausz. Attosecond science. Nature Phys 3, 381–387 (2007).
[3] M. F. Kling and M.J.J. Vrakking. Attosecond electron dynamics. Annu. Rev. Phys. Chem.
59, 463-492 (2008).
[4] G. Lukas, C. Claudio, and Ursula K. Attosecond science: Recent highlights and future trends.
Annu. Rev. Phys. Chem. 63, 447-469 (2012).
[5] M.C. Chen, C. Mancuso, C. Hernández-Garcı́a, F. Dollar, B. Galloway, D. Popmintchev,
P.C. Huang, B. Walker, L. Plaja, and et al Jaroń-Becker A. Generation of bright isolated
attosecond soft X-ray pulses driven by multicycle midinfrared lasers. Proc. Natl. Acad. Sci.
USA 111, E2361–E2367 (2014) .
[6] K. Zhao, Q. Zhang, M. Chini, Y. Wu, X. Wang, and Z. Chang. Tailoring a 67 attosecond
pulse through advantageous phase-mismatch,. Opt. Lett. 37, 3891-3893 (2012).
[7] A. Einstein. Über einen die Erzeugung und Verwandlung des Lichtes betreffenden heuristischen Gesichtspunkt. Ann. der Phys. 17, 132–148 (1905).
[8] M. Göppert-Mayer. Über Elementarakte mit zwei Quantensprüngen. Ann. der Phys. 401,
273–294 (1931).
[9] W. Kaiser and C. G. B. Garrett. Two-photon excitation in CaF2 :Eu2+ . Phys. Rev. Lett. 7,229 (1961).
[10] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich. Generation of optical harmonics.
Phys. Rev. Lett. 7, 118 (1961).
[11] L. Gallmann, I. Jordan, H.J. Wörner, L. Castiglioni, M. Hengsberger, J. Osterwalder, C. Arrell, M. Chergui, E. Liberatore, U. Rothlisberger, and U. Keller. Photoemission and photoionization time delays and rates. Struct. Dyn. 4, 061502 (2017).
[12] I. Shchatsinin. Free clusters and free molecules in strong, shaped laser fields. PhD thesis,
Max-Born-Institut fur Nichtlineare Optik und Kurzzeitspektroskopie Berlin, Germany, 20(2009).
[13] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman. Free-free transitions
following six-photon ionization of Xenon atoms. Phys. Rev. Lett 42, 1127 (1979).

117[14] E. Cormier and P. Lambropoulos. Above-threshold ionization spectrum of hydrogen usingb-spline functions. J. Phys. B: At. Mol. Opt. Phys. 30, 77 (1997).
[15] P. Agostini. The simple Rahman’s theory. Cambridge University Press, (2004).
[16] C. J. Joachain, N. J. Kylstra, and R. M. Potvliege. Atoms in intense laser fields. Cambridge
University Press (2011).
[17] N. B. Delone and V. P. Krainov. AC stark shift of atomic energy levels. Physics-Uspekhi 42669-687 (1999).
[18] P. Kruit, J. Kimman, H. G. Muller, and M. J van der Wiel. Electron spectra from multiphoton
ionization of Xenon at 1064, 532, and 355 nm. Phys. Rev. A 28, 248 (1983).
[19] H.G. Muller, A. Tip, and M.J. Van derWiel. Ponderomotive force and AC stark shift in
multiphoton ionisation. J. Phys. B: Atom. Mol. Phys. 16, L679 (1983).
[20] L.V. Keldysh. Ionization in the Field of a Strong Electromagnetic Wave. Sov. Phys. JETP
20, 1307 (1965).
[21] N.B. Delone and V.P. Krainov. Tunnel ionization of complex atoms and of atomic ions in analternating electromagnetic field. Sov. Phys. JETP. 64 (1986).
[22] M.V. Ammosov, N.B. Delone, and V.P. Krainov. Energy and angular electron spectra forthe tunnel ionization of atoms by strong low-frequency radiation. J. Opt. Soc. Am. B 8,
1207-1211 (1991).
[23] A.M. Perelomov, V.S. Popov, and M.V. Terent’ev. Ionization of Atoms in an Alternating
Electric Field. Sov. Phys. JETP 23, (1966).
[24] A.M. Perelomov, V.S. Popov, and M.V. Terent’ev. Ionization of Atoms in an Alternating
Electric Field:II. Sov. Phys. JETP 24, (1967).
[25] G.L. Yudin and M.Y. Ivanov. Nonadiabatic tunnel ionization: Looking inside a laser cycle.
Phys. Rev. A 64, 013409 (2001).
[26] I. Dreissigacker and M. Lein. Quantitative theory for the lateral momentum distribution afterstrong-field ionization. Chemical Physics 414, (2013).
[27] F.H.M. Faisal. Multiple absorption of laser photons by atoms. J. Phys. B: Atom. Mol. Phys.
6, L89 (1973).
[28] H.R. Reiss. Effect of an intense electromagnetic field on a weakly bound system. Phys. Rev.
A. 22, 1786 (1980).
[29] V. P. Krainov. Ionization rates and energy and angular distributions at the barrier-suppression
ionization of complex atoms and atomic ions. J. Opt. Soc. Am. B 14, 425-431 (1997).
[30] I. Barth and O. Smirnova. Nonadiabatic tunneling in circularly polarized laser fields. ii.
derivation of formulas. Phys. Rev. A 87, 013433 (2013).
[31] S. Camp, K.J. Schafer, and M.B. Gaarde. Interplay between resonant enhancement andquantum path dynamics in harmonic generation in helium. Phys. Rev. A 92, 013404 (2015).

118[32] L. Fechner, N. Camus, A. Krupp, J. Ullrich, T. Pfeifer, and R. Moshammer. Creation andsurvival of autoionizing states in strong laser fields. Phys. Rev. A 92, 051403 (2015).
[33] M. Fushitani, C.N. Liu, A. Matsuda, T. Endo, Y. Toida, M. Nagasono, T. Togashi,
M. Yabashi, T. Ishikawa, Y. Hikosaka, T. Morishita, and A. Hishikawa. Femtosecond twophoton rabi oscillations in excited He driven by ultrashort intense laser felds. Nature Photon
10, 102–105 (2016).
[34] X. Gao, G. Patwardhan, S. Schrauth, D. Zhu, T. Popmintchev, H.C. Kapteyn, M.M. Murnane, D.A. Romanov, R.J. Levis, and A.L. Gaeta. Picosecond ionization dynamics in femtosecond filaments at high pressures. Phys. Rev. A 95, 013412 (2017).
[35] N.A. Hart, J. Strohaber, A.A. Kolomenskii, G.G. Paulus, D. Bauer, and H.A. Schuessler.
Selective strong-field enhancement and suppression of ionization with short laser pulses. Phys.
Rev. A 93, 063426 (2016).
[36] I.A. Ivanov, C.H. Nam, and K.T. Kim. Photoionization in the presence of circularly polarized
fundamental and odd-order harmonic fields. Phys. Rev. A 95, 053401 (2017).
[37] S. Larimian, S. Erattupuzha, C. Lemell, S. Yoshida, S. Nagele, R. Maurer, A. Baltuška,
J. Burgdörfer, M. Kitzler, and X. Xie. Coincidence spectroscopy of high-lying Rydberg states
produced in strong laser fields. Phys. Rev. A 94, 033401 (2016).
[38] Q. Li, X.M. Tong, T. Morishita, C. Jin, H. Wei, and C.D. Lin. Rydberg states in the strong
field ionization of hydrogen by 800, 1200 and 1600 nm lasers. J. Phys. B 47, 204019 (2014).
[39] Q. Li, X.M. Tong, T. Morishita, H. Wei, and C.D. Lin. Fine structures in the intensity
dependence of excitation and ionization probabilities of hydrogen atoms in intense 800-nm
laser pulses. Phys. Rev. A 89, 023421 (2014).
[40] P.F. O’Mahony A. Galstyan B. Piraux, F. Mota-Furtado and Yu. V. Popov. Excitation of
Rydberg wave packets in the tunneling regime. Phys. Rev. A 96, 043403 (2017).
[41] M. Li, P. Zhang, S. Luo, Y. Zhou, Q. Zhang, P. Lan, and P. Lu. Selective enhancement ofresonant multiphoton ionization with strong laser fields. Phys. Rev. A 92, 063404 (2015).
[42] H. Lv, W. Zuo, L. Zhao, H. Xu, M. Jin, D. Ding, S. Hu, and J. Chen. Comparative study on
atomic and molecular Rydberg-state excitation in strong infrared laser fields. Phys. Rev. A
93, 033415 (2016).
[43] E. E. Serebryannikov and A. M. Zheltikov. Strong-Field Photoionization as Excited-State
Tunneling. Phys. Rev. Lett 116, 123901 (2016).
[44] K. Krajewska, I.I. Fabrikant, and A.F. Starace. Threshold effects in strong-field ionization:
Energy shifts and Rydberg structures. Phys. Rev. A 86, 053410 (2012).
[45] D.G. Arbó, K.I. Dimitriou, E. Persson, and J. Burgdörfer. Sub-poissonian angular momentum
distribution near threshold in atomic ionization by short laser pulses. Phys. Rev. A 78, 013406(2008).
[46] J. Venzke, R. Reiff, Z. Xue, A. Jaro ń Becker, and A. Becker. Angular momentum distribution
in Rydberg states excited by a strong laser pulse. Phys. Rev. A 98, 043434 (2018).

119[47] D.B. Milosevic. Possibility of introducing spin into attoscience with spin-polarized electrons
produced by a bichromatic circularly polarized laser field. Phys. Rev. A 93, 051402(R) (2016).
[48] C.A. Mancuso, D.D. Hickstein, J.L. Ellis K.M. Dorney, E. Hasovic, R. Knut, P. Grychtol,
C. Gentry, M. Gopalakrishnan, D. Zusin, F.J. Dollar, X.M. Tong, D.B. Milosevic, W. Becker,
H.C. Kapteyn, and M.M. Murnane. Controlling electron-ion rescattering in two-color circularly polarized femtosecond laser fields. Phys. Rev. A 93, 053406 (2016).
[49] M.M. Liu, Y. Shao, M. Han, P. Ge, Y. Deng, C. Wu, Q. Gong, and Y. Liu. Energy- andmomentum-resolved photoelectron spin polarization in multiphoton ionization of Xe by circularly polarized fields. Phys. Rev. Lett. 120, 043201 (2018).
[50] E. Arimondo and G. Orriols. Nonabsorbing atomic coherences by coherent two-photon transitions in a three-level optical pumping. Nuovo Cimento 17, 333–338 (1976).
[51] R. Whitley. Double optical resonance. University of Rochester, Rochester, N.Y. (1977).
[52] A. Ludwig, J. Maurer, B. W. Mayer, C. R. Phillips, L. Gallmann, and U. Keller. Breakdownof the dipole approximation in strong-field ionization. Phys. Rev. Lett. 113, 243001 (2014).
[53] R. Reiff, T. Joyce, A. Jaroń-Becker, and A Becker. Single-active electron calculations ofhigh-order harmonic generation from valence shells in atoms for quantitative comparison
with TDDFT calculations. J. Phys. Commun. 4, 065011 (2020).
[54] V. Hernandez, J.E. Roman, and V. Vidal. SLEPc: A scalable and flexible toolkit for thesolution of eigenvalue problems. ACM Trans. Math. Software 31.3.351–362 (2005).
[55] T. Sauer. Numerical Analysis. Pearson Education, Inc., Boston, 2012.
[56] J. Venzke. Theoretical studies on imaging of electron motion and excitation on ultrafast timescales. PhD thesis, University of Colorado at Boulder, 2020.
[57] S. Balay, W. Gropp, L.C. McInnes, and B.F Smith. PETSc, the portable, extensible toolkit
for scientific computation. Argonne National Laboratory 2, 17 (1998).
[58] L.D. Dalcin, R.R. Paz, P.A. Kler, and A. Cosimo. Parallel distributed computing using
python. Advances in Water Resources, 34.9,1124-1139(2011).
[59] F. He, C. Ruiz, and A. Becker. Absorbing boundaries in numerical solutions of the timedependent schrödinger equation on a grid using exterior complex scaling. Phys. Rev. A. 75,
053407 (2007.
[60] B. Fornberg. Generation of finite difference formulas on arbitrarily spaced grids. Math. Comp.
51, 699-706 (1988), .
[61] A.N Pfeiffer, C. Cirelli, M. Smolarski, D. Dimitrovski, M. Abu-samha, L.B. Madsen, and
U. Keller. Attoclock reveals natural coordinates of the laser-induced tunnelling current flow
in atoms. Nat. Phys. 8, 76-80 (2012).
[62] A. N. Pfeiffer, C. Cirelli, A. S. Landsman, M. Smolarski, D. Dimitrovski, L. B. Madsen, and
U. Keller. Probing the longitudinal momentum spread of the electron wave packet at thetunnel exit. Phys. Rev. Lett. 109, 083002 (2012).

120[63] C. Hofmann, A.S. Landsman, C. Cirelli, A.N. Pfeiffer, and U. Keller. Comparison of different
approaches to the longitudinal momentum spread after tunnel ionization. J. Phys. B: At.
Mol. Opt. Phys. 46, 125601 (2013).
[64] I. Barth and O. Smirnova. Nonadiabatic tunneling in circularly polarized laser fields. ii.
derivation of formulas. Phys. Rev. A. 87, 013433 (2013).
[65] M. Li, M.M. Liu, J.W. Geng, M. Han, X. Sun, Y. Shao, Y. Deng, C. Wu, L.Y. Peng, Q. Gong,
and Y. Liu. Experimental verification of the nonadiabatic effect in strong-field ionization withelliptical polarization. Phys. Rev. A. 95, 053425 (2017).
[66] R. Aav and K.A Mishra. The breaking of symmetry leads to chirality in cucurbituril-type
hosts. Symmetry 10, 98 (2018).
[67] Z. Xiao, W. Quan, S. Yu, X. Lai, X. Liu, Z. Wei, and J. Chen. Nonadiabatic strong field
ionization of noble gas atoms in elliptically polarized laser pulses. Opt. Express 30, 1487314885 (2022).
[68] M. Li, J.W. Geng, M. Han, M.M. Liu, L.Y. Peng, Q. Gong, and Y. Liu. Subcycle nonadiabatic
strong-field tunneling ionization. Phys. Rev. A. 93, 013402 (2016).
[69] T.M. Yan and D. Bauer. Sub-barrier Coulomb effects on the interference pattern in tunnelingionization photoelectron spectra. Phys. Rev. A. 86, 053403 (2012).
[70] M. Li, J.W. Geng, H. Liu, Y. Deng, C. Wu, L.Y. Peng, Q. Gong, and Y. Liu. Classicalquantum correspondence for above-threshold ionization. Phys. Rev. Lett. 112, 113002 (2014).
[71] N.I. Shvetsov-Shilovski, M. Lein, L.B. Madsen, E. Räsänen, C. Lemell, J. Burgdörfer, D. G.
Arbó, and K. Tőkési. Semiclassical two-step model for strong-field ionization. Phys. Rev. A.
94, 013415 (2016).
[72] D. Trabert, N. Anders, S. Brennecke, M. S. Schöffler, T. Jahnke, L. Ph. H. Schmidt, M. Kunitski, M. Lein, R. Dörner, and S. Eckart. Nonadiabatic strong field ionization of atomic
hydrogen. Phys. Rev. Lett. 127, 273201 (2021).
[73] S. Yu, X. Lai, S. Xu, Y. Wang, L. Hua, W. Quan, and X. Liu. Nonadiabatic effect on temporal
double-slit interference structures in an orthogonal two-color laser field. Phys. Rev. A. 105,
013116 (2022).
[74] L. Yupeng, W. Xie, M. Li, C. Cao, Y. Zhou, and P. Lu. Nonadiabatic tunneling ionizationof atoms in few-cycle elliptically polarized laser pulses. J. Phys. B: At. Mol. Opt. Phys. 56
105601 (2023).
[75] N.I. Shvetsov-Shilovski, D. Dimitrovski, and L.B. Madsen. Ionization in elliptically polarized pulses: Multielectron polarization effects and asymmetry of photoelectron momentum
distributions. Phys. Rev. A 85, 023428 (2012).
[76] N.I. Shvetsov-Shilovski, S.P. Goreslavski, S.V. Popruzhenko, and W. Becker. Capture into
Rydberg states and momentum distributions of ionized electrons. Laser. Phys. 19, 1550
(2009).

121[77] Y.N. Qin, M. Li, Y. Feng, S. Luo, Y. Zhou, and P. Lu. Extracting the phase distributionof the electron wave packet ionized by an elliptically polarized laser pulse. Front. Phys. 16,32502 (2021).
[78] G. Casella, C.P. Robert, and M.T. Wells. Generalized accept-reject sampling schemes. Lecture
Notes-Monograph Series 45, 342–347 (2004).
[79] S. Chelkowski and A. Bandrauk. Sensitivity of spatial photoelectron distributions to theabsolute phase of an ultrashort intense laser pulse. Phys. Rev. A 65, 061802 (2002).
[80] J. Venzke, T. Joyce, Z. Xue, A. Becker, and A. Jaron-Becker. Central frequency of few-cycle
laser pulses in strong-field processes. Phys. Rev. A 98, 063409 (2018).
[81] Z. Xue A. Jaron-Becker J. Venzke, R. Reiff and A. Becker. Angular momentum distribution
in Rydberg states excited by a strong laser pulse. Phys. Rev. A 98, 043434 (2018).
[82] H. Eichmann, A. Egbert, S. Nolte, C. Momma, B. Wellegehausen, W. Becker, S. Long, andJ.K. McIver. Polarization-dependent high-order two-color mixing. Phys. Rev. A. 51, R3414(1995).
[83] S. Long, W. Becker, and J.K. McIver. Model calculations of polarization-dependent two-color
high-harmonic generation. Phys. Rev. A. 52, 2262 (1995).
[84] A. Fleischer, O. Kfir, T. Diskin, P. Sidorenko, and O. Cohen. Spin angular momentum andtunable polarization in high-harmonic generation. Nat. Photon. 8, 543 (2014).
[85] E. Pisanty, S. Sukiasyan, and M. Ivanov. Spin conservation in high-order-harmonic generation
using bicircular fields. Phys. Rev. A. 90, 043829 (2014).
[86] O. Kfir, P. Grychtol, E. Turgut, R. Knut, D. Zusin, D. Popmintchev, T. Popmintchev,
H. Nembach, J.M. Shaw, A. Fleischer, H. Kapteyn, M. Murnane, and O. Cohen. Spin conservation in high-order-harmonic generation using bicircular fields. Nat. Photon. 9, 99–105(2015).
[87] T. Fan, P. Grychtol, R. Knut, C. Hernandez-Garcia, D.D. Hickstein, D. Zusin, C. Gentry,
F.J. Dollar, C.A. Mancuso, C.W. Hogle, O. Kfir, D. Legut, J.L. Ellis K. Carva, K.M. Dorney,
C. Chen, O.G. Shpyrko, E.E. Fullerton, O. Cohen, P.M. Oppeneer, D.B. Milosevic, A. Becker,
A. Jaron-Becker, T. Popmintchev, M.M. Murnane, and H.C. Kapteyn. Bright circularly
polarized soft x-ray high harmonics for x-ray magnetic circular dichroism. Proc. Natl. Acad.
Sci. U.S.A. 112, 14206–14211 (2015).
[88] D.D. Hickstein, F.J. Dollar, P. Grychtol, R. Knut J. Ellis, C. Hernandez-Garcia, D. Zusin,
C. Gentry, J. Shaw, T. Fan, K. Dorney, A. Becker, A. Jaron-Becker, H.C. Kapteyn, and M.M.
Murnane. Non-collinear generation of angularly isolated circularly polarized high harmonics.
Nat. Photon. 9, 743 (2015).
[89] K.M. Dorney, F. Tingting, Q.D. Nguyen, J.L. Ellis, D.D. Hickstein, N. Brooks, D. Zusin,
C. Gentry, H.C. Kapteyn C. Hernández-Garcı́a, and M.M. Murnane. Bright, single helicity,
high harmonics driven by mid-infrared bicircular laser fields. Opt. Express 29, 38119-38128(2021).

122[90] D.B. Milosevic. High-order harmonic generation by a bichromatic elliptically polarized field:
conservation of angular momentum. J. Phys. B: At. Mol. Opt. Phys. 48, 171001 (2015).
[91] L. Medisauskas, J. Wragg, H. van der Hart, and M.Y. Ivanov. Generating isolated ellipticallypolarized attosecond pulses using bichromatic counterrotating circularly polarized laser fields.
Phys. Rev. Lett. 115, 153001 (2015).
[92] D.B. Milosevic. Generation of elliptically polarized attosecond pulse trains. Opt. Lett. 40,
2381 (2015).
[93] C. Chen et al. Tomographic reconstruction of circularly polarized high-harmonic fields: 3d
attosecond metrology. Sci. Adv. 2, e1501333 (2016).
[94] N. Lin D. Baykusheva, M.S. Ahsan and H.J. Wörner. Bicircular high-harmonic spectroscopy
reveals dynamical symmetries of atoms and molecules. Phys. Rev. Lett. 116, 123001 (2016).
[95] D.D. Hickstein T. Popmintchev A. Meier M.M. Murnane H.C. Kapteyn I.J. Sola A. JaronBecker C. Hernandez-Garcia, C.G. Durfee and A. Becker. Schemes for generation of isolated
attosecond pulses of pure circular polarization. Phys. Rev. A. 93, 043855 (2016).
[96] X. Liu, X. Zhu, L. Li, Y. Li, Q. Zhang, P. Lan., and P. Lu. Selection rules of high-orderharmonic generation: Symmetries of molecules and laser fields. Phys. Rev. A. 94, 033410(2016).
[97] F. Mauger, A.D. Bandrauk, and T. Uzer. Circularly polarized molecular high harmonic
generation using a bicircular laser. J. Phys. B: At. Mol. Opt. Phys. 49, 10LT01 (2016).
[98] A.D. Bandrauk, F. Mauger, and K.J. Yuan. Circularly polarized harmonic generation byintense bicircular laser pulses: electron recollision dynamics and frequency dependent helicity.
J. Phys. B: At. Mol. Opt. Phys. 49, 23LT01 (2016).
[99] D.M. Reich and L.B. Madsen. Illuminating molecular symmetries with bicircular high-orderharmonic generation. Phys. Rev. Lett. 117, 133902 (2016).
[100] S. Odzaka, E. Hasovic, and D.B. Milosevic. High-order harmonic generation in polyatomic
molecules induced by a bicircular laser field. Phys. Rev. A. 94, 033419 (2016).
[101] K.M. Dorney, J.L. Ellis, C. Hernández-Garcı́a, D.D. Hickstein, C.A. Mancuso, N. Brooks,
T. Fan, G. Fan, D. Zusin, C. Gentry, P. Grychtol, H.C. Kapteyn, and M.M. Murnane. Helicityselective enhancement and polarization control of attosecond high harmonic waveforms driven
by bichromatic circularly polarized laser fields. Phys. Rev. Lett. 119, 063201 (2017).
[102] A. Fleischer, E. Bordo, O. Kfir, P. Sidorenko, and O. Cohen. Polarization-fan high-orderharmonics. J. Phys. B 50, 034001 (2017).
[103] N. Zhavoronkov and M. Ivanov. Extended ellipticity control for attosecond pulses by highharmonic generation. Opt. Lett. 42, 4720 (2017).
[104] E. Pisanty and A. Jimenez-Galan. Extended ellipticity control for attosecond pulses by highharmonic generation. Phys. Rev. A 96, 063401 (2017).
[105] D. Baykusheva, S. Brennecke, M. Lein, and H.J. Wörner. Signatures of electronic structure
in bicircular high-harmonic spectroscopy. Phys. Rev. Lett. 119, 203201 (2017). .

123[106] G. Lerner, T. Diskin, O. Neufeld, O. Kfir, and O. Cohen. Selective suppression of high-orderharmonics within phase-matched spectral regions. Opt. Lett. 42, 1349 (2017). .
[107] D. Ayuso, P. Decleva, S. Patchkovskii, and O. Smirnova. Chiral dichroism in bi-ellipticalhigh-order harmonic generation. J. Phys. B 51, 06LT01 (2018).
[108] G. Dixit, A. Jimenez-Galan, L. Medisauskas, and M. Ivanov. Control of the helicity of highorder harmonic radiation using bichromatic circularly polarized laser fields. Phys. Rev. A 98,
053402 (2018).
[109] L. Barreau, K. Veyrinas, V. Gruson, S.J. Weber, T. Auguste, J.F. Hergott, F. Lepetit,
B. Carre, J.C. Houver, D. Dowek, and P. Salieres. Evidence of depolarization and ellipticity of high harmonics driven by ultrashort bichromatic circularly polarized fields. Nat.
Commun. 9, 4727 (2018).
[110] J.T. Huang P.Y. Huang C.H. Lu L. Rego D.D. Hickstein J.L. Ellis A. Jaron-Becker A. Becker
S.D. Yang C.G. Durfee L. Plaja H.C. Kapteyn M.M. Murnane A.H. Kung P.C. Huang, C.
Hernandez-Garcia and M.C. Chen. Polarization control of isolated high-harmonic pulses.
Nat. Photon. 12, 349 (2018).
[111] D. Ayuso F. Morales S. Patchkovskii M. Schloz E. Pisanty O. Smirnova A. Jimenez-Galan,
N. Zhavoronkov and M. Ivanov. Control of attosecond light polarization in two-color bicircularfields. Phys. Rev. A 97, 023409 (2018).
[112] D.A. Telnov J. Heslar and S.I Chu. Controlling electron quantum paths for generation of
circularly polarized high-order harmonics by H+
2 subject to tailored (ω, 2ω) counter-rotating
laser fields. Phys. Rev. A 97, 043419 (2018).
[113] B. Börning W. Paufler and S. Fritzsche. Tailored orbital angular momentum in high-orderharmonic generation with bicircular laguerre-gaussian beams. Phys. Rev. A 98, 011401(R)(2018).
[114] G.R. Jia M.Z. Li, Y. Xu and X.B. Bian. Controlling polarization of high-order harmonic
generation by molecular alignment in a bicircular laser field. Phys. Rev. A 100, 033410(2019).
[115] L.B. Madsen N.L. Manakov A.V. Meremianin J.M. Ngoko Djiokap, S.X. Hu and A.F. Starace.
Electron vortices in photoionization by circularly polarized attosecond pulses. Phys. Rev.
Lett. 115, 113004 (2015).
[116] S. Chelkowski K.J. Yuan and A.D. Bandrauk. Photoelectron momentum distributions ofmolecules in bichromatic circularly polarized attosecond UV laser fields. Phys. Rev. A 93,
053425 (2016).
[117] D.B. Milosevic and W. Becker. Improved strong-field approximation and quantum-orbit
theory: Application to ionization by a bicircular laser field. Phys. Rev. A 93, 063418 (2016).
[118] D.D. Hickstein J.L. Chaloupka X.M. Tong J.L. Ellis H.C. Kapteyn C.A. Mancuso, K.M. Dorney and M.M. Murnane. Observation of ionization enhancement in two-color circularly polarized laser fields. Phys. Rev. A 96, 023402 (2017).

124[119] D. Pengel, S. Kerbstadt, D. Johannmeyer, L. Englert, T. Bayer, , and M. Wollenhaupt.
Electron vortices in femtosecond multiphoton ionization. Phys. Rev. Lett. 118, 053003 (2017).
[120] M. Busuladzic, A. Gazibegovic-Busuladzic, and D.B. Milosevic. Phys. Rev. A 95, 033411(2017).
[121] V.H. Hoang, V.H. Le, C.D. Lin, and A.-T. Le. Retrieval of target structure information fromlaser-induced photoelectrons by few-cycle bicircular laser fields. Phys. Rev. A 95, 031402(R)(2017).
[122] K. Lin, X. Jia, Z. Yu, F. He, J. Ma, H. Li, X. Gong, Q. Song, Q. Ji, W. Zhang, H. Li, P. Lu,
H. Zeng, J. Chen, , and J. Wu. Comparison study of strong-field ionization of molecules andatoms by bicircular two-color femtosecond laser pulses. Phys. Rev. Lett. 119, 203202 (2017).
[123] M. Abu-samha and L.B. Madsen. Probing atomic and molecular targets by intense bicircularcounter-rotating laser fields. J. Phys. B 51, 135401 (2018).
[124] M. Busuladzic, A. Cerkic, A. Gazibegovic-Busuladzic, E. Hasovic, and D.B. Milosevic.
Molecular-orientation-dependent interference and plateau structures in strong-field ionizationof a diatomic molecule by a corotating bichromatic elliptically polarized laser field. Phys. Rev.
A 98, 013413 (2018).
[125] S. Eckart, M. Kunitski, I. Ivanov, M. Richter, K. Fehre, A. Hartung, J. Rist, K. Henrichs, D. Trabert, N. Schlott, L. Ph.H. Schmidt, T. Jahnke, M.S. Schöffler, A. Kheifets,
and R. Dörner. Subcycle interference upon tunnel ionization by counter-rotating two-color
fields. Phys. Rev. A 97, 041402(R) (2018).
[126] M. Li, W.-C. Jiang, H. Xie, S. Luo, Y. Zhou, and P. Lu. Strong-field photoelectron holographyof atoms by bicircular two-color laser pulses. Phys. Rev. A 97, 023415 (2018).
[127] M. Han, P. Ge, Y. Shao, Q. Gong, and Y. Liu. Attoclock photoelectron interferometry with
two-color corotating circular fields to probe the phase and the amplitude of emitting wave
packets. Phys. Rev. Lett. 120, 073202 (2018).
[128] S. Eckart, K. Fehre, N. Eicke, A. Hartung, J. Rist, D. Trabert, N. Strenger, A. Pier, L.Ph.H.
Schmidt, T. Jahnke, M.S. Schöffler, M. Lein, M. Kunitski, and R. Dörner. Direct experimental
access to the nonadiabatic initial momentum offset upon tunnel ionization. Phys. Rev. Lett.
121, 163202 (2018).
[129] N. Eicke and M. Lein. Attoclock with counter-rotating bicircular laser fields. Phys. Rev. A
99, 031402(R) (2019).
[130] P. Ge, M. Han, Y. Deng, Q. Gong, and Y. Liu. Universal description of the attoclock with
two-color corotating circular fields. Phys. Rev. Lett. 122, 013201 (2019).
[131] S. Kerbstadt, K. Eickhoff, T. Bayer, and M. Wollenhaupt. Odd electron wave packets from
cycloidal ultrashort laser fields. Nat. Commun. 10, 658 (2019).
[132] M. Abu-samha and L.B. Madsen. Pulse-length effects in strong-field ionization of atoms by
co-rotating and counter-rotating bicircular laser pulses. Phys. Rev. A 100, 043415 (2019).

125[133] M.Y. Ma, J.P. Wang, W.Q. Jing, Z. Guan, Z.H. Jiao, G.L Wang, J.H. Chen, , and S.F. Zhao.
Controlling the atomic-orbital-resolved photoionization for neon atoms by counter-rotating
circularly polarized attosecond pulses. Opt. Express 29, 33245-33256 (2021).
[134] R.R. Wang, M.Y. Ma, L.C. Wen, Z. Guan, Z.Q. Yang, Z.H. Jiao, G.L. Wang, and S.F.
Zhao. Comparative study of electron vortices in photoionization of molecules and atoms bycounter-rotating circularly polarized laser pulses. J. Opt. Soc. Am. B 40, 1749-1755 (2023).
[135] C.A. Mancuso, K.M. Dorney, D.D. Hickstein, J.L. Chaloupka, J.L. Ellis, F.J. Dollar, R. Knut,
P. Grychtol, D. Zusin, C. Gentry, M. Gopalakrishnan, H.C. Kapteyn, and M.M. Murnane.
Controlling nonsequential double ionization in two-color circularly polarized femtosecond laser
fields. Phys. Rev. Lett. 117, 133201 (2016).
[136] S. Eckart, M. Richter, M. Kunitski, A. Hartung, J. Rist, K. Henrichs, N. Schlott, T. Bauer
H. Kang, H. Sann, L.Ph.H. Schmidt, M. Schöffler, T. Jahnke, and R. Dörner. Nonsequential
double ionization by counterrotating circularly polarized two-color laser fields. Phys. Rev.
Lett. 117, 133202 (2016).
[137] S. Ben, P.Y. Guo, X.F. Pan, T.T. Xu, K.L. Song, and X.S. Liu. Recollision induced excitationionization with counter-rotating two-color circularly polarized laser field. Chem. Phys. Lett.
679, 38 (2017).
[138] J.M. Ngoko Djiokap, A.V. Meremianin, N.L. Manakov, S.X. Hu, L.B. Madsen, and A.F.
Starace. Kinematical vortices in double photoionization of helium by attosecond pulses.
Phys. Rev. A 96, 013405 (2017).
[139] W.W. Yu, S. Ben, L. Yi, Y.S. Liu, J. Guo, and X.S. Liu. Enhancement of nonsequential
double ionization in counter-rotating two-color circularly polarized laser fields. Chem. Phys.
Lett. 706, 62 (2018).
[140] C. Huang, M. Zhong, and Z. Wu. Intensity-dependent two-electron emission dynamics innonsequential double ionization by counter-rotating two-color circularly polarized laser fields.
Opt. Express 26, 26045 (2018).
[141] X. Ma, Y. Zhou, Y. Chen, M. Li, Y. Li, Q. Zhang, and P. Lu. Timing the release of thecorrelated electrons in strong-field nonsequential double ionization by circularly polarized
two-color laser fields. Opt. Express 27, 1825 (2019).
[142] M. Peng, L.H. Bai, and Z. Guo. Influence of relative phase on the nonsequential double
ionization process of CO2 molecules by counter-rotating two-color circularly polarized laser
fields. Commun. Theor. Phys. 73, 075501 (2021).
[143] Z. Chen, J.S., X. Zeng, X. Huang, Y. Li, and C. Huang. Electron angular correlation in
nonsequential double ionization of molecules by counter-rotating two-color circularly polarized
fields. Opt. Express 29, 29576-29586 (2021).
[144] J. Su, Z. Liu, J. Liao, X. Huang, Y. Li, and C. Huang. Electron correlation and recollision dynamics in nonsequential double ionization by counter-rotating two-color elliptically polarizedlaser fields. Opt. Express 30, 24898-24908 (2022).
[145] K.J. Yuan and A.D. Bandrauk. Attosecond-magnetic-field-pulse generation by electronic
currents in bichromatic circularly polarized uv laser fields. Phys. Rev. A 92, 063401 (2015).

126[146] G. Buica. Circular dichroism in angular distribution of electron-hydrogen scattering in a
two-color bicircular laser field. Phys. Rev. A 98, 053427 (2018).
[147] J. Guo, K.J. Yuan, H. Lu, and A.D. Bandrauk. Spatiotemporal evolution of ultrafast
magnetic-field generation in molecules with intense bichromatic circularly polarized uv laser
pulses. Phys. Rev. A 99, 053416 (2019).
[148] I. Barth and O. Smirnova. Nonadiabatic tunneling in circularly polarized laser fields: Physicalpicture and calculations. Phys. Rev. A 84, 063415 (2011).
[149] A.K. Kazansky, A.V. Grigorieva, and N.M. Kabachnik. Dichroism in short-pulse two-color
XUV plus IR multiphoton ionization of atoms. Phys. Rev. A 85, 053409 (2012).
[150] T. Herath, L. Yan, S.K. Lee, and W. Li. Strong-field ionization rate depends on the sign ofthe magnetic quantum number. Phys. Rev. Lett. 109, 043004 (2012).
[151] J.H. Bauer, F. Mota-Furtado, P.F. O’Mahony, B. Piraux, and K. Warda. Ionization andexcitation of the excited hydrogen atom in strong circularly polarized laser fields. Phys. Rev.
A 90, 063402 (2014).
[152] I. Barth and M. Lein. Numerical verification of the theory of nonadiabatic tunnel ionizationin strong circularly polarized laser fields. J. Phys. B 47, 204016 (2014).
[153] C.H.R. Ooi, W.L. Ho, and A.D. Bandrauk. Photoelectron angular distributions of excitedatoms in intense laser fields. Phys. Rev. A 90, 013417 (2014).
[154] A. Hartung, F. Morales, M. Kunitski, K. Henrichs, A. Laucke, M. Richter, T. Jahnke,
A. Kalinin, M. Schöffler, L.Ph.H. Schmidt, M. Ivanov, O. Smirnova, and R. Dörner. Electron
spin polarization in strong-field ionization of Xenon atoms. Nat. Photon. 10, 526 (2016).
[155] N. Douguet, A.N. Grum-Grzhimailo, E.V. Gryzlova, E.I. Staroselskaya, J. Venzke, and
K. Bartschat. Photoelectron angular distributions in bichromatic atomic ionization inducedby circularly polarized VUV femtosecond pulses. Phys. Rev. A 93, 033402 (2016).
[156] J. Wätzel and J. Berakdar. Discerning on a sub-optical-wavelength the attosecond time delays
in electron emission from magnetic sublevels by optical vortices. Phys. Rev. A 94, 033414(2016).
[157] J.P. Wang and F. He. Tunneling ionization of neon atoms carrying different orbital angularmomenta in strong laser fields. Phys. Rev. A 95, 043420 (2017).
[158] Q. Zhang, G. Basnayake, A. Winney, Y.F. Lin, D. Debrah, S.K. Lee, and W. Li. Orbitalresolved nonadiabatic tunneling ionization. Phys. Rev. A 96, 023422 (2017).
[159] S. Eckart, M. Kunitski, M. Richter, A. Hartung, J. Rist, F. Trinter, K. Fehre, N. Schlott,
K. Henrichs, L.Ph.H. Schmidt, T. Jahnke, M. Schöffler, K. Liu, I. Barth, J. Kaushal,
F. Morales, M. Ivanov, O. Smirnova, and R. Dörner. Ultrafast preparation and detectionof ring currents in single atoms. Nat. Phys. 14, 701 (2018).
[160] K. Liu, H. Ni, K. Renziehausen, J.M. Rost, and I. Barth. Deformation of atomic p± orbitals
in strong elliptically polarized laser fields: Ionization time drifts and spatial photoelectron
separation. Phys. Rev. Lett. 121, 203201 (2018).

127[161] D. Trabert, A. Hartung, S. Eckart, F. Trinter, A. Kalinin, M. Schöffler, L.Ph.H. Schmidt,
T. Jahnke, M. Kunitski, and R. Dörner. Spin and angular momentum in strong-field ionization. Phys. Rev. Lett. 120, 043202 (2018).
[162] A.H.N.C. De Silva, T. Moon, K.L. Romans, B.P. Acharya, S. Dubey, K. Foster, O. Russ,
C. Rischbieter, N. Douguet, K. Bartschat, and D Fischer. Circular dichroism in atomic
resonance-enhanced few-photon ionization. Phys. Rev. A 103, 053125 (2021).
[163] S. Walker, L. Kolanz, J. Venzke, and A. Becker. Enhanced ionization of counter-rotating
electrons via doorway states in ultrashort circularly polarized laser pulses. Phys. Rev. A 103,
L061101 (2021).
[164] S. Walker, L. Kolanz, J. Venzke, and A. Becker. Selectivity in electron emission induced byultrashort circularly polarized laser pulses. Phys. Rev. Research 3, 043051 (2021).
[165] A.N. Artemyev, E. Kutscher, B.M. Lagutin, and P.V. Demekhin. Theoretical study of spinpolarization in multiphoton ionization of Xe. J. Chem. Phys. 158 15 (2023).
[166] T. Nubbemeyer, K. Gorling, A. Saenz, U. Eichmann, and W. Sandner. Strong-field tunnelingwithout ionization. Phys. Rev. Lett. 101, 233001 (2008).
[167] U. Eichmann, A. Saenz, S. Eilzer, T. Nubbemeyer, and W. Sandner. Observing Rydberg
atoms to survive intense laser fields. Phys. Rev. Lett. 110, 203002 (2013).
[168] M. Chini, X. Wang, Y. Cheng, H. Wang, Y. Wu, E. Cunningham, P.C. Li, J. Heslar, D.A.
Telnow, S.I Chu, and Z. Chang. Coherent phase-matched VUV generation by field-controlled
bound states. Nat. Photon. 8, 437 (2014).
[169] H. Zimmermann, J. Buller, S. Eilzer, and U. Eichmann. Strong-field excitation of helium:
Bound state distribution and spin effects. Phys. Rev. Lett. 114, 123003 (2015).
[170] R.R. Freeman, P.H. Bucksbaum, H. Milchberg, S. Darack, D. Schumacher, and M.E. Geusic.
Above-threshold ionization with subpicosecond laser pulses. Phys. Rev. Lett. 59, 1092 (1987).
[171] M.D. Perry, A. Szoke, and K.C. Kulander. Resonantly enhanced above-threshold ionizationof helium. Phys. Rev. Lett. 63, 1058 (1989).
[172] P. Agostini, P. Breger, A. L’Huillier, H.G. Muller, G. Petite, A. Antonetti, and A. Migus.
Giant stark shifts in multiphoton ionization. Phys. Rev. Lett. 63, 2208 (1989).
[173] W.H. Xiong, J.-Z. Jin, L.-Y. Peng, and Q. Gong. Numerical observation of two sets of
low-order harmonics near the ionization threshold. Phys. Rev. A 96, 023418 (2017).
[174] J. Venzke, Y. Gebre, A. Becker, and A. Jaron-Becker. Pathways to excitation of atoms withbicircular laser pulses. Phys. Rev. A 101, 053425 (2020).
[175] H.R. Gray, R.M. Whitley, and C.R. Stroud. Coherent trapping of atomic populations. Opt.
Lett. 3, 218 (1978).
[176] M.V. Fedorov, M.M. Tehranchi, and S.M. Fedorov. Interference stabilization of Rydberg
atoms: numerical calculations and physical models. J. Phys. B: At. Mol. Opt. Phys. 29, 2907
(1996).

128[177] D.G. Arbo, K.I. Dimitriou, E. Persson, and J. Burgdörfer. Sub-poissonian angular momentum
distribution near threshold in atomic ionization by short laser pulses. Phys. Rev. A 78, 013406(2008).
[178] C.A. Mancuso, K.M. Dorney, D.D. Hickstein, J.L. Chaloupka, X.M. Tong, J.L. Ellis, H.C.
Kapteyn, and M.M. Murnane. Observation of ionization enhancement in two-color circularly
polarized laser fields. Phys. Rev. A 96, 023402 (2017).
[179] Y. Gebre, J. Venzke, A. Jaron-Becker, and A. Becker. Transitions between Rydberg statesin two-color corotating circularly polarized laser pulses. Phys. Rev. A 103, 013101 (2021).
[180] E. Arimondo and G. Orriols. Nonabsorbing atomic coherences by coherent two-photon transitions in a three-level optical pumping. Lett. Nuovo Cimento 17, 333–338 (1976).
[181] H. R. Gray, R. M. Whitley, and C. R. Stroud. Coherent trapping of atomic populations. Opt.
Lett. 3, 218-220 (1978).
[182] A. Imamoǧlu and S.E. Harris. Lasers without inversion: interference of dressed lifetimebroadened states. Opt. Lett. 14, 1344-1346 (1989).
[183] U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann. Population transfer between
molecular vibrational levels by stimulated Raman scattering with partially overlapping laser
fields. A new concept and experimental results. J. Chem. Phys. 92, 5363–5376 (1990).
[184] P. Marte, P. Zoller, and J. L. Hall. Coherent atomic mirrors and beam splitters by adiabatic
passage in multilevel systems. Phys. Rev. A 44, R4118 (1991).
[185] K.J. Boller, A. Imamoğlu, and S.E. Harris. Observation of electromagnetically induced transparency. Phys. Rev. Lett 66, 2593 (1991).
[186] J. E. Field, K. H. Hahn, and S. E. Harris. Observation of electromagnetically induced transparency in collisionally broadened lead vapor. Phys. Rev. Lett 67, 3062 (1991).
[187] M.V. Fedorov and A.M. Movsesian. Field-induced effects of narrowing of photoelectron spectra and stabilisation of Rydberg atoms. J. Phys. B: At. Mol. Opt. Phys. 21, L155 (1998).
[188] M.V. Fedorov, M.M. Tehranchi, and S.M. Fedorov. Interference stabilization of Rydberg
atoms: numerical calculations and physical models. J. Phys. B: At. Mol. Opt. Phys. 29, 2907
(1996).
[189] A. Talebpour, C.Y. Chien, and S.L. Chin. Population trapping in rare gases. J. Phys. B: At.
Mol. Opt. Phys. 29, 5725 (1996).
[190] A. Talebpour, Y. Liang, and S. L. Chin. Population trapping in the CO molecule. J. Phys.
B: At. Mol. Opt. Phys. 29, 3435 (1996).
[191] D. Peng, B. Wu, P. Fu, B. Wang, J. Gong, and Z.C. Yan. Sensitive frequency dependenceof the carrier-envelope phase effect on bound-bound transitions: An interference perspective.
Phys. Rev. A 82, 053407 (2010).
[192] Z. Zhai, D. Peng, X. Zhao, F. Guo, Y. Yang, P. Fu, J. Chen, Z.C. Yan, and B. Wang. Carrierenvelope-phase effect on laser-driven bound-bound transitions in the high-frequency region.
Phys. Rev. A 86, 043432 (2012).

129[193] S. Patsch, D.M. Reich, J.M. Raimond, M. Brune, S. Gleyzes, and C.P. Koch. Fast andaccurate circularization of a Rydberg atom. Phys. Rev. A 97, 053418 (2018).
[194] J. Solanpää and E. Räsänen. Control of Rydberg-state population with realistic femtosecond
laser pulses. Phys. Rev. A 98, 053422 (2018).
[195] P. Lambropoulos. Spin-orbit coupling and photoelectron polarization in multiphoton ionization of atoms. Phys. Rev. Lett 30, 413 (1973).
[196] E.H.A Granneman, M. Klewer, K.J. Nygaard, and M.J. Van der Wiel. Polarization effects inresonant two-photon ionization of Caesium. J. Phys. B: At. Mol. Phys. 9, L87 (1976).
[197] Y. Li, P. Lan, H. Xie, M. He, X. Zhu, Q. Zhang, and P. Lu. Nonadiabatic tunnel ionizationin strong circularly polarized laser fields: counterintuitive angular shifts in the photoelectron
momentum distribution. Opt. Express 23, 28801-28807 (2015).
[198] D.B. Milošević, G.G. Paulus, D. Bauer, and W. Becker. Above-threshold ionization by fewcycle pulses. J. Phys. B: At. Mol. Opt. Phys. 39, R203 (2006).
[199] J. Heinz, R. Reiff, T. Joyce, A. Becker, A. Jaroń-Becker, and M. Kolesik. Nonlinear polarization and ionization in O2 : metastable electronic state model. Opt. Express 28, 25022-25036(2020).
[200] R. Reiff. Atomic characteristics and high-order harmonic spectra: extension of ab-initio
numerical calculations to larger. PhD thesis, University of Colorado at Boulder, 2021.
[201] B.H. Bransden and C.J. Joachain. Physics of Atoms and Molecules. Pearson Education.
Prentice Hall (2003).
[202] F.W.R Bader. An introduction to the electronic structure of atoms and molecules. McMaster
University / Hamilton, Ontario.

ProQuest Number: 30574915
INFORMATION TO ALL USERS
The quality and completeness of this reproduction is dependent on the quality
and completeness of the copy made available to ProQuest.

Distributed by ProQuest LLC ( 2023 ).
Copyright of the Dissertation is held by the Author unless otherwise noted.

This work may be used in accordance with the terms of the Creative Commons licenseor other rights statement, as indicated in the copyright statement or in the metadata
associated with this work. Unless otherwise specified in the copyright statement
or the metadata, all rights are reserved by the copyright holder.

This work is protected against unauthorized copying under Title 17,
United States Code and other applicable copyright laws.

Microform Edition where available © ProQuest LLC. No reproduction or digitizationof the Microform Edition is authorized without permission of ProQuest LLC.

ProQuest LLC
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346 USA